[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                       
======


<<categories>> Mathematics | Arts and crafts of Tcl-Tk programming