## Thursday, October 29, 2015

### Integer factoring

Factoring of large integers is a great challenge in computational mathematics.
When the numbers are very large, no efficient, non-quantum integer factorization algorithm is known; an effort by several researchers concluded in 2009, factoring a 232-digit number (RSA-768), utilizing hundreds of machines over a span of two years. Wikipedia
However, there are efficient and simple algorithms for cell size numbers (about 19 digit numbers in 64 bit systems) that immediately delivers the prime factorization. One such algorithm is Pollard rho method.

: func ( numb n -- m ) dup um* 1 0 d+ rot um/mod drop ;

: pollard1 ( numb start -- pfac flag )
\ numb is an odd composite

dup dup loc{ numb start alpha beta } 0
begin drop numb alpha func to alpha
numb dup beta func func to beta
alpha beta - abs        \ |alpha-beta|
numb ugcd dup numb =    \ gcd flag
if false exit
then dup 1 = 0=         \ gcd<>1
until true ;

: sqrtf \ m -- n             \ floor

0 d>f fsqrt f>s ;

: sqrtfslow \ m -- n

dup 4 u< if dup
if drop 1 then exit then
dup 1- >r 1 rshift
begin r@ over 0 swap um/mod nip
over + 1 rshift tuck u> 0=
until r> drop ;

: sqrtc \ m -- n             \ ceiling

1- sqrtf 1+ ;

: pollard2 \ numb -- prime numb>1
dup 1 and 0= if drop 2 exit then
dup prime if exit then
dup sqrtf dup * over =
if sqrtf exit then -1 2    \ numb 100 0
do dup i pollard1          \ numb pfac flag
if leave                \ numb pfac
then drop               \ numb
loop nip ;                 \ pfac

: pollard \ n -- p1 ... pk

dup pollard2 2dup =        \ n q f=
if drop exit               \ n
then dup >r
0 swap um/mod nip recurse  \ q n/q
r> recurse ;

4267640728818788929 pollard .s 1234567891 3456789019  ok

12345678987654321 pollard cr .s
333667 37 3 3 333667 37 3 3  ok

Sometimes the greatest prime factor is of interest and I guess there is no other way to find it than factoring and sorting the stack.

: lower \ m1 ... mn m n+1 -- m1 ... m ... mi n+1

\ lower m on stack until next number not is greater
dup 2 =
if drop 2dup u>
if swap
then 2 exit
then >r 2dup u> 0=
if r> exit
then r> rot >r
1- recurse 1+ r> swap ;

: sort \ m1 ... mn n -- s1 ... sn n o(n²)

dup 3 <
if dup 2 =
if drop 2dup u>
if swap
then 2
then exit
then dup >r
1- recurse roll
r> lower ;

But to use sort one have to know the number of factors to sort.

: pollard# \ numb -- p1 ... pk k
depth >r
pollard depth r> - 1+ ;

12345678987654321 pollard# cr .s
333667 37 3 3 333667 37 3 3 8  ok
sort .s 3 3 3 3 37 37 333667 333667 8  ok
over . 333667  ok