Jacobi symbol computation error

[andrewwhitehead]

For some inputs with >= 4 limbs, the result of Jacobi symbol computation may be incorrect. This seems to be an inherent issue with using the optimized binary GCD for this purpose. When the top limbs of the compacted representation match but the full-precision inputs do not, the numerator may become negative mid-batch computation. The terms are made positive after the batch completes, but there is no good way to correct the Jacobi symbol sign at this point.

There aren't a lot of constant-time alternatives here besides just using classic binary GCD (slow for large inputs). At the moment the variable-time implementation has the same issue, but that's easier to switch to another method, GMP uses Lehmer's algorithm I believe.

#[test]
fn jacobi_mismatch() {
    use crate::{Odd, U64, U128, U256, Uint};

    fn check<const LIMBS: usize>(a: Uint<LIMBS>, b: Odd<Uint<LIMBS>>) {
        let (gcd_1, jacobi_neg_1) = b.classic_bingcd_(&a);
        let (gcd_2, jacobi_neg_2) = b
            .optimized_bingcd_::<{ U64::BITS }, { U64::LIMBS }, { U128::LIMBS }>(
                &a,
                U64::BITS - 2,
            );
        assert_eq!(gcd_1, gcd_2);
        assert_eq!(jacobi_neg_1, jacobi_neg_2);
    }

    check(
        U256::from_be_hex("0000000000000002108DEAFCB180F023912BEED0186CEEAD593A8507B7DA4E9B"),
        Odd::<U256>::from_be_hex(
            "0000000000000002108DEAFCB180F023912BEED0186CEEED593A8507B7DA4E9B",
        ),
    );
}

[tarcieri]

I think crypto-primes already has a vartime Jacobi symbol implementation. Maybe we could upstream it? Otherwise Lehmer's algorithm looks interesting /cc @fjarri

What are the actual use cases of a constant-time Jacobi symbol calculation? Could the compacted representation somehow be expanded to have enough information to disambiguate these cases without it negatively impacting other usages?

[andrewwhitehead]

I believe the constant-time Jacobi symbol is used in EC point decompression for some curves to choose between which sqrt to evaluate, especially if it's faster than the modular sqrt operation.

I don't think you can just expand the compacted representation, you'd need to look at all the intermediate bits and then it's just classic bingcd (which does work).

There are a few options for a vartime version, including the standard Euclidean method and Sorenson's k-ary method (see "Efficient Algorithms for Computing the Jacobi Symbol"). For constant-time it would be nice, and faster, to use safegcd but that method has potential patent issues (although no application seems to have been filed). There's also Aranha's paper which extends that method to the Kronecker symbol, but it has the same patent issue once you get into the optimized batching version.

Incidentally I've already rewritten the GCD implementations to be much faster (including vartime xgcd and BoxedUint xgcd), but that seems to be blocked on this issue.

[fjarri]

The implementation in crypto-primes assumes one of the arguments is limb-sized, but I guess with no constant-time requirements it's easy to extend to both arguments being large.

[kayabaNerve]

Using the classic binary GCD, while unfortunate, is great simply as a fix as there's already the code implementing it. For the special-case one number is prime (the Legendre symbol), Euler's criterion is a single modpow and sufficiently simple/efficient (even if suboptimal). I know that isn't the most helpful comment, I just wanted to raise the anecdote in case it ends up being helpful for the general case or justifying for its own dedicated method (though that'd require a Prime as we have an Odd, which likely isn't feasible...).