[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 ====== <> Mathematics | Arts and crafts of Tcl-Tk programming