primeCheckLucas

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:

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

Dependencies:

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)
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:

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