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:
- http://en.wikipedia.org/wiki/Jacobi_symbol
- http://2000clicks.com/mathhelp/NumberTh27JacobiSymbolAlgorithm.aspx
Dependencies:
- none ??
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
}
}
