RosettaCodeData/Task/Elliptic-Curve-Digital-Signature-Algorithm/Phix/elliptic-curve-digital-signature-algorithm.phix

requires(64)
enum X, Y            -- rational ec point
enum A, B, N, G, R  -- elliptic curve parameters
                    -- also signature pair(A,B)

constant mxN = 1073741789   -- maximum modulus
constant mxr = 1073807325   -- max order G = mxN + 65536
constant inf = -2147483647  -- symbolic infinity

sequence e = {0,0,0,{0,0},0}    -- single global curve
constant zerO = {inf,0}         -- point at infinity zerO

bool inverr -- impossible inverse mod N

function exgcd(atom v, u)
-- return mod(v^-1, u)
atom q, t, r = 0, s = 1
    if v<0 then v += u end if

    while v do
        q = floor(u/v)
        t = u-q*v
        u = v
        v = t
        t = r-q*s
        r = s
        s = t
    end while

    if u!=1 then
        printf(1," impossible inverse mod N, gcd = %d\n",{u})
        inverr = true
    end if

    return r
end function

function modn(atom a)
-- return mod(a, N)
    a = mod(a,e[N])
    if a<0 then a += e[N] end if
    return a
end function

function modr(atom a)
-- return mod(a, r)
    a = mod(a,e[R])
    if a<0 then a += e[R] end if
    return a
end function

function disc()
-- return the discriminant of E
    atom a = e[A], b = e[B],
         c = 4*modn(a*modn(a*a))
    return modn(-16*(c+27*modn(b*b)))
end function

function isO(sequence p)
-- return true if P = zerO
    return (p[X]=inf and p[Y]=0)
end function

function ison(sequence p)
-- return true if P is on curve E
    atom r = 0, s = 0
    if not isO(p) then
        r = modn(e[B]+p[X]*modn(e[A]+p[X]*p[X]))
        s = modn(p[Y]*p[Y])
    end if
    return (r=s)
end function

procedure pprint(string f, sequence p)
-- print point P with prefix f
    if isO(p) then
        printf(1,"%s (0)\n",{f})
    else
        atom y = p[Y]
        if y>e[N]-y then y -= e[N] end if
        printf(1,"%s (%d,%d)\n",{f,p[X],y})
    end if
end procedure

function padd(sequence p, q)
-- full ec point addition
atom la, t

    if isO(p) then return q end if
    if isO(q) then return p end if

    if p[X]!=q[X] then --                   R := P + Q
        t = p[Y]-q[Y]
        la = modn(t*exgcd(p[X]-q[X], e[N]))

    else --                                 P = Q, R := 2P
        if (p[Y]=q[Y]) and (p[Y]!=0) then
            t = modn(3*modn(p[X]*p[X])+e[A])
            la = modn(t*exgcd(2*p[Y], e[N]))
        else
            return zerO --                  P = -Q, R := O
        end if
    end if

    t = modn(la*la-p[X]-q[X])
    sequence r = deep_copy(zerO)
    r[Y] = modn(la*(p[X]-t)-p[Y])
    r[X] = t
    if inverr then r = zerO end if
    return r
end function

function pmul(sequence p, atom k)
-- R:= multiple kP
    sequence s = zerO, q = p

    while k do
        if and_bits(k,1) then
            s = padd(s, q)
        end if
        if inverr then s = zerO; exit end if
        k = floor(k/2)
        q = padd(q, q)
    end while
    return s
end function

function ellinit(sequence i)
-- initialize elliptic curve
atom a = i[1], b = i[2]
    inverr = false
    e[N] = i[3]

    if (e[N]<5) or (e[N]>mxN) then return 0 end if

    e[A] = modn(a)
    e[B] = modn(b)
    e[G][X] = modn(i[4])
    e[G][Y] = modn(i[5])
    e[R] = i[6]

    if (e[R]<5) or (e[R]>mxr) then return 0 end if

    printf(1,"E: y^2 = x^3 + %dx + %d (mod %d)\n",{a,b,e[N]})
    pprint("base point G", e[G])
    printf(1,"order(G, E) = %d\n",{e[R]})

    return -1
end function

function signature(atom s, f)
-- signature primitive
atom c, d, u, u1
sequence V

    printf(1,"signature computation\n")
    while true do
        while true do
