## Version 2 of primeCheckLucas

Updated 2017-04-13 00:22:00 by RKzn

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:

1. fips183-6 - http://csrc.nist.gov/publications/fips/fips186-3/fips_186-3.pdf
2. sci.crypt discussion of lucas testing: https://groups.google.com/forum/?fromgroups#!msg/sci.crypt/uGko8m1lGQM/IhQTSXw6EzAJ

Dependencies:

1. modmath::div2 - in modular math if you need to divide an odd number by 2 you can add the modulus first to make the divisor even (ie A/2 mod M can be rewritten as (A+M)/2 mod M when needed)
2. jacobi - a routine to calculate the jacobi symbol

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]```

 Category Mathematics Arts and Crafts of Tcl-Tk Programming