1 #include "sympow.h"
2 #define DEBUG (FALSE || GLOBAL_DEBUG)
3
4 typedef struct {int x,y; int i;} AP_pt;
5
6 #define MUL(x,y,z,p) z=(((llint) x)*((llint) y))%p
7 #define INV(x,i,p)\
8 {int rr=0,ss=1,cc,aa=p,bb=x; while (1)\
9 {cc=aa/bb; aa-=bb*cc;\
10 if (aa==0) {if (bb>0) i=ss; else i=-ss; break;} rr-=cc*ss;\
11 cc=bb/aa; bb-=aa*cc;\
12 if (bb==0) {if (aa>0) i=rr; else i=-rr; break;} ss-=cc*rr;}}
13 #define NORMALISE(x,p) while ((x)>=(p)) (x)-=(p); while ((x)<0) (x)+=(p)
14
15 #define FROM_LINE(P,Q,R,s,p)\
16 {int T,x1; AP_pt V;\
17 MUL(s,s,T,p); V.x=T-(P.x+Q.x); NORMALISE(V.x,p);\
18 x1=P.x-V.x; MUL(x1,s,T,p); V.y=T-P.y; NORMALISE(V.y,p); *(R)=V;}
19
ec_AP_pt_double(int a4,AP_pt P,AP_pt * R,int p)20 static int ec_AP_pt_double(int a4,AP_pt P,AP_pt *R,int p)
21 {int T,n,i,s,tt; /* NORMALISE(P.x,p); NORMALISE(P.y,p); */
22 if (P.y==0) return TRUE; MUL(P.x,P.x,T,p); n=T+T+T+a4; NORMALISE(n,p);
23 tt=P.y+P.y; INV(tt,i,p); MUL(n,i,s,p); FROM_LINE(P,P,R,s,p); return FALSE;}
24
ec_AP_pt_add(int a4,AP_pt P,AP_pt Q,AP_pt * R,int p)25 static int ec_AP_pt_add(int a4,AP_pt P,AP_pt Q,AP_pt *R,int p)
26 {int Ld,Ln,i,s; /* NORMALISE(P.x,p); NORMALISE(P.y,p); etc. */
27 if (P.x==Q.x)
28 {if (P.y==Q.y) return ec_AP_pt_double(a4,P,R,p); return TRUE;}
29 Ld=P.x-Q.x; Ln=P.y-Q.y; /* allow Ld to be negative */
30 INV(Ld,i,p); MUL(Ln,i,s,p); FROM_LINE(P,Q,R,s,p); return FALSE;}
31
32 #define ADDI(P,Q,R,p,i) {int s; MUL(P.y-Q.y,i,s,p); FROM_LINE(P,Q,R,s,p);}
33 #define DOUBLEI(a4,P,R,p,i)\
34 {int TT,nn,s; MUL(P.x,P.x,TT,p); nn=TT+TT+TT+a4;\
35 MUL(nn,i,s,p); FROM_LINE(P,P,R,s,p);}
36
ec_AP_pt_mult(int a4,AP_pt P,int n,AP_pt * R,int p)37 static int ec_AP_pt_mult(int a4,AP_pt P,int n,AP_pt *R,int p)
38 {int m=n; AP_pt Q=P,T,U; int z=TRUE; /* P is not O */
39 if (n==0) return TRUE; if (n==1) {*R=P; return FALSE;}
40 if (m&1) {U=Q; z=FALSE;}
41 while (m)
42 {m=m>>1; if (ec_AP_pt_double(a4,Q,&T,p)) break; Q=T;
43 if (m&1)
44 {if (z) {U=Q; z=FALSE;} else {z=ec_AP_pt_add(a4,Q,U,&T,p); if (!z) U=T;}}}
45 if (!z) *R=U; return z;}
46
ec_AP_pt_order_from_mult(int a4,AP_pt P,int m,int pr)47 static int ec_AP_pt_order_from_mult(int a4,AP_pt P,int m,int pr)
48 {int i,e,k,t; int o=1,n=m,p,pe; AP_pt Q=P,R; LIST f;
49 f.p[0]=0; ifactor(n,&f,1,1000); NORMALISE(Q.x,pr); NORMALISE(Q.y,pr);
50 for (i=0;f.p[i]!=0;i++)
51 {p=f.p[i]; e=f.e[i]; for (t=0;t<e;t++) n=n/p;
52 if (ec_AP_pt_mult(a4,Q,n,&R,pr)) continue;
53 for (k=1;k<e;k++) {if (ec_AP_pt_mult(a4,R,p,&R,pr)) break;}
54 e=k; pe=1; for (t=0;t<e;t++) pe=pe*p; o=o*pe;
55 if (ec_AP_pt_mult(a4,Q,pe,&R,pr)) break; Q=R;}
56 return o;}
57
58 #define TYPE AP_pt
59 #define NAME ec_AP_qsort
60 #define CMP(A,B) (A->x>B->x)
61 #define SWAP_VARS(type,A,B) {type _ttemmp; _ttemmp=A; A=B; B=_ttemmp;}
62 #define SWAP_DATA(A,B) {SWAP_VARS(int,A->x,B->x);\
63 SWAP_VARS(int,A->y,B->y); SWAP_VARS(int,A->i,B->i);}
64
NAME(TYPE * ARR,int n)65 static void NAME(TYPE *ARR,int n)
66 {TYPE *left,*right,*ip,*jp,*last,*mid; int l;
67 if (DEBUG) printf("ec_AP_qsort %i\n",n);
68 if (n<=1) return; left=ARR;
69 if (n==2) {ip=left+1; if (CMP(left,ip)) SWAP_DATA(left,ip); return;}
70 right=left+(n-1);
71 if (n<=5)
72 {for (ip=left;ip<right;ip++)
73 for (jp=ip+1;jp<=right;jp++) if (!CMP(jp,ip)) SWAP_DATA(ip,jp); return;}
74 mid=left+(n/2);
75 if (!CMP(left,mid))
76 {if (!CMP(mid,right)) SWAP_DATA(left,mid)
77 else if (!CMP(left,right)) SWAP_DATA(left,right)}
78 else
79 {if (CMP(mid,right)) SWAP_DATA(left,mid)
80 else if (CMP(left,right)) SWAP_DATA(left,right)}
81 l=0; last=left;
82 for (ip=left+1;ip<=right;ip++)
83 if (!CMP(ip,left)) {l++; last++; SWAP_DATA(last,ip)}
84 if (left!=last) SWAP_DATA(left,last)
85 NAME(left,l); NAME(last+1,n-l-1);}
86
87 #define HEAP_SEARCH(n,A,sz,k)\
88 {int LL=0,HH=sz,w;\
89 while (LL<HH) {w=(LL+HH)>>1; if (A[w].x<n) LL=w+1; else HH=w;}\
90 if (DEBUG>=2) printf("heap_search %i %i %i\n",HH,A[HH].x,n);\
91 if ((HH<sz) && (A[HH].x==n)) k=HH; else k=-1;}
92 #define DONE(x) {free(ARR); return x;}
93 #define DDONE(x) {if (mo==0) mo=(x); else {free(ARR); return gcd(mo,x);}}
94 #define DONEBS(x) {free(B); free(S); DONE(x);}
95 #define MONTGOMERY_INVERSE_BABY_STEPS TRUE
96
ec_AP_bsgs(int a4,AP_pt P,int p,int L)97 static llint ec_AP_bsgs(int a4,AP_pt P,int p,int L)
98 {int i,j,k,bs,ns,gs,G,n; AP_pt *ARR,BS=P,GS,T,U; int *B,*S,u,t,x,mo=0;
99 if (DEBUG) printf("ec_AP_bsgs a4:%i x:%i y:%i p:%i\n",a4,P.x,P.y,p);
100 ns=1+2*(int) Sqrt(4.0*(double) p); bs=(int) ((1.0+Sqrt((double) ns))/2.0);
101 ARR=malloc(bs*sizeof(AP_pt)); ARR[0]=BS; ARR[0].i=1;
102 if (!MONTGOMERY_INVERSE_BABY_STEPS)
103 {for (i=1;i<bs;i++)
104 {if (ec_AP_pt_add(a4,BS,ARR[i-1],&(ARR[i]),p)) DONE(i+1);
105 if (ARR[i].y==0) DONE(2*i+2); ARR[i].i=i+1;}}
106 else
107 {for (i=1;i<=3;i++)
108 {if (ec_AP_pt_add(a4,BS,ARR[i-1],&(ARR[i]),p)) DONE(i+1);}
109 if (DEBUG) printf("Done with first baby steps, Total: %i\n",bs);
110 B=malloc((bs/2)*sizeof(int)); S=malloc((bs/2)*sizeof(int)); n=3;
111 while (n<bs-1)
112 {x=ARR[n].x; G=minimum(n,bs-n-1);
113 if (DEBUG) printf("Baby steps %i %i %i\n",n,ARR[n].x,ARR[n].y);
114 for (i=0;i<G;i++)
115 {B[i]=ARR[i].x-x; if (B[i]==0) DONEBS(i+n+2);}
116 if (G!=bs-n-1)
117 {if (ARR[n].y==0) DONEBS(2*n+2); B[n]=ARR[n].y+ARR[n].y; G++;}
118 if (DEBUG) printf("Doing Montgomery %i\n",G);
119 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);
120 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;
121 if (DEBUG) printf("Montgomery done %i\n",bs-n-1);
122 for (i=0;i<G;i++) ADDI(ARR[i],ARR[n],&(ARR[n+i+1]),p,B[i]);
123 if (n+n+1<bs) DOUBLEI(a4,ARR[n],&(ARR[n+n+1]),p,B[n]); n=n+n+1;}
124 free(B); free(S);
125 for (i=0;i<bs;i++) {if (ARR[i].y==0) DONE(2*i+2); ARR[i].i=i+1;}
126 }
127 if (DEBUG) printf("%i baby steps\n",bs);
128 if (DEBUG>=1) for (i=0;i<bs;i++) printf("bs:%i %i %i\n",i,ARR[i].x,ARR[i].y);
129 gs=2*bs; if (ec_AP_pt_double(a4,ARR[bs-1],&GS,p)) DONE(gs);
130 if (ec_AP_pt_mult(a4,P,L,&T,p)) DONE(L); ec_AP_qsort(ARR,bs);
131 if (DEBUG>=1) printf("T.x:%i T.y:%i\n",T.x,T.y);
132 HEAP_SEARCH(T.x,ARR,bs,k);
133 if (k!=-1) {i=ARR[k].i; if (T.y==ARR[k].y) k=L-i; else k=L+i; if (k) DONE(k);}
134 for (j=1;j*(2*bs-1)<=ns;j++)
135 {if (ec_AP_pt_add(a4,T,GS,&U,p)) {DDONE(L+gs*j); j++; T=GS;}
136 else T=U; HEAP_SEARCH(T.x,ARR,bs,k);
137 if (k!=-1) /* idea of Geoff Bailey; avoid computing exact order */
138 {i=ARR[k].i; if (T.y==ARR[k].y) k=L-i; else k=L+i;
139 if (DEBUG) printf("FOUND %i %i %i %i\n",i,k,j,gs); DDONE(k+gs*j);}}
140 ASSERT(mo!=0); free(ARR); return -mo;}
141
ec_bsgs_ap(QD C4,QD C6,llint pr)142 llint ec_bsgs_ap(QD C4,QD C6,llint pr)
143 {int aa,bb; if (DEBUG) printf("ec_bsgs_ap %lli\n",pr);
144 if (pr>(1<<29)) return ec_bsgs_ap_large(C4,C6,pr);
145 aa=(int) ((-27*QD_modll(C4,pr))%pr); bb=(int) ((-54*QD_modll(C6,pr))%pr);
146 return ec_bsgs_ap_AB(aa,bb,pr);}
147
ec_bsgs_ap_AB(int aa,int bb,llint pr)148 llint ec_bsgs_ap_AB(int aa,int bb,llint pr)
149 {int a4,L,U,WIDTH,v,S,k,ko,x,m,o; AP_pt P;
150 S=(int) Floor(2.0*Sqrt(pr)); L=(int) pr+1-S; U=(int) pr+1+S; WIDTH=U-L+1;
151 if (DEBUG) printf("a4:%i a6:%i\n",aa,bb);
152 if (DEBUG) printf("L:%i U:%i W:%i\n",L,U,WIDTH);
153 v=bb; NORMALISE(aa,pr); NORMALISE(v,pr); k=kron(v,(int) pr); ko=-k; x=0;
154 while (1)
155 {while (!k || (k==ko))
156 {x++; v=(int) ((((llint) aa+(llint) (x*x))*(llint) x+(llint) bb)%pr);
157 k=kron(v,pr);}
158 ko=k; if (DEBUG) printf("x:%i k:%i v:%i\n",x,k,v);
159 MUL(x,v,P.x,pr); MUL(v,v,P.y,pr); MUL(aa,P.y,a4,pr);
160 NORMALISE(P.x,pr); NORMALISE(P.y,pr); NORMALISE(a4,pr);
161 if (DEBUG) printf("a4:%i Px:%i Py:%i\n",a4,P.x,P.y);
162 m=ec_AP_bsgs(a4,P,pr,L); if (DEBUG) printf("multiple of order is %i\n",m);
163 if (m>0) {if (m<WIDTH) continue; o=ec_AP_pt_order_from_mult(a4,P,m,pr);}
164 else o=-m; if (DEBUG) printf("ord is %i\n",o);
165 if (o<WIDTH) continue; m=o*(L/o); if (m!=L) m+=o; m-=((int) pr+1);
166 if (k==1) m=-m; return (llint) m;}}
167