ECProcedures.sagews

Elliptic Curve Procedures

import numpy
import itertools

############################################################################################################################
# For all these procedures we are working with elliptic curve (EC) groups over Z_p where p is a prime number larger than 3.#
# An Elliptic Curve group will be represented by a triple [A,B,p] where                                                    #
#      (1) p is a prime number larger than 3                                                                               #
#      (2) A is the coefficient of x and                                                                                   #
#      (3) B is the constant term in  y^2 = x^3 + A*x + B mod p.                                                           #
#                                                                                                                          #
# A point in the group will be represented as a pair [x,y] where                                                           #
#      (1)  x is the x-coordinate of a solution                                                                            #
#      (2)  y is the y-coordinate of a solution                                                                            #
#      (3)  The identity will be represented by []                                                                         #
############################################################################################################################

def TonSh (a, Prime):
    if kronecker(a, Prime) == -1:
        pass
        print("{0} does not have a square root mod {1}".format(a, Prime))
        return None
    elif Prime % 4 == 3:
        x = power_mod(ZZ(a), ZZ((Prime+1)/4),ZZ(Prime))
        return(x)
    else:
        #################################################################################
        # Determine e so that Prime-1 = (2^e)*q for some odd number q                   #
        #################################################################################

        e = ZZ(0)
        q = ZZ(Prime - 1)
#        for j in range(1, Prime):
        for j in itertools.count(1):
            if q % 2 == 0:
                e = ZZ(e + 1)
                q = ZZ(q / 2)
            else:
                break
        for i in range(1, 101):
            n = i
            if kronecker(n, Prime) == -1:
                break
        z = power_mod(ZZ(n),ZZ(q),ZZ(Prime))
        y = ZZ(z)
        r = ZZ(e)
        x = power_mod(ZZ(a),ZZ( (q-1)/2),ZZ( Prime))
        b = (a * x ** 2) % Prime
        x = (a * x) % Prime
        #for i in range(1, e + 1):
        for i in itertools.count(1):
            if b == 1:
                break
            else:
                for m in itertools.count(1):
                    t = power_mod(ZZ(b), ZZ(2^m) ,  ZZ(Prime))
                    if t == 1:
                        mm = m
                        break
                t = power_mod(ZZ(y), ZZ(2^(r - mm - 1)),ZZ(Prime))
                y = power_mod(ZZ(t), ZZ(2),ZZ(Prime))
                r = mm
                x = (x * t) % Prime
                b = (b * y) % Prime
        return(x)

ECSearch is a procedure that searches for a point on a given elliptic curve with first coordinate between "lowerbound" and "upperbound". The lowerbound and the upperbound are positive random numbers.

import itertools
def ECSearch (lowerbound, upperbound, Group):
    p = Group[2]
    a = Group[0]
    b = Group[1]
    if (4 * a ** 3 + 27 * b ** 2) % p == 0:
        print("This is not an elliptic curve. ")
    else:
        for j in itertools.count(lowerbound):
            if j==lowerbound-1 or j>upperbound or j>upperbound:
                return "not found"
            j=j%p
            #print "test "+str(j)
            #print (kronecker(j ** 3 + a * j + b, p) )
            if kronecker(j ** 3 + a * j + b, p) == 1:
                x = (j ** 3 + a * j + b) % p
                var('z')
                #print ("Modsqrt({0}, {1})".format(x,p))
                #L = solve_mod(z ** 2 == x, p)
                L=TonSh(x,p)
                y = L
                #y = L[0][0]
                return([j,y])

EXAMPLE: Find a point on the elliptic curve $y^2=x^3+2x+14\;mod\; 1291$ such that the first coordinate of the point is between (234,789).

P=ECSearch(234,789,[2,14,1291])
P
Q=ECSearch(345,891,[2,14,1291])
Q

[236, 1131] [346, 749]

import sympy.ntheory
def FermatGreat (N):
    if sympy.ntheory.isprime(N):
        if N % 4 == 1:
            Z = two_squares(N)
            #print(Z)
            return(Z)
        else:
            print("{0} is prime, but {1} mod 4 = 3").format(N, N)
    else:
        print("{0} is not prime\n").format(N)

FermatGreat(17)

(1, 4)

ECAdd is a procedure that computes the sum of two points on a given elliptic curve using the elliptic curve operation.

def ECAdd (Point1, Point2, Group):
    a = Group[0]
    b = Group[1]
    p = Group[2]
    if(Point1!=[]):
        x1 = Point1[0]
        y1 = Point1[1]
    if(Point2!=[]):
        x2 = Point2[0]
        y2 = Point2[1]
    if (4 * a ** 3 + 27 * b ** 2) % p == 0:
        print("This is not an elliptic curve. ")
    elif Point1 != [] and y1 ** 2 % p != (x1 ** 3 + a * x1 + b) % p:
        print("Point 1 is not on the elliptic curve. \n")
    elif Point2 != [] and y2 ** 2 % p != (x2 ** 3 + a * x2 + b) % p:
        print("Point 2 is not on the elliptic curve. \n")
    else:
        if Point1 == []:
            R = Point2
        elif Point2 == []:
            R = Point1
        else:
            if x1 == x2 and 0 == (y1 + y2) % p:
                R = []
            elif x1 == x2 and y1 == y2:
                R = ECDouble(Point1, Group)
            else:
                s = ((y1 - y2) / (x1 - x2)) % p
                x = (s ** 2 - x1 - x2) % p
                y = (s * (x1 - x) - y1) % p
                R = [x,y]
    return(R)

