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

cairo_vm/math_utils/
mod.rs
