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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177

///!
///! A set of mathematical utility functions.
///! TODO: T might be a signed integer! Everything in here only considers unsigned really.
///!
use crate::prelude::*;

mod ct_poly;
pub mod ct_util;
pub mod poly;

use ct_poly::*;
use poly::*;

/// polynomial subtraction, calculates a - b mod modulo
#[cfg_attr(feature = "use_attributes", unsafe_hacspec)]
pub fn sub_poly<T: Numeric + Copy>(a: &Seq<T>, b: &Seq<T>, modulo: T) -> Seq<T> {
    let result = Seq::from_native_slice(&poly_sub(&a.b, &b.b, modulo));
    make_positive(&result, modulo)
}

/// Polynomial Addition, calculates a + b mod modulo
#[cfg_attr(feature = "use_attributes", unsafe_hacspec)]
pub fn add_poly<T: Numeric + Copy>(a: &Seq<T>, b: &Seq<T>, modulo: T) -> Seq<T> {
    let result = Seq::from_native_slice(&poly_add(&a.b, &b.b, modulo));
    make_positive(&result, modulo)
}

/// Simple polynomial multiplication for two fixed size polynomials O(n²) with `a * b mod n`
pub fn mul_poly<T: Numeric + Copy>(a: &Seq<T>, b: &Seq<T>, n: T) -> Seq<T> {
    Seq::from_native_slice(&poly_mul(&a.b, &b.b, n))
}

/// Euclidean polynomial division, calculates `a/b` in `R_n`.
/// Returns `Ok(quotient, remainder)` or `Err("Can't divide these two polynomials")`
pub fn div_poly<T: Integer + Copy>(
    a: &Seq<T>,
    b: &Seq<T>,
    n: T,
) -> Result<(Seq<T>, Seq<T>), &'static str> {
    let (q, r) = poly_div(&a.b, &b.b, n);
    Ok((Seq::from_native_slice(&q), Seq::from_native_slice(&r)))
}

/// Scalar division in `R_p`.
/// Returns `a / scalar mod p`.
#[cfg_attr(feature = "use_attributes", not_hacspec)]
pub fn div_scalar<T: Integer + Copy>(a: &Seq<T>, scalar: T, p: T) -> Seq<T> {
    Seq::from_native_slice(&scalar_div(&a.b, scalar, p))
}

/// Returns degree of polynomial, e.g. for  3x² + 2x + 1 -> 2
#[cfg_attr(feature = "use_attributes", unsafe_hacspec)]
pub fn degree_poly<T: Integer + Copy>(poly: &Seq<T>) -> usize {
    deg(&poly.b)
}

/// Euclidean algorithm to compute the inverse of x in yℤ\[x\]
#[cfg_attr(feature = "use_attributes", unsafe_hacspec)]
pub fn extended_euclid<T: Integer + Copy>(
    x: &Seq<T>,
    irr: &Seq<T>,
    n: T,
) -> Result<Seq<T>, &'static str> {
    let result = extended_euclid_internal(&x.b, &irr.b, n)?;
    Ok(Seq::from_native_slice(&result))
}

/// Returns number of coefficient != 0, e.g. for  -3x⁵ + 3x² + 2x + 1 -> 4
pub fn weight<T: Integer + Copy>(poly: &Seq<T>) -> usize {
    let tmp = Seq::from_seq(poly);
    let mut weight = 0;
    for i in 0..tmp.len() {
        if !tmp[i].equal(T::default()) {
            weight = weight + 1;
        }
    }
    weight
}

/// makes coefficients positiv, e.g. -3 mod 4 = 1
pub fn make_positive<T: Numeric + Copy>(poly: &Seq<T>, q: T) -> Seq<T> {
    Seq::from_native_slice(&make_positive_internal(&poly.b, q))
}

/// Extended euclidean algorithm to compute the inverse of x in ℤ/n
///
/// **Panics** if x is not invertible.
///
#[inline]
#[cfg_attr(feature = "use_attributes", not_hacspec)]
pub(crate) fn extended_euclid_invert<T: Integer + Copy>(x: T, n: T, signed: bool) -> T {
    let mut t = T::ZERO();
    let mut r = n;
    let mut new_t = T::ONE();
    let mut new_r = x;

    while !new_r.equal(T::ZERO()) {
        let q: T = r.divide(new_r);

        let tmp = new_r;
        // XXX: a little hacky
        let tmp_prod = q * new_r;
        let mut tmp_r = r;
        while tmp_r.less_than(tmp_prod) {
            tmp_r = tmp_r + n;
        }
        new_r = tmp_r - tmp_prod;
        r = tmp;

        let tmp = new_t;
        // XXX: a little hacky
        let tmp_prod = q * new_t;
        let mut tmp_t = t;
        while tmp_t.less_than(tmp_prod) {
            tmp_t = tmp_t + n;
        }
        new_t = tmp_t - tmp_prod;
        t = tmp;
    }

    if r.greater_than(T::ONE()) && !x.equal(T::ZERO()) {
        panic!("{:x?} is not invertible in ℤ/{:x?}", x, n);
    }
    if t.less_than(T::ZERO()) {
        if signed {
            t = t.absolute()
        } else {
            t = t + n
        };
    };

    t
}

/// Makes poly to an element of R_modulo \ irr
#[cfg_attr(feature = "use_attributes", unsafe_hacspec)]
pub fn poly_to_ring<T: Integer + Copy>(irr: &Seq<T>, poly: &Seq<T>, modulus: T) -> (Seq<T>, bool) {
    match div_poly(&poly, &irr, modulus) {
        Ok(pre) => {
            let pre = pre.1;
            (make_positive(&pre, modulus), true)
        }
        Err(_) => (Seq::new(1usize), false),
    }
}

/// Polynomial multiplication of two size fixed polynomials in R_modulo \ irr
#[cfg_attr(feature = "use_attributes", unsafe_hacspec)]
pub fn mul_poly_irr<T: Integer + Copy>(a: &Seq<T>, b: &Seq<T>, irr: &Seq<T>, modulo: T) -> Seq<T> {
    let mut result: Seq<T> = Seq::new(a.len());
    for i in 0..a.len() {
        if a[i].equal(T::default()) {
            continue;
        }
        for j in 0..b.len() {
            if b[j].equal(T::default()) {
                continue;
            }

            if i + j > a.len() - 2 {
                // modulo irr | factor is the coefficient
                let factor = a[i] * b[j];
                result[1 + ((i + j) % (a.len() - 1))] =
                    result[1 + (i + j) % (a.len() - 1)] - factor * irr[1];
                result[(i + j) % (a.len() - 1)] = result[(i + j) % (a.len() - 1)] - factor * irr[0];
                continue;
            }
            result[i + j] = result[i + j] + a[i] * b[j];
        }
    }
    if modulo.greater_than(T::default()) {
        for i in 0..result.len() {
            result[i] = result[i].modulo(modulo);
        }
    }
    make_positive(&result, modulo)
}