Making Jacobi symbol function faster

General Usage Performance
11

The function produced errors:

Failed: jacobi(2, 5) = 1, expected -1
Failed: jacobi(7, 9) = -1, expected 1

which I did not know how to make AI to fix for me. I managed Claude to spit out a following code for BigInt:

function jacobi(n::BigInt, k::BigInt)

    # Get the last limb (word) of a BigInt
    get_limb(n::BigInt, i::Int) = unsafe_load(n.d, 1 + i)

    iseven(k) && throw("Argument k=$k should be odd.")
    n = mod(n, k)
    t = 1
    k_1 = true  # k is odd by requirement
    last_limb = get_limb(k, 0)
    k_2 = (last_limb & 2 == 2)  # second bit
    k_3 = (last_limb & 4 == 4)  # third bit
    
    while !iszero(n)
        z = trailing_zeros(n)
        n >>= z
        # After right shift, check bits directly on the last limb
        last_limb = get_limb(n, 0)
        n_2 = (last_limb & 2 == 2)  # second bit
        
        # For k ≡ 3,5 (mod 8) when z is odd
        if isodd(z) && (k_2 ⊻ k_3)
            t = -t
        end
        
        # For quadratic reciprocity when both are ≡ 3 (mod 4)
        if k_2 && n_2
            t = -t
        end
        
        # Save current n's bits for next iteration's k
        k_2 = n_2
        k_3 = (last_limb & 4 == 4)
        n, k = k, n
        n = mod(n, k)
    end
    return isone(k) ? t : 0
end

which is also around 4x faster.

2 Likes