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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
use std::{u64, i64};
use std::cmp::{max, min};
use arithmetic::mod_abs;
pub fn gcd(a: u64, b: u64) -> u64 {
let mut bigger = max(a,b);
let mut smaller = min(a,b);
let mut remainder = bigger % smaller;
while remainder != 0 {
bigger = smaller;
smaller = remainder;
remainder = bigger % smaller
}
smaller
}
pub fn split_odd_even(n: u64) -> (u64, u64) {
let mut even_power = 0;
let mut n_shift = n;
loop {
if n_shift & 1 == 0 {
even_power += 1;
n_shift = n_shift.rotate_right(1)
} else {
return (n_shift, even_power)
}
}
}
pub fn mut_even_power(n: &mut u64) -> u64 {
let mut even_power = 0;
loop {
if *n & 1 == 0 {
even_power += 1;
*n = n.rotate_right(1);
} else {
return even_power;
}
}
}
fn jacobi_even_power(pow: u64, d:u64) -> i8 {
if d % 8 == 1 || d % 8 == 7 || pow % 2 == 0 {
1
} else {
-1
}
}
fn jacobi_neg(d: u64) -> i8 {
if d % 4 == 1 {
1
} else {
-1
}
}
fn jacobi_flip(n: u64, d: u64) -> bool {
if n % 4 == 1 || d % 4 == 1 {
false
} else {
true
}
}
pub fn jacobi_symbol(n: i64, d: u64) -> i8 {
let mut num = mod_abs(n, d);
let mut den = d;
let mut jacobi: i8 = 1;
loop{
num = num % den;
if gcd(num, den) != 1 {
return 0
}
let num_even_pow = mut_even_power(&mut num);
jacobi *= jacobi_even_power(num_even_pow, d);
if num == 1 {
return jacobi
}
if num == d-1 {
return jacobi * jacobi_neg(den)
}
if jacobi_flip(num, den) {
jacobi *= -1;
}
std::mem::swap(&mut num, &mut den);
}
}