Exponentiation (jl)


function msb(x)
        local i=0
        t=copy(x)
        while t>0 
          i+=1
          t>>=1 
        end
        i
end

function bad_pow(x,n)
    n==0 && return oneunit(x)
    n<0  && return bad_pow(inv(x),-n)
    x*bad_pow(x,n-1)
end

function better_pow(x,n)
    c=oneunit(x)
    n==0 && return c
    n<0  && return better_pow(inv(x),-n)
    for i in msb(n):-1:0
        c *= c
        (((n>>i) & 1)==1) && ( c*= x)
    end
    c
end

function my_gcd(x,y)
    x>y || return my_gcd(y,x)
    y==0 && return x
    @show (x,y)
    my_gcd(y,mod(x,y))
end

function my_xgcd(x,y)
    local U=[oneunit(x),zero(x),x]
    local V=[zero(y),oneunit(y),y]
    while V[3]!=zero(x)
        q=U[3]÷V[3]
        U,V=V,U-V*q
        @show U,V
    end
    U
end

function SPprime(n,a)
    t,s = n-1,0
    while (t%2==0)
        s +=  1
        t >>= 1
    end
    b=powermod(a,t,n)
    ((b==1) | (b==n-1)) && return true
    for j in 1:s-1
        b=powermod(b,2,n)
        b==n-1 && return true
    end
    false
end

function Miller(n)
    W = Integer(min(round(2*log(n)^2),n-1))
    for a in 2:W
        SPprime(n,a) || return false
    end
    true
end

function NextPrime(n)
    (n<=2)    && return 2
    (n%2==0)  && return NextPrime(n+1)
    Miller(n) && return n
    NextPrime(n+2)
end


### DEFINE Zp

import Base.+, Base.-, Base.*,Base.^,Base.inv,Base.==,Base.oneunit,Base.zero

struct ZpElement <: Number
    value::Integer
    prime::Integer
end

function oneunit(x::ZpElement)
    ZpElement(1,x.prime)
end

function zero(x::ZpElement)
    ZpElement(0,x.prime)
end

function +(x::ZpElement, y::ZpElement)
    ZpElement(mod(x.value+y.value,x.prime),x.prime)
end

function *(t::BigInt, y::ZpElement)
    ZpElement(mod(t*y.value,y.prime),y.prime)
end

function -(y::ZpElement)
   -1*y
end


function -(x::ZpElement,y::ZpElement)
   x+-y
end

function ==(x::ZpElement,y::ZpElement)
	(x.value==y.value) & (x.prime==y.prime)
end

function *(x::ZpElement,y::ZpElement)
    ZpElement(mod(x.value*y.value,x.prime),x.prime)
end


function inv(x::ZpElement)
    v=my_xgcd(x.value,x.prime)
    return ZpElement(mod(v[1],x.prime),x.prime)
end

function ^(x::ZpElement,i::BigInt)
    better_pow(x,i)
end

function CRT( a::Vector{T}, b::Vector{T} ) where T <: Integer
    ap,bp,r=my_xgcd(a[2],b[2])
    @assert r≡oneunit(r) "Must be coprime"
    [ mod(a[1]*bp*b[2]+b[1]*ap*a[2],a[2]*b[2]), a[2]*b[2] ]
end

function CRT( x::Vector{Vector{T}} ) where T <: Integer
    reduce(CRT,x)
end