Prime Factors without using / Div or Mod.

[Richard]

The aim here was to find prime factors without using / \ or Mod.
It can then run on the simplest hardware or microcontroller.
Has anyone seen any similar / \ Div free code ?

The number here is represented in the form; n = a * b + c.
It searches for c = 0, when a and b will be factors.

This compact code can run fast in the CPU registers.
Comparisons are fast because all registers are the same length.
It produces a list of prime factors increasing in size.

This simple algorithm could be improved. It is inefficient because:
1. It goes around the loop very many times, so it does not like really big numbers.
2. It tests all even candidates, when it could do a two step.
3. After finding a prime factor, p, it begins again from scratch, not at p.

Here is the code inside a loop that exercises it for a block of numbers.

Code: Select all

' Prime Factors without using DIV or MOD
' n = ( a * b ) + c, when c = 0, a and b are factors
Dim As Ulong a, b, c, n, i = 2^20
Print " Prime Factors."
For n = i - 9 To i + 9
 
 Print n; " = ";
 
 a = n ' n = a * b + c
 b = 1 ' trial factor
 c = 0 ' mod, remainder
 
 Do 
 a -= 1
 c += b
 If c >= a Then ' carry
 b += 1
 c -= a
 End If 
 
 If c = 0 Then ' b is a factor
 Print b; " * ";
 b = 1
 End If
 
 Loop While b < a ' square root
 
 Print a * b + c ' terminal factor may be 1
 
Next n
Print " Done."
Sleep

[srvaldez]

interesting Richard, but it would probably be too slow for integers larger than 32-bit, unless there's some way to make it massively parallel

[coderJeff]

Richard, wow, that's compact.
Clever. I never would have thought of Loop While b < a ' square root to check bounds on the square root. Good tip!

2. It tests all even candidates, when it could do a two step.

yeah, pulling out the 2's first will give a number that has only odd factors allow to do steps of 2.

Code: Select all

	'' factor out all the 2's
	while( (n and 1) = 0 )
		print " 2 * ";
		n shr= 1
	wend

Has anyone seen any similar / \ Div free code ?

In my own code years ago. Though I was trying find the 2 factors of n closest to sqr(n) using only add/sub/cmp in the inner loop, so the flow is different. I visualized the problem as the graph n=xy, and tracked a trial product for comparisons so it involved a multiply instruction in the set-up. Plus I started at sqr(n) so also involved an integer approximation of the square root. After looking at your code I feel like I should be able recheck the math to get rid of more comparisons.

Code: Select all

function SquareRoot( byval n as ulong ) as ulong
	dim as ulong x = 0ul
	dim as ulong b = &h8000ul
	while( b )
		dim as ulong d = (x or b) 
		dim as ulong dd = d * d
		if( dd <= n ) then
			x = d
		end if
		b shr= 1
	wend
	return x	  
end function 

sub Factor( byval n as ulong, byref x as ulong, byref y as ulong )
	dim c as ulong = any
	x = SquareRoot( n )
	y = x
	c = x * y
	while( c <> n )
		if( c < n ) then
			x += 1
			c += y
		end if
		if( c > n ) then
			 y -= 1
			 c -= x
		end if
	wend
end sub 

dim n as ulong = (2ul * 2 * 3 * 5 * 5 * 13 * 13 * 29 * 101)
dim as ulong x, y

Factor( n, x, y )
print n & " = " & x & " * " & y
print n & " = " & (x * y) 

sleep

[Richard]

I'm on holiday exploring the wide algorithm space. Obviously the linear search will run slow for big numbers, but I am not homesick for speed, so I will not abandon the holiday immediately. Also; George Polya wrote. “It is better to solve one problem in 5 different ways than to solve 5 different problems”.

This is an absolute minimum solution that clearly demonstrates the algorithm. I have more complex optimised code that runs faster, but the changes hide the incredibly simple geometry of the algorithm. The fact that it works as well as it does surprised me.

The square root search limit "trick" is neat in that it works without recomputation for the decreasing residue as the smaller factors are removed.

I am annoyed by the spurious ultimate factor of one, which only occurs if the penultimate factor was two. There are a few solutions to eliminate the factor of 1, all of which waste clock cycles, increase complexity and confuse the programmer. So I will let that sleeping dog lie, until I can first cast out twos and go double stepping in the search.

One dilemma is application of the no Div or Mod rule. Should it cover all, or only the inner loop. The use of Div or Mod, to efficiently restart the algorithm after finding a factor, may or may not be acceptable in some final forms. As a purist on holiday, I imagine building the algorithm from logic chips in hardware, which precludes Div or Mod anywhere.

[jj2007]

For your inspiration: we had some fun six years ago

[Richard]

Richard, wow, that's compact.
Clever. I never would have thought of Loop While b < a ' square root to check bounds on the square root. Good tip!After looking at your code I feel like I should be able recheck the math to get rid of more comparisons.

Here is my code to find the rectangle closest to a square. It eliminates the square root.