--          u = rand(e[R]-1)
            u = 571533488       -- match FreeBASIC output
--          u = 605163545       -- match C output
            V = pmul(e[G], u)
            c = modr(V[X])
            if c!=0 then exit end if
        end while

        u1 = exgcd(u, e[R])
        d = modr(u1*(f+modr(s*c)))
        if d!=0 then exit end if
    end while
    printf(1,"one-time u = %d\n",u)
    pprint("V = uG", V)
    return {c,d}
end function

function verify(sequence W, atom f, sequence sg)
-- verification primitive
atom c = sg[A], d = sg[B],
     t, c1, h1, h2, h
sequence V, V2

   --domain check
    t = (c>0) and (c<e[R])
    t = t and (d>0) and (d<e[R])
    if not t then return 0 end if

    printf(1,"\nsignature verification\n")
    h = exgcd(d, e[R])
    h1 = modr(f*h)
    h2 = modr(c*h)
    printf(1,"h1,h2 = %d,%d\n",{h1,h2})
    V = pmul(e[G], h1)
    V2 = pmul(W, h2)
    pprint("h1G", V)
    pprint("h2W", V2)
    V = padd(V, V2)
    pprint("+ =", V)
    if isO(V) then return 0 end if
    c1 = modr(V[X])
    printf(1,"c' = %d\n",c1)

    return (c1=c)
end function

procedure errmsg()
    printf(1,"invalid parameter set\n")
    printf(1,"_____________________\n")
end procedure

procedure ec_dsa(atom f, d)
-- digital signature on message hash f, error bit d
atom i, s, t
sequence sg, W

   --parameter check
    t = (disc()=0)
    t = t or isO(e[G])
    W = pmul(e[G], e[R])
    t = t or not isO(W)
    t = t or not ison(e[G])
    if t then errmsg() return end if

    puts(1,"\nkey generation\n")
--  s = rand(e[R] - 1)
    s = 509100772       -- match FreeBASIC output
--  s = 1343570         -- match C output
    W = pmul(e[G], s)
    printf(1,"private key s = %d\n",{s})
    pprint("public key W = sG", W)

   --next highest power of 2 - 1
    t = e[R]
    i = 1
    while i<32 do
        t = or_bits(t,floor(t/power(2,i)))
        i *= 2
    end while
    while f>t do
        f = floor(f/2)
    end while
    printf(1,"\naligned hash %x\n\n",{f})

    sg = signature(s, f)
    if inverr then errmsg() return end if
    printf(1,"signature c,d = %d,%d\n",sg)

    if d>0 then
        while d>t do
            d = floor(d/2)
        end while
        f = xor_bits(f,d)
        printf(1,"corrupted hash %x\n",{f})
    end if

    t = verify(W, f, sg)
    if inverr then errmsg() return end if
    if t then
        printf(1,"Valid\n_____\n")
    else
        printf(1,"invalid\n_______\n")
    end if
end procedure

--Test vectors: elliptic curve domain parameters,
--short Weierstrass model y^2 = x^3 + ax + b (mod N)

constant tests = {
--                  a,   b,  modulus N, base point G, order(G, E), cofactor
                  {355, 671, 1073741789, 13693, 10088, 1073807281},
                  {  0,   7,   67096021,  6580,   779,   16769911}, --   4
                  { -3,   1,     877073,     0,     1,     878159},
                  {  0,  14,      22651,    63,    30,        151}, -- 151
                  {  3,   2,          5,     2,     1,          5},

                    --ecdsa may fail if...
                    --the base point is of composite order
                  {  0,   7,   67096021,  2402,  6067,   33539822}, --   2
                    --the given order is a multiple of the true order
                  {  0,   7,   67096021,  6580,   779,   67079644}, --   1
                    --the modulus is not prime (deceptive example)
                  {  0,   7,     877069,     3, 97123,     877069},
                    --fails if the modulus divides the discriminant
                  { 39, 387,      22651,    95,    27,      22651}}

--Digital signature on message hash f,
--set d > 0 to simulate corrupted data
atom f = #789ABCDE,
     d = 0

--for i=1 to length(tests) do
for i=1 to 1 do
    if not ellinit(tests[i]) then exit end if
    ec_dsa(f, d)
end for