Assignment 7.sagews

import itertools
import numpy
import itertools
import math

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)

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])

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]])
    
    
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)

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);

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

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)

def ECInverse (Point, Group):
    if Point == []:
        return(Point)
    else:
        p = Group[2]
        x = Point[0]
        y = Point[1]
        return([x,(p - y) % p])
def findptOrder(point,group):
    E = EllipticCurve(GF(group[2]),[group[0],group[1]])
    Ep = E(point[0],point[1])
    return Ep.additive_order()

def listproduct (primelist):
    x = primelist
    k = len(primelist)
    #print x
    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 ):
            try:
                #print result / p
                result = ZZ(result / p)
                test = ECTimes(pt, result, Group)
                if test != ECIdentity:
                    result = result * p
                    ord = ord + 1
            except:
                pass
        if 0 < ord:
            orderlist.append([p,ord])
    #print orderlist
    returnlist = []
    return orderlist#(convert(orderlist, list))

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


def powerHPSonEC (y,g,og,p,a,Group):
    gp = ECTimes(g, ZZ(og / p), Group)
    newog = og
    newy = y
    c = 0
    partx = 0
    for i in range(0, a ):
        #print (y, ECInverse(ECTimes(g, partx, Group), Group))
        newy = ECAdd(y, ECInverse(ECTimes(g, partx, Group), Group), Group)
        newog = ZZ(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)
    #print K
    for j in range(0, k):
        #print (K,j)
        p = K[j][0] #op(1, op(j, K))
        powerp = K[j][1] #op(2, op(j, K))
        #print powerp
        partx = powerHPSonEC(y, g, ordg, p, powerp, Group)
        #print partx
        if j == 0:
            X = partx
            Ord = p ** powerp
            #print "firstif"
            #print X
            #print Ord
        else:
            #print [X,partx]
            #print [Ord,p^powerp]
            X = crt([X,partx], [Ord,p ^ powerp])
            Ord = Ord * (p ^ powerp)
        #print X
    return(X)

p=72282850616325306432153954925097782772391298102482728458113409101633304290542793749999999
A=1
B=0
g=[17801318456920812504189779759683982034207488925189310604312386489063003522183777461745820,13493528240040346907430282491247753306238353560185327103875081990944182571181305682391385]
q=36141425308162653216076977462548891386195649051241364229056704550816652145271396875000000
b=[1,19157523638784122540010643100163963974332197057266595252451199031758721152257560625189409]
R=[298757533787476246810285089213201709898407152028354406106,2951975782515941124925610013860677667037257530745608706186687445372449787917550124573688]
r=26256091035727029708701170901895407906552904172593915065416857031520110161809746688938699
G=[A,B,p]


E = EllipticCurve(GF(p),[A,B])
Eord = E.order()
print "Group Order:", Eord
factor(Eord)



w=49792922297912707801714181535533618316401192004725734351
k=2^6 * 3^7 * 5^11 * 7^20 * 11^3



E=ECTimes(b,w,G)
W=ECTimes(g,w,G)



z=HPSonEC(E,W,G,[[2,7],[3,7],[5,11],[7,20],[11,3],[211,1],[49792922297912707801714181535533618316401192004725734351,1]])
z
c=R[0]/w


Message = "Sell today!"
m = ASCIIPad(Message,p)[0]

y=ZZ(mod(R[0],q))

s=ZZ(mod((m - (w*z*c))/r,q))
print "Signature:",(y,s)

Group Order: 72282850616325306432153954925097782772391298102482728458113409101633304290542793750000000 2^7 * 3^7 * 5^11 * 7^20 * 11^3 * 49792922297912707801714181535533618316401192004725734351 181458647335093656933513281250000 Signature: (298757533787476246810285089213201709898407152028354406106, 18070712654081326608038488731274445693097824525620750136953758746692073075057726446961767)

u = ZZ(mod((m/s),q))

v = ZZ(mod(y/s,q))

S=ECAdd(ECTimes(g,u,G),ECTimes(b,v,G),G)
S[0]
y

298757533787476246810285089213201709898407152028354406106 298757533787476246810285089213201709898407152028354406106