Code: Select all

' Given n, find the rectangle x * y, closest to a square, using add & sub only
'   n = x * y + r;   when r = 0, x and y are integer factors

Dim As Ulong n = 24 ' 1000 * 1024 ' 149 ' 147 ' 120 '

Dim As Ulong x, y, r, u, v ' width, height, residue, last solution found
x = n   ' a wide slab
y = 1   '   of one layer
r = 0   ' with no residue

Do while x >= y                     ' will stop after the square
    if r = 0 then u = x : v = y     ' latest valid solution found
    x -= 1                          ' reduce width
    r += y                          ' accumulate residue until
    If r >= x Then r -= x : y += 1  ' allocate residue to build next layer
Loop

Print  u * v; " = "; u; " * ";  v
Print " Done."
Sleep

Edit:
Terse means brief, or using very few words.
A terse reply or command may seem rude or unfriendly.
Curt means brief or terse, especially to the point of being rude.

This code is curt. Like a politician, it is rude.
It will always find the right solution in the end, but only after trying every possible alternative.

[coderJeff]

Here is my code to find the rectangle closest to a square. It eliminates the square root.

Very cool. Parallelizing and cpu cycles aside, I find elegance in your minimalistic solution. I don't know how to create code / mathematical proofs, but if I look at the rules for each step of the iteration, I can see that it works for any value of N without underflow / overflow, etc. And I can see that if one iteration (or several) were put on a logic chip then the outputs of one chip could feed directly in to the next. Neato.

[Richard]

Here is a faster version of the prime factors code, with two stepping, and without back tracking.

Code: Select all

' Prime Factors of a number, without using MOD to test within the inner loop
'
' n = ( a * b ) + c, when c = 0, a and b are factors of n
'
' Features; This code shifts out factors of two, before testing for odd candidates
' It only uses a slow Div after finding and removing a prime factor
'
' WARNING: This code does appear to work, but is unproven
' my fear is that it may yield a false prime, that is actually composite
'
'-----------------------------------------------------------------------
Dim As Ulong n, i = 2^24 ' 2^31
Print " Prime Factors."
For n = i - 9 To i + 9 ' test a block of integers either side of i
 If n = 0 Then
 Print " Zero."
 Exit For
 End If
 
 '-------------------------------------------------------------------
 ' prime factor search algorithm starts here
 
 Dim As Ulong a, b, c, check = 1 ' will check product of factors as found
 Print n; " = "; 
 
 ' first cast out all factors of two
 a = n
 Do Until ( a And 1 ) ' until odd
 a Shr= 1
 Print "2 * ";
 check *= 2
 Loop
 ' a is now quite odd
 
 ' the problem has now been reduced to odd factors only
 c = a ' c will be the remainder after the a * b rectangle is removed
 b = 3 ' three will be the first factor to test
 a \= b ' fast precompute a, which has a 50% chance of a becoming even
 If a Then a -= 1 ' so force a
 a Or= 1 ' to be odd again
 c -= a * b ' remainder will be in c
 
 '----------------------------------------------------------
 ' n = a * b + c ; when b is a prime factor, c is zero ; 
 ' we now manipulate a rectangle of width a, and height b, remainder c
 
 ' main loop employs fast primitive operations on registers inside CPU
 Do ' removes factors, restarting search with previous factor found
 
 Do While c = 0 ' remove all repetitions of prime factor b
 Print b; " * "; ' b is a prime factor
 check *= b
 
 c = a ' remainder for later
 a \= b ' a has a 50% chance of becoming even
 If a Then a -= 1 ' must avoid underflow while
 a Or= 1 ' forcing a to be odd again
 c -= a * b ' remainder, test again for repeat of factor b
 Loop
 
 If a <= b Then Exit Do ' has passed the square root, so is final prime
 
 ' advance to next odd, to test if possible prime factor
 a -= 2 ' demolish two columns of height b, accumulate the
 c += b + b ' debris with the remainder in c
 
 If c >= ( a + a ) Then ' enough to make two more layers
 b += 2 ' jump up to next odd factor
 c -= a + a ' balance the remainder account
 End If ' b is a prime factor when remainder c is zero
 
 Loop
 Print Culng( a * b + c ); ' terminal prime, but may be 1
 
 '----------------------------------------------------------
 ' verify product of all factors
 check *= a * b + c
 If check <> n Then Print " WARNING Error "; : Sleep
 Print ' end line of factors
 
Next n
Print " Done."
Sleep

[Munair]

I see "Dim As Ulong a, b, c, check = 1" in your for-loop. Is that with performance in mind?

[Richard]

@Munair
I don't think performance is that important at the moment. I was trying to understand the two-stepping algorithm well enough to write code that might work. The code needed is now more obvious. There are enough comments remaining that I might be able to follow the logic again next week.

Those variables are needed for the algorithm, and check = 1 is needed for each new n during testing, so when I throw out the exercise loop, the Dim statement will remain with the code that uses it.

