secp256k1 #979

jacobi_divsteps.py

N = 62

def count_trailing_zeros(x):
    assert x != 0
    ans = 0
    while x & 1 == 0:
        ans += 1
        x >>= 1
    return ans

def update_fg(f, g, t):
    u, v, q, r = t
    cf, cg = u * f + v * g, q * f + r * g
    assert cf % 2**N == 0
    assert cg % 2**N == 0
    return cf >> N, cg >> N

def divsteps_N_matrix(eta, f, g, jac):
    u, v, q, r = 1, 0, 0, 1
    i = N
    while True:
        z = min(i, count_trailing_zeros(g))
        eta, g, u, v = eta - z, g >> z, u << z, v << z
        i -= z
        if z & 1 and (f % 8 == 3 or f % 8 == 5):
            jac = -jac
        if i == 0:
            break
        assert (g & 1) == 1
        if eta < 0:
            eta, f, g, u, v, q, r = -eta, g, f, q, r, u, v
            if f % 4 == 3 and g % 4 == 3:
                jac = -jac
            limit = min(i, eta + 1, 6)
            assert limit > 0 and limit <= N
            m = (1 << limit) - 1
            w = (f * g * (f * f - 2)) & m
        else:
            limit = min(i, eta + 1, 4)
            assert limit > 0 and limit <= N
            m = (1 << limit) - 1
            w = f + (((f + 1) & 4) << 1)
            w = (-w * g) & m
        g, q, r = g + f * w, q + u * w, r + v * w
        assert g % (2**limit) == 0
    return eta, (u, v, q, r), jac

def fastjac_divsteps(x, M):
    assert x > 0 and M & 1
    jac = 1
    eta, f, g = -1, M, x
    while f != 1:
        assert f & 1
        eta, t, jac = divsteps_N_matrix(eta, f, g, jac)
        f, g = update_fg(f, g, t)
    return jac

def slowjac(a, n):
    jac = 1
    while True:
        assert n & 1 == 1
        if a == 1:
            return jac
        if a == n - 1:
            if n & 3 == 3:
                jac = -jac
            return jac
        while a & 1 == 0:
            if n & 7 in (3, 5):
                jac = -jac
            a >>= 1
        if a & 3 == 3 and n & 3 == 3:
            jac = -jac
        a, n = n % a, a
    return jac

if __name__ == '__main__':
    from math import gcd
    from random import randint
    for i in range(10000):
        while True:
            a = randint(1, 10**6)
            M = 2 * randint(1, 10**6) + 1
            if M > 1 and gcd(a, M) == 1:
                break
        assert fastjac_divsteps(a, M) == slowjac(a, M)