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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

//! Helper functions related to [quadratic reciprocity](https://en.wikipedia.org/wiki/Quadratic_reciprocity).
//!
//! This library calculates the [Legendre](http://mathworld.wolfram.com/LegendreSymbol.html)
//! and [Jacobi symbol](http://mathworld.wolfram.com/JacobiSymbol.html) using the
//! [standard algorithm](https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol).
//!
//! Since the Jacobi symbol is a generalization of the Legendre symbol,
//! only a `jacobi` function is provided.
//!

extern crate num;
use num::Integer;
use num::traits::ToPrimitive;


/// Supplementary law
fn two_over(n: usize) -> (i8) {
    if n % 8 == 1 || n % 8 == 7 { 1 } else { -1 }
}

/// Legendre's version of quadratic reciprocity law
/// Returns the sign change needed to keep track of if flipping
fn reciprocity(num: usize, den: usize) -> (i8) {
    if num % 4 == 3 && den % 4 == 3 { -1 } else { 1 }
}

/// Returns the value of the Jacobi symbol for `(a \ n)`
/// Where `a` is an integer and `n` is an odd positive integer
///
/// If the intergers passed in do not match these assumptions,
/// the return value will be `None`. Otherwise the algorithm can be assumed to
/// always succeed.
///
/// # Example
///
/// ```rust
/// # extern crate quadratic;
/// // compute the Jacobi symbol of (2 \ 5)
/// let symb = quadratic::jacobi(2, 5);
/// assert_eq!(Some(-1), symb);
/// ```

pub fn jacobi(a: isize, n: isize) -> Option<i8> {
    // jacobi symbol is only defined for odd positive moduli
    if n.is_even() || n <= 0 {
        return None;
    }

    // Raise a mod n, then start the unsigned algorithm
    let mut acc = 1;
    let mut num = a.mod_floor(&n).to_usize().unwrap();
    let mut den = n as usize;
    loop {
        // reduce numerator
        num = num % den;
        if num == 0 {
            return Some(0);
        }

        // extract factors of two from numerator
        while num.is_even() {
            acc *= two_over(den);
            num /= 2;
        }
        // if numerator is 1 => this sub-symbol is 1
        if num == 1 {
            return Some(acc);
        }
        // shared factors => one sub-symbol is zero
        if num.gcd(&den) > 1 {
            return Some(0);
        }
        // num and den are now odd co-prime, use reciprocity law:
        acc *= reciprocity(num, den);
        let tmp = num;
        num = den;
        den = tmp;
    }
}

#[cfg(test)]
mod tests {
    use jacobi;

    // legendre tests
    #[test]
    fn minus_one_over_p() {
        // 1 mod 4 => 1
        assert_eq!(Some(1), jacobi(-1, 5));
        assert_eq!(Some(1), jacobi(-1, 13));
        // 3 mod 4 => -1
        assert_eq!(Some(-1), jacobi(-1, 3));
        assert_eq!(Some(-1), jacobi(-1, 7));
    }
    #[test]
    fn two_over_p() {
        assert_eq!(Some(-1), jacobi(2, 3));
        assert_eq!(Some(-1), jacobi(2, 5));
        assert_eq!(Some(1), jacobi(2, 7));
        assert_eq!(Some(1), jacobi(2, 17)); // 17 = 1 mod 8
    }
    #[test]
    fn three_over_p() {
        assert_eq!(Some(0), jacobi(3, 3));
        assert_eq!(Some(-1), jacobi(3, 5));
        assert_eq!(Some(-1), jacobi(3, 7));
    }
    #[test]
    fn periodicity() {
        assert_eq!(jacobi(3,5), jacobi(-2,5));
        assert_eq!(jacobi(-1,5), jacobi(4,5));
        assert_eq!(jacobi(11,7), jacobi(4,7));
        assert_eq!(jacobi(-3,7), jacobi(4,7));
        assert_eq!(jacobi(10,7), jacobi(3,7));
    }

    #[test]
    fn jacobi_simple() { // 45 = 3*3*5
        assert_eq!(Some(-1), jacobi(2, 45)); // -1 * -1 * -1
        assert_eq!(Some(0), jacobi(3, 45)); // 0 * 0 * -1
        assert_eq!(Some(-1), jacobi(7, 45)); // -1 * -1 * -1
        assert_eq!(Some(1), jacobi(2, 15)); // (2\5) * (2\3)
        assert_eq!(Some(-1), jacobi(1001, 9907)); // wikepedia example
    }

    #[test]
    fn even_moduli_fails() {
        assert_eq!(None, jacobi(2, 4));
    }
    #[test]
    fn negative_moduli_fails() {
        assert_eq!(None, jacobi(2, -3));
    }
}