Since the prime factors come from three separate places in the code, a function() will need to return an array of prime factors. The three variables used will be local to that function.

[coderJeff]

@Richard, I messed around with your original post to add the stepping by 2 and keep the spirit of not using mod or div. The 'check' variable is a good idea for testing, but I didn't add it here.

Code: Select all

function factor( byval n as ulong ) as ulong

    '' Prime Factors without using DIV or MOD
    '' n = ( a * b ) + c, when c = 0, a and b are factors
    '' as each factor is found, divide n by b and find
    '' the factors for the remaining n
    ''
    '' return the number of prime factors (excluding 1)

    dim As ulong a, b, c, count = 0

    '' only print '*' symbol if needed
    #macro printFactor( x )
        if( count > 0 ) then
            print " *";
        end if
        print " " & x;
    #endmacro

    '' original number
    print n; " =";

    '' special cases - terminate if 0 or 1
    '' 0 causes inifite loop when checking for 2's
    '' 1 won't be printed when we fall past the main loop
    if( n <= 1 ) then
        printFactor( n )
        print
        return count   '' return number of prime factors (0)
    end if

    '' factor out all the 2's
    while( (n and 1) = 0 )
        n shr= 1
        printFactor( 2 )
        count += 1
    wend

    '' by now n is odd, so all factors must be odd
    a = n  '' n = a * b + c
    b = 1  '' trial factor
    c = 0  '' mod, remainder
     
    '' only while square root not exceeded
    while( b < a )
        a -= 2
        c += b
        c += b

        '' check c >= (a + a),
        '' but avoid 32-bit overflow of (a+a) for n > 2^31
        if( c >= a ) then
            if( (c - a) >= a ) then '' carry
                b += 2
                c -= a
                c -= a
            end if
        end if

        '' c = 0? then b is a factor?
        if c = 0 Then
            printFactor( b )
            count += 1
            '' dividing n=a*b+c by factor b
            '' gives n'=a because c=0
            n = a
            b = 1
            '' continue iterations with new value of n to
            '' find the remaining factors
        end If
    wend

    '' show the last not yet printed factor?
    if( n > 1 ) then
        printFactor( n )
        count += 1
    end if
    print
    return count
end function


factor( 0 )
factor( 1 )
factor( 2 )
factor( 3 )
factor( 4 )
factor( 5 )
factor( 8 )
factor( 8 * 3 )
factor( 3 * 7 )
factor( 2 * 7 * 11 )
factor( 3 * 7 * 11 )

dim as ulong i = 2^20
For n as integer = i - 5 To i + 5
    factor( n )
next

sleep

[Richard]

@coderJeff
As a purist, I really like it.
I also like it when people drop rocks in my wheel ruts, to help me escape onto alternative tracks.

[dodicat]

Hi Richard.
I see you use integer\division, is this allowed now?

[Richard]

@dodicat.
My aim was to eliminate the slow Div or Mod test from the inside loop. I have done that. Now I am looking at ways to speed up the resulting envelope by eliminating all silly re-computation.

The slowest part of the process was the duplicated counting up again to the last factor found. That first repeated loop starting at one, is actually just a very slow algorithm for Div 3. If instead of starting at one (coderJeff's purist code) we started at three, the first loop would be 1/3 as far, in steps of 3, so nine times faster.

For that reason I allow very limited Div to be used, once at the start of the outside loop, then again once after each factor is found. Since factors of two were pre-stripped, the worst case Div count will be for numbers of the form 3^m. For 32 bit Ulong, that takes a worst case maximum of;
1 + ( Log( 2^32 ) / Log( 3 ) ) = 21.

That limited use of Div, a maximum of 21 times, removes all that duplication and makes the whole process fast and practical. That was the code that I posted earlier using Div.

If I used a Sub Div(a,b), written using primitives, would that be an acceptable solution ?

[dodicat]

I'll work on the one div I use here.
I have several ways, all very slow.
Mod is not really needed anyway, only one value of mod is required, (zero), so that is easily got around.
This is an old routine partially revived:

Code: Select all

Function getprimefactors(number As Ulongint) As long

      #macro fill(x)
      Print x;" "; 'normally fill an array in here
      check*=x
      #endmacro
      
      #macro Eliminate(x)
      while true
      k=number\x   '<<< ----------------  one div only
      if k*x=number then
            number=k
            fill(x)
      else
            exit while
      end if
      wend
      #endmacro
      
      Print number;" ",
      Dim As Ulongint original=number,check=1,divisor,k
            Eliminate(2ull)
            Eliminate(3ull)
      While divisor*divisor<=number
            divisor+=6 
            Eliminate((divisor-1))
            Eliminate((divisor+1))
      Wend
      If number>1 Then fill(number)
      print
      Return check=original
End Function

Dim As Double t=Timer
For n As Ulongint=2^31-9 To 2^31+9
      If getprimefactors(n)=0 Then Print "error"
Next n
Print
Print Timer-t;" seconds"
Sleep