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 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
                #if y>p/2:
                #    y=p/2-(y-p/2)
                #y = L[0][0]
                return([j,y])

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)

def ECAdd (Point1, Point2, Group):
    
    a = Group[0]
    b = Group[1]
    p = Group[2]
    '''
    #This could be an optimization for code readability
    E = EllipticCurve(GF(p),[a,b])
    p1 = E(Point1[0],Point1[1])
    p2 = E(Point2[0],Point2[1])
    toReturn = p1+p2
    return [toReturn[0],toReturn[1]]
    '''
    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 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 ECInverse (Point, Group):
    if Point == []:
        return(Point)
    else:
        p = Group[2]
        x = Point[0]
        y = Point[1]
        return([x,(p - y) % p])

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

def findptOrder(point,group):
    E = EllipticCurve(GF(group[2]),[group[0],group[1]])
    Ep = E(point[0],point[1])
    return Ep.additive_order()

def Dirichlet (q):
    j=0

    #for j in numpy.arange(0.10e1, infinity + 0.10e1, 0.10e1):
    while True:
        j+=1
        p = 2 * j * q + 1
        if is_prime(p) == True:
            return(p)

#Constructing El Gamal key resistant against Hellman Pohlig Silver attack
#Choose a large prime number q
q=next_prime(283743928792183019830289149837287463847938240218302198308921094832985743875834);
#Use the Dirichlet procedure to find a prime number p such that (p-1) has q as one of its prime factors
p=Dirichlet(q);
p
A=3476283764
B=0
G=[A,B,p]
k=factor(p+1);
k

#Constructing a signature key that is resistant to the HPSonP attacks

q1 = next_prime(37493278489327498732487398498393879480898493287489732974921094832985743875834)
p1 = Dirichlet(q1)
fL1=list(factor(p1-1))
for test in fL1:
    print test
c11 = 1
c21 = 100
c31 = 374932784893274987324873984983938794808984932874897329749210948329857438758
x1 = crt([c11,c21,c31],[fL1[0][0],fL1[1][0],fL1[2][0]])
x1
g1=findLargePrimitiveRoot(10038748923784,p1)
g1
b1 = power_mod(g1,x1,p1)
b1

q2 = next_prime(89324793827492374982748732402938332894798230912830912921094832985743875834)
p2 = Dirichlet(q2)
fL2=list(factor(p2-1))
c12 = 1
c22 = 30
c32 = 8932479382749237498274873240293833289479823091283091292109483298574387616
x2 = crt([c12,c22,c32],[fL2[0][0],fL2[1][0],fL2[2][0]])
x2
g2=findLargePrimitiveRoot(10038748923784,p2)
g2
b2 = power_mod(g2,x2,p2)
b2

q3 = next_prime(898328932749832794872478923748923749873289473892743829749238742189723879294832985743875834)
p3 = Dirichlet(q3)
fL3=list(factor(p3-1))
c13 = 1
c23 = 50
c33 = 898328932749832794872478923748923749873289473892743829749238742189723879294832985743875859
x3 = crt([c13,c23,c33],[fL3[0][0],fL3[1][0],fL3[2][0]])
x3