640 lines
16 KiB
Text
640 lines
16 KiB
Text
// ------------------------------------------------------------
|
|
// Fermat numbers with factorization
|
|
//
|
|
// Using FutureBasic 7.0.36
|
|
// November 2025, R.W.
|
|
//
|
|
// This code is larger than I originally imagined but this is
|
|
// simply because my original coding idea was to use
|
|
// Pollard-Rho to do the factorization until I realize, it
|
|
// can be a very long time for factors of F[9] and above using
|
|
// PR and it may still returns nothing. So I did a pseudo ECM
|
|
// reasonably small but works faster.
|
|
// ------------------------------------------------------------
|
|
|
|
include "gmp.incl"
|
|
|
|
CFTimeInterval tim
|
|
|
|
// Global mpz_t bigints reused for speed
|
|
mpz_t fm // work
|
|
mpz_t num, den, inv, t1, t2, g
|
|
|
|
// Limit is set at 15 since anything greater than 14 most
|
|
// likely will exceed the 'under 5 minutes' Rosetta task limit
|
|
mpz_t gF(15) // Fermat numbers F0..F15
|
|
|
|
mpz_t gFactors(31) // factor buffer
|
|
int gFactorCount
|
|
|
|
long gFN // next Fermat index for factoring
|
|
int gDidFactorHeader // print header once
|
|
|
|
// -----------------------------------------------------------
|
|
// Helper: CFString "F0", "F1", ...
|
|
// -----------------------------------------------------------
|
|
local fn FermatLabel( n as long ) as CFStringRef
|
|
CFStringRef out
|
|
out = @"\nF"
|
|
out = concat( out, str( n ) )
|
|
end fn = out
|
|
|
|
// -----------------------------------------------------------
|
|
// Helper: Abbreviate big numeric CFString: head...tail
|
|
// -----------------------------------------------------------
|
|
local fn AbbrevString( s as CFStringRef, head as int, ¬
|
|
tail as int ) as CFStringRef
|
|
|
|
if s == 0 then RETURN @""
|
|
CFStringRef a, b, out
|
|
int nLen
|
|
nLen = len( s )
|
|
if nLen <= head + tail then RETURN s
|
|
a = left( s, head )
|
|
b = right( s, tail )
|
|
out = @""
|
|
out = concat( a, @"...", b )
|
|
end fn = out
|
|
|
|
// -----------------------------------------------------------
|
|
// Helper: Get mpz_t result and abbreviate if > 53 digits
|
|
// -----------------------------------------------------------
|
|
local fn Abbrev( ts as mpz_t ) as CFStringRef
|
|
CFStringRef result
|
|
result = fn mpz_cf2( ts ) // FB wrapper mpz_t to CFString
|
|
int nLen
|
|
nLen = len( result )
|
|
if nLen <= 53 then RETURN result
|
|
result = fn AbbrevString( result, 25, 25 )
|
|
end fn = result
|
|
|
|
// -----------------------------------------------------------
|
|
// Helper: Is n probably prime?
|
|
// -----------------------------------------------------------
|
|
local fn IsProbablePrime( n as mpz_t ) as boolean
|
|
int r
|
|
r = fn mpz_probab_prime_p( n, 25 ) // 25 reps should be sufficient
|
|
// r = 0 (composite), 1 (probable prime), 2 (definite prime)
|
|
if r > 0 then RETURN _true
|
|
end fn = _false
|
|
|
|
// -----------------------------------------------------------
|
|
// Helper: Simple 32-bit primality test for Stage 1 ECM
|
|
// -----------------------------------------------------------
|
|
local fn IsSmallPrime32( p as UInt32 ) as boolean
|
|
if p < 2 then return _false
|
|
if p = 2 then return _true
|
|
if (p AND 1) == 0 then return _false
|
|
|
|
UInt32 d
|
|
d = 3
|
|
while d * d <= p
|
|
if (p MOD d) == 0 then return _false
|
|
d = d + 2
|
|
wend
|
|
|
|
end fn = _true
|
|
|
|
// Clear global mpz
|
|
void local fn clear_MPZ
|
|
mpz_clear( num )
|
|
mpz_clear( den )
|
|
mpz_clear( inv )
|
|
mpz_clear( t1 )
|
|
mpz_clear( t2 )
|
|
mpz_clear( g )
|
|
end fn
|
|
|
|
// Init global mpz
|
|
void local fn Init_MPZ
|
|
mpz_init( num )
|
|
mpz_init( den )
|
|
mpz_init( inv )
|
|
mpz_init( t1 )
|
|
mpz_init( t2 )
|
|
mpz_init( g )
|
|
end fn
|
|
|
|
// ---------------------------------------------------------------
|
|
// ECM_PointDouble
|
|
// Point doubling: P = (x,y) → 2P on y^2 = x^3 + a x + b (MOD n)
|
|
// Returns _true if a nontrivial factor found (in 'factor'),
|
|
// else _false.
|
|
// On success, 'factor' holds GCD > 1, < n, and x,y are undefined
|
|
// On normal completion, x,y are updated to 2P.
|
|
// ---------------------------------------------------------------
|
|
|
|
local fn ECM_PointDouble( n as mpz_t, a as mpz_t, ¬
|
|
x as mpz_t, y as mpz_t, factor as mpz_t ) as boolean
|
|
|
|
fn Init_MPZ
|
|
|
|
mpz_mul( t1, x, x ) // num = 3 * x^2 + a (MOD n)
|
|
mpz_mod( t1, t1, n ) // t1 = x^2
|
|
mpz_mul_ui( num, t1, 3 ) // num = 3*x^2
|
|
mpz_add( num, num, a )
|
|
mpz_mod( num, num, n )
|
|
|
|
mpz_mul_ui( den, y, 2 ) // den = 2*y (MOD n)
|
|
mpz_mod( den, den, n )
|
|
|
|
if fn mpz_invert( inv, den, n ) == 0 // try to invert den
|
|
// inversion failed -> GCD(den, n) may give a factor?
|
|
// do you feel lucky, punk?
|
|
mpz_GCD( g, den, n )
|
|
if (fn mpz_cmp_ui( g, 1 ) > 0) && (fn mpz_cmp( g, n ) < 0)
|
|
mpz_set( factor, g )
|
|
fn Clear_MPZ
|
|
return _true
|
|
end if
|
|
// GCD == 1 or GCD == n -> useless curve (sigh)
|
|
fn Clear_MPZ
|
|
return _false
|
|
end if
|
|
|
|
mpz_mul( t1, num, inv ) // lambda = num * inv (MOD n)
|
|
mpz_mod( t1, t1, n ) // t1 = lambda
|
|
|
|
// x3 = lambda^2 - 2 * x (MOD n)
|
|
mpz_mul( t2, t1, t1 ) // lambda^2
|
|
mpz_mod( t2, t2, n )
|
|
mpz_mul_ui( num, x, 2 ) // 2x
|
|
mpz_mod( num, num, n )
|
|
mpz_sub( t2, t2, num )
|
|
mpz_mod( t2, t2, n ) // t2 = x3
|
|
|
|
// y3 = lambda * (x - x3) - y (MOD n)
|
|
mpz_sub( num, x, t2 ) // x - x3
|
|
mpz_mod( num, num, n )
|
|
mpz_mul( num, t1, num ) // lambda * (x - x3)
|
|
mpz_mod( num, num, n )
|
|
mpz_sub( num, num, y ) // lambda * (x - x3) - y
|
|
mpz_mod( num, num, n ) // num = y3
|
|
|
|
mpz_set( x, t2 ) // update P ← (x3, y3)
|
|
mpz_set( y, num )
|
|
|
|
fn Clear_MPZ
|
|
end fn = _false
|
|
|
|
// -----------------------------------------------------------
|
|
// ECM_PointAdd
|
|
// Point addition: P1 = (x1,y1) ← P1 + P2, with P2 = (x2,y2)
|
|
// Curve: y^2 = x^3 + a x + b (MOD n)
|
|
// Returns _true if a factor found (in 'factor'), else _false.
|
|
// -----------------------------------------------------------
|
|
local fn ECM_PointAdd( n as mpz_t, a as mpz_t, x1 as mpz_t, ¬
|
|
y1 as mpz_t, x2 as mpz_t, y2 as mpz_t, factor as mpz_t ) as boolean
|
|
|
|
fn Init_MPZ
|
|
|
|
mpz_sub( t1, n, y2 ) // t1 = -y2 MOD n
|
|
mpz_mod( t1, t1, n )
|
|
// If x1 == x2 and y1 == -y2, point is infinity
|
|
// For our purposes, treat this as a useless -> just return _false.
|
|
if (fn mpz_cmp( x1, x2 ) == 0) && (fn mpz_cmp( y1, t1 )== 0)
|
|
fn Clear_MPZ
|
|
return _false
|
|
end if
|
|
|
|
mpz_sub( num, y2, y1 ) // num = y2 - y1 (MOD n)
|
|
mpz_mod( num, num, n )
|
|
|
|
mpz_sub( den, x2, x1 ) // den = x2 - x1 (MOD n)
|
|
mpz_mod( den, den, n )
|
|
|
|
// try to invert den
|
|
if fn mpz_invert( inv, den, n ) == 0
|
|
mpz_GCD( g, den, n )
|
|
if (fn mpz_cmp_ui( g, 1 ) > 0) && (fn mpz_cmp( g, n ) < 0)
|
|
mpz_set( factor, g )
|
|
fn Clear_MPZ
|
|
return _true
|
|
end if
|
|
fn Clear_MPZ
|
|
return _false
|
|
end if
|
|
|
|
mpz_mul( t1, num, inv ) // lambda = num * inv (MOD n)
|
|
mpz_mod( t1, t1, n ) // t1 = lambda
|
|
|
|
// x3 = lambda^2 - x1 - x2 (MOD n)
|
|
mpz_mul( t2, t1, t1 ) // lambda^2
|
|
mpz_mod( t2, t2, n )
|
|
mpz_sub( t2, t2, x1 )
|
|
mpz_sub( t2, t2, x2 )
|
|
mpz_mod( t2, t2, n ) // t2 = x3
|
|
|
|
// y3 = lambda * (x1 - x3) - y1 (MOD n)
|
|
mpz_sub( num, x1, t2 ) // x1 - x3
|
|
mpz_mod( num, num, n )
|
|
mpz_mul( num, t1, num )
|
|
mpz_mod( num, num, n )
|
|
mpz_sub( num, num, y1 )
|
|
mpz_mod( num, num, n ) // num = y3
|
|
|
|
mpz_set( x1, t2 )
|
|
mpz_set( y1, num )
|
|
|
|
fn Clear_MPZ
|
|
|
|
end fn = _false
|
|
|
|
// ------------------------------------------------------------
|
|
// ECM_ScalarMul
|
|
// Scalar multiply P = (x,y) by k: P ← k * P
|
|
// Uses binary double-and-add. Returns _true if a factor found.
|
|
// ------------------------------------------------------------
|
|
|
|
local fn ECM_ScalarMul( n as mpz_t, a as mpz_t, x as mpz_t, ¬
|
|
y as mpz_t, k as UInt64, factor as mpz_t ) as boolean
|
|
|
|
if k <= 1 then return _false
|
|
|
|
mpz_t rx, ry, qx, qy
|
|
mpz_init( rx )
|
|
mpz_init( ry )
|
|
mpz_init( qx )
|
|
mpz_init( qy )
|
|
|
|
int hasR
|
|
hasR = _false // R (result) is "not set" yet
|
|
|
|
mpz_set( qx, x ) // Q = P (base)
|
|
mpz_set( qy, y )
|
|
|
|
// find highest bit of k
|
|
int sbit
|
|
sbit = 63
|
|
while (sbit > 0) && ((k AND (1ULL << sbit)) = 0)
|
|
sbit -- // decrement bit
|
|
wend
|
|
|
|
// process bits from MSB down to 0
|
|
while sbit >= 0
|
|
if hasR
|
|
// R = 2R
|
|
if fn ECM_PointDouble( n, a, rx, ry, factor )
|
|
// factor found
|
|
mpz_clear( rx )
|
|
mpz_clear( ry )
|
|
mpz_clear( qx )
|
|
mpz_clear( qy )
|
|
return _true
|
|
end if
|
|
end if
|
|
|
|
if (k AND (1ULL << sbit)) != 0
|
|
if hasR = _false
|
|
// first 1 bit: R = Q
|
|
mpz_set( rx, qx )
|
|
mpz_set( ry, qy )
|
|
hasR = _true
|
|
else
|
|
// R = R + Q
|
|
if fn ECM_PointAdd( n, a, rx, ry, qx, qy, factor )
|
|
mpz_clear( rx )
|
|
mpz_clear( ry )
|
|
mpz_clear( qx )
|
|
mpz_clear( qy )
|
|
return _true
|
|
end if
|
|
end if
|
|
end if
|
|
|
|
sbit -- // decrement bit
|
|
wend
|
|
|
|
if hasR
|
|
mpz_set( x, rx )
|
|
mpz_set( y, ry )
|
|
end if
|
|
|
|
mpz_clear( rx )
|
|
mpz_clear( ry )
|
|
mpz_clear( qx )
|
|
mpz_clear( qy )
|
|
end fn = _false
|
|
|
|
|
|
// ----------------------------------------------------------------
|
|
// ECM_Stage1Curve
|
|
// Stage-1 ECM on one curve determined by 'sigma'
|
|
// Curve: y^2 = x^3 + a x + b (MOD n) with (x,y) chosen from sigma.
|
|
// B1: first-stage bound, e.g. 10000.
|
|
// Returns _true and sets 'factor' if a factor was
|
|
// found, else _false.
|
|
// ----------------------------------------------------------------
|
|
|
|
local fn ECM_Stage1Curve( n as mpz_t, factor as mpz_t, ¬
|
|
B1 as UInt32, sigma as UInt32 ) as boolean
|
|
|
|
mpz_t a, b, x, y, tmp
|
|
mpz_init( a )
|
|
mpz_init( b )
|
|
mpz_init( x )
|
|
mpz_init( y )
|
|
mpz_init( tmp )
|
|
|
|
mpz_set_ui( x, sigma ) // x = sigma (MOD n)
|
|
mpz_mod( x, x, n )
|
|
|
|
mpz_set_ui( y, sigma + 1 ) // y = sigma + 1 (MOD n)
|
|
mpz_mod( y, y, n )
|
|
|
|
mpz_set_ui( a, sigma + 2 ) // a = sigma + 2 (MOD n)
|
|
mpz_mod( a, a, n )
|
|
|
|
// b = y^2 - x^3 - a*x (MOD n)
|
|
mpz_mul( b, y, y ) // y^2
|
|
mpz_mod( b, b, n )
|
|
|
|
mpz_mul( tmp, x, x ) // x^2
|
|
mpz_mod( tmp, tmp, n )
|
|
mpz_mul( tmp, tmp, x ) // x^3
|
|
mpz_mod( tmp, tmp, n )
|
|
mpz_sub( b, b, tmp ) // y^2 - x^3
|
|
|
|
mpz_mul( tmp, a, x ) // a * x
|
|
mpz_mod( tmp, tmp, n )
|
|
mpz_sub( b, b, tmp ) // y^2 - x^3 - a*x
|
|
mpz_mod( b, b, n )
|
|
|
|
// Multiply P by prime powers p^e with p <= B1
|
|
UInt32 p, pe
|
|
p = 2
|
|
|
|
while p <= B1
|
|
if fn IsSmallPrime32( p )
|
|
// pe = p^e <= B1
|
|
pe = p
|
|
while (pe * p) <= B1
|
|
pe = pe * p
|
|
wend
|
|
|
|
if fn ECM_ScalarMul( n, a, x, y, pe, factor )
|
|
// factor found
|
|
mpz_clear( a )
|
|
mpz_clear( b )
|
|
mpz_clear( x )
|
|
mpz_clear( y )
|
|
mpz_clear( tmp )
|
|
return _true
|
|
end if
|
|
end if
|
|
|
|
if p == 2 then p = 3 else p = p + 2
|
|
wend
|
|
|
|
mpz_clear( a )
|
|
mpz_clear( b )
|
|
mpz_clear( x )
|
|
mpz_clear( y )
|
|
mpz_clear( tmp )
|
|
end fn = _false
|
|
|
|
// -----------------------------------------------------------
|
|
// ECM_FindFactor
|
|
// Try several Stage-1 ECM curves with increasing sigma.
|
|
// Returns _true and sets 'factor' if a factor is found.
|
|
// -----------------------------------------------------------
|
|
|
|
local fn ECM_FindFactor( n as mpz_t, factor as mpz_t ) as boolean
|
|
UInt32 B1, sigma
|
|
int curves
|
|
int i
|
|
B1 = 10000 // bound (bigger = slower)
|
|
curves = 20 // number of different curves to try
|
|
sigma = 2
|
|
i = 0
|
|
while i < curves
|
|
if fn ECM_Stage1Curve( n, factor, B1, sigma )
|
|
return _true
|
|
end if
|
|
sigma ++
|
|
i ++
|
|
wend
|
|
end fn = _false
|
|
|
|
// Safely divide mpz_t
|
|
|
|
local fn SafeMpzDiv( z as mpz_t, x as mpz_t, y as mpz_t ) as boolean
|
|
if ( fn mpz_sgn( y ) == 0 ) // div by zero?
|
|
mpz_set_ui( z, 0)
|
|
return _false
|
|
end if
|
|
// Call GMP division (truncate toward zero)
|
|
mpz_tdiv_q( z, x, y )
|
|
end fn = _true
|
|
|
|
// -----------------------------------------------------------
|
|
// Recursive factorization using trial division + ECM Stage 1
|
|
// Results go into global gFactors() / gFactorCount
|
|
// -----------------------------------------------------------
|
|
void local fn FactorRec( n as mpz_t )
|
|
// if n == 1 → nothing
|
|
if fn mpz_cmp_ui( n, 1 ) = 0 then RETURN
|
|
// if prime, store it
|
|
if fn IsProbablePrime( n )
|
|
mpz_set( gFactors(gFactorCount), n )
|
|
gFactorCount = gFactorCount + 1
|
|
RETURN
|
|
end if
|
|
|
|
mpz_t d, q
|
|
mpz_init( d )
|
|
mpz_init( q )
|
|
|
|
Boolean got = _False
|
|
|
|
got = fn ECM_FindFactor( n, d )
|
|
|
|
if got == _false
|
|
// ECM did not find a nontrivial factor; treat n as a "final" factor
|
|
mpz_set( gFactors(gFactorCount), n )
|
|
gFactorCount = gFactorCount + 1
|
|
mpz_clear( d )
|
|
mpz_clear( q )
|
|
RETURN
|
|
end if
|
|
|
|
if fn SafeMpzDiv( q, n, d) // q = n / d
|
|
// recurse on d and q
|
|
fn FactorRec( d )
|
|
fn FactorRec( q )
|
|
end if
|
|
|
|
mpz_clear( d )
|
|
mpz_clear( q )
|
|
end fn
|
|
|
|
|
|
// -----------------------------------------------------------
|
|
// Init Fermat numbers F0..F9
|
|
// F_n = 2^(2^n) + 1
|
|
// -----------------------------------------------------------
|
|
local fn InitFermat
|
|
int i, expPow
|
|
// init working bigint
|
|
mpz_init( fm )
|
|
// init Fermat array
|
|
i = 0
|
|
while i <= 15
|
|
mpz_init( gF(i) )
|
|
i = i + 1
|
|
wend
|
|
// init factor buffer mpz_t
|
|
i = 0
|
|
while i <= 30
|
|
mpz_init( gFactors(i) )
|
|
i = i + 1
|
|
wend
|
|
// compute F0..F15 once
|
|
i = 0
|
|
while i <= 15
|
|
expPow = 1
|
|
// expPow = 2^i
|
|
int k
|
|
k = 0
|
|
while k < i
|
|
expPow = expPow * 2
|
|
k = k + 1
|
|
wend
|
|
// fm = 2^(2^i) + 1
|
|
mpz_ui_pow_ui( fm, 2, expPow )
|
|
mpz_add_ui( fm, fm, 1 )
|
|
mpz_set( gF(i), fm )
|
|
i ++
|
|
wend
|
|
gFN = 0
|
|
gDidFactorHeader = 0
|
|
end fn
|
|
|
|
// -----------------------------------------------------------
|
|
// Print first 10 Fermat numbers
|
|
// -----------------------------------------------------------
|
|
local fn Main
|
|
CFStringRef lineStr, fStr
|
|
int i
|
|
print @"\nThe first 10 Fermat numbers are:"
|
|
tim = fn CACurrentMediaTime
|
|
|
|
i = 0
|
|
while i <= 9
|
|
fStr = fn Abbrev( gF(i) )
|
|
lineStr = @""
|
|
lineStr = concat( lineStr, fn FermatLabel( i ) )
|
|
lineStr = concat( lineStr, @" = " )
|
|
lineStr = concat( lineStr, fStr )
|
|
print lineStr;
|
|
i ++
|
|
wend
|
|
|
|
printf @"\nElapsed time: %.3f ms", (fn CACurrentMediaTime - tim) * 1000
|
|
print @"\nClick [Next] to display factors of Fermat numbers.\n"
|
|
end fn
|
|
|
|
// -----------------------------------------------------------
|
|
// Show next Fermat factors (one F_n per click)
|
|
// Calculation time displayed if greater than 1 sec
|
|
// -----------------------------------------------------------
|
|
void local fn Factor_More
|
|
long td
|
|
CFStringRef lineStr, factStr
|
|
|
|
if gDidFactorHeader == 0
|
|
print @"Factors of Fermat numbers + product if multiple factors found:"
|
|
gDidFactorHeader = 1
|
|
end if
|
|
|
|
tim = fn CACurrentMediaTime
|
|
|
|
// factor current Fermat number gF(gFN)
|
|
gFactorCount = 0
|
|
|
|
mpz_t nCopy
|
|
mpz_init( nCopy )
|
|
mpz_set( nCopy, gF(gFN) )
|
|
fn FactorRec( nCopy )
|
|
mpz_clear( nCopy )
|
|
// Fk = [f1, f2, ...]
|
|
int j
|
|
lineStr = @""
|
|
lineStr = concat( lineStr, fn FermatLabel( gFN ) )
|
|
lineStr = concat( lineStr, @" = [" )
|
|
j = 0
|
|
while j < gFactorCount
|
|
if j > 0
|
|
lineStr = concat( lineStr, @", " )
|
|
end if
|
|
factStr = fn Abbrev( gFactors(j) )
|
|
lineStr = concat( lineStr, factStr )
|
|
j = j + 1
|
|
wend
|
|
|
|
lineStr = concat( lineStr, @"]" )
|
|
print lineStr;
|
|
// if >= 2 factors, second line = product line
|
|
if gFactorCount >= 2
|
|
lineStr = @""
|
|
lineStr = concat( lineStr, fn FermatLabel( gFN ) )
|
|
lineStr = concat( lineStr, @" = " )
|
|
CFStringRef fAbbrev
|
|
fAbbrev = fn Abbrev( gF(gFN) )
|
|
lineStr = @"\n "
|
|
lineStr = concat( lineStr, fAbbrev )
|
|
print lineStr;
|
|
end if
|
|
|
|
td = fn CACurrentMediaTime - tim
|
|
if td > 1
|
|
lineStr = fn StringWithFormat (@" (%.3f secs)", (fn CACurrentMediaTime - tim))
|
|
print @lineStr;
|
|
end if
|
|
gFN = gFN + 1
|
|
// Rosetta task sets a limit to under 5 minutes
|
|
// so kill the next button
|
|
if gFN > 14 then button 2, NO
|
|
end fn
|
|
|
|
// ------------------------------------------------------------
|
|
// Dialog handler
|
|
// ------------------------------------------------------------
|
|
void local fn DoDialog( ev as long, tag as long, wnd as long, obj as CFTypeRef )
|
|
select ( ev )
|
|
case _btnClick
|
|
select ( tag )
|
|
case 2
|
|
fn Factor_More
|
|
end select
|
|
|
|
case _windowShouldClose
|
|
// clear any uncleared mpz
|
|
int i
|
|
mpz_clear( fm )
|
|
i = 0
|
|
while i <= 15
|
|
mpz_clear( gF(i) )
|
|
i ++
|
|
wend
|
|
i = 0
|
|
while i <= 30
|
|
mpz_clear( gFactors(i) )
|
|
i ++
|
|
wend
|
|
end
|
|
end select
|
|
end fn
|
|
|
|
on dialog fn DoDialog
|
|
|
|
// ------------------------------------------------------------
|
|
// Main
|
|
// ------------------------------------------------------------
|
|
window 1, @"Fermat numbers", (0,0,700,540)
|
|
button 2, YES,, @"Next", (580,20,100,20),,, 1
|
|
|
|
fn InitFermat
|
|
fn Main
|
|
|
|
HandleEvents
|