1 /* Even Faster Duursma-Lee char 2 Tate pairing based on eta_T pairing */
2 /* cl /O2 /GX dl2.cpp ec2.cpp gf2m4x.cpp gf2m.cpp big.cpp miracl.lib */
3 /* Half sized loop so nearly twice as fast! */
4 /* 14th March 2005 */
5
6 #include <iostream>
7 #include <ctime>
8
9 #include "gf2m.h"
10 #include "gf2m4x.h"
11 #include "ec2.h"
12
13 // set TYPE = 1 if B=0 && (M=1 or 7 mod 8), else TYPE = 2
14 // set TYPE = 1 if B=1 && (M=3 or 5 mod 8), else TYPE = 2
15
16 // some example curves to play with
17
18 //#define M 283
19 //#define T 249
20 //#define U 219
21 //#define V 27
22 //#define B 1
23 //#define TYPE 1
24 //#define CF 1
25
26 //#define M 367
27 //#define T 21
28 //#define U 0
29 //#define V 0
30 //#define B 1
31 //#define TYPE 2
32 //#define CF 1
33
34 //#define M 379
35 //#define T 317
36 //#define U 315
37 //#define V 283
38 //#define B 1
39 //#define TYPE 1
40 //#define CF 1
41
42 #define M 1223
43 #define T 255
44 #define U 0
45 #define V 0
46 #define B 0
47 #define TYPE 1
48 #define CF 5
49
50 //#define M 271
51 //#define T 207
52 //#define U 175
53 //#define V 111
54 //#define B 0
55 //#define TYPE 1
56 //#define CF 487805
57
58 //#define M 353
59 //#define B 1
60 //#define T 215
61 //#define U 0
62 //#define V 0
63 //#define TYPE 2
64 //#define CF 1
65
66 //#define M 271
67 //#define U 0
68 //#define V 0
69 //#define T 201
70 //#define B 0
71 //#define TYPE 1
72 //#define CF 487805
73
74 #define IMOD4 ((M+1)/2)%4
75
76 //#define XX (IMOD4%2)
77 //#define YY (IMOD4/2)
78 //#define NXX (1-XX)
79
80 using namespace std;
81
82 Miracl precision(40,0);
83
84 //
85 // Extract ECn point in internal GF2m format
86 //
87
extract(EC2 & A,GF2m & x,GF2m & y)88 void extract(EC2& A,GF2m& x,GF2m& y)
89 {
90 x=(A.get_point())->X;
91 y=(A.get_point())->Y;
92 }
93
94 //
95 // Tate Pairing - note miller -> miller variable
96 // Loop unrolled x2 for speed
97 //
98
tate(EC2 & P,EC2 & Q)99 GF2m4x tate(EC2& P,EC2& Q)
100 {
101 GF2m xp,yp,xq,yq,t;
102 GF2m4x miller,w,u,u0,u1,v,f,res;
103 int i,m=M;
104
105 normalise(P); normalise(Q);
106 extract(P,xp,yp);
107 extract(Q,xq,yq);
108
109 // first calculate the contribution of adding P or -P to 2^[(m+1)/2].P
110 //
111 // Note that 2^[(m+1)/2].Point(x,y) = Point(x^2+1,x^2+y^2) or something similar....
112 // Line slope is x or x+1 (thanks Steven!)
113 //
114 // Then the truncated loop, four flavours...
115
116 #if IMOD4 == 1
117 // (X=1)
118 // (Y=0)
119 t=xp; // 0 (X+1)
120 f.set(t*(xp+xq+1)+yq+yp+B,t+xq+1,t+xq,0); // 0 (Y)
121
122
123 miller=1;
124 for (i=0;i<(m-3)/2;i+=2)
125 {
126
127 t=xp+1; xp=sqrt(xp); yp=sqrt(yp); // 1 (X)
128 u0.set(t*(xp+xq+1)+yp+yq,t+xq+1,t+xq,0); // 1 0 (X) ((X+1)*(xp+1)+Y)
129 xq*=xq; yq*=yq;
130
131 t=xp+1; xp=sqrt(xp); yp=sqrt(yp);
132 u1.set(t*(xp+xq+1)+yp+yq,t+xq+1,t+xq,0);
133 xq*=xq; yq*=yq;
134
135 u=mul(u0,u1);
136 miller*=u;
137 }
138
139 // final step
140
141 t=xp+1; xp=sqrt(xp); yp=sqrt(yp);
142 u.set(t*(xp+xq+1)+yp+yq,t+xq+1,t+xq,0);
143 miller*=u;
144
145 #endif
146
147 #if IMOD4 == 0
148 // (X=0)
149 // (Y=0)
150 t=xp+1; // 1 (X+1)
151 f.set(t*(xq+xp+1)+yq+yp+B,t+xq+1,t+xq,0); // 0 (Y)
152 miller=1;
153
154 for (i=0;i<(m-1)/2;i+=2)
155 {
156 // loop is unrolled x 2
157 t=xp; xp=sqrt(xp); yp=sqrt(yp); // 0 (X)
158 u0.set(t*(xp+xq)+yp+yq+xp+1,t+xq+1,t+xq,0); // 0 xp+1 (X) ((X+1)*(xp+1)+Y
159 xq*=xq; yq*=yq;
160
161 t=xp; xp=sqrt(xp); yp=sqrt(yp);
162 u1.set(t*(xp+xq)+yp+yq+xp+1,t+xq+1,t+xq,0);
163 xq*=xq; yq*=yq;
164
165 u=mul(u0,u1);
166 miller*=u;
167 }
168
169 #endif
170
171 #if IMOD4 == 2
172 // (X=0) // (Y=1)
173 t=xp+1; // 1 (X+1)
174 f.set(t*(xq+xp+1)+yq+yp+B+1,t+xq+1,t+xq,0); // 1 (Y)
175 miller=1;
176 for (i=0;i<(m-1)/2;i+=2)
177 {
178
179 t=xp; xp=sqrt(xp); yp=sqrt(yp); // 0 (X)
180 u0.set(t*(xp+xq)+yp+yq+xp,t+xq+1,t+xq,0); // 0 xp+0 (X) ((X+1)*(xp+1)+Y)
181 xq*=xq; yq*=yq;
182
183 t=xp; xp=sqrt(xp); yp=sqrt(yp);
184 u1.set(t*(xp+xq)+yp+yq+xp,t+xq+1,t+xq,0);
185 xq*=xq; yq*=yq;
186
187 u=mul(u0,u1);
188 miller*=u;
189 }
190
191 #endif
192
193 #if IMOD4 == 3
194 // (X=1) // (Y=1)
195 t=xp; // 0 (X+1)
196 f.set(t*(xq+xp+1)+yq+yp+B+1,t+xq+1,t+xq,0); // 1 (Y)
197
198 miller=1;
199 for (i=0;i<(m-3)/2;i+=2)
200 {
201
202 t=xp+1; xp=sqrt(xp); yp=sqrt(yp); // 1 (X)
203 u0.set(t*(xp+xq+1)+yp+yq+1,t+xq+1,t+xq,0); // 1 1 (X) ((X+1)*(xp+1)+Y)
204 xq*=xq; yq*=yq;
205
206 t=xp+1; xp=sqrt(xp); yp=sqrt(yp);
207 u1.set(t*(xp+xq+1)+yp+yq+1,t+xq+1,t+xq,0);
208 xq*=xq; yq*=yq;
209
210 u=mul(u0,u1);
211 miller*=u;
212 }
213
214 // final step
215
216 t=xp+1; xp=sqrt(xp); yp=sqrt(yp);
217 u.set(t*(xp+xq+1)+yp+yq+1,t+xq+1,t+xq,0);
218 miller*=u;
219
220 #endif
221
222 miller*=f;
223
224 // raising to the power (2^m-2^[m+1)/2]+1)(2^[(m+1)/2]+1)(2^(2m)-1) (TYPE 2)
225 // or (2^m+2^[(m+1)/2]+1)(2^[(m+1)/2]-1)(2^(2m)-1) (TYPE 1)
226 // 6 Frobenius, 4 big field muls...
227
228 u=v=w=miller;
229 for (i=0;i<(m+1)/2;i++) u*=u;
230
231 #if TYPE == 1
232
233 u.powq();
234 w.powq();
235 v=w;
236 w.powq();
237 res=w;
238 w.powq();
239 w*=u;
240 w*=miller;
241 res*=v;
242 u.powq();
243 u.powq();
244 res*=u;
245
246 #else
247
248 u.powq();
249 v.powq();
250 w=u*v;
251 v.powq();
252 w*=v;
253 v.powq();
254 u.powq();
255 u.powq();
256 res=v*u;
257 res*=miller;
258
259 #endif
260
261 res/=w;
262 return res;
263 }
264
main()265 int main()
266 {
267 EC2 P,Q,W;
268 Big bx,s,r,x,y,order;
269 GF2m4x res;
270 time_t seed;
271 int i;
272 miracl *mip=&precision;
273
274 time(&seed);
275 irand((long)seed);
276
277 if (!ecurve2(-M,T,U,V,(Big)1,(Big)B,TRUE,MR_PROJECTIVE)) // -M indicates Super-Singular
278 {
279 cout << "Problem with the curve" << endl;
280 return 0;
281 }
282
283 // Curve order = 2^M+2^[(M+1)/2]+1 or 2^M-2^[(M+1)/2]+1 is nearly prime
284
285 cout << "IMOD4= " << IMOD4 << endl;
286 cout << "M%8= " << M%8 << endl;
287
288 forever
289 {
290 bx=rand(M,2);
291 if (P.set(bx,bx)) break;
292 }
293
294 forever
295 {
296 bx=rand(M,2);
297 if (Q.set(bx,bx)) break;
298 }
299
300 /* for (int i=0;i<10000;i++) */
301
302 P*=CF; // cofactor multiplication
303 Q*=CF;
304
305 // order=pow((Big)2,M)-pow((Big)2,(M+1)/2)+1;
306 // P*=order;
307 // cout << "P= " << P << endl;
308 // exit(0);
309 /*
310 mip->IOBASE=16;
311 cout << "P= " << P << endl;
312 cout << "Q= " << Q << endl;
313
314 Big ddd=pow((Big)2,32);
315 Big sx,sy;
316
317
318 P.get(x,y);
319 sx=x;
320 sy=y;
321
322 while (x>0)
323 {
324 cout << "0x" << hex << x%ddd << ",";
325 x/=ddd;
326 }
327 cout << endl;
328
329 while (y>0)
330 {
331 cout << "0x" << hex << y%ddd << ",";
332 y/=ddd;
333 }
334 cout << endl;
335
336 ddd=256;
337 x=sx;
338 y=sy;
339
340 while (x>0)
341 {
342 cout << "0x" << hex << x%ddd << ",";
343 x/=ddd;
344 }
345 cout << endl;
346
347 while (y>0)
348 {
349 cout << "0x" << hex << y%ddd << ",";
350 y/=ddd;
351 }
352 cout << endl;
353
354
355 Q.get(x,y);
356 sx=x;
357 sy=y;
358 ddd=pow((Big)2,32);
359
360 while (x>0)
361 {
362 cout << "0x" << hex << x%ddd << ",";
363 x/=ddd;
364 }
365 cout << endl;
366
367 while (y>0)
368 {
369 cout << "0x" << hex << y%ddd << ",";
370 y/=ddd;
371 }
372 cout << endl;
373
374 ddd=256;
375 x=sx;
376 y=sy;
377
378 while (x>0)
379 {
380 cout << "0x" << hex << x%ddd << ",";
381 x/=ddd;
382 }
383 cout << endl;
384
385 while (y>0)
386 {
387 cout << "0x" << hex << y%ddd << ",";
388 y/=ddd;
389 }
390 cout << endl;
391
392
393
394 exit(0);
395 */
396 //for (i=0;i<1000;i++)
397 res=tate(P,Q);
398
399 s=rand(256,2);
400 r=rand(256,2);
401 res=pow(res,s); res=pow(res,r);
402 //mip->IOBASE=16;
403 cout << "e(P,Q)^sr= " << res << endl;
404
405 P*=s;
406 Q*=r;
407
408 res=tate(Q,P);
409
410 cout << "e(sP,rQ)= " << res << endl;
411
412 //Big q=pow((Big)2,M)-pow((Big)2,(M+1)/2)+1;
413
414 //cout << pow(res,q) << endl;
415
416 return 0;
417 }
418