EXAMPLE: Find the sum of the points (236, 1131) and (346, 749) on the elliptic curve $y^2=x^3+2x+14\;mod\; 1291$.

ECAdd([236,1131],[346,749],[2,14,1291])

[1109, 774]

def ECDouble (Point, Group):
    a = Group[0]
    b = Group[1]
    p = Group[2]
    if Point != []:
        x1 = Point[0]
        y1 = Point[1]
        if (4 * a ** 3 + 27 * b ** 2) % p == 0:
            print("This is not an elliptic curve. ")
        elif Point != [] and y1 ** 2 % p != (x1 ** 3 + a * x1 + b) % p:
            print("The point to double is not on the elliptic curve. ")
        elif y1 == 0:
            R = []
        else:
            s = mod((3 * x1 ** 2 + a) / y1 / 2, p)
            x = (s ** 2 - 2 * x1) % p
            y = (s * (x1 - x) - y1) % p
            R = [x,y]
    else:
        R = []
    return(R)

EXAMPLE: Find the sum of the point (6,9) to itself i.e. $2\cdot(6,9)$ on the curve $y^2=x^3+x+2\;mod\; 13$.

ECDouble([6,9],[1,2,13])

[2, 8]

ECInverse is a procedure that finds the inverse of a point on a given elliptic curve using the elliptic curve operation.

def ECInverse (Point, Group):
    if Point == []:
        return(Point)
    else:
        p = Group[2]
        x = Point[0]
        y = Point[1]
        return([x,(p - y) % p])

EXAMPLE: Find an inverse of the point (6,9) on the curve $y^2=x^3+x+2\;mod\; 13$.

ECInverse([6,9],[1,2,13])

[6, 4]

ECTimes is a procedure that adds a point k-times (scalar) to itself for a given elliptic curve using the elliptic curve operation.

def ECTimes (Point, scalar, Group):
    ECIdentity=[]
    if Point == ECIdentity or scalar == 0:
        return(ECIdentity)
    else:
        E=EllipticCurve(GF(Group[2]),[Group[0],Group[1]])
        EPoint = E(Point[0],Point[1])
        #print type(EPoint)
        cgret = scalar*EPoint
        #print cgret[0]
        if(cgret[0]==0 and cgret[1]==1 and cgret[2]==0):
            return ECIdentity
        return([cgret[0],cgret[1]])

EXAMPLE: Multiply the point (1,11) by 71 (add the point (1,11) to itself 71 times) on the curve $y^2=x^3+x+2\;mod\; 13$.

ECTimes([1,11],71,[1,2,13])
G=[3,2,next_prime(11696876878768768767867687687687687678768768766876879878111)]
G ####### Elliptic Curve Group
g=ECSearch(2378, 892749287373,G)
g ####### base element
x=34239473274983742
x ####### private key
b=ECTimes(g,x,G)

[1, 2] [3, 2, 11696876878768768767867687687687687678768768766876879878131] [2378, 1882214856783904319500075249934504476880494552940632700484] 34239473274983742

def listproduct (primelist):
    x = primelist
    k = len(primelist)
    print(k)
    result = 1
    for i in range(0, k):
        result = result * x[i][0] ^ x[i][1] #op(1, op(i, x)) ** op(2, op(i, x))
    return(result)

def listptorder (pt, Group, factoredGpOrder):
    ECIdentity = []
    x = factoredGpOrder
    k = len(factoredGpOrder)
    orderlist = set([])
    result = listproduct(factoredGpOrder)
    for i in range(0, k ):
        p =  x[i][0] #op(1, op(i, x))
        a =  x[i][1] #op(2, op(i, x))
        ord = 0
        for j in range(0, a ):
            result = result / p
            test = ECTimes(pt, result, Group)
            if test != ECIdentity:
                result = result * p
                ord = ord + 1
        if 0 < ord:
            orderlist = orderlist | set([[p,ord]])
    return(convert(orderlist, list))


def primeHPSonEC (y,g,Group,p):
    newg = []#ECIdentity
    for j in range(0, p ):
        if y == newg:
            break
        newg = ECAdd(newg, g, Group)
    return(j)


def powerHPSonEC (y,g,og,p,a,Group):
    gp = ECTimes(g, og / p, Group)
    newog = og
    newy = y
    c = 0
    partx = 0
    for i in range(0, a ):
        newy = ECAdd(y, ECInverse(ECTimes(g, partx, Group), Group), Group)
        newog = newog / p
        ypower = ECTimes(newy, newog, Group)
        c = primeHPSonEC(ypower, gp, Group, p)
        partx = partx + c * p ** i
    return(partx)

