SharedSage Procedures / LUC Group Iterated Encryption.sagewsOpen in CoCalc

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)


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):
    #print "the padded number is"
    #print 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);


import itertools
def IsIn(pt, lucgp):
    if(((pt[0]^2-lucgp[0]*(pt[1]^2)) % lucgp[1]) == (4 % lucgp[1])):
        return True
    else:
        return False
    
def LUCIdentity():
    return [2,0]

def LUCInverse(point, lucasgroup):
    return([point[0],lucasgroup[1]-point[1]])

def LUCTimes(pt1, pt2, lucgp):
    if(not IsIn(pt1,lucgp)):
        print "Point 1 is not in the group"
        raise Exception("Point is not in group")
    if(not IsIn(pt2,lucgp)):
        print "Point 2 is not in the group"
        raise Exception("Point is not in group")
    else:
        luctimesx=((lucgp[1]+1)/2)*(pt1[0]*pt2[0]+lucgp[0]*pt1[1]*pt2[1]) % lucgp[1];
        luctimesy=((lucgp[1]+1)/2)*(pt1[0]*pt2[1]+pt1[1]*pt2[0]) % lucgp[1];
        return([luctimesx,luctimesy])

def LUCSquare(pt,lucgp):
    lucsqx=(pt[0]^2-2) % lucgp[1];
    lucsqy=(pt[0]*pt[1]) % lucgp[1];
    lucsqresult=[lucsqx,lucsqy]
    return lucsqresult

def LUCexp(Point,scalar,Group):
    if Point == LUCIdentity() or scalar == 0:
        return LUCIdentity()
    else:
        m = scalar
        pt = Point
        x = LUCIdentity()
        for j in itertools.count(1):
            if m%2==0:
                m = m/2
            else:
                m=(m-1)/2
                x = LUCTimes(x,pt,Group)
            if m==0:
                return x
            else:
                pt = LUCSquare(pt, Group)

def LUCSearch(lower,higher,lucgp):
    ##for j in xrange(lower,higher+1):
    for j in itertools.count(lower):
        if j >= higher:
            print "No point found in range"
            break
        luctester = (( j ^ 2 - 4 ) / lucgp[0]) % lucgp[1]
        #print luctester
        if(kronecker(luctester,lucgp[1])==1):
            foundx=j
            foundy=TonSh(luctester,lucgp[1])
            return [foundx,foundy]
        
def GroupOrder(group):
    d = group[0]
    N = group[1]
    return(N-kronecker(d,N))
def listproduct(primelist):
    x = primelist
    k = len(primelist)
    result = 1
    for i in xrange(k):
        result = result * x[i][0]^x[i][1]
    return result


def ptorder(pt, Group, factoredGpOrder):
    x = factoredGpOrder
    k = len(factoredGpOrder)
    result = listproduct(factoredGpOrder)
    for i in xrange(k):
        p = x[i][0]
        a = x[i][1]
        for j in xrange(a):
            result = result/p
            test = LUCexp(pt,result,Group)
            if (test!=LUCIdentity()):
                result = result*p
    return result

def listptorder(pt, Group, factoredGpOrder):
    x = factoredGpOrder
    k = len(factoredGpOrder)
    orderList = [];
    result = listproduct(factoredGpOrder)
    for i in xrange(k):
        p = x[i][0]
        a = x[i][1]
        ord =0
        for j in xrange(a):
            result = result/p
            test = LUCexp(pt,result,Group)
            if (test!=LUCIdentity()):
                result = result*p
                ord = ord+1
        if ord>0:
            orderlist.append([p,ord])
    return orderlist



    return [2,0]

def LUCNumberembed(xvalue,Gp,tolerance):
    if ((xvalue+1)*tolerance-1 < Gp[1]):
        pt=LUCSearch(xvalue*tolerance, (xvalue+1)*tolerance-1,Gp)
        #print pt
        return pt
    else:
        print "The embedding interval is too large for the group"

def LUCEmbed(Message,Gp,tolerance):
    AMessage = ASCIIPad(Message,QQ(Gp[1]/(tolerance+1)).floor())
    #print AMessage
    numPackets = len(AMessage)
    pt = [0]*numPackets
    for j in xrange(numPackets):
        N =AMessage[j]
        pt[j]=LUCNumberembed(N,Gp,tolerance)
    return pt


def LUCUnembed(ptlist,tolerance):
    k = len(ptlist)
    X=[0]*k
    for j in xrange(k):
        X[j]=QQ(ptlist[j][0]/tolerance).floor()
    #print X
    outStr=""
    for numInX in X:
        #print numInX
        outStr=outStr+ASCIIDepad(numInX)
    return outStr

p = 589427896723589236699649629649894277947212612604461
q = 535370272895504328886755262985028294645348311475249496914313 
R = p*q 
gammaR = lcm(lcm(p-1, q-1),lcm( p+1, q+1))
n = 1835588100642725811278316299646156385857628660058608776628872255463681001129596115865795858135978837517124882698892069345209251312769967638417909171551219459094621692202465346003400944724101186482909006146466321559204641
gcd(n, gammaR)
m = (1/n) % gammaR
n
print 
m 
print 
gammaR
1 1835588100642725811278316299646156385857628660058608776628872255463681001129596115865795858135978837517124882698892069345209251312769967638417909171551219459094621692202465346003400944724101186482909006146466321559204641 897546185455421193873195606403069779343597886021115530378295123322640162170476797485920169861351799715204572901206966125965517545802671764900882146517656036884602699958419591274867168097559891567931329661425223976514761 2074572616871462504295931902539412428803323602058731016624602959353732538193423888608074452837794138587865805219862448125803731909942996945278174776030455832494512328815350440129662697992840904878516008732568284442271320
message = "Nice system!"
mess = ASCIIPad(message, R)
M =  mess[0]
DD = (M^2-4) % R
G = [DD, R]
pt = [M, 1]
ept = LUCexp(pt, n, G)
message; M; DD; G; pt; ept
'Nice system!' 178205199201132215221215216201209133 31757093022315213918069141249722233075207263075986889676020571202611685 [31757093022315213918069141249722233075207263075986889676020571202611685, 315562173921131111550753275258954314886929684614433857059068225329589007516900986203500713113234143069878550293] [178205199201132215221215216201209133, 1] [83998001190803516381268591520945089727506133337704160937617393083672591255727337006232243124900847240919478370, 180083374451388718720354846216799073863947458877212034881497682919879387959145785788125908599174273371748135361]
e = ept[0]
encryption = ept[0]
for j in xrange(24563453455):
    dd = (e^2-4)%R
    G1 = [dd,R]
    pt1 = [e,1]
    e1 = LUCexp(pt1,n,G1)
    if e1[0]==encryption:
        A = ASCIIDepad(pt1[0])
        print A
        break
    e = e1[0]
Nice system!