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;


//classical euclidean algorithm implemented non-recusively
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
}

//returns the odd part of n and the exponent of the even power
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)
		}
	}
}


//transforms n into its odd part and returns the exponent of its even part
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;
		}
	}
}

//Jacobi symbol of 2 to the power of n
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;

		//can this step be done only once at the beginning?
		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);
	}
}