1 #include "sympow.h"
2 #define DEBUG (FALSE || GLOBAL_DEBUG)
3
4 typedef struct {double x,y; int i;} ap_pt;
5
6 #define QD_split(a,hi,lo)\
7 {double xt=134217729.0*a,t1; t1=xt-a; hi=xt-t1; lo=a-hi;}
8 #define QD_prod(a,b,d,e)\
9 {double ah,al,bh,bl,u1,u2,u3; d=a*b; QD_split(a,ah,al); QD_split(b,bh,bl);\
10 u1=ah*bh; u2=u1-d; u3=ah*bl; u1=u2+u3; u2=al*bh; u3=u1+u2; u2=al*bl; e=u2+u3;}
11
12 #define BIG_FLOAT 6755399441055744.0 /* 2^52+2^51 */
13 #define ROUND(x,r) {double y=(x),w1; w1=y+BIG_FLOAT; r=w1-BIG_FLOAT;}
14
15 #define MUL(x,y,z,p)\
16 {double zhi,zlo,thi,tlo,fl,v1,v2; QD_prod((x),(y),zhi,zlo); ROUND(zhi/p,fl);\
17 QD_prod(p,fl,thi,tlo); v1=zhi-thi; v2=zlo-tlo; z=v1+v2;}
18 #define INV(x,i,p)\
19 {double rr=0.0,ss=1.0,cc,aa=p,bb=x; while (1)\
20 {ROUND(aa/bb,cc); aa-=bb*cc;\
21 if (aa==0.0) {if (bb>0.0) i=ss; else i=-ss; break;} rr-=cc*ss;\
22 ROUND(bb/aa,cc); bb-=aa*cc;\
23 if (bb==0.0) {if (aa>0.0) i=rr; else i=-rr; break;} ss-=cc*rr;}}
24 #define NORMALISE(x,p) while ((x)>=(p)) (x)-=(p); while ((x)<0) (x)+=(p)
25
26 #define FROM_LINE(P,Q,R,s,p)\
27 {double T,x1; ap_pt V;\
28 MUL(s,s,T,p); V.x=T-(P.x+Q.x); NORMALISE(V.x,p);\
29 x1=P.x-V.x; MUL(x1,s,T,p); V.y=T-P.y; NORMALISE(V.y,p); *(R)=V;}
30
ec_ap_large_pt_double(double a4,ap_pt P,ap_pt * R,double p)31 static int ec_ap_large_pt_double(double a4,ap_pt P,ap_pt *R,double p)
32 {double T,n,i,s,tt;
33 if (P.y==0.0) return TRUE; MUL(P.x,P.x,T,p); n=T+T+T+a4; /*NORMALISE(n,p);*/
34 tt=P.y+P.y; INV(tt,i,p); MUL(n,i,s,p); FROM_LINE(P,P,R,s,p); return FALSE;}
35
ec_ap_large_pt_add(double a4,ap_pt P,ap_pt Q,ap_pt * R,double p)36 static int ec_ap_large_pt_add(double a4,ap_pt P,ap_pt Q,ap_pt *R,double p)
37 {double Ld,Ln,i,s;
38 if (P.x==Q.x)
39 {if (P.y==Q.y) return ec_ap_large_pt_double(a4,P,R,p); return TRUE;}
40 Ld=P.x-Q.x; Ln=P.y-Q.y; /* allow Ld to be negative */
41 INV(Ld,i,p); MUL(Ln,i,s,p); FROM_LINE(P,Q,R,s,p); return FALSE;}
42
43 #define ADDI(P,Q,R,p,i) {double s; MUL(P.y-Q.y,i,s,p); FROM_LINE(P,Q,R,s,p);}
44 #define DOUBLEI(a4,P,R,p,i)\
45 {double TT,nn,s; MUL(P.x,P.x,TT,p); nn=TT+TT+TT+a4;\
46 MUL(nn,i,s,p); FROM_LINE(P,P,R,s,p);}
47
ec_ap_large_pt_mult(double a4,ap_pt P,llint n,ap_pt * R,double p)48 static int ec_ap_large_pt_mult(double a4,ap_pt P,llint n,ap_pt *R,double p)
49 {llint m=n; ap_pt Q=P,T,U; int z=TRUE; /* P is not O */
50 if (n==0) return TRUE; if (n==1) {*R=P; return FALSE;}
51 if (m&1) {U=Q; z=FALSE;}
52 while (m)
53 {m=m>>1;
54 if (ec_ap_large_pt_double(a4,Q,&T,p)) break; Q=T;
55 if (m&1)
56 {if (z) {U=Q; z=FALSE;}
57 else {z=ec_ap_large_pt_add(a4,Q,U,&T,p); if (!z) U=T;}}}
58 if (!z) *R=U; return z;}
59
ec_ap_large_order_from_mult(double a4,ap_pt P,llint m,double pr)60 static llint ec_ap_large_order_from_mult(double a4,ap_pt P,llint m,double pr)
61 {int i,e,k,t; llint o=1,n=m,p,pe; ap_pt Q=P,R; LIST f;
62 f.p[0]=0; ifactor(n,&f,1,1000);
63 for (i=0;f.p[i]!=0;i++)
64 {p=f.p[i]; e=f.e[i]; for (t=0;t<e;t++) n=n/p;
65 if (ec_ap_large_pt_mult(a4,Q,n,&R,pr)) continue;
66 for (k=1;k<e;k++) {if (ec_ap_large_pt_mult(a4,R,p,&R,pr)) break;}
67 e=k; pe=1; for (t=0;t<e;t++) pe=pe*p; o=o*pe;
68 if (ec_ap_large_pt_mult(a4,Q,pe,&R,pr)) break; Q=R;}
69 return o;}
70
71 #define TYPE ap_pt
72 #define NAME ec_ap_large_qsort
73 #define CMP(A,B) (A->x>B->x)
74 #define SWAP_VARS(type,A,B) {type _ttemmp; _ttemmp=A; A=B; B=_ttemmp;}
75 #define SWAP_DATA(A,B) {SWAP_VARS(double,A->x,B->x);\
76 SWAP_VARS(double,A->y,B->y); SWAP_VARS(int,A->i,B->i);}
77
NAME(TYPE * ARR,int n)78 static void NAME(TYPE *ARR,int n)
79 {TYPE *left,*right,*ip,*jp,*last,*mid; int l;
80 if (DEBUG) printf("ec_ap_large_qsort %i\n",n);
81 if (n<=1) return; left=ARR;
82 if (n==2) {ip=left+1; if (CMP(left,ip)) SWAP_DATA(left,ip); return;}
83 right=left+(n-1);
84 if (n<=5)
85 {for (ip=left;ip<right;ip++)
86 for (jp=ip+1;jp<=right;jp++) if (!CMP(jp,ip)) SWAP_DATA(ip,jp); return;}
87 mid=left+(n/2);
88 if (!CMP(left,mid))
89 {if (!CMP(mid,right)) SWAP_DATA(left,mid)
90 else if (!CMP(left,right)) SWAP_DATA(left,right)}
91 else
92 {if (CMP(mid,right)) SWAP_DATA(left,mid)
93 else if (CMP(left,right)) SWAP_DATA(left,right)}
94 l=0; last=left;
95 for (ip=left+1;ip<=right;ip++)
96 if (!CMP(ip,left)) {l++; last++; SWAP_DATA(last,ip)}
97 if (left!=last) SWAP_DATA(left,last)
98 NAME(left,l); NAME(last+1,n-l-1);}
99
100 #define HEAP_SEARCH(n,A,sz,k)\
101 {int LL=0,HH=sz,w;\
102 while (LL<HH) {w=(LL+HH)>>1; if (A[w].x<n) LL=w+1; else HH=w;}\
103 if (DEBUG) printf("heap_search %i %f %f\n",HH,A[HH].x,n);\
104 if ((HH<sz) && (A[HH].x==n)) k=HH; else k=-1;}
105 #define DONE(x) {free(ARR); return ((llint) (x));}
106 #define DONEBS(x) {free(B); free(S); DONE((x));}
107 #define MONTGOMERY_INVERSE_BABY_STEPS TRUE
108
ec_ap_large_bsgs(double a4,ap_pt P,double p,llint L)109 static llint ec_ap_large_bsgs(double a4,ap_pt P,double p,llint L)
110 {int i,j,k,bs,ns,gs,G,n; ap_pt *ARR,BS=P,GS,T,U; double *B,*S,u,t,x; llint K;
111 if (DEBUG) printf("ec_ap_large_bsgs a4:%f x:%f y:%f p:%f\n",a4,P.x,P.y,p);
112 ns=1+2*(int) Sqrt(4.0*p); bs=(int) ((1.0+Sqrt((double) ns))/2.0);
113 ARR=malloc(bs*sizeof(ap_pt));
114 NORMALISE(BS.y,p); NORMALISE(BS.x,p); ARR[0]=BS; ARR[0].i=1;
115 if (!MONTGOMERY_INVERSE_BABY_STEPS)
116 {for (i=1;i<bs;i++)
117 {if (ec_ap_large_pt_add(a4,BS,ARR[i-1],&(ARR[i]),p)) DONE(i+1);
118 if (ARR[i].y==0.0) DONE(2*i+2); ARR[i].i=i+1;}}
119 else
120 {for (i=1;i<=3;i++)
121 {if (ec_ap_large_pt_add(a4,BS,ARR[i-1],&(ARR[i]),p)) DONE(i+1);}
122 if (DEBUG) printf("Done with first baby steps, Total: %i\n",bs);
123 B=malloc((bs/2)*sizeof(double)); S=malloc((bs/2)*sizeof(double)); n=3;
124 while (n<bs-1)
125 {x=ARR[n].x; G=minimum(n,bs-n-1);
126 if (DEBUG) printf("Baby steps %i %f %f\n",n,ARR[n].x,ARR[n].y);
127 for (i=0;i<G;i++)
128 {B[i]=ARR[i].x-x; if (B[i]==0.0) DONEBS(i+n+2);}
129 if (G!=bs-n-1)
130 {if (ARR[n].y==0.0) DONEBS(2*n+2); B[n]=ARR[n].y+ARR[n].y; G++;}
131 if (DEBUG) printf("Doing Montgomery %i\n",G);
132 S[0]=B[0]; for (i=1;i<G;i++) MUL(S[i-1],B[i],S[i],p); INV(S[G-1],u,p);
133 for (i=G-1;i>0;i--) {MUL(u,S[i-1],t,p); MUL(u,B[i],u,p); B[i]=t;} B[0]=u;
134 if (DEBUG) printf("Montgomery done %i\n",bs-n-1);
135 for (i=0;i<G;i++) ADDI(ARR[i],ARR[n],&(ARR[n+i+1]),p,B[i]);
136 if (n+n+1<bs) DOUBLEI(a4,ARR[n],&(ARR[n+n+1]),p,B[n]); n=n+n+1;}
137 free(B); free(S);
138 for (i=0;i<bs;i++) {if (ARR[i].y==0.0) DONE(2*i+2); ARR[i].i=i+1;}}
139 if (DEBUG) printf("%i baby steps\n",bs);
140 if (DEBUG>=1) for (i=0;i<bs;i++) printf("bs:%i %f %f\n",i,ARR[i].x,ARR[i].y);
141 gs=2*bs; if (ec_ap_large_pt_double(a4,ARR[bs-1],&GS,p)) DONE(gs);
142 if (ec_ap_large_pt_mult(a4,P,L,&T,p)) DONE(L); ec_ap_large_qsort(ARR,bs);
143 if (DEBUG>=1) printf("T.x:%f T.y:%f\n",T.x,T.y);
144 HEAP_SEARCH(T.x,ARR,bs,k);
145 if (k!=-1) {i=ARR[k].i; if (T.y==ARR[k].y) K=L-i; else K=L+i; if (K) DONE(K);}
146 for (j=1;;j++)
147 {if (ec_ap_large_pt_add(a4,T,GS,&U,p)) DONE((L+(gs*j)));
148 T=U; HEAP_SEARCH(T.x,ARR,bs,k);
149 if (k!=-1) /* Don't use the GB trick, as it probably isn't faster */
150 {i=ARR[k].i; if (T.y==ARR[k].y) K=L-i; else K=L+i;
151 if (DEBUG) printf("FOUND %i %lli %i %i\n",i,K,j,gs); DONE((K+(gs*j)));}}
152 errorit("ec_bsgs_bsgs error"); return -1;}
153
ec_bsgs_ap_large(QD C4,QD C6,llint pr)154 llint ec_bsgs_ap_large(QD C4,QD C6,llint pr)
155 {llint aa,bb,v,U,L,S,WIDTH,m,o; QD AA,BB; int x; double V;
156 int k,ko; double a4,p=(double) pr; ap_pt P;
157 if (DEBUG) printf("ec_bsgs_ap_large %lli\n",pr);
158 if (pr<500) errorit("p<500 in ec_bsgs_ap\n");
159 QD_mul1(wmax,C4,-27.0,AA); aa=QD_modll(AA,p);
160 QD_mul1(wmax,C6,-54.0,BB); bb=QD_modll(BB,p); v=bb;
161 S=(llint) Floor(Sqrt(4*p)); L=pr+1-S; U=pr+1+S; WIDTH=U-L+1;
162 if (DEBUG) printf("a4:%lli a6:%lli\n",aa,bb);
163 if (DEBUG) printf("L:%lli U:%lli W:%lli\n",L,U,WIDTH);
164 k=kronll(v,pr); ko=-k; x=0;
165 while (1)
166 {while (!k || (k==ko)) {x++; v=((aa+x*x)*x+bb)%pr; k=kronll(v,pr);}
167 ko=k; if (DEBUG) printf("x:%i k:%i v:%lli\n",x,k,v); P.x=(double) ((x*v)%pr);
168 V=(double) v; MUL(V,V,P.y,p); MUL(P.y,(double) aa,a4,p); if (a4+a4>p) a4-=p;
169 if (DEBUG) printf("a4:%f Px:%f Py:%f\n",a4,P.x,P.y);
170 m=ec_ap_large_bsgs(a4,P,p,L); if (DEBUG) printf("mult of order is %lli\n",m);
171 o=ec_ap_large_order_from_mult(a4,P,m,p); if (DEBUG) printf("ord %lli\n",o);
172 if (o<WIDTH) continue; m=o*(L/o); if (m!=L) m+=o; m-=(pr+1);
173 if (k==1) m=-m; return m;}}
174