diff --git a/.gitignore b/.gitignore index 2ed75c3..78aba10 100644 --- a/.gitignore +++ b/.gitignore @@ -246,3 +246,4 @@ TSWLatexianTemp* *.spl .ipynb_checkpoints/ +.vscode/ diff --git a/ctnt2018/ECM.py b/ctnt2018/ECM.py index 975dbbe..dc01c5f 100644 --- a/ctnt2018/ECM.py +++ b/ctnt2018/ECM.py @@ -1,118 +1,15 @@ # %load ec3.py from math import factorial, gcd, log +import pickle import numpy as np +import gmpy2 -P=np.zeros(10000) -P[0]=1 -P[1]=1 -for i in range(2,10000): - if P[i]==0: - j=2 - while i*j<1000: - P[i*j]=1 - j=j+1 -Primes10000=[i for i,x in enumerate(P) if x==0 ] - - -def mexp(a,x,N): - m,s=1,a - while x>0: - if x % 2 ==1: - m=((m*s) % N) - s=((s*s) % N) - x=x//2 - return m - -def euclid(u,v): - if v==0: - raise ArithmeticError('Division by Zero') - x0,x1=u,v - a0,a1=1,0 - b0,b1=0,1 - while x1!=0: - q=x0//x1 - x2=x0-q*x1 - a2=a0-q*a1 - b2=b0-q*b1 - x0,a0,b0=x1,a1,b1 - x1,a1,b1=x2,a2,b2 - if x0<0: - return -x0,-a0,-b0 - else: - return x0,a0,b0 - -def mod_inv(u,N): - d,a,b=euclid(u,N) - if d==1: - return a - else: - raise ArithmeticError('Common factor is '+str(d)) - -def two_p(x,y,a,b,N): - Lu=(3*x**2+a) % N - # print(Lu) - Lb=mod_inv(2*y,N) - # print(Lb) - L=Lu*Lb % N - x_two=(L*L-2*x) % N - y_two=(L*(x-x_two)-y) %N - return x_two,y_two - -def sum_p(x1,y1,x2,y2,a,b,N): - Lu=(y2-y1) % N - Lb=mod_inv(x2-x1,N) - L=(Lu*Lb) % N - x_sum=(L*L-x1-x2) %N - y_sum=(L*(x1-x_sum)-y1) %N - return x_sum,y_sum - -def exp_p(x,y,a,b,m,N): - sx,sy=x,y - first=True - while m>0: - if m%2==1: - if first: - xm,ym=sx,sy - first=False - else: - xm,ym=sum_p(xm,ym,sx,sy,a,b,N) - sx,sy=two_p(sx,sy,a,b,N) - m=m//2 - return xm,ym - -#def mexp(a,x,N): - # m,s=1,a - # while x>0: - # if x % 2 ==1: - # m=((m*s) % N) - # s=((s*s) % N) - # x=x//2 - # return m - -def ecm_trial(N,arange=50,krange=30): - for a in range(-arange,arange): - xm,ym=0,1 - print(a) - for k in range(2,krange): - try: - xm,ym=exp_p(xm,ym,a,1,k,N) - except ArithmeticError: - print('try the following: a=',a,' and k=',k) - break - - - -#N=149185656432189838133 -#ecm_trial(N,arange=20,krange=10000) -N=2**128+1 -#ecm_trial(N,arange=100,krange=10000) -#xm,ym=exp_p(0,1,-91,1,factorial(7883),N) - - +with open('primes10000.pickle','rb') as f: + Primes10000=pickle.load(f) + def ecm(N,arange=100,krange=10000): for a in range(-arange,arange): x,y=0,1 - print(a) for B in Primes10000: k=B**(int(log(10000)/log(B))) sx,sy,t=x,y,k @@ -124,7 +21,7 @@ def ecm(N,arange=100,krange=10000): first=False else: # x1=xm,y1=ym, x2=sx,y2=sy,a=a,b=1,N=N - d,u,v=euclid(sx-xm,N) + d,u,v=gmpy2.gcdext(sx-xm,N) if d>1: if d==N: break @@ -135,7 +32,7 @@ def ecm(N,arange=100,krange=10000): ym=(L*(xm-x_sum)-ym) % N xm=x_sum # sx=x,sy=y,a=a,b=1,N - d,u,v=euclid(2*sy,N) + d,u,v=gmpy2.gcdext(2*sy,N) if d>1: if d==N: break diff --git a/ctnt2018/MakePrimes10000.py b/ctnt2018/MakePrimes10000.py new file mode 100644 index 0000000..485d6d9 --- /dev/null +++ b/ctnt2018/MakePrimes10000.py @@ -0,0 +1,18 @@ +import numpy as np +import pickle + + + +P=np.zeros(10000) +P[0]=1 +P[1]=1 +for i in range(2,10000): + if P[i]==0: + j=2 + while i*j<1000: + P[i*j]=1 + j=j+1 + +Primes10000=[i for i,x in enumerate(P) if x==0] +with open('primes10000.pickle','wb') as f: + pickle.dump(Primes10000,f) diff --git a/ctnt2018/old_ecm.py b/ctnt2018/old_ecm.py new file mode 100644 index 0000000..16d8de3 --- /dev/null +++ b/ctnt2018/old_ecm.py @@ -0,0 +1,94 @@ +def mexp(a,x,N): + m,s=1,a + while x>0: + if x % 2 ==1: + m=((m*s) % N) + s=((s*s) % N) + x=x//2 + return m + +def euclid(u,v): + if v==0: + raise ArithmeticError('Division by Zero') + x0,x1=u,v + a0,a1=1,0 + b0,b1=0,1 + while x1!=0: + q=x0//x1 + x2=x0-q*x1 + a2=a0-q*a1 + b2=b0-q*b1 + x0,a0,b0=x1,a1,b1 + x1,a1,b1=x2,a2,b2 + if x0<0: + return -x0,-a0,-b0 + else: + return x0,a0,b0 + +def mod_inv(u,N): + d,a,b=euclid(u,N) + if d==1: + return a + else: + raise ArithmeticError('Common factor is '+str(d)) + +def two_p(x,y,a,b,N): + Lu=(3*x**2+a) % N + # print(Lu) + Lb=mod_inv(2*y,N) + # print(Lb) + L=Lu*Lb % N + x_two=(L*L-2*x) % N + y_two=(L*(x-x_two)-y) %N + return x_two,y_two + +def sum_p(x1,y1,x2,y2,a,b,N): + Lu=(y2-y1) % N + Lb=mod_inv(x2-x1,N) + L=(Lu*Lb) % N + x_sum=(L*L-x1-x2) %N + y_sum=(L*(x1-x_sum)-y1) %N + return x_sum,y_sum + +def exp_p(x,y,a,b,m,N): + sx,sy=x,y + first=True + while m>0: + if m%2==1: + if first: + xm,ym=sx,sy + first=False + else: + xm,ym=sum_p(xm,ym,sx,sy,a,b,N) + sx,sy=two_p(sx,sy,a,b,N) + m=m//2 + return xm,ym + +#def mexp(a,x,N): + # m,s=1,a + # while x>0: + # if x % 2 ==1: + # m=((m*s) % N) + # s=((s*s) % N) + # x=x//2 + # return m + +def ecm_trial(N,arange=50,krange=30): + for a in range(-arange,arange): + xm,ym=0,1 + print(a) + for k in range(2,krange): + try: + xm,ym=exp_p(xm,ym,a,1,k,N) + except ArithmeticError: + print('try the following: a=',a,' and k=',k) + break + + + +#N=149185656432189838133 +#ecm_trial(N,arange=20,krange=10000) +N=2**128+1 +#ecm_trial(N,arange=100,krange=10000) +#xm,ym=exp_p(0,1,-91,1,factorial(7883),N) + \ No newline at end of file diff --git a/primes10000.pickle b/primes10000.pickle new file mode 100644 index 0000000..2d12955 Binary files /dev/null and b/primes10000.pickle differ