888 lines
23 KiB
Text
888 lines
23 KiB
Text
// ------------------------------------------------------------
|
||
// Home primes
|
||
//
|
||
// Using FutureBasic 7.0.37
|
||
// November 2025, R.W.
|
||
//
|
||
// This code is a lot longer than needed. Simply because I
|
||
// refused to use Pollard Rho (it's very slow and unreliable).
|
||
// So I used pseudo ECM (ECF) for factorization. Much faster.
|
||
//
|
||
// Output includes Rosetta Task requiring HP65 to be included
|
||
// it also has the "impossible" HP49 (just click [next]..)
|
||
// ------------------------------------------------------------
|
||
|
||
include "gmp.incl"
|
||
|
||
CFTimeInterval tim
|
||
|
||
long gN // the next starting integer for [Next] batches
|
||
// Global mpz_t bigints reused for speed
|
||
mpz_t num, den, inv, t1, t2, g
|
||
|
||
mpz_t pp, tmpN, tmpFactor
|
||
mpz_t fm
|
||
mpz_t gFactors(31)
|
||
int gFactorCount
|
||
|
||
|
||
// -----------------------------------------------------------
|
||
// 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 > 43 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 <= 43 then RETURN result
|
||
result = fn AbbrevString( result, 20, 20 )
|
||
end fn = result
|
||
|
||
// -----------------------------------------------------------
|
||
// Helper: Convert mpz_t to CFString
|
||
// -----------------------------------------------------------
|
||
local fn mpz_to_CFString( m as mpz_t ) as CFStringRef
|
||
CFStringRef result = @""
|
||
result = fn mpz_cf2( m ) //FB wrapper mpz_t to CFString base 10
|
||
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
|
||
|
||
// ------------------------------------------------------------
|
||
// Small prime test for Int input
|
||
// Uses mpz_probab_prime_p with gmp
|
||
// ------------------------------------------------------------
|
||
local fn IsPrimeInt( n as long ) as Boolean
|
||
if n < 2 then RETURN _false
|
||
fn mpz_set_si( pp, n )
|
||
int r
|
||
r = fn mpz_probab_prime_p( pp, 15 )
|
||
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 ECM Stage 1
|
||
// Results go into global gFactors() / gFactorCount
|
||
// trustProbPrime = _true → use IsProbablePrime shortcut
|
||
// trustProbPrime = _false → NEVER trust IsProbablePrime; rely on ECM
|
||
// -----------------------------------------------------------
|
||
void local fn FactorRec( n as mpz_t, trustProbPrime as boolean )
|
||
// if n == 1 → nothing
|
||
if fn mpz_cmp_ui( n, 1 ) = 0 then RETURN
|
||
|
||
// If allowed, trust mpz_probab_prime_p
|
||
if trustProbPrime
|
||
if fn IsProbablePrime( n )
|
||
mpz_set( gFactors( gFactorCount ), n )
|
||
gFactorCount = gFactorCount + 1
|
||
RETURN
|
||
end if
|
||
end if
|
||
|
||
mpz_t d, q
|
||
mpz_init( d )
|
||
mpz_init( q )
|
||
|
||
Boolean got
|
||
got = fn ECM_FindFactor( n, d )
|
||
|
||
if got == _false
|
||
// ECM did not find a nontrivial factor; treat n as 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 with the same mode
|
||
fn FactorRec( d, trustProbPrime )
|
||
fn FactorRec( q, trustProbPrime )
|
||
end if
|
||
|
||
mpz_clear( d )
|
||
mpz_clear( q )
|
||
end fn
|
||
|
||
// ------------------------------------------------------------
|
||
// FactorStringUnderscore
|
||
// Factor j into primes and build "p1_p2_p3" as CFStringRef
|
||
// forceFull = _false → use IsProbablePrime shortcut (fast)
|
||
// forceFull = _true → ignore IsProbablePrime, full ECM
|
||
// ------------------------------------------------------------
|
||
local fn FactorStringUnderscore( j as mpz_t, forceFull as boolean ) as CFStringRef
|
||
CFStringRef outStr = @""
|
||
|
||
// Work copy of j in fm, factor |j|
|
||
fn mpz_set( fm, j )
|
||
fn mpz_abs( fm, fm )
|
||
|
||
// If fm < 2, nothing to do
|
||
if fn mpz_cmp_ui( fm, 2 ) < 0 then RETURN outStr
|
||
|
||
// Reset factor buffer
|
||
gFactorCount = 0
|
||
|
||
// Fill gFactors[0..gFactorCount-1] with prime factors (not sorted)
|
||
// trustProbPrime = NOT forceFull
|
||
Boolean trustProbPrime
|
||
trustProbPrime = _true
|
||
if forceFull then trustProbPrime = _false
|
||
|
||
fn FactorRec( fm, trustProbPrime )
|
||
|
||
// If FactorRec somehow didn't produce anything, bail out
|
||
if gFactorCount <= 0 then RETURN outStr
|
||
|
||
// sort factors in ascending order (same as before)
|
||
int i, jidx
|
||
mpz_t tmp
|
||
mpz_init( tmp )
|
||
|
||
i = 0
|
||
while i < gFactorCount - 1
|
||
jidx = i + 1
|
||
while jidx < gFactorCount
|
||
if fn mpz_cmp( gFactors(i), gFactors(jidx) ) > 0
|
||
fn mpz_set( tmp, gFactors(i) )
|
||
fn mpz_set( gFactors(i), gFactors(jidx) )
|
||
fn mpz_set( gFactors(jidx), tmp )
|
||
end if
|
||
jidx++
|
||
wend
|
||
i++
|
||
wend
|
||
|
||
mpz_clear( tmp )
|
||
|
||
// build "p1_p2_p3_..." string
|
||
outStr = @""
|
||
i = 0
|
||
while i < gFactorCount
|
||
if i > 0
|
||
outStr = concat( outStr, @"_" )
|
||
end if
|
||
outStr = concat( outStr, fn mpz_to_CFString( gFactors(i) ) )
|
||
i++
|
||
wend
|
||
|
||
end fn = outStr
|
||
|
||
// ------------------------------------------------------------
|
||
// FactorConcat (ECM-based)
|
||
// Factor j into primes using FactorRec (ECM + primality test),
|
||
// collect factors in ascending order, and concatenate them
|
||
// ------------------------------------------------------------
|
||
local fn FactorConcat( j as mpz_t ) as CFStringRef
|
||
CFStringRef outStr = @""
|
||
|
||
// Work copy of j in fm, factor |j|
|
||
fn mpz_set( fm, j )
|
||
fn mpz_abs( fm, fm )
|
||
|
||
// If fm < 2, nothing to do
|
||
if fn mpz_cmp_ui( fm, 2 ) < 0 then RETURN outStr
|
||
|
||
// Reset factor buffer
|
||
gFactorCount = 0
|
||
|
||
// Fill gFactors[0..gFactorCount-1] with prime factors (not sorted)
|
||
fn FactorRec( fm, _true ) // for chain math, trust probable primes
|
||
|
||
// If FactorRec somehow didn't produce anything, bail out
|
||
if gFactorCount <= 0 then RETURN outStr
|
||
|
||
// sort factors in ascending order
|
||
int i, jidx
|
||
mpz_t tmp
|
||
mpz_init( tmp )
|
||
|
||
i = 0
|
||
while i < gFactorCount - 1
|
||
jidx = i + 1
|
||
while jidx < gFactorCount
|
||
if fn mpz_cmp( gFactors(i), gFactors(jidx) ) > 0
|
||
// swap gFactors(i) and gFactors(jidx)
|
||
fn mpz_set( tmp, gFactors(i) )
|
||
fn mpz_set( gFactors(i), gFactors(jidx) )
|
||
fn mpz_set( gFactors(jidx), tmp )
|
||
end if
|
||
jidx++
|
||
wend
|
||
i++
|
||
wend
|
||
|
||
mpz_clear( tmp )
|
||
|
||
// build concatenated decimal string of prime factors
|
||
outStr = @""
|
||
i = 0
|
||
while i < gFactorCount
|
||
outStr = concat( outStr, fn mpz_to_CFString( gFactors(i) ) )
|
||
i++
|
||
wend
|
||
|
||
end fn = outStr
|
||
|
||
|
||
// Max chain length we’re willing to handle
|
||
_kMaxChain = 32
|
||
|
||
// ------------------------------------------------------------
|
||
// Home_Prime_Chain
|
||
// Build Home Prime chain for integer n
|
||
// - Uses FactorConcat (ECM-backed) for factoring
|
||
// - Stops when value is prime OR when no progress is made
|
||
// ------------------------------------------------------------
|
||
local fn Home_Prime_Chain( n as long ) as CFStringRef
|
||
CFStringRef result = @""
|
||
CFStringRef pad = @" "
|
||
CFStringRef chain
|
||
CFStringRef nStr, iterStr, tmpCF
|
||
|
||
// Guard: n < 2 → empty
|
||
if n < 2 then RETURN result
|
||
|
||
// If n is prime: "HPn = n"
|
||
if fn IsPrimeInt( n )
|
||
nStr = mid( str( n ), 1 )
|
||
chain = concat( @"HP", nStr )
|
||
chain = concat( chain, left( pad, 10 - len( chain ) ), @"= " )
|
||
chain = concat( chain, nStr )
|
||
result = chain
|
||
RETURN result
|
||
end if
|
||
|
||
// ----------------------------------------
|
||
// Otherwise, build the chain using GMP
|
||
// h[0..iterCount] hold the successive values
|
||
// ----------------------------------------
|
||
_kMaxSteps = 64
|
||
mpz_t h( _kMaxSteps - 1 )
|
||
int i
|
||
|
||
// initialize mpz_t array
|
||
i = 0
|
||
while i < _kMaxSteps
|
||
fn mpz_init( h(i) )
|
||
i++
|
||
wend
|
||
|
||
int iterCount
|
||
iterCount = 0
|
||
|
||
CFStringRef concatStr
|
||
CFStringRef oneDigit
|
||
|
||
// Track previous value to detect "no progress" stalls
|
||
mpz_t prevN
|
||
fn mpz_init( prevN )
|
||
int stalled
|
||
stalled = _false
|
||
|
||
// h[0] = initial n
|
||
fn mpz_set_si( h(0), n )
|
||
fn mpz_set_si( tmpN, n ) // current value in the chain
|
||
|
||
// iterate until prime, stall, or max steps
|
||
while iterCount < _kMaxSteps
|
||
int isPrime
|
||
isPrime = fn mpz_probab_prime_p( tmpN, 15 )
|
||
if isPrime > 0
|
||
// Probable prime: normal termination
|
||
exit while
|
||
end if
|
||
|
||
// Save previous value so we can detect no progress
|
||
fn mpz_set( prevN, tmpN )
|
||
|
||
// Not prime: factor tmpN, build concatenation string of prime factors
|
||
concatStr = fn FactorConcat( tmpN ) // decimal string of concatenated factors
|
||
|
||
// Build new tmpN from decimal digits of concatStr
|
||
fn mpz_set_ui( tmpN, 0 )
|
||
|
||
int L, dpos, digit
|
||
L = len( concatStr )
|
||
dpos = 1
|
||
|
||
while dpos <= L
|
||
// CFString mid() is 0-based: start at dpos-1
|
||
oneDigit = mid( concatStr, dpos - 1, 1 )
|
||
digit = IntVal( oneDigit ) // CFStringRef → int digit (0–9)
|
||
fn mpz_mul_ui( tmpN, tmpN, 10 )
|
||
fn mpz_add_ui( tmpN, tmpN, digit )
|
||
dpos++
|
||
wend
|
||
|
||
// Check for "no progress": new tmpN == prevN ⇒ factoring failed / stalled
|
||
if fn mpz_cmp( tmpN, prevN ) == 0
|
||
stalled = _true
|
||
exit while
|
||
end if
|
||
|
||
// Store this term in the chain
|
||
iterCount++
|
||
if iterCount < _kMaxSteps
|
||
fn mpz_set( h(iterCount), tmpN )
|
||
end if
|
||
wend
|
||
|
||
// Safety clamp if we bailed by hitting the step cap
|
||
if iterCount >= _kMaxSteps
|
||
iterCount = _kMaxSteps - 1
|
||
stalled = _true
|
||
end if
|
||
|
||
// iterCount is the index of the final value in h[]
|
||
CFStringRef finalStr
|
||
finalStr = fn Abbrev( h(iterCount) )
|
||
|
||
// ----------------------------------------------
|
||
// Build the chain in multi-line readable format
|
||
// ----------------------------------------------
|
||
chain = @""
|
||
|
||
// First line: HPn(iterCount)
|
||
pad = @" "
|
||
nStr = mid( str( n ), 1 )
|
||
iterStr = mid( str( iterCount ), 1 )
|
||
|
||
chain = concat( chain, @"HP" )
|
||
chain = concat( chain, nStr )
|
||
if iterCount > 1
|
||
chain = concat( chain, @" (" )
|
||
chain = concat( chain, iterStr )
|
||
chain = concat( chain, @")" )
|
||
end if
|
||
chain = concat( chain, left( pad, 10 - len( chain ) ), @"= " )
|
||
|
||
// ----------------------------------------------------
|
||
// Print each step on its own line, as factor strings:
|
||
// stepIndex = 1 → factors of h(0) (original n)
|
||
// stepIndex = 2 → factors of h(1), etc.
|
||
// remaining = iterCount - stepIndex
|
||
// ----------------------------------------------------
|
||
int stepIndex
|
||
stepIndex = 1
|
||
while stepIndex <= iterCount - 1
|
||
int remaining
|
||
remaining = iterCount - stepIndex
|
||
|
||
// Factor h(stepIndex-1) and show "p1_p2_p3"
|
||
// For intermediate steps we can use the fast mode (trust probable primes).
|
||
tmpCF = fn FactorStringUnderscore( h(stepIndex - 1), _false )
|
||
iterStr = mid( str( remaining ), 1 )
|
||
|
||
if stepIndex == 1
|
||
pad = @""
|
||
else
|
||
pad = @" "
|
||
end if
|
||
|
||
chain = concat( chain, pad )
|
||
tmpCF = fn AbbrevString( tmpCF, 20, 20 )
|
||
chain = concat( chain, tmpCF )
|
||
chain = concat( chain, @" (" )
|
||
chain = concat( chain, iterStr )
|
||
chain = concat( chain, @")" )
|
||
chain = concat( chain, @"\n" )
|
||
|
||
stepIndex++
|
||
wend
|
||
|
||
// Final line: just show the final value (probable prime)
|
||
if stepIndex > 1
|
||
chain = concat( chain, @" " )
|
||
end if
|
||
|
||
//CFStringRef finalStr
|
||
finalStr = fn Abbrev( h(iterCount) )
|
||
chain = concat( chain, finalStr )
|
||
if stalled
|
||
chain = concat( chain, @" (stalled – factoring limit)" )
|
||
end if
|
||
|
||
result = chain
|
||
|
||
// clear h[] and prevN
|
||
i = 0
|
||
while i < _kMaxSteps
|
||
fn mpz_clear( h(i) )
|
||
i++
|
||
wend
|
||
fn mpz_clear( prevN )
|
||
|
||
end fn = result
|
||
|
||
// ------------------------------------------------------------
|
||
// More_HP
|
||
// show the next 20 Home Primes
|
||
// ------------------------------------------------------------
|
||
void local fn More_HP
|
||
CFStringRef pline
|
||
long n
|
||
int count
|
||
|
||
tim = fn CACurrentMediaTime
|
||
|
||
count = 0
|
||
n = gN
|
||
while count < 20
|
||
pline = fn Home_Prime_Chain( n )
|
||
print @pline
|
||
n++
|
||
count++
|
||
wend
|
||
|
||
gN = n
|
||
if gN > 160 then button 2, NO
|
||
printf @"\nElapsed time: %.3f secs\n", ( fn CACurrentMediaTime - tim )
|
||
end fn
|
||
|
||
// ------------------------------------------------------------
|
||
// Main — compute 2 to 20 Home primes
|
||
// ------------------------------------------------------------
|
||
local fn Main
|
||
CFStringRef pline
|
||
long k
|
||
// init globals
|
||
fn mpz_init( pp )
|
||
fn mpz_init( tmpN )
|
||
fn mpz_init( tmpFactor )
|
||
|
||
int i
|
||
i = 0
|
||
while i <= 30
|
||
mpz_init( gFactors(i) )
|
||
i++
|
||
wend
|
||
|
||
print @"Home Prime for integer 2 through 20:\n"
|
||
|
||
tim = fn CACurrentMediaTime
|
||
|
||
k = 2
|
||
while k <= 20
|
||
pline = fn Home_Prime_Chain( k )
|
||
print @pline
|
||
k++
|
||
wend
|
||
|
||
printf @"\nElapsed time: %.3f secs\n", ( fn CACurrentMediaTime - tim )
|
||
print @"Click [Next] for the next 20 Home primes.\n"
|
||
|
||
// Prepare next batch start
|
||
gN = 21
|
||
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 More_HP
|
||
end select
|
||
|
||
case _windowShouldClose
|
||
// clear any uncleared mpz
|
||
fn mpz_clear( pp )
|
||
fn mpz_clear( tmpN )
|
||
fn mpz_clear( tmpFactor )
|
||
end
|
||
end select
|
||
end fn
|
||
|
||
on dialog fn DoDialog
|
||
|
||
// ------------------------------------------------------------
|
||
// Main
|
||
// ------------------------------------------------------------
|
||
window 1, @"Home primes", (0,0,800,600)
|
||
button 2, YES,, @"Next", (680,20,100,20),,, 1
|
||
|
||
// init globals
|
||
fn mpz_init( pp )
|
||
fn mpz_init( tmpN )
|
||
fn mpz_init( tmpFactor )
|
||
|
||
fn Main
|
||
|
||
HandleEvents
|