Jacobi symbol
The Jacobi symbol is a multiplicative function that generalizes the Legendre symbol. Specifically, the Jacobi symbol (a | n) equals the product of the Legendre symbols (a | p_i)^(k_i), where n = p_1^(k_1)*p_2^(k_2)*...*p_i^(k_i) and the Legendre symbol (a | p) denotes the value of a ^ ((p-1)/2) (mod p)
- (a | p) ≡ 1 if a is a square (mod p)
- (a | p) ≡ -1 if a is not a square (mod p)
- (a | p) ≡ 0 if a ≡ 0
If n is prime, then the Jacobi symbol (a | n) equals the Legendre symbol (a | n).
- Task
Calculate the Jacobi symbol (a | n).
- Reference
Factor
The jacobi
word already exists in the math.extras
vocabulary. See the implementation here.
Go
The big.Jacobi function in the standard library (for 'big integers') returns the Jacobi symbol for given values of 'a' and 'n'.
This translates the Lua code in the above referenced Wikipedia article to Go (for 8 byte integers) and checks that it gives the same answers for a small table of values - which it does. <lang go>package main
import (
"fmt" "log" "math/big"
)
func jacobi(a, n uint64) int {
if n%2 == 0 { log.Fatal("'n' must be a positive odd integer") } a %= n result := 1 for a != 0 { for a%2 == 0 { a /= 2 nn := n % 8 if nn == 3 || nn == 5 { result = -result } } a, n = n, a if a%4 == 3 && n%4 == 3 { result = -result } a %= n } if n == 1 { return result } return 0
}
func main() {
fmt.Println("Using hand-coded version:") fmt.Println("n/a 0 1 2 3 4 5 6 7 8 9") fmt.Println("---------------------------------") for n := uint64(1); n <= 17; n += 2 { fmt.Printf("%2d ", n) for a := uint64(0); a <= 9; a++ { fmt.Printf(" % d", jacobi(a, n)) } fmt.Println() }
ba, bn := new(big.Int), new(big.Int) fmt.Println("\nUsing standard library function:") fmt.Println("n/a 0 1 2 3 4 5 6 7 8 9") fmt.Println("---------------------------------") for n := uint64(1); n <= 17; n += 2 { fmt.Printf("%2d ", n) for a := uint64(0); a <= 9; a++ { ba.SetUint64(a) bn.SetUint64(n) fmt.Printf(" % d", big.Jacobi(ba, bn)) } fmt.Println() }
}</lang>
- Output:
Using hand-coded version: n/a 0 1 2 3 4 5 6 7 8 9 --------------------------------- 1 1 1 1 1 1 1 1 1 1 1 3 0 1 -1 0 1 -1 0 1 -1 0 5 0 1 -1 -1 1 0 1 -1 -1 1 7 0 1 1 -1 1 -1 -1 0 1 1 9 0 1 1 0 1 1 0 1 1 0 11 0 1 -1 1 1 1 -1 -1 -1 1 13 0 1 -1 1 1 -1 -1 -1 -1 1 15 0 1 1 0 1 0 0 -1 1 0 17 0 1 1 -1 1 -1 -1 -1 1 1 Using standard library function: n/a 0 1 2 3 4 5 6 7 8 9 --------------------------------- 1 1 1 1 1 1 1 1 1 1 1 3 0 1 -1 0 1 -1 0 1 -1 0 5 0 1 -1 -1 1 0 1 -1 -1 1 7 0 1 1 -1 1 -1 -1 0 1 1 9 0 1 1 0 1 1 0 1 1 0 11 0 1 -1 1 1 1 -1 -1 -1 1 13 0 1 -1 1 1 -1 -1 -1 -1 1 15 0 1 1 0 1 0 0 -1 1 0 17 0 1 1 -1 1 -1 -1 -1 1 1
Julia
<lang julia>function jacobi(a, n)
a %= n result = 1 while a != 0 while iseven(a) a ÷= 2 ((n % 8) in [3, 5]) && (result *= -1) end a, n = n, a (a % 4 == 3) && (n % 4 == 3) && (result *= -1) a %= n end return n == 1 ? result : 0
end
print(" Table of jacobi(a, n) for a 1 to 12, n 1 to 31\n 1 2 3 4 5 6 7 8",
" 9 10 11 12\nn\n_____________________________________________________")
for n in 1:2:31
print("\n", rpad(n, 3)) for a in 1:11 print(lpad(jacobi(a, n), 4)) end
end
</lang>
- Output:
Table of jacobi(a, n) for a 1 to 12, n 1 to 31 1 2 3 4 5 6 7 8 9 10 11 12 n _____________________________________________________ 1 1 1 1 1 1 1 1 1 1 1 1 3 1 -1 0 1 -1 0 1 -1 0 1 -1 5 1 -1 -1 1 0 1 -1 -1 1 0 1 7 1 1 -1 1 -1 -1 0 1 1 -1 1 9 1 1 0 1 1 0 1 1 0 1 1 11 1 -1 1 1 1 -1 -1 -1 1 -1 0 13 1 -1 1 1 -1 -1 -1 -1 1 1 -1 15 1 1 0 1 0 0 -1 1 0 0 -1 17 1 1 -1 1 -1 -1 -1 1 1 -1 -1 19 1 -1 -1 1 1 1 1 -1 1 -1 1 21 1 -1 0 1 1 0 0 -1 0 -1 -1 23 1 1 1 1 -1 1 -1 1 1 -1 -1 25 1 1 1 1 0 1 1 1 1 0 1 27 1 -1 0 1 -1 0 1 -1 0 1 -1 29 1 -1 -1 1 1 1 1 -1 1 -1 -1 31 1 1 -1 1 1 -1 1 1 1 1 -1
Perl 6
<lang perl6># Jacobi function sub infix:<J> (Int $k is copy, Int $n is copy where * % 2) {
$k %= $n; my $jacobi = 1; while $k { while $k %% 2 { $k div= 2; $jacobi *= -1 if $n % 8 == 3 | 5; } ($k, $n) = $n, $k; $jacobi *= -1 if 3 == $n%4 & $k%4; $k %= $n; } $n == 1 ?? $jacobi !! 0
}
- Testing
my $maxa = 30; my $maxn = 29;
say 'n\k ', (1..$maxa).fmt: '%3d'; say ' ', '-' x 4 * $maxa; for 1,*+2 … $maxn -> $n {
print $n.fmt: '%3d'; for 1..$maxa -> $k { print ($k J $n).fmt: '%4d'; } print "\n";
}</lang>
- Output:
n\k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ------------------------------------------------------------------------------------------------------------------------ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 5 1 -1 -1 1 0 1 -1 -1 1 0 1 -1 -1 1 0 1 -1 -1 1 0 1 -1 -1 1 0 1 -1 -1 1 0 7 1 1 -1 1 -1 -1 0 1 1 -1 1 -1 -1 0 1 1 -1 1 -1 -1 0 1 1 -1 1 -1 -1 0 1 1 9 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 11 1 -1 1 1 1 -1 -1 -1 1 -1 0 1 -1 1 1 1 -1 -1 -1 1 -1 0 1 -1 1 1 1 -1 -1 -1 13 1 -1 1 1 -1 -1 -1 -1 1 1 -1 1 0 1 -1 1 1 -1 -1 -1 -1 1 1 -1 1 0 1 -1 1 1 15 1 1 0 1 0 0 -1 1 0 0 -1 0 -1 -1 0 1 1 0 1 0 0 -1 1 0 0 -1 0 -1 -1 0 17 1 1 -1 1 -1 -1 -1 1 1 -1 -1 -1 1 -1 1 1 0 1 1 -1 1 -1 -1 -1 1 1 -1 -1 -1 1 19 1 -1 -1 1 1 1 1 -1 1 -1 1 -1 -1 -1 -1 1 1 -1 0 1 -1 -1 1 1 1 1 -1 1 -1 1 21 1 -1 0 1 1 0 0 -1 0 -1 -1 0 -1 0 0 1 1 0 -1 1 0 1 -1 0 1 1 0 0 -1 0 23 1 1 1 1 -1 1 -1 1 1 -1 -1 1 1 -1 -1 1 -1 1 -1 -1 -1 -1 0 1 1 1 1 -1 1 -1 25 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 27 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 29 1 -1 -1 1 1 1 1 -1 1 -1 -1 -1 1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 1 -1 -1 1 0 1
Python
<lang python>def jacobi(a, n):
a %= n result = 1 while a != 0: while a % 2 == 0: a /= 2 n_mod_8 = n % 8 if n_mod_8 in (3, 5): result = -result a, n = n, a if a % 4 == 3 and n % 4 == 3: result = -result a %= n if n == 1: return result else: return 0</lang>
REXX
A little extra code was added to make a prettier grid.
<lang rexx>/*REXX pgm computes/displays the Jacobi symbol, the # of rows & columns can be specified*/
parse arg rows cols . /*obtain optional arguments from the CL*/
if rows= | rows=="," then rows= 17 /*Not specified? Then use the default.*/
if cols= | cols=="," then cols= 16 /* " " " " " " */
call hdrs /*display the (two) headers to the term*/
do r=1 by 2 to rows; _= right(r, 3) /*build odd (numbered) rows of a table.*/ do c=0 to cols /* [↓] build a column for a table row.*/ _= _ ! right(jacobi(c, r), 2); != '│' /*reset grid end char.*/ end /*c*/ say _ '║'; != '║' /*display a table row; reset grid glyph*/ end /*r*/
say translate(@.2, '╩╧╝', "╬╤╗") /*display the bottom of the grid border*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ hdrs: @.1= 'n/a ║'; do c=0 to cols; @.1= @.1 || right(c, 3)" "; end
L= length(@.1); @.1= left(@.1, L - 1) ; say @.1 @.2= '════╬'; do c=0 to cols; @.2= @.2 || "════╤" ; end L= length(@.2); @.2= left(@.2, L - 1)"╗" ; say @.2 != '║' ; return /*define an external grid border glyph.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/ jacobi: procedure; parse arg a,n; er= '***error***'; $ = 1 /*define result.*/
if n//2==0 then do; say er n " must be a positive odd integer."; exit 13 end a= a // n /*obtain A modulus N */ do while a\==0 /*perform while A isn't zero. */ do while a//2==0; a= a % 2 /*divide A (as a integer) by 2 */ if n//8==3 | n//8==5 then $= -$ /*use N mod 8 */ end /*while a//2==0*/ parse value a n with n a /*swap values of variables: A N */ if a//4==3 & n//4==3 then $= -$ a= a // n end /*while a\==0*/ if n==1 then return $ return 0</lang>
- output when using the default inputs:
n/a ║ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ════╬════╤════╤════╤════╤════╤════╤════╤════╤════╤════╤════╤════╤════╤════╤════╤════╤════╗ 1 ║ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 │ 1 ║ 3 ║ 0 │ 1 │ -1 │ 0 │ 1 │ -1 │ 0 │ 1 │ -1 │ 0 │ 1 │ -1 │ 0 │ 1 │ -1 │ 0 │ 1 ║ 5 ║ 0 │ 1 │ -1 │ -1 │ 1 │ 0 │ 1 │ -1 │ -1 │ 1 │ 0 │ 1 │ -1 │ -1 │ 1 │ 0 │ 1 ║ 7 ║ 0 │ 1 │ 1 │ -1 │ 1 │ -1 │ -1 │ 0 │ 1 │ 1 │ -1 │ 1 │ -1 │ -1 │ 0 │ 1 │ 1 ║ 9 ║ 0 │ 1 │ 1 │ 0 │ 1 │ 1 │ 0 │ 1 │ 1 │ 0 │ 1 │ 1 │ 0 │ 1 │ 1 │ 0 │ 1 ║ 11 ║ 0 │ 1 │ -1 │ 1 │ 1 │ 1 │ -1 │ -1 │ -1 │ 1 │ -1 │ 0 │ 1 │ -1 │ 1 │ 1 │ 1 ║ 13 ║ 0 │ 1 │ -1 │ 1 │ 1 │ -1 │ -1 │ -1 │ -1 │ 1 │ 1 │ -1 │ 1 │ 0 │ 1 │ -1 │ 1 ║ 15 ║ 0 │ 1 │ 1 │ 0 │ 1 │ 0 │ 0 │ -1 │ 1 │ 0 │ 0 │ -1 │ 0 │ -1 │ -1 │ 0 │ 1 ║ 17 ║ 0 │ 1 │ 1 │ -1 │ 1 │ -1 │ -1 │ -1 │ 1 │ 1 │ -1 │ -1 │ -1 │ 1 │ -1 │ 1 │ 1 ║ ════╩════╧════╧════╧════╧════╧════╧════╧════╧════╧════╧════╧════╧════╧════╧════╧════╧════╝
Scheme
<lang scheme>(define jacobi (lambda (a n) (let ((a-mod-n (modulo a n))) (if (zero? a-mod-n) (if (= n 1) 1 0) (if (even? a-mod-n) (case (modulo n 8) ((3 5) (- (jacobi (/ a-mod-n 2) n))) ((1 7) (jacobi (/ a-mod-n 2) n))) (if (and (= (modulo a-mod-n 4) 3) (= (modulo n 4) 3)) (- (jacobi n a-mod-n)) (jacobi n a-mod-n)))))))</lang>
zkl
<lang zkl>fcn jacobi(a,n){
if(n.isEven or n<1) throw(Exception.ValueError("'n' must be a positive odd integer")); a=a%n; result,t := 1,0; while(a!=0){ while(a.isEven){
a/=2; n_mod_8:=n%8; if(n_mod_8==3 or n_mod_8==5) result=-result;
} t,a,n = a,n,t; if(a%4==3 and n%4==3) result=-result; a=a%n; } if(n==1) result else 0
}</lang> <lang zkl>println("Using hand-coded version:"); println("n/a 0 1 2 3 4 5 6 7 8 9"); println("---------------------------------"); foreach n in ([1..17,2]){
print("%2d ".fmt(n)); foreach a in (10){ print(" % d".fmt(jacobi(a,n))) } println();
}</lang>
GNU Multiple Precision Arithmetic Library
<lang zkl>var [const] BI=Import.lib("zklBigNum"); // libGMP println("\nUsing BigInt library function:"); println("n/a 0 1 2 3 4 5 6 7 8 9"); println("---------------------------------"); foreach n in ([1..17,2]){
print("%2d ".fmt(n)); foreach a in (10){ print(" % d".fmt(BI(a).jacobi(n))) } println();
}</lang>
- Output:
Using hand-coded version: n/a 0 1 2 3 4 5 6 7 8 9 --------------------------------- 1 1 1 1 1 1 1 1 1 1 1 3 0 1 -1 0 1 -1 0 1 -1 0 5 0 1 -1 -1 1 0 1 -1 -1 1 7 0 1 1 -1 1 -1 -1 0 1 1 9 0 1 1 0 1 1 0 1 1 0 11 0 1 -1 1 1 1 -1 -1 -1 1 13 0 1 -1 1 1 -1 -1 -1 -1 1 15 0 1 1 0 1 0 0 -1 1 0 17 0 1 1 -1 1 -1 -1 -1 1 1 Using BigInt library function: n/a 0 1 2 3 4 5 6 7 8 9 --------------------------------- 1 1 1 1 1 1 1 1 1 1 1 3 0 1 -1 0 1 -1 0 1 -1 0 5 0 1 -1 -1 1 0 1 -1 -1 1 7 0 1 1 -1 1 -1 -1 0 1 1 9 0 1 1 0 1 1 0 1 1 0 11 0 1 -1 1 1 1 -1 -1 -1 1 13 0 1 -1 1 1 -1 -1 -1 -1 1 15 0 1 1 0 1 0 0 -1 1 0 17 0 1 1 -1 1 -1 -1 -1 1 1