Assignment 6.sagews

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

#Encryption Group is:
A_Encryption = 32423423
B_Encryption = 0
P_Encryption = 111111756372541867431527
G_Encryption = [32423423, 0, 111111756372541867431527]
#Public Encryption point b:
Apb_Encryption = [77012625572753817011635, 107177820655281658062820]
Bpb_Encryption = [63567155336615021877739, 36007178347937208049611]
Cpb_Encryption = [104225718616102818923406, 48710727653703914714853]
#Base point g for encryption is:
gpoint_Encryption = [133213, 42904742913851518927444]
#Notice that B=0. So if p mod 4 is 3, then the 
P_Encryption % 4

3

#So the size of our group is going to be p+1
gorder_Encryption=P_Encryption+1
factor(gorder_Encryption)

2^3 * 3^2 * 13 * 17 * 89 * 211 * 839 * 466637 * 949777

#So this group is likely vulnerable to HPS
#Recall that each secret key is x such that b=g^x=ECTimes(g,x,G)
#Let's write down what we need to use to find x_A
#x_a=HPSonEC(bA,g,G,factorization)
x_a=HPSonEC(Apb_Encryption,gpoint_Encryption,G_Encryption,[[2,3],[3,2],[13,1],[17,1],[89,1],[211,1],[839,1],[466637,1],[949777,1]])
x_a

432123

x_a=432123
ECTimes(gpoint_Encryption,x_a,G_Encryption)

[77012625572753817011635, 107177820655281658062820]

#x_a=432123
x_b=HPSonEC(Bpb_Encryption,gpoint_Encryption,G_Encryption,[[2,3],[3,2],[13,1],[17,1],[89,1],[211,1],[839,1],[466637,1],[949777,1]])
x_b

675348

#x_b=675348
x_c=HPSonEC(Cpb_Encryption,gpoint_Encryption,G_Encryption,[[2,3],[3,2],[13,1],[17,1],[89,1],[211,1],[839,1],[466637,1],[949777,1]])
x_c

876549

#signature group: G_s=[0,B,p]
#so let's check to see if p mod 3=2
#It is!
#now we factor(p+1)=2*3^5*51*620176035713316419537
#Also note that they all use the same base point g
#factor(order(g))
#but g is prime! It is the big prime from before! Yay!
#Cauchy's Theorem
#If (G,^) is any finite group, and g \in G, then order(g) divides into |G|
#signature b from Alice=sbA, etc
#again, 

#
A_Signature = 0
B_Signature = 98766778
P_Signature = 12357627687623542975694261
G_Signature = [0, 98766778, 12357627687623542975694261]
#
Apb_Signature = [5246859784562186496635824, 5457612259949488195005309]
Bpb_Signature = [5278780115068892215096216, 6857888556017625502158912]
Cpb_Signature = [11500587084934239979269346, 4630177618473717713077113]
#
T=45
gpoint_Signature = [2539662136869871783297850, 1414894921445613162605806]
gorder_Signature = 620176035713316419537
#
M_Signature = [522228913141586484647, 496418275430782450147]
y=522228913141586484647
s=496418275430782450147
M1 = [64822910729849959453845, 54324913307643112006572]
M2 = [40965317570056334177487, 76589794174623293315486]
#
#Check that p mod 3 = 2
P_Signature % 3

2

#It is! So the order of G is p+1
factor(P_Signature+1)

2 * 3^5 * 41 * 620176035713316419537

#Now for the signatures
sx_a=HPSonEC(Apb_Signature,gpoint_Signature,G_Signature,[[2,1],[3,5],[41,1],[620176035713316419537,1]])
sx_a

897672

sx_b=HPSonEC(Bpb_Signature,gpoint_Signature,G_Signature,[[2,1],[3,5],[41,1],[620176035713316419537,1]])
sx_b

213543

sx_c=HPSonEC(Cpb_Signature,gpoint_Signature,G_Signature,[[2,1],[3,5],[41,1],[620176035713316419537,1]])
sx_c

432198

#Now we have all the secret keys for encryption!

#Let's unembed M
r=HPSonEC(M2,gpoint_Encryption,G_Encryption,[[2,3],[3,2],[13,1],[17,1],[89,1],[211,1],[839,1],[466637,1],[949777,1]])
r

1234

#Now we have the r they used for encryption, we can use this to find the inverse s's
sA= ECTimes(Apb_Encryption,r,G_Encryption)
sB= ECTimes(Bpb_Encryption,r,G_Encryption)
sC= ECTimes(Cpb_Encryption,r,G_Encryption)
siA=ECInverse(sA,G_Encryption)
siB=ECInverse(sB,G_Encryption)
siC=ECInverse(sC,G_Encryption)

Mt1=ECAdd(siA,M1,G_Encryption)
Mt2=ECAdd(siB,M1,G_Encryption)
Mt3=ECAdd(siC,M1,G_Encryption)
m1=ECUnembed([Mt1],T)
#m2=ECUnembed([Mt2],T)
#m3=ECUnembed([Mt3],T)
m1
#m2
#m3

'Tokyo'

#m1 was the only one that worked, so Alice must be the recipient
M=ASCIIPad(m1,gorder_Encryption)
M

[184211207221211]

#Let's now check the signature on the message. These are my notes, put here for easy checking
#
#Verification:
#(1) Compute u=y*s^-1 mod q
#(2) Compute v=M*s^-1 mod q
#(3) Compute P=g*V+ (elliptic curve addition) u*b^-1
#Note that P is a point on the elliptic curve E
#P also has coordinates (x_1,y_1)
#The signature is valid if y=x_1 mod q
#Otherwise the signature is NOT valid
#
#So:
U= y*(s^-1) % gorder_Signature
U
z = 184211207221211
V = z*(s^-1) % gorder_Signature
V
p_1=ECTimes(gpoint_Signature,V,G_Signature)
p_2=ECTimes(ECInverse(Cpb_Signature,G_Signature),U,G_Signature)
P=ECAdd(p_1,p_2,G_Signature)
P
p=12093334749287098451036610 % gorder_Signature
p
y

364163352104018133215 537581219848672452424 [12093334749287098451036610, 411433098304526179989068] 522228913141586484647 522228913141586484647

#Because p and y match, we know that Carol sent the message
#Next, we know that r=(M-x*y)/s mod q, so for the random number for the signature, we can:
r_Signature=(z-sx_c*y)*(s^-1) % gorder_Signature
r_Signature

432554325

#Now let's check and see if it works, checking it against our given s value and s calculated with our r
#s=(M-x*y)*r^-1 mod q
s
s_check=(z-sx_c*y)*(r_Signature^-1) % gorder_Signature
s_check

496418275430782450147 496418275430782450147

#Because these are the same, we are done.