cairo_vm/math_utils/
mod.rs

1mod is_prime;
2
3pub use is_prime::is_prime;
4
5use core::cmp::min;
6
7use crate::stdlib::{boxed::Box, ops::Shr, prelude::Vec};
8use crate::types::errors::math_errors::MathError;
9use crate::utils::CAIRO_PRIME;
10use crate::Felt252;
11use lazy_static::lazy_static;
12use num_bigint::{BigInt, BigUint, RandBigInt, ToBigInt};
13use num_integer::Integer;
14use num_traits::{One, Signed, Zero};
15use rand::{rngs::SmallRng, SeedableRng};
16use starknet_types_core::felt::NonZeroFelt;
17
18lazy_static! {
19    pub static ref SIGNED_FELT_MAX: BigUint = (&*CAIRO_PRIME).shr(1_u32);
20    static ref POWERS_OF_TWO: Vec<NonZeroFelt> =
21        core::iter::successors(Some(Felt252::ONE), |x| Some(x * Felt252::TWO))
22            .take(252)
23            .map(|x| x.try_into().unwrap())
24            .collect::<Vec<_>>();
25}
26
27pub const STWO_PRIME: u64 = (1 << 31) - 1;
28const STWO_PRIME_U128: u128 = STWO_PRIME as u128;
29const MASK_36: u64 = (1 << 36) - 1;
30const MASK_8: u64 = (1 << 8) - 1;
31
32/// Returns the `n`th (up to the `251`th power) power of 2 as a [`Felt252`]
33/// in constant time.
34/// It silently returns `1` if the input is out of bounds.
35pub fn pow2_const(n: u32) -> Felt252 {
36    // If the conversion fails then it's out of range and we compute the power as usual
37    POWERS_OF_TWO
38        .get(n as usize)
39        .unwrap_or(&POWERS_OF_TWO[0])
40        .into()
41}
42
43/// Returns the `n`th (up to the `251`th power) power of 2 as a [`&stark_felt::NonZeroFelt`]
44/// in constant time.
45/// It silently returns `1` if the input is out of bounds.
46pub fn pow2_const_nz(n: u32) -> &'static NonZeroFelt {
47    // If the conversion fails then it's out of range and we compute the power as usual
48    POWERS_OF_TWO.get(n as usize).unwrap_or(&POWERS_OF_TWO[0])
49}
50
51/// Converts [`Felt252`] into a [`BigInt`] number in the range: `(- FIELD / 2, FIELD / 2)`.
52///
53/// # Examples
54///
55/// ```
56/// # use cairo_vm::{Felt252, math_utils::signed_felt};
57/// # use num_bigint::BigInt;
58/// let positive = Felt252::from(5);
59/// assert_eq!(signed_felt(positive), BigInt::from(5));
60///
61/// let negative = Felt252::MAX;
62/// assert_eq!(signed_felt(negative), BigInt::from(-1));
63/// ```
64pub fn signed_felt(felt: Felt252) -> BigInt {
65    let biguint = felt.to_biguint();
66    if biguint > *SIGNED_FELT_MAX {
67        BigInt::from_biguint(num_bigint::Sign::Minus, &*CAIRO_PRIME - &biguint)
68    } else {
69        biguint.to_bigint().expect("cannot fail")
70    }
71}
72
73/// QM31 utility function, used specifically for Stwo.
74/// QM31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
75/// Reads four u64 coordinates from a single Felt252.
76/// STWO_PRIME fits in 36 bits, hence each coordinate can be represented by 36 bits and a QM31
77/// element can be stored in the first 144 bits of a Felt252.
78/// Returns an error if the input has over 144 bits or any coordinate is unreduced.
79fn qm31_packed_reduced_read_coordinates(felt: Felt252) -> Result<[u64; 4], MathError> {
80    let limbs = felt.to_le_digits();
81    if limbs[3] != 0 || limbs[2] >= 1 << 16 {
82        return Err(MathError::QM31UnreducedError(Box::new(felt)));
83    }
84    let coordinates = [
85        (limbs[0] & MASK_36),
86        ((limbs[0] >> 36) + ((limbs[1] & MASK_8) << 28)),
87        ((limbs[1] >> 8) & MASK_36),
88        ((limbs[1] >> 44) + (limbs[2] << 20)),
89    ];
90    for x in coordinates.iter() {
91        if *x >= STWO_PRIME {
92            return Err(MathError::QM31UnreducedError(Box::new(felt)));
93        }
94    }
95    Ok(coordinates)
96}
97
98/// QM31 utility function, used specifically for Stwo.
99/// QM31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
100/// Reduces four u64 coordinates and packs them into a single Felt252.
101/// STWO_PRIME fits in 36 bits, hence each coordinate can be represented by 36 bits and a QM31
102/// element can be stored in the first 144 bits of a Felt252.
103pub fn qm31_coordinates_to_packed_reduced(coordinates: [u64; 4]) -> Felt252 {
104    let bytes_part1 = ((coordinates[0] % STWO_PRIME) as u128
105        + (((coordinates[1] % STWO_PRIME) as u128) << 36))
106        .to_le_bytes();
107    let bytes_part2 = ((coordinates[2] % STWO_PRIME) as u128
108        + (((coordinates[3] % STWO_PRIME) as u128) << 36))
109        .to_le_bytes();
110    let mut result_bytes = [0u8; 32];
111    result_bytes[0..9].copy_from_slice(&bytes_part1[0..9]);
112    result_bytes[9..18].copy_from_slice(&bytes_part2[0..9]);
113    Felt252::from_bytes_le(&result_bytes)
114}
115
116/// QM31 utility function, used specifically for Stwo.
117/// QM31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
118/// Computes the addition of two QM31 elements in reduced form.
119/// Returns an error if either operand is not reduced.
120pub fn qm31_packed_reduced_add(felt1: Felt252, felt2: Felt252) -> Result<Felt252, MathError> {
121    let coordinates1 = qm31_packed_reduced_read_coordinates(felt1)?;
122    let coordinates2 = qm31_packed_reduced_read_coordinates(felt2)?;
123    let result_unreduced_coordinates = [
124        coordinates1[0] + coordinates2[0],
125        coordinates1[1] + coordinates2[1],
126        coordinates1[2] + coordinates2[2],
127        coordinates1[3] + coordinates2[3],
128    ];
129    Ok(qm31_coordinates_to_packed_reduced(
130        result_unreduced_coordinates,
131    ))
132}
133
134/// QM31 utility function, used specifically for Stwo.
135/// QM31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
136/// Computes the negative of a QM31 element in reduced form.
137/// Returns an error if the input is not reduced.
138pub fn qm31_packed_reduced_neg(felt: Felt252) -> Result<Felt252, MathError> {
139    let coordinates = qm31_packed_reduced_read_coordinates(felt)?;
140    Ok(qm31_coordinates_to_packed_reduced([
141        STWO_PRIME - coordinates[0],
142        STWO_PRIME - coordinates[1],
143        STWO_PRIME - coordinates[2],
144        STWO_PRIME - coordinates[3],
145    ]))
146}
147
148/// QM31 utility function, used specifically for Stwo.
149/// QM31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
150/// Computes the subtraction of two QM31 elements in reduced form.
151/// Returns an error if either operand is not reduced.
152pub fn qm31_packed_reduced_sub(felt1: Felt252, felt2: Felt252) -> Result<Felt252, MathError> {
153    let coordinates1 = qm31_packed_reduced_read_coordinates(felt1)?;
154    let coordinates2 = qm31_packed_reduced_read_coordinates(felt2)?;
155    let result_unreduced_coordinates = [
156        STWO_PRIME + coordinates1[0] - coordinates2[0],
157        STWO_PRIME + coordinates1[1] - coordinates2[1],
158        STWO_PRIME + coordinates1[2] - coordinates2[2],
159        STWO_PRIME + coordinates1[3] - coordinates2[3],
160    ];
161    Ok(qm31_coordinates_to_packed_reduced(
162        result_unreduced_coordinates,
163    ))
164}
165
166/// QM31 utility function, used specifically for Stwo.
167/// QM31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
168/// Computes the multiplication of two QM31 elements in reduced form.
169/// Returns an error if either operand is not reduced.
170pub fn qm31_packed_reduced_mul(felt1: Felt252, felt2: Felt252) -> Result<Felt252, MathError> {
171    let coordinates1_u64 = qm31_packed_reduced_read_coordinates(felt1)?;
172    let coordinates2_u64 = qm31_packed_reduced_read_coordinates(felt2)?;
173    let coordinates1 = coordinates1_u64.map(u128::from);
174    let coordinates2 = coordinates2_u64.map(u128::from);
175
176    let result_coordinates = [
177        ((5 * STWO_PRIME_U128 * STWO_PRIME_U128 + coordinates1[0] * coordinates2[0]
178            - coordinates1[1] * coordinates2[1]
179            + 2 * coordinates1[2] * coordinates2[2]
180            - 2 * coordinates1[3] * coordinates2[3]
181            - coordinates1[2] * coordinates2[3]
182            - coordinates1[3] * coordinates2[2])
183            % STWO_PRIME_U128) as u64,
184        ((STWO_PRIME_U128 * STWO_PRIME_U128
185            + coordinates1[0] * coordinates2[1]
186            + coordinates1[1] * coordinates2[0]
187            + 2 * (coordinates1[2] * coordinates2[3] + coordinates1[3] * coordinates2[2])
188            + coordinates1[2] * coordinates2[2]
189            - coordinates1[3] * coordinates2[3])
190            % STWO_PRIME_U128) as u64,
191        2 * STWO_PRIME * STWO_PRIME + coordinates1_u64[0] * coordinates2_u64[2]
192            - coordinates1_u64[1] * coordinates2_u64[3]
193            + coordinates1_u64[2] * coordinates2_u64[0]
194            - coordinates1_u64[3] * coordinates2_u64[1],
195        coordinates1_u64[0] * coordinates2_u64[3]
196            + coordinates1_u64[1] * coordinates2_u64[2]
197            + coordinates1_u64[2] * coordinates2_u64[1]
198            + coordinates1_u64[3] * coordinates2_u64[0],
199    ];
200    Ok(qm31_coordinates_to_packed_reduced(result_coordinates))
201}
202
203/// M31 utility function, used specifically for Stwo.
204/// M31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
205/// Computes the inverse in the M31 field using Fermat's little theorem, i.e., returns
206/// `v^(STWO_PRIME-2) modulo STWO_PRIME`, which is the inverse of v unless v % STWO_PRIME == 0.
207pub fn pow2147483645(v: u64) -> u64 {
208    let t0 = (sqn(v, 2) * v) % STWO_PRIME;
209    let t1 = (sqn(t0, 1) * t0) % STWO_PRIME;
210    let t2 = (sqn(t1, 3) * t0) % STWO_PRIME;
211    let t3 = (sqn(t2, 1) * t0) % STWO_PRIME;
212    let t4 = (sqn(t3, 8) * t3) % STWO_PRIME;
213    let t5 = (sqn(t4, 8) * t3) % STWO_PRIME;
214    (sqn(t5, 7) * t2) % STWO_PRIME
215}
216
217/// M31 utility function, used specifically for Stwo.
218/// M31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
219/// Computes `v^(2^n) modulo STWO_PRIME`.
220fn sqn(v: u64, n: usize) -> u64 {
221    let mut u = v;
222    for _ in 0..n {
223        u = (u * u) % STWO_PRIME;
224    }
225    u
226}
227
228/// QM31 utility function, used specifically for Stwo.
229/// QM31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
230/// Computes the inverse of a QM31 element in reduced form.
231/// Returns an error if the denominator is zero or either operand is not reduced.
232pub fn qm31_packed_reduced_inv(felt: Felt252) -> Result<Felt252, MathError> {
233    if felt.is_zero() {
234        return Err(MathError::DividedByZero);
235    }
236    let coordinates = qm31_packed_reduced_read_coordinates(felt)?;
237
238    let b2_r = (coordinates[2] * coordinates[2] + STWO_PRIME * STWO_PRIME
239        - coordinates[3] * coordinates[3])
240        % STWO_PRIME;
241    let b2_i = (2 * coordinates[2] * coordinates[3]) % STWO_PRIME;
242
243    let denom_r = (coordinates[0] * coordinates[0] + STWO_PRIME * STWO_PRIME
244        - coordinates[1] * coordinates[1]
245        + 2 * STWO_PRIME
246        - 2 * b2_r
247        + b2_i)
248        % STWO_PRIME;
249    let denom_i =
250        (2 * coordinates[0] * coordinates[1] + 3 * STWO_PRIME - 2 * b2_i - b2_r) % STWO_PRIME;
251
252    let denom_norm_squared = (denom_r * denom_r + denom_i * denom_i) % STWO_PRIME;
253    let denom_norm_inverse_squared = pow2147483645(denom_norm_squared);
254
255    let denom_inverse_r = (denom_r * denom_norm_inverse_squared) % STWO_PRIME;
256    let denom_inverse_i = ((STWO_PRIME - denom_i) * denom_norm_inverse_squared) % STWO_PRIME;
257
258    Ok(qm31_coordinates_to_packed_reduced([
259        coordinates[0] * denom_inverse_r + STWO_PRIME * STWO_PRIME
260            - coordinates[1] * denom_inverse_i,
261        coordinates[0] * denom_inverse_i + coordinates[1] * denom_inverse_r,
262        coordinates[3] * denom_inverse_i + STWO_PRIME * STWO_PRIME
263            - coordinates[2] * denom_inverse_r,
264        2 * STWO_PRIME * STWO_PRIME
265            - coordinates[2] * denom_inverse_i
266            - coordinates[3] * denom_inverse_r,
267    ]))
268}
269
270/// QM31 utility function, used specifically for Stwo.
271/// QM31 operations are to be relocated into https://github.com/lambdaclass/lambdaworks.
272/// Computes the division of two QM31 elements in reduced form.
273/// Returns an error if the input is zero.
274pub fn qm31_packed_reduced_div(felt1: Felt252, felt2: Felt252) -> Result<Felt252, MathError> {
275    let felt2_inv = qm31_packed_reduced_inv(felt2)?;
276    qm31_packed_reduced_mul(felt1, felt2_inv)
277}
278
279///Returns the integer square root of the nonnegative integer n.
280///This is the floor of the exact square root of n.
281///Unlike math.sqrt(), this function doesn't have rounding error issues.
282pub fn isqrt(n: &BigUint) -> Result<BigUint, MathError> {
283    /*    # The following algorithm was copied from
284    # https://stackoverflow.com/questions/15390807/integer-square-root-in-python.
285    x = n
286    y = (x + 1) // 2
287    while y < x:
288        x = y
289        y = (x + n // x) // 2
290    assert x**2 <= n < (x + 1) ** 2
291    return x*/
292
293    let mut x = n.clone();
294    //n.shr(1) = n.div_floor(2)
295    let mut y = (&x + 1_u32).shr(1_u32);
296
297    while y < x {
298        x = y;
299        y = (&x + n.div_floor(&x)).shr(1_u32);
300    }
301
302    if !(&BigUint::pow(&x, 2_u32) <= n && n < &BigUint::pow(&(&x + 1_u32), 2_u32)) {
303        return Err(MathError::FailedToGetSqrt(Box::new(n.clone())));
304    };
305    Ok(x)
306}
307
308/// Performs integer division between x and y; fails if x is not divisible by y.
309pub fn safe_div(x: &Felt252, y: &Felt252) -> Result<Felt252, MathError> {
310    let (q, r) = x.div_rem(&y.try_into().map_err(|_| MathError::DividedByZero)?);
311
312    if !r.is_zero() {
313        Err(MathError::SafeDivFail(Box::new((*x, *y))))
314    } else {
315        Ok(q)
316    }
317}
318
319/// Performs integer division between x and y; fails if x is not divisible by y.
320pub fn safe_div_bigint(x: &BigInt, y: &BigInt) -> Result<BigInt, MathError> {
321    if y.is_zero() {
322        return Err(MathError::DividedByZero);
323    }
324
325    let (q, r) = x.div_mod_floor(y);
326
327    if !r.is_zero() {
328        return Err(MathError::SafeDivFailBigInt(Box::new((
329            x.clone(),
330            y.clone(),
331        ))));
332    }
333
334    Ok(q)
335}
336
337/// Performs integer division between x and y; fails if x is not divisible by y.
338pub fn safe_div_usize(x: usize, y: usize) -> Result<usize, MathError> {
339    if y.is_zero() {
340        return Err(MathError::DividedByZero);
341    }
342
343    let (q, r) = x.div_mod_floor(&y);
344
345    if !r.is_zero() {
346        return Err(MathError::SafeDivFailUsize(Box::new((x, y))));
347    }
348
349    Ok(q)
350}
351
352///Returns num_a^-1 mod p
353pub(crate) fn mul_inv(num_a: &BigInt, p: &BigInt) -> BigInt {
354    if num_a.is_zero() {
355        return BigInt::zero();
356    }
357    let mut a = num_a.abs();
358    let x_sign = num_a.signum();
359    let mut b = p.abs();
360    let (mut x, mut r) = (BigInt::one(), BigInt::zero());
361    let (mut c, mut q);
362    while !b.is_zero() {
363        (q, c) = a.div_mod_floor(&b);
364        x -= &q * &r;
365        (a, b, r, x) = (b, c, x, r)
366    }
367
368    x * x_sign
369}
370
371///Returns x, y, g such that g = x*a + y*b = gcd(a, b).
372fn igcdex(num_a: &BigInt, num_b: &BigInt) -> (BigInt, BigInt, BigInt) {
373    match (num_a, num_b) {
374        (a, b) if a.is_zero() && b.is_zero() => (BigInt::zero(), BigInt::one(), BigInt::zero()),
375        (a, _) if a.is_zero() => (BigInt::zero(), num_b.signum(), num_b.abs()),
376        (_, b) if b.is_zero() => (num_a.signum(), BigInt::zero(), num_a.abs()),
377        _ => {
378            let mut a = num_a.abs();
379            let x_sign = num_a.signum();
380            let mut b = num_b.abs();
381            let y_sign = num_b.signum();
382            let (mut x, mut y, mut r, mut s) =
383                (BigInt::one(), BigInt::zero(), BigInt::zero(), BigInt::one());
384            let (mut c, mut q);
385            while !b.is_zero() {
386                (q, c) = a.div_mod_floor(&b);
387                x -= &q * &r;
388                y -= &q * &s;
389                (a, b, r, s, x, y) = (b, c, x, y, r, s)
390            }
391            (x * x_sign, y * y_sign, a)
392        }
393    }
394}
395
396///Finds a nonnegative integer x < p such that (m * x) % p == n.
397pub fn div_mod(n: &BigInt, m: &BigInt, p: &BigInt) -> Result<BigInt, MathError> {
398    let (a, _, c) = igcdex(m, p);
399    if !c.is_one() {
400        return Err(MathError::DivModIgcdexNotZero(Box::new((
401            n.clone(),
402            m.clone(),
403            p.clone(),
404        ))));
405    }
406    Ok((n * a).mod_floor(p))
407}
408
409pub(crate) fn div_mod_unsigned(
410    n: &BigUint,
411    m: &BigUint,
412    p: &BigUint,
413) -> Result<BigUint, MathError> {
414    // BigUint to BigInt conversion cannot fail & div_mod will always return a positive value if all values are positive so we can safely unwrap here
415    div_mod(
416        &n.to_bigint().unwrap(),
417        &m.to_bigint().unwrap(),
418        &p.to_bigint().unwrap(),
419    )
420    .map(|i| i.to_biguint().unwrap())
421}
422
423pub fn ec_add(
424    point_a: (BigInt, BigInt),
425    point_b: (BigInt, BigInt),
426    prime: &BigInt,
427) -> Result<(BigInt, BigInt), MathError> {
428    let m = line_slope(&point_a, &point_b, prime)?;
429    let x = (m.clone() * m.clone() - point_a.0.clone() - point_b.0).mod_floor(prime);
430    let y = (m * (point_a.0 - x.clone()) - point_a.1).mod_floor(prime);
431    Ok((x, y))
432}
433
434/// Computes the slope of the line connecting the two given EC points over the field GF(p).
435/// Assumes the points are given in affine form (x, y) and have different x coordinates.
436pub fn line_slope(
437    point_a: &(BigInt, BigInt),
438    point_b: &(BigInt, BigInt),
439    prime: &BigInt,
440) -> Result<BigInt, MathError> {
441    debug_assert!(!(&point_a.0 - &point_b.0).is_multiple_of(prime));
442    div_mod(
443        &(&point_a.1 - &point_b.1),
444        &(&point_a.0 - &point_b.0),
445        prime,
446    )
447}
448
449///  Doubles a point on an elliptic curve with the equation y^2 = x^3 + alpha*x + beta mod p.
450/// Assumes the point is given in affine form (x, y) and has y != 0.
451pub fn ec_double(
452    point: (BigInt, BigInt),
453    alpha: &BigInt,
454    prime: &BigInt,
455) -> Result<(BigInt, BigInt), MathError> {
456    let m = ec_double_slope(&point, alpha, prime)?;
457    let x = ((&m * &m) - (2_i32 * &point.0)).mod_floor(prime);
458    let y = (m * (point.0 - &x) - point.1).mod_floor(prime);
459    Ok((x, y))
460}
461/// Computes the slope of an elliptic curve with the equation y^2 = x^3 + alpha*x + beta mod p, at
462/// the given point.
463/// Assumes the point is given in affine form (x, y) and has y != 0.
464pub fn ec_double_slope(
465    point: &(BigInt, BigInt),
466    alpha: &BigInt,
467    prime: &BigInt,
468) -> Result<BigInt, MathError> {
469    debug_assert!(!point.1.is_multiple_of(prime));
470    div_mod(
471        &(3_i32 * &point.0 * &point.0 + alpha),
472        &(2_i32 * &point.1),
473        prime,
474    )
475}
476
477// Adapted from sympy _sqrt_prime_power with k == 1
478pub fn sqrt_prime_power(a: &BigUint, p: &BigUint) -> Option<BigUint> {
479    if p.is_zero() || !is_prime(p) {
480        return None;
481    }
482    let two = BigUint::from(2_u32);
483    let a = a.mod_floor(p);
484    if p == &two {
485        return Some(a);
486    }
487    if !(a < two || (a.modpow(&(p - 1_u32).div_floor(&two), p)).is_one()) {
488        return None;
489    };
490
491    if p.mod_floor(&BigUint::from(4_u32)) == 3_u32.into() {
492        let res = a.modpow(&(p + 1_u32).div_floor(&BigUint::from(4_u32)), p);
493        return Some(min(res.clone(), p - res));
494    };
495
496    if p.mod_floor(&BigUint::from(8_u32)) == 5_u32.into() {
497        let sign = a.modpow(&(p - 1_u32).div_floor(&BigUint::from(4_u32)), p);
498        if sign.is_one() {
499            let res = a.modpow(&(p + 3_u32).div_floor(&BigUint::from(8_u32)), p);
500            return Some(min(res.clone(), p - res));
501        } else {
502            let b = (4_u32 * &a).modpow(&(p - 5_u32).div_floor(&BigUint::from(8_u32)), p);
503            let x = (2_u32 * &a * b).mod_floor(p);
504            if x.modpow(&two, p) == a {
505                return Some(x);
506            }
507        }
508    };
509
510    Some(sqrt_tonelli_shanks(&a, p))
511}
512
513fn sqrt_tonelli_shanks(n: &BigUint, prime: &BigUint) -> BigUint {
514    // Based on Tonelli-Shanks' algorithm for finding square roots
515    // and sympy's library implementation of said algorithm.
516    if n.is_zero() || n.is_one() {
517        return n.clone();
518    }
519    let s = (prime - 1_u32).trailing_zeros().unwrap_or_default();
520    let t = prime >> s;
521    let a = n.modpow(&t, prime);
522    // Rng is not critical here so its safe to use a seeded value
523    let mut rng = SmallRng::seed_from_u64(11480028852697973135);
524    let mut d;
525    loop {
526        d = RandBigInt::gen_biguint_range(&mut rng, &BigUint::from(2_u32), &(prime - 1_u32));
527        let r = legendre_symbol(&d, prime);
528        if r == -1 {
529            break;
530        };
531    }
532    d = d.modpow(&t, prime);
533    let mut m = BigUint::zero();
534    let mut exponent = BigUint::one() << (s - 1);
535    let mut adm;
536    for i in 0..s as u32 {
537        adm = &a * &d.modpow(&m, prime);
538        adm = adm.modpow(&exponent, prime);
539        exponent >>= 1;
540        if adm == (prime - 1_u32) {
541            m += BigUint::from(1_u32) << i;
542        }
543    }
544    let root_1 =
545        (n.modpow(&((t + 1_u32) >> 1), prime) * d.modpow(&(m >> 1), prime)).mod_floor(prime);
546    let root_2 = prime - &root_1;
547    if root_1 < root_2 {
548        root_1
549    } else {
550        root_2
551    }
552}
553
554/* Disclaimer: Some asumptions have been taken based on the functions that rely on this function, make sure these are true before calling this function individually
555Adpted from sympy implementation, asuming:
556    - p is an odd prime number
557    - a.mod_floor(p) == a
558Returns the Legendre symbol `(a / p)`.
559
560    For an integer ``a`` and an odd prime ``p``, the Legendre symbol is
561    defined as
562
563    .. math ::
564        \genfrac(){}{}{a}{p} = \begin{cases}
565             0 & \text{if } p \text{ divides } a\\
566             1 & \text{if } a \text{ is a quadratic residue modulo } p\\
567            -1 & \text{if } a \text{ is a quadratic nonresidue modulo } p
568        \end{cases}
569*/
570fn legendre_symbol(a: &BigUint, p: &BigUint) -> i8 {
571    if a.is_zero() {
572        return 0;
573    };
574    if is_quad_residue(a, p).unwrap_or_default() {
575        1
576    } else {
577        -1
578    }
579}
580
581// Ported from sympy implementation
582// Simplified as a & p are nonnegative
583// Asumes p is a prime number
584pub(crate) fn is_quad_residue(a: &BigUint, p: &BigUint) -> Result<bool, MathError> {
585    if p.is_zero() {
586        return Err(MathError::IsQuadResidueZeroPrime);
587    }
588    let a = if a >= p { a.mod_floor(p) } else { a.clone() };
589    if a < BigUint::from(2_u8) || p < &BigUint::from(3_u8) {
590        return Ok(true);
591    }
592    Ok(
593        a.modpow(&(p - BigUint::one()).div_floor(&BigUint::from(2_u8)), p)
594            .is_one(),
595    )
596}
597
598#[cfg(test)]
599mod tests {
600    use super::*;
601    use crate::utils::test_utils::*;
602    use crate::utils::CAIRO_PRIME;
603    use assert_matches::assert_matches;
604
605    use num_traits::Num;
606    use rand::Rng;
607
608    #[cfg(feature = "std")]
609    use num_prime::RandPrime;
610
611    #[cfg(feature = "std")]
612    use proptest::prelude::*;
613
614    // Only used in proptest for now
615    #[cfg(feature = "std")]
616    use num_bigint::Sign;
617
618    #[cfg(target_arch = "wasm32")]
619    use wasm_bindgen_test::*;
620
621    #[test]
622    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
623    fn calculate_divmod_a() {
624        let a = bigint_str!(
625            "11260647941622813594563746375280766662237311019551239924981511729608487775604310196863705127454617186486639011517352066501847110680463498585797912894788"
626        );
627        let b = bigint_str!(
628            "4020711254448367604954374443741161860304516084891705811279711044808359405970"
629        );
630        assert_eq!(
631            bigint_str!(
632                "2904750555256547440469454488220756360634457312540595732507835416669695939476"
633            ),
634            div_mod(
635                &a,
636                &b,
637                &BigInt::from_str_radix(&crate::utils::PRIME_STR[2..], 16)
638                    .expect("Couldn't parse prime")
639            )
640            .unwrap()
641        );
642    }
643
644    #[test]
645    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
646    fn calculate_divmod_b() {
647        let a = bigint_str!(
648            "29642372811668969595956851264770043260610851505766181624574941701711520154703788233010819515917136995474951116158286220089597404329949295479559895970988"
649        );
650        let b = bigint_str!(
651            "3443173965374276972000139705137775968422921151703548011275075734291405722262"
652        );
653        assert_eq!(
654            bigint_str!(
655                "3601388548860259779932034493250169083811722919049731683411013070523752439691"
656            ),
657            div_mod(
658                &a,
659                &b,
660                &BigInt::from_str_radix(&crate::utils::PRIME_STR[2..], 16)
661                    .expect("Couldn't parse prime")
662            )
663            .unwrap()
664        );
665    }
666
667    #[test]
668    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
669    fn calculate_divmod_c() {
670        let a = bigint_str!(
671            "1208267356464811040667664150251401430616174694388968865551115897173431833224432165394286799069453655049199580362994484548890574931604445970825506916876"
672        );
673        let b = bigint_str!(
674            "1809792356889571967986805709823554331258072667897598829955472663737669990418"
675        );
676        assert_eq!(
677            bigint_str!(
678                "1545825591488572374291664030703937603499513742109806697511239542787093258962"
679            ),
680            div_mod(
681                &a,
682                &b,
683                &BigInt::from_str_radix(&crate::utils::PRIME_STR[2..], 16)
684                    .expect("Couldn't parse prime")
685            )
686            .unwrap()
687        );
688    }
689
690    #[test]
691    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
692    fn compute_safe_div() {
693        let x = Felt252::from(26);
694        let y = Felt252::from(13);
695        assert_matches!(safe_div(&x, &y), Ok(i) if i == Felt252::from(2));
696    }
697
698    #[test]
699    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
700    fn compute_safe_div_non_divisor() {
701        let x = Felt252::from(25);
702        let y = Felt252::from(4);
703        let result = safe_div(&x, &y);
704        assert_matches!(
705            result,
706            Err(MathError::SafeDivFail(bx)) if *bx == (Felt252::from(25), Felt252::from(4)));
707    }
708
709    #[test]
710    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
711    fn compute_safe_div_by_zero() {
712        let x = Felt252::from(25);
713        let y = Felt252::ZERO;
714        let result = safe_div(&x, &y);
715        assert_matches!(result, Err(MathError::DividedByZero));
716    }
717
718    #[test]
719    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
720    fn compute_safe_div_usize() {
721        assert_matches!(safe_div_usize(26, 13), Ok(2));
722    }
723
724    #[test]
725    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
726    fn compute_safe_div_usize_non_divisor() {
727        assert_matches!(
728            safe_div_usize(25, 4),
729            Err(MathError::SafeDivFailUsize(bx)) if *bx == (25, 4)
730        );
731    }
732
733    #[test]
734    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
735    fn compute_safe_div_usize_by_zero() {
736        assert_matches!(safe_div_usize(25, 0), Err(MathError::DividedByZero));
737    }
738
739    #[test]
740    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
741    fn compute_line_slope_for_valid_points() {
742        let point_a = (
743            bigint_str!(
744                "3139037544796708144595053687182055617920475701120786241351436619796497072089"
745            ),
746            bigint_str!(
747                "2119589567875935397690285099786081818522144748339117565577200220779667999801"
748            ),
749        );
750        let point_b = (
751            bigint_str!(
752                "3324833730090626974525872402899302150520188025637965566623476530814354734325"
753            ),
754            bigint_str!(
755                "3147007486456030910661996439995670279305852583596209647900952752170983517249"
756            ),
757        );
758        let prime = (*CAIRO_PRIME).clone().into();
759        assert_eq!(
760            bigint_str!(
761                "992545364708437554384321881954558327331693627531977596999212637460266617010"
762            ),
763            line_slope(&point_a, &point_b, &prime).unwrap()
764        );
765    }
766
767    #[test]
768    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
769    fn compute_double_slope_for_valid_point_a() {
770        let point = (
771            bigint_str!(
772                "3143372541908290873737380228370996772020829254218248561772745122290262847573"
773            ),
774            bigint_str!(
775                "1721586982687138486000069852568887984211460575851774005637537867145702861131"
776            ),
777        );
778        let prime = (*CAIRO_PRIME).clone().into();
779        let alpha = bigint!(1);
780        assert_eq!(
781            bigint_str!(
782                "3601388548860259779932034493250169083811722919049731683411013070523752439691"
783            ),
784            ec_double_slope(&point, &alpha, &prime).unwrap()
785        );
786    }
787
788    #[test]
789    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
790    fn compute_double_slope_for_valid_point_b() {
791        let point = (
792            bigint_str!(
793                "1937407885261715145522756206040455121546447384489085099828343908348117672673"
794            ),
795            bigint_str!(
796                "2010355627224183802477187221870580930152258042445852905639855522404179702985"
797            ),
798        );
799        let prime = (*CAIRO_PRIME).clone().into();
800        let alpha = bigint!(1);
801        assert_eq!(
802            bigint_str!(
803                "2904750555256547440469454488220756360634457312540595732507835416669695939476"
804            ),
805            ec_double_slope(&point, &alpha, &prime).unwrap()
806        );
807    }
808
809    #[test]
810    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
811    fn calculate_ec_double_for_valid_point_a() {
812        let point = (
813            bigint_str!(
814                "1937407885261715145522756206040455121546447384489085099828343908348117672673"
815            ),
816            bigint_str!(
817                "2010355627224183802477187221870580930152258042445852905639855522404179702985"
818            ),
819        );
820        let prime = (*CAIRO_PRIME).clone().into();
821        let alpha = bigint!(1);
822        assert_eq!(
823            (
824                bigint_str!(
825                    "58460926014232092148191979591712815229424797874927791614218178721848875644"
826                ),
827                bigint_str!(
828                    "1065613861227134732854284722490492186040898336012372352512913425790457998694"
829                )
830            ),
831            ec_double(point, &alpha, &prime).unwrap()
832        );
833    }
834
835    #[test]
836    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
837    fn calculate_ec_double_for_valid_point_b() {
838        let point = (
839            bigint_str!(
840                "3143372541908290873737380228370996772020829254218248561772745122290262847573"
841            ),
842            bigint_str!(
843                "1721586982687138486000069852568887984211460575851774005637537867145702861131"
844            ),
845        );
846        let prime = (*CAIRO_PRIME).clone().into();
847        let alpha = bigint!(1);
848        assert_eq!(
849            (
850                bigint_str!(
851                    "1937407885261715145522756206040455121546447384489085099828343908348117672673"
852                ),
853                bigint_str!(
854                    "2010355627224183802477187221870580930152258042445852905639855522404179702985"
855                )
856            ),
857            ec_double(point, &alpha, &prime).unwrap()
858        );
859    }
860
861    #[test]
862    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
863    fn calculate_ec_double_for_valid_point_c() {
864        let point = (
865            bigint_str!(
866                "634630432210960355305430036410971013200846091773294855689580772209984122075"
867            ),
868            bigint_str!(
869                "904896178444785983993402854911777165629036333948799414977736331868834995209"
870            ),
871        );
872        let prime = (*CAIRO_PRIME).clone().into();
873        let alpha = bigint!(1);
874        assert_eq!(
875            (
876                bigint_str!(
877                    "3143372541908290873737380228370996772020829254218248561772745122290262847573"
878                ),
879                bigint_str!(
880                    "1721586982687138486000069852568887984211460575851774005637537867145702861131"
881                )
882            ),
883            ec_double(point, &alpha, &prime).unwrap()
884        );
885    }
886
887    #[test]
888    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
889    fn calculate_ec_add_for_valid_points_a() {
890        let point_a = (
891            bigint_str!(
892                "1183418161532233795704555250127335895546712857142554564893196731153957537489"
893            ),
894            bigint_str!(
895                "1938007580204102038458825306058547644691739966277761828724036384003180924526"
896            ),
897        );
898        let point_b = (
899            bigint_str!(
900                "1977703130303461992863803129734853218488251484396280000763960303272760326570"
901            ),
902            bigint_str!(
903                "2565191853811572867032277464238286011368568368717965689023024980325333517459"
904            ),
905        );
906        let prime = (*CAIRO_PRIME).clone().into();
907        assert_eq!(
908            (
909                bigint_str!(
910                    "1977874238339000383330315148209250828062304908491266318460063803060754089297"
911                ),
912                bigint_str!(
913                    "2969386888251099938335087541720168257053975603483053253007176033556822156706"
914                )
915            ),
916            ec_add(point_a, point_b, &prime).unwrap()
917        );
918    }
919
920    #[test]
921    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
922    fn calculate_ec_add_for_valid_points_b() {
923        let point_a = (
924            bigint_str!(
925                "3139037544796708144595053687182055617920475701120786241351436619796497072089"
926            ),
927            bigint_str!(
928                "2119589567875935397690285099786081818522144748339117565577200220779667999801"
929            ),
930        );
931        let point_b = (
932            bigint_str!(
933                "3324833730090626974525872402899302150520188025637965566623476530814354734325"
934            ),
935            bigint_str!(
936                "3147007486456030910661996439995670279305852583596209647900952752170983517249"
937            ),
938        );
939        let prime = (*CAIRO_PRIME).clone().into();
940        assert_eq!(
941            (
942                bigint_str!(
943                    "1183418161532233795704555250127335895546712857142554564893196731153957537489"
944                ),
945                bigint_str!(
946                    "1938007580204102038458825306058547644691739966277761828724036384003180924526"
947                )
948            ),
949            ec_add(point_a, point_b, &prime).unwrap()
950        );
951    }
952
953    #[test]
954    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
955    fn calculate_ec_add_for_valid_points_c() {
956        let point_a = (
957            bigint_str!(
958                "1183418161532233795704555250127335895546712857142554564893196731153957537489"
959            ),
960            bigint_str!(
961                "1938007580204102038458825306058547644691739966277761828724036384003180924526"
962            ),
963        );
964        let point_b = (
965            bigint_str!(
966                "1977703130303461992863803129734853218488251484396280000763960303272760326570"
967            ),
968            bigint_str!(
969                "2565191853811572867032277464238286011368568368717965689023024980325333517459"
970            ),
971        );
972        let prime = (*CAIRO_PRIME).clone().into();
973        assert_eq!(
974            (
975                bigint_str!(
976                    "1977874238339000383330315148209250828062304908491266318460063803060754089297"
977                ),
978                bigint_str!(
979                    "2969386888251099938335087541720168257053975603483053253007176033556822156706"
980                )
981            ),
982            ec_add(point_a, point_b, &prime).unwrap()
983        );
984    }
985
986    #[test]
987    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
988    fn calculate_isqrt_a() {
989        let n = biguint!(81);
990        assert_matches!(isqrt(&n), Ok(x) if x == biguint!(9));
991    }
992
993    #[test]
994    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
995    fn calculate_isqrt_b() {
996        let n = biguint_str!("4573659632505831259480");
997        assert_matches!(isqrt(&BigUint::pow(&n, 2_u32)), Ok(num) if num == n);
998    }
999
1000    #[test]
1001    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1002    fn calculate_isqrt_c() {
1003        let n = biguint_str!(
1004            "3618502788666131213697322783095070105623107215331596699973092056135872020481"
1005        );
1006        assert_matches!(isqrt(&BigUint::pow(&n, 2_u32)), Ok(inner) if inner == n);
1007    }
1008
1009    #[test]
1010    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1011    fn calculate_isqrt_zero() {
1012        let n = BigUint::zero();
1013        assert_matches!(isqrt(&n), Ok(inner) if inner.is_zero());
1014    }
1015
1016    #[test]
1017    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1018    fn safe_div_bigint_by_zero() {
1019        let x = BigInt::one();
1020        let y = BigInt::zero();
1021        assert_matches!(safe_div_bigint(&x, &y), Err(MathError::DividedByZero))
1022    }
1023
1024    #[test]
1025    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1026    fn test_sqrt_prime_power() {
1027        let n: BigUint = 25_u32.into();
1028        let p: BigUint = 18446744069414584321_u128.into();
1029        assert_eq!(sqrt_prime_power(&n, &p), Some(5_u32.into()));
1030    }
1031
1032    #[test]
1033    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1034    fn test_sqrt_prime_power_p_is_zero() {
1035        let n = BigUint::one();
1036        let p: BigUint = BigUint::zero();
1037        assert_eq!(sqrt_prime_power(&n, &p), None);
1038    }
1039
1040    #[test]
1041    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1042    fn test_sqrt_prime_power_non_prime() {
1043        let p: BigUint = BigUint::from_bytes_be(&[
1044            69, 15, 232, 82, 215, 167, 38, 143, 173, 94, 133, 111, 1, 2, 182, 229, 110, 113, 76, 0,
1045            47, 110, 148, 109, 6, 133, 27, 190, 158, 197, 168, 219, 165, 254, 81, 53, 25, 34,
1046        ]);
1047        let n = BigUint::from_bytes_be(&[
1048            9, 13, 22, 191, 87, 62, 157, 83, 157, 85, 93, 105, 230, 187, 32, 101, 51, 181, 49, 202,
1049            203, 195, 76, 193, 149, 78, 109, 146, 240, 126, 182, 115, 161, 238, 30, 118, 157, 252,
1050        ]);
1051
1052        assert_eq!(sqrt_prime_power(&n, &p), None);
1053    }
1054
1055    #[test]
1056    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1057    fn test_sqrt_prime_power_none() {
1058        let n: BigUint = 10_u32.into();
1059        let p: BigUint = 602_u32.into();
1060        assert_eq!(sqrt_prime_power(&n, &p), None);
1061    }
1062
1063    #[test]
1064    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1065    fn test_sqrt_prime_power_prime_two() {
1066        let n: BigUint = 25_u32.into();
1067        let p: BigUint = 2_u32.into();
1068        assert_eq!(sqrt_prime_power(&n, &p), Some(BigUint::one()));
1069    }
1070
1071    #[test]
1072    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1073    fn test_sqrt_prime_power_prime_mod_8_is_5_sign_not_one() {
1074        let n: BigUint = 676_u32.into();
1075        let p: BigUint = 9956234341095173_u64.into();
1076        assert_eq!(
1077            sqrt_prime_power(&n, &p),
1078            Some(BigUint::from(9956234341095147_u64))
1079        );
1080    }
1081
1082    #[test]
1083    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1084    fn test_sqrt_prime_power_prime_mod_8_is_5_sign_is_one() {
1085        let n: BigUint = 130283432663_u64.into();
1086        let p: BigUint = 743900351477_u64.into();
1087        assert_eq!(
1088            sqrt_prime_power(&n, &p),
1089            Some(BigUint::from(123538694848_u64))
1090        );
1091    }
1092
1093    #[test]
1094    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1095    fn test_legendre_symbol_zero() {
1096        assert!(legendre_symbol(&BigUint::zero(), &BigUint::one()).is_zero())
1097    }
1098
1099    #[test]
1100    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1101    fn test_is_quad_residue_prime_zero() {
1102        assert_eq!(
1103            is_quad_residue(&BigUint::one(), &BigUint::zero()),
1104            Err(MathError::IsQuadResidueZeroPrime)
1105        )
1106    }
1107
1108    #[test]
1109    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1110    fn test_is_quad_residue_prime_a_one_true() {
1111        assert_eq!(is_quad_residue(&BigUint::one(), &BigUint::one()), Ok(true))
1112    }
1113
1114    #[test]
1115    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1116    fn mul_inv_0_is_0() {
1117        let p = &(*CAIRO_PRIME).clone().into();
1118        let x = &BigInt::zero();
1119        let x_inv = mul_inv(x, p);
1120
1121        assert_eq!(x_inv, BigInt::zero());
1122    }
1123
1124    #[test]
1125    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1126    fn igcdex_1_1() {
1127        assert_eq!(
1128            igcdex(&BigInt::one(), &BigInt::one()),
1129            (BigInt::zero(), BigInt::one(), BigInt::one())
1130        )
1131    }
1132
1133    #[test]
1134    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1135    fn igcdex_0_0() {
1136        assert_eq!(
1137            igcdex(&BigInt::zero(), &BigInt::zero()),
1138            (BigInt::zero(), BigInt::one(), BigInt::zero())
1139        )
1140    }
1141
1142    #[test]
1143    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1144    fn igcdex_1_0() {
1145        assert_eq!(
1146            igcdex(&BigInt::one(), &BigInt::zero()),
1147            (BigInt::one(), BigInt::zero(), BigInt::one())
1148        )
1149    }
1150
1151    #[test]
1152    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1153    fn igcdex_4_6() {
1154        assert_eq!(
1155            igcdex(&BigInt::from(4), &BigInt::from(6)),
1156            (BigInt::from(-1), BigInt::one(), BigInt::from(2))
1157        )
1158    }
1159
1160    #[test]
1161    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1162    fn qm31_packed_reduced_read_coordinates_over_144_bits() {
1163        let mut felt_bytes = [0u8; 32];
1164        felt_bytes[18] = 1;
1165        let felt = Felt252::from_bytes_le(&felt_bytes);
1166        assert_matches!(
1167            qm31_packed_reduced_read_coordinates(felt),
1168            Err(MathError::QM31UnreducedError(bx)) if *bx == felt
1169        );
1170    }
1171
1172    #[test]
1173    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1174    fn qm31_packed_reduced_read_coordinates_unreduced() {
1175        let mut felt_bytes = [0u8; 32];
1176        felt_bytes[0] = 0xff;
1177        felt_bytes[1] = 0xff;
1178        felt_bytes[2] = 0xff;
1179        felt_bytes[3] = (1 << 7) - 1;
1180        let felt = Felt252::from_bytes_le(&felt_bytes);
1181        assert_matches!(
1182            qm31_packed_reduced_read_coordinates(felt),
1183            Err(MathError::QM31UnreducedError(bx)) if *bx == felt
1184        );
1185    }
1186
1187    #[test]
1188    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1189    fn test_qm31_packed_reduced_add() {
1190        let x_coordinates = [1414213562, 1732050807, 1618033988, 1234567890];
1191        let y_coordinates = [1234567890, 1414213562, 1732050807, 1618033988];
1192        let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1193        let y = qm31_coordinates_to_packed_reduced(y_coordinates);
1194        let res = qm31_packed_reduced_add(x, y).unwrap();
1195        let res_coordinates = qm31_packed_reduced_read_coordinates(res);
1196        assert_eq!(
1197            res_coordinates,
1198            Ok([
1199                (1414213562 + 1234567890) % STWO_PRIME,
1200                (1732050807 + 1414213562) % STWO_PRIME,
1201                (1618033988 + 1732050807) % STWO_PRIME,
1202                (1234567890 + 1618033988) % STWO_PRIME,
1203            ])
1204        );
1205    }
1206
1207    #[test]
1208    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1209    fn test_qm31_packed_reduced_neg() {
1210        let x_coordinates = [1749652895, 834624081, 1930174752, 2063872165];
1211        let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1212        let res = qm31_packed_reduced_neg(x).unwrap();
1213        let res_coordinates = qm31_packed_reduced_read_coordinates(res);
1214        assert_eq!(
1215            res_coordinates,
1216            Ok([
1217                STWO_PRIME - x_coordinates[0],
1218                STWO_PRIME - x_coordinates[1],
1219                STWO_PRIME - x_coordinates[2],
1220                STWO_PRIME - x_coordinates[3]
1221            ])
1222        );
1223    }
1224
1225    #[test]
1226    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1227    fn test_qm31_packed_reduced_sub() {
1228        let x_coordinates = [
1229            (1414213562 + 1234567890) % STWO_PRIME,
1230            (1732050807 + 1414213562) % STWO_PRIME,
1231            (1618033988 + 1732050807) % STWO_PRIME,
1232            (1234567890 + 1618033988) % STWO_PRIME,
1233        ];
1234        let y_coordinates = [1414213562, 1732050807, 1618033988, 1234567890];
1235        let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1236        let y = qm31_coordinates_to_packed_reduced(y_coordinates);
1237        let res = qm31_packed_reduced_sub(x, y).unwrap();
1238        let res_coordinates = qm31_packed_reduced_read_coordinates(res);
1239        assert_eq!(
1240            res_coordinates,
1241            Ok([1234567890, 1414213562, 1732050807, 1618033988])
1242        );
1243    }
1244
1245    #[test]
1246    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1247    fn test_qm31_packed_reduced_mul() {
1248        let x_coordinates = [1414213562, 1732050807, 1618033988, 1234567890];
1249        let y_coordinates = [1259921049, 1442249570, 1847759065, 2094551481];
1250        let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1251        let y = qm31_coordinates_to_packed_reduced(y_coordinates);
1252        let res = qm31_packed_reduced_mul(x, y).unwrap();
1253        let res_coordinates = qm31_packed_reduced_read_coordinates(res);
1254        assert_eq!(
1255            res_coordinates,
1256            Ok([947980980, 1510986506, 623360030, 1260310989])
1257        );
1258    }
1259
1260    #[test]
1261    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1262    fn test_qm31_packed_reduced_inv() {
1263        let x_coordinates = [1259921049, 1442249570, 1847759065, 2094551481];
1264        let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1265        let res = qm31_packed_reduced_inv(x).unwrap();
1266        assert_eq!(qm31_packed_reduced_mul(x, res), Ok(Felt252::from(1)));
1267
1268        let x_coordinates = [1, 2, 3, 4];
1269        let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1270        let res = qm31_packed_reduced_inv(x).unwrap();
1271        assert_eq!(qm31_packed_reduced_mul(x, res), Ok(Felt252::from(1)));
1272
1273        let x_coordinates = [1749652895, 834624081, 1930174752, 2063872165];
1274        let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1275        let res = qm31_packed_reduced_inv(x).unwrap();
1276        assert_eq!(qm31_packed_reduced_mul(x, res), Ok(Felt252::from(1)));
1277    }
1278
1279    // TODO: Refactor using proptest and separating particular cases
1280    #[test]
1281    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1282    fn test_qm31_packed_reduced_inv_extensive() {
1283        let mut rng = SmallRng::seed_from_u64(11480028852697973135);
1284        #[derive(Clone, Copy)]
1285        enum Configuration {
1286            Zero,
1287            One,
1288            MinusOne,
1289            Random,
1290        }
1291        let configurations = [
1292            Configuration::Zero,
1293            Configuration::One,
1294            Configuration::MinusOne,
1295            Configuration::Random,
1296        ];
1297        let mut cartesian_product = vec![];
1298        for &a in &configurations {
1299            for &b in &configurations {
1300                for &c in &configurations {
1301                    for &d in &configurations {
1302                        cartesian_product.push([a, b, c, d]);
1303                    }
1304                }
1305            }
1306        }
1307
1308        for test_case in cartesian_product {
1309            let x_coordinates: [u64; 4] = test_case
1310                .iter()
1311                .map(|&x| match x {
1312                    Configuration::Zero => 0,
1313                    Configuration::One => 1,
1314                    Configuration::MinusOne => STWO_PRIME - 1,
1315                    Configuration::Random => rng.gen_range(0..STWO_PRIME),
1316                })
1317                .collect::<Vec<u64>>()
1318                .try_into()
1319                .unwrap();
1320            if x_coordinates == [0, 0, 0, 0] {
1321                continue;
1322            }
1323            let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1324            let res = qm31_packed_reduced_inv(x).unwrap();
1325            assert_eq!(qm31_packed_reduced_mul(x, res), Ok(Felt252::from(1)));
1326        }
1327    }
1328
1329    #[test]
1330    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test)]
1331    fn test_qm31_packed_reduced_div() {
1332        let x_coordinates = [1259921049, 1442249570, 1847759065, 2094551481];
1333        let y_coordinates = [1414213562, 1732050807, 1618033988, 1234567890];
1334        let xy_coordinates = [947980980, 1510986506, 623360030, 1260310989];
1335        let x = qm31_coordinates_to_packed_reduced(x_coordinates);
1336        let y = qm31_coordinates_to_packed_reduced(y_coordinates);
1337        let xy = qm31_coordinates_to_packed_reduced(xy_coordinates);
1338
1339        let res = qm31_packed_reduced_div(xy, y).unwrap();
1340        assert_eq!(res, x);
1341
1342        let res = qm31_packed_reduced_div(xy, x).unwrap();
1343        assert_eq!(res, y);
1344    }
1345
1346    #[cfg(feature = "std")]
1347    proptest! {
1348        #[test]
1349        fn pow2_const_in_range_returns_power_of_2(x in 0..=251u32) {
1350            prop_assert_eq!(pow2_const(x), Felt252::TWO.pow(x));
1351        }
1352
1353        #[test]
1354        fn pow2_const_oob_returns_1(x in 252u32..) {
1355            prop_assert_eq!(pow2_const(x), Felt252::ONE);
1356        }
1357
1358        #[test]
1359        fn pow2_const_nz_in_range_returns_power_of_2(x in 0..=251u32) {
1360            prop_assert_eq!(Felt252::from(pow2_const_nz(x)), Felt252::TWO.pow(x));
1361        }
1362
1363        #[test]
1364        fn pow2_const_nz_oob_returns_1(x in 252u32..) {
1365            prop_assert_eq!(Felt252::from(pow2_const_nz(x)), Felt252::ONE);
1366        }
1367
1368        #[test]
1369        // Test for sqrt_prime_power_ of a quadratic residue. Result should be the minimum root.
1370        fn sqrt_prime_power_using_random_prime(ref x in any::<[u8; 38]>(), ref y in any::<u64>()) {
1371            let mut rng = SmallRng::seed_from_u64(*y);
1372            let x = &BigUint::from_bytes_be(x);
1373            // Generate a prime here instead of relying on y, otherwise y may never be a prime number
1374            let p : &BigUint = &RandPrime::gen_prime(&mut rng, 384,  None);
1375            let x_sq = x * x;
1376            if let Some(sqrt) = sqrt_prime_power(&x_sq, p) {
1377                if &sqrt != x {
1378                    prop_assert_eq!(&(p - sqrt), x);
1379                } else {
1380                prop_assert_eq!(&sqrt, x);
1381                }
1382            }
1383        }
1384
1385        #[test]
1386        fn mul_inv_x_by_x_is_1(ref x in any::<[u8; 32]>()) {
1387            let p = &(*CAIRO_PRIME).clone().into();
1388            let pos_x = &BigInt::from_bytes_be(Sign::Plus, x);
1389            let neg_x = &BigInt::from_bytes_be(Sign::Minus, x);
1390            let pos_x_inv = mul_inv(pos_x, p);
1391            let neg_x_inv = mul_inv(neg_x, p);
1392
1393            prop_assert_eq!((pos_x * pos_x_inv).mod_floor(p), BigInt::one());
1394            prop_assert_eq!((neg_x * neg_x_inv).mod_floor(p), BigInt::one());
1395        }
1396    }
1397}
cairo_vm/math_utils/mod.rs

cairo_vm/math_utils/
mod.rs