def HPSonEC (y,g,Group,factoredGpOrder):
    X = 0
    Ord = 1
    K = listptorder(g, Group, factoredGpOrder)
    k = len(K)
    ordg = listproduct(K)
    for j in range(1, k + 1):
        p = K[j][0]#op(1, op(j, K))
        powerp = K[j][1]#op(2, op(j, K))
        partx = powerHPSonEC(y, g, ordg, p, powerp, Group)
        if j == 1:
            X = partx
            Ord = p ** powerp
        else:
            X = crt([X,partx], [Ord,p ** powerp])
            Ord = Ord * p ** powerp
    return(X)

def ECEmbed (Message, gp, tolerance):
    p = ZZ(math.floor(gp[2] / (tolerance + 1)))
    M = ASCIIPad(Message, p)
    packets = len(M)
    #print packets
    pointlist = [0]*packets
    for j in range(0, packets ):
        N = M[j]
#        N:=(op(0,op(1,M))[j]);
        pointlist[j] = ECSearch(tolerance * N, tolerance * (N + 1) - 1, gp)
        #print pointlist
    return(pointlist)

def ECUnembed (pointlist, tolerance):
    #print pointlist
    k = len(pointlist)
    paddedasciilist=[0]*k
    for j in range(0, k ):
        #print type(pointlist[j][0]/tolerance)
        #print pointlist
        #print paddedasciilist
        #print "test this "+str((pointlist[j][0]/tolerance))
        #.floor() returns correct while math.floor() does not work
        #print (pointlist[j][0])
        pointlist[j][0]=ZZ(pointlist[j][0])
        toType = ZZ(QQ((pointlist[j][0])/tolerance).floor())
#        typed = type(toType)
#        print typed
        paddedasciilist[j] = ((pointlist[j][0])/tolerance).floor()
    returnStr = ""
    for paddedItem in paddedasciilist:
        #print paddedItem
        buffer = ASCIIDepad(paddedItem)
        returnStr+= buffer
        #print buffer
    return returnStr


import math

def ASCIIPad(Mess,Mod):
    K = []
    for letter in Mess:
        K.append(ord(letter))
    L = Mod.ndigits()
    l = len(K)
    y = ZZ(math.floor(L/3))
    count = 0
    padded = []
    buffer = ""
    for numChar in K:
        numChar+=100
        buffer+=str(numChar)
        count+=1
        if count==y:
            padded.append(ZZ(buffer))
            buffer = ""
            count = 0
    if len(buffer)>0:
        padded.append(ZZ(buffer))
    return padded


def ASCIIDepad(Number):
    N = "";
    #Number = ZZ(Number[0])
    n = Number.ndigits() % 3;
    if (n > 0):
        print("This is not a padded ASCII string\n");
    else:
        L = [((Number - (Number % (1000^i)))/1000^i)%1000 - 100 for i in range(Number.ndigits()/3)];
        for i in range(Number.ndigits()/3):
            #print L[i]
            N = chr(L[i]) + N;

    return(N);

ECSearch is a procedure that searches for a point on a given elliptic curve with first coordinate between "lowerbound" and "upperbound". The lowerbound and the upperbound are positive random numbers.

import itertools
def ECSearch (lowerbound, upperbound, Group):
    p = Group[2]
    a = Group[0]
    b = Group[1]
    if (4 * a ** 3 + 27 * b ** 2) % p == 0:
        print("This is not an elliptic curve. ")
    else:
        for j in itertools.count(lowerbound):
            if j==lowerbound-1 or j>upperbound or j>upperbound:
                return("not found")
            j=j%p
            #print "test "+str(j)
            #print (kronecker(j ** 3 + a * j + b, p) )
            if kronecker(j ** 3 + a * j + b, p) == 1:
                x = (j ** 3 + a * j + b) % p
                var('z')
                #print("stuck here")
                #print ("Modsqrt({0}, {1})".format(x,p))
                #L = solve_mod(z ** 2 == x, p)
                L=TonSh(x,p)
                y = L
                #y = L[0][0]
                return([j,y])

#pr=next_prime(17416141711513452221)
#Ec=[21,3,pr]
#pt=[1,5]
#pk=834756
#rn = 324765348769
#G1 = EllipticCurve(GF(pr,'a'),[21,3])
#p1 = G1(1,5)
#pk = pk*p1
#rp = rn*p1
#Go = G1.order()
##fgo = factor(Go)
##print(pr, pk, Go, fgo)
#Message = "Deposit $50 in account number 23456789"
#EMessage = ECEmbed(Message,G1,50)
#EMessage
##AM = ASCIIPad(Message,Go)
##AM
##D=ECSearch(29,31,E)
##Q=ECSearch(34,36,E)
##D
##Q
##

2398456394587698346003