rwm -- I thought the primes page was getting too crowded so I added this as a new page. My intent on adding this was to help teach the math and provide a simple Tcl implementation of the algorithm. (It might belong in a cookbook of useful recipes??)
Warning - This is a probabalistic test and there are several know "lucas pseudo-primes" (Numbers that are not flagged as composite, but are in fact composite.)
References:
Dependencies:
binaryModExpo - a modular exponentiation function (I like the one at the bottom of RSA in Tcl)
set debug(primeCheckLucas) 0
set modes(fips186) 0
proc primeCheckLucas {N} {
variable modes
variable debug
## implemented by rwm 9/24/12 from ANSI X9.31
# find the first D in the sequence {5 -7 9 -11 13 -15}
# for which the Jacobi symbol (D/N)=-1.
# Set P=1 and Q=(1-D)/4
# If U(n-1) = 0 mod N then N is probably prime.
# To calculate Uk where K= N+1 from the lucas sequence perform the steps
#
if {$debug(primeCheckLucas)} {
puts stderr "${::STATUS}: primeCheckLucas [hex $N]"
}
for {set i 0} {1} {incr i} {
set D [expr {$i&1? (-2*$i-5) : (5+2*$i)}]
if {$debug(primeCheckLucas)} {
puts stderr "${::STATUS}: primeCheckLucas testing D=$D"
}
if {[jacobi $D $N] == 0} {
if {$modes(fips186)} {
if {$debug(primeCheckLucas)} {
puts stderr "${::STATUS}: number is composite since Jacobi($D,N)==0"
}
return 0 ;# for any D id the Jacobi is 0 N is composite
}
}
if {[jacobi $D $N] == -1} break
}
if {$debug(primeCheckLucas)} {
puts stderr "${::STATUS}: primeCheckLucas D=$D ([hex $D]) satisfies jacobi = -1"
}
## Now we have found a D that satisfies jacobi(D,N)=-1
eval {;# steps
set P 1
set Q [expr {(1-$D)/4}]
## if (UK % N) == 0 N is probably prime
## UK is U @ K=N+1
## or U(n+1) which I'll call Uk
#...
## following steps
eval {;#1
set D [expr {($P*$P)-4*$Q}]
}
eval {;#2 ??
set K [expr {$N+1}]
# and an associated binary expansion over r {$K&(1<<$r)} or {($K>>$r)&1}
set Kbinary [bin $K]
set Klength [string length $Kbinary]
set r [expr {$Klength-1}]
if {$debug(primeCheckLucas)} {
puts stderr "${::STATUS}: primeCheckLucas K expands as: $Kbinary"
puts stderr "${::STATUS}: primeCheckLucas K is $Klength bits long"
puts stderr "${::STATUS}: primeCheckLucas r=$r"
}
}
eval {;#3
set U 1
set V $P
}
set p $N ;#<-------------- RWM GUESS
eval {;#4
for {set i [expr {$r-1}]} {$i >= 0} {incr i -1} {
set nextU [expr {($U*$V)%$p}] ;# <---------- p ???
set nextV [modmath::div2 [expr {($V*$V)+($D*$U*$U)}] $p];# <---------- p ???
set U $nextU
set V $nextV
if {$K&(1<<$i)} {
set nextU [modmath::div2 [expr {$P*$U+$V}] $p]
set nextV [modmath::div2 [expr {$P*$V+$D*$U}] $p]
if {$modes(fips186)} {
## -- fips 186 is different
set nextU [modmath::div2 [expr {$U+$V}] $p]
set nextV [modmath::div2 [expr {$V+$D*$U}] $p]
}
set U $nextU
set V $nextV
}
if {$debug(primeCheckLucas)} {
puts stderr "${::STATUS}: primeCheckLucas i=$i (U,V)=($U,$V) Ki=[expr {($K&(1<<$i))?1:0}]"
}
}
}
eval {;#5
set Uk $U
set Vk $V
}
}
if {$modes(fips186)} {# fips 186 is different
if {$U == 0} {
return 1;# probably prime
} else {
return 0;# known composite
}
}
# so do final check on Uk
if {(($U*$K) % $N) == 0} {
return 1;# probably prime
} else {
return 0;# known composite
}
}Here are the 1st 20 pseudo-primes I get by comparing the above code to small-prime trial division and 27 roundes of MR testing starting at N=3 and stepping by 2:
5 11 323 377 1159 1829 3827 5459 5777 9071 9179 10877 11419 11663 13919 14839 16109 16211 18407 18971
RKzn 2017-04-13 I have not checked any code, but I note that "5" and "11" from the above list are actual prime numbers, not 'pseudo-primes'. The remaining are not prime numbers, as they may be factored in two, or three, largish prime factors (using maxima):
% factor ([323, 377, 1159, 1829, 3827, 5459, 5777, 9071, 9179,
10877, 11419, 11663, 13919, 14839, 16109, 16211, 18407, 18971]);
% [17 19, 13 29, 19 61, 31 59, 43 89, 53 103, 53 109, 47 193, 67 137,
73 149, 19 601, 107 109, 31 449, 11 19 71, 89 181, 13 29 43, 79 233, 61 311]