149 lines
3.8 KiB
Text
149 lines
3.8 KiB
Text
// -----------------------------------------------
|
|
// Wagstaff primes
|
|
// Using FutureBasic 7.0.36
|
|
//
|
|
// November 2025, R.W
|
|
//
|
|
// -----------------------------------------------
|
|
|
|
include "gmp.incl"
|
|
|
|
_kPrimeReps = 25
|
|
_kNeedFirst = 10
|
|
_kNeedMore = 14
|
|
|
|
// mpz_sizeinbase sometimes would return one digit too many especially
|
|
// when the value is close to the edge of power of 10. This helper
|
|
// accounts for that and adjusts.
|
|
|
|
local fn mpz_digits10( z as mpz_t ) as int
|
|
int d
|
|
mpz_t pow10
|
|
d = fn mpz_sizeinbase( z, 10 )
|
|
if d <= 1 then RETURN d
|
|
mpz_init( pow10 )
|
|
// pow10 = 10^(d-1)
|
|
fn mpz_ui_pow_ui( pow10, 10, d - 1 )
|
|
// If z < 10^(d-1) then mpz_sizeinbase returned one too many
|
|
if fn mpz_cmp( z, pow10 ) < 0 then d = d - 1 //if so -> decrement
|
|
mpz_clear( pow10 )
|
|
end fn = d
|
|
|
|
local fn IsSmallPrime( n as UInt32 ) as boolean
|
|
if n < 2 then RETURN _false
|
|
if n = 2 then RETURN _true
|
|
if (n AND 1) == 0 then RETURN _false
|
|
mpz_t t
|
|
mpz_init( t )
|
|
fn mpz_set_ui( t, n )
|
|
int r : r = fn mpz_probab_prime_p( t, _kPrimeReps )
|
|
mpz_clear( t )
|
|
if r > 0 then RETURN _true
|
|
end fn = _false
|
|
|
|
// Find the next upstream prime
|
|
|
|
local fn NextPrimeAfter( startN as UInt32 ) as UInt32
|
|
UInt32 p
|
|
if startN < 2 then RETURN 2
|
|
p = startN + 1
|
|
if (p AND 1) = 0 then p = p + 1
|
|
while _true
|
|
if fn IsSmallPrime( p ) then RETURN p
|
|
p = p + 2
|
|
wend
|
|
end fn = 0
|
|
|
|
// Abbreviate long strings
|
|
|
|
local fn AbbrevLongStr( longStr as CFStringRef ) as CFStringRef
|
|
CFStringRef outStr, leftStr,rightStr
|
|
CFIndex strlen
|
|
outStr = @""
|
|
strlen = len(longStr)
|
|
if strlen < 60 then return longStr
|
|
leftStr = left(longStr, 27 )
|
|
rightStr = right(longStr, 28 )
|
|
outStr = concat( leftStr, @"...", rightStr)
|
|
end fn = outStr
|
|
|
|
// WagStaff_P = (2^p + 1) / 3 → outW
|
|
|
|
local fn WagstaffFromExponent( p as UInt32, outW as mpz_t ) as boolean
|
|
mpz_t N
|
|
mpz_init( N )
|
|
fn mpz_ui_pow_ui( N, 2, p )
|
|
fn mpz_add_ui( N, N, 1 )
|
|
unsigned long remd
|
|
remd = fn mpz_fdiv_q_ui( outW, N, 3 )
|
|
mpz_clear( N )
|
|
if remd != 0 then RETURN _false
|
|
end fn = _true
|
|
|
|
// Use GMP probable prime function with 25 reps
|
|
// sufficient to determine a large prime
|
|
|
|
local fn IsWagstaffPrime( p as UInt32, outW as mpz_t ) as boolean
|
|
if fn WagstaffFromExponent( p, outW ) = _false then RETURN _false
|
|
int res : res = fn mpz_probab_prime_p( outW, _kPrimeReps )
|
|
// 0 = composite, 1 = probable, 3 = positive
|
|
if res > 0 then RETURN _true
|
|
end fn = _false
|
|
|
|
// MAIN
|
|
//
|
|
local fn Main
|
|
UInt32 p
|
|
UInt32 foundFirst, foundMore
|
|
foundFirst = 0 : foundMore = 0
|
|
CFStringRef s2
|
|
|
|
CFTimeInterval t0, t1
|
|
t0 = fn CACurrentMediaTime
|
|
|
|
// Rosetta requirement for 10 Wagstaff primes
|
|
print @"First 10 Wagstaff primes:\n"
|
|
p = fn NextPrimeAfter( 1 )
|
|
while foundFirst < _kNeedFirst
|
|
mpz_t W
|
|
mpz_init( W )
|
|
if fn IsWagstaffPrime( p, W )
|
|
foundFirst = foundFirst + 1
|
|
// FB has a handy wrapper to obtain CFString from mpz_t output
|
|
// toolbox fn mpz_cf2( mpz_t op ) = CFStringRef // mpz_t → CFString (base 10)
|
|
CFStringRef s = fn mpz_cf2( W )
|
|
printf @"%4lu: %@", p,s
|
|
end if
|
|
mpz_clear( W )
|
|
p = fn NextPrimeAfter( p )
|
|
wend
|
|
t1 = fn CACurrentMediaTime
|
|
printf @"\nElapsed time: %.3f ms\n\n", (t1 - t0) * 1000.0
|
|
|
|
// Rosetta additional task
|
|
|
|
print @"Wagstaff next 14:\n"
|
|
|
|
CFTimeInterval t2 = fn CACurrentMediaTime
|
|
foundMore = 0
|
|
while foundMore < _kNeedMore
|
|
mpz_t W2
|
|
mpz_init( W2 )
|
|
if fn IsWagstaffPrime( p, W2 )
|
|
foundMore = foundMore + 1
|
|
int digits : digits = fn mpz_digits10( W2 )
|
|
s2 = fn mpz_cf2( W2 ) //mpz_cf2 is an FB wrapper
|
|
s2 = fn AbbrevLongStr(s2)
|
|
printf @"%4lu: %@ (%d digits)", p, s2, digits
|
|
end if
|
|
mpz_clear( W2 )
|
|
p = fn NextPrimeAfter( p )
|
|
wend
|
|
|
|
CFTimeInterval t3 = fn CACurrentMediaTime
|
|
printf @"\nElapsed time: %.3f secs\n", (t3 - t2)
|
|
end fn
|
|
|
|
Window 1,@"Wagstaff Primes",(0,0,640,540)
|
|
fn Main
|
|
HandleEvents
|