[rwm] -- I needed this when implementing a [primeCheckLucas] function. 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??) References: 1. http://en.wikipedia.org/wiki/Jacobi_symbol 2. http://2000clicks.com/mathhelp/NumberTh27JacobiSymbolAlgorithm.aspx Dependencies: 1. binaryModExpo - a modular exponentiation function (I like the one at the bottom of [RSA in Tcl]) ====== set debug(jacobi) 0 proc jacobi {a b} { variable debug ## coded by rmelton 9/25/12 # ref: http://2000clicks.com/mathhelp/NumberTh27JacobiSymbolAlgorithm.aspx if {\$debug(jacobi)} { puts stderr "\${::STATUS}: jacobi \$a \$b" } if {\$b<=0 || (\$b&1)==0} { if {\$debug(jacobi)} { puts stderr "\${::STATUS}: jacobi=0 because b is negative, or even" } return 0; } set j 1 if {\$a<0} { set a [expr {0-\$a}] if {(\$b & 3) == 3} { if {\$debug(jacobi)} { puts stderr "\${::STATUS}: use identity: Jacobi(-a,b) = -Jacobi(a,b) if b≡3 (mod 4)" } set j [expr {0-\$j}] } else { if {\$debug(jacobi)} { puts stderr "\${::STATUS}: use identity: Jacobi(-a,b) = Jacobi(a,b) if b≡1 (mod 4)" } } } while {\$a != 0} { while {(\$a&1) == 0} { ##/* Process factors of 2: Jacobi(2,b)=-1 if b=3,5 (mod 8) */ set a [expr {\$a>>1}] if {((\$b & 7)==3) || ((\$b & 7)==5)} { set j [expr {0-\$j}] if {\$debug(jacobi)} { puts stderr "\${::STATUS}: removed factors of 2 by identity: Jacobi(2a,b) = -Jacobi(a,b) if b≡3 (mod 8)" } } else { if {\$debug(jacobi)} { puts stderr "\${::STATUS}: removed factors of 2 by identity: Jacobi(2a,b) = Jacobi(a,b) if b≡1 (mod 8)" } } } ##/* Quadratic reciprocity: Jacobi(a,b)=-Jacobi(b,a) if a=3,b=3 (mod 4) */ lassign [list \$a \$b] b a if {((\$a & 3)==3) && ((\$b & 3)==3)} { set j [expr {0-\$j}] if {\$debug(jacobi)} { puts stderr "\${::STATUS}: use Identity Jacobi(a,b)=-Jacobi(b,a) if a=3,b=3 (mod 4)" } } else { if {\$debug(jacobi)} { puts stderr "\${::STATUS}: use Identity Jacobi(a,b)=Jacobi(b,a) if a≡1 or b≡1 (mod 4)" } } if {\$debug(jacobi)} { puts stderr "\${::STATUS}: use Identity Jacobi(a Mod b,b) = Jacobi(a,b)" } set a [expr {\$a % \$b}] if {\$debug(jacobi)} { puts stderr "\${::STATUS}: Jacobi has been reduced to-> \$j*Jacobi(\$a,\$b)" } } if {\$b==1} { if {\$debug(jacobi)} { puts stderr "\${::STATUS}: jacobi=\$j (because we reduced to \$j*Jacobi(0,1)=1)" } return \$j } else { if {\$debug(jacobi)} { puts stderr "\${::STATUS}: jacobi=0 (because we reduced to \$j*Jacobi(0,1)=0 n!=1)" } return 0 } } ====== <> Mathematics | Arts and crafts of Tcl-Tk programming