1 /* Author: Michael Scott */
2 /* Date: Dec 2007        */
3 /* Even Faster Duursma-Lee char 2 Tate pairing based on eta_T pairing */
4 /* See MIRACL dl2.cpp for more readable C++ version */
5 /* cl /O2 etat271.c miracl.lib  */
6 /* 8-bit version */
7 /* Half sized loop so nearly twice as fast! */
8 
9 /* MIRACL mirdef.h
10  * For Atmel AVR (e.g. ATmega128L) set up mirdef.h as follows
11 
12 #define MR_LITTLE_ENDIAN
13 #define MIRACL 8
14 #define mr_utype char
15 #define MR_IBITS 16
16 #define MR_LBITS 32
17 #define mr_unsign32 unsigned long
18 #define mr_dltype int
19 #define MR_STATIC 34
20 #define MR_ALWAYS_BINARY
21 #define MR_NOASM
22 #define MR_STRIPPED_DOWN
23 #define MR_GENERIC_MT
24 #define MAXBASE ((mr_small)1<<(MIRACL-1))
25 #define MR_BITSINCHAR 8
26 #define MR_NOKOBLITZ
27 #define MR_NO_STANDARD_IO
28 #define MR_NO_FILE_IO
29 #define MR_SIMPLE_BASE
30 #define MR_SIMPLE_IO
31 #define MR_AVR
32 #define SP271
33 
34 */
35 
36 /* use this mirdef.h to mimic 8-bit implementation on a PC
37 #define MR_LITTLE_ENDIAN
38 #define MIRACL 8
39 #define mr_utype char
40 #define MR_IBITS 32
41 #define MR_LBITS 32
42 #define mr_unsign32 unsigned int
43 #define mr_dltype short
44 #define MR_STATIC 34
45 #define MR_ALWAYS_BINARY
46 #define MR_NOASM
47 #define MR_STRIPPED_DOWN
48 #define MR_GENERIC_MT
49 #define MAXBASE ((mr_small)1<<(MIRACL-1))
50 #define MR_BITSINCHAR 8
51 #define MR_NOKOBLITZ
52 
53 */
54 
55 /* rem build using this batch file for PC
56 rem Compile MIRACL modules
57 cl /c /O2 /W3 mrcore.c
58 cl /c /O2 /W3 mrarth0.c
59 cl /c /O2 /W3 mrarth1.c
60 cl /c /O2 /W3 mrio1.c
61 cl /c /O2 /W3 mrbits.c
62 cl /c /O2 /W3 mrgf2m.c
63 cl /c /O2 /W3 mrec2m.c
64 
65 rem
66 rem Create library 'miracl.lib'
67 del miracl.lib
68 
69 lib /OUT:miracl.lib mrio1.obj
70 lib /OUT:miracl.lib miracl.lib mrbits.obj
71 lib /OUT:miracl.lib miracl.lib mrarth0.obj mrarth1.obj mrcore.obj
72 lib /OUT:miracl.lib miracl.lib mrec2m.obj mrgf2m.obj
73 del mr*.obj
74 
75 cl /O2 etat271.c miracl.lib
76 
77 On the ARM use a header like
78 
79 #define MR_LITTLE_ENDIAN
80 #define MIRACL 32
81 #define mr_utype int
82 #define MR_IBITS 32
83 #define MR_LBITS 32
84 #define mr_unsign32 unsigned int
85 #define mr_dltype long long
86 #define MR_STATIC 9
87 #define MR_ALWAYS_BINARY
88 #define MR_NOASM
89 #define MR_STRIPPED_DOWN
90 #define MR_GENERIC_MT
91 #define MAXBASE ((mr_small)1<<(MIRACL-1))
92 #define MR_BITSINCHAR 8
93 #define MR_NOKOBLITZ
94 
95 
96 /* define one curve or the other.. */
97 
98 #include <stdio.h>
99 #include <string.h>
100 #include "miracl.h"
101 
102 #define M 271
103 #define T 207
104 #define U 175
105 #define V 111
106 
107 #define B 0
108 #define TYPE 1
109 
110 /* points P and Q from ROM */
111 
112 /* WORDS = number of words needs to store GF(2^m) = size of bigs */
113 
114 /* elements of GF(2^m) are stored in bigs */
115 /* elements of the quartic extension field GF(2^{4m}) are stored as an array of 4 bigs */
116 /* = {a,b,c,d} = d.X^3+c.X^2+b.X+a */
117 
118 /* fast inlined addition code */
119 
120 #if MIRACL==64
121 
122 #define WORDS 5
123 #define NPW 16 /* nibbles per word */
124 #define ROMSZ 20
125 
126 static const mr_small rom[]={
127 0x591B401498D66271,0xA16F0C4E5357F2F6,0xD76AEF912696E510,0x75C041258C778D1D,0x10B1,
128 0x80DC7F385B9C26BF,0x2B65C2A7BAF3B9FD,0x6A84C19620F8D8B9,0x6D0DB856E16E7097,0x7C02,
129 0x4EDF428FD0EE2151,0x8A4509E6D6013138,0xBB5FBE66F7C468E7,0xA2740AF91652325E,0x2C67,
130 0x329B869A3E833026,0xB3716EC7D5F80608,0x3EE35C892B03AE59,0x5AF93E7449ABB134,0x48FB
131 };
132 
fincr2(big a,big c)133 void fincr2(big a,big c)
134 {
135     mr_small *aa,*cc;
136     aa=a->w; cc=c->w;
137 
138     cc[0]^=aa[0];
139     cc[1]^=aa[1];
140     cc[2]^=aa[2];
141     cc[3]^=aa[3];
142     cc[4]^=aa[4];
143 
144     c->len=WORDS;
145     if (cc[4]==0) mr_lzero(c);
146 }
147 
fadd2(big a,big b,big c)148 void fadd2(big a,big b,big c)
149 {
150     mr_small *aa,*bb,*cc;
151     aa=a->w; bb=b->w; cc=c->w;
152 
153     cc[0]=aa[0]^bb[0];
154     cc[1]=aa[1]^bb[1];
155     cc[2]=aa[2]^bb[2];
156     cc[3]=aa[3]^bb[3];
157     cc[4]=aa[4]^bb[4];
158 
159     c->len=WORDS;
160     if (cc[4]==0) mr_lzero(c);
161 }
162 
163 /* fast inlined copy code - replaces copy(.) */
164 
fcopy2(big a,big b)165 void fcopy2(big a,big b)
166 {
167     mr_small *aa,*bb;
168     aa=a->w; bb=b->w;
169 
170     bb[0]=aa[0];
171     bb[1]=aa[1];
172     bb[2]=aa[2];
173     bb[3]=aa[3];
174     bb[4]=aa[4];
175 
176     b->len=a->len;
177 }
178 
179 
180 #endif
181 
182 #if MIRACL==32
183 
184 #define WORDS 9
185 #define NPW 8 /* nibbles per word */
186 #define ROMSZ 36
187 
188 static const mr_small rom[]={
189 0x98D66271,0x591B4014,0x5357F2F6,0xA16F0C4E,0x2696E510,0xD76AEF91,0x8C778D1D,0x75C04125,0x10B1,
190 0x5B9C26BF,0x80DC7F38,0xBAF3B9FD,0x2B65C2A7,0x20F8D8B9,0x6A84C196,0xE16E7097,0x6D0DB856,0x7C02,
191 0xD0EE2151,0x4EDF428F,0xD6013138,0x8A4509E6,0xF7C468E7,0xBB5FBE66,0x1652325E,0xA2740AF9,0x2C67,
192 0x3E833026,0x329B869A,0xD5F80608,0xB3716EC7,0x2B03AE59,0x3EE35C89,0x49ABB134,0x5AF93E74,0x48FB
193 };
194 
fincr2(big a,big c)195 void fincr2(big a,big c)
196 {
197     mr_small *aa,*cc;
198     aa=a->w; cc=c->w;
199 
200     cc[0]^=aa[0];
201     cc[1]^=aa[1];
202     cc[2]^=aa[2];
203     cc[3]^=aa[3];
204     cc[4]^=aa[4];
205     cc[5]^=aa[5];
206     cc[6]^=aa[6];
207     cc[7]^=aa[7];
208     cc[8]^=aa[8];
209 
210 
211     c->len=WORDS;
212     if (cc[8]==0) mr_lzero(c);
213 }
214 
fadd2(big a,big b,big c)215 void fadd2(big a,big b,big c)
216 {
217     mr_small *aa,*bb,*cc;
218     aa=a->w; bb=b->w; cc=c->w;
219 
220     cc[0]=aa[0]^bb[0];
221     cc[1]=aa[1]^bb[1];
222     cc[2]=aa[2]^bb[2];
223     cc[3]=aa[3]^bb[3];
224     cc[4]=aa[4]^bb[4];
225     cc[5]=aa[5]^bb[5];
226     cc[6]=aa[6]^bb[6];
227     cc[7]=aa[7]^bb[7];
228     cc[8]=aa[8]^bb[8];
229 
230     c->len=WORDS;
231     if (cc[8]==0) mr_lzero(c);
232 }
233 
234 /* fast inlined copy code - replaces copy(.) */
235 
fcopy2(big a,big b)236 void fcopy2(big a,big b)
237 {
238     mr_small *aa,*bb;
239     aa=a->w; bb=b->w;
240 
241     bb[0]=aa[0];
242     bb[1]=aa[1];
243     bb[2]=aa[2];
244     bb[3]=aa[3];
245     bb[4]=aa[4];
246     bb[5]=aa[5];
247     bb[6]=aa[6];
248     bb[7]=aa[7];
249     bb[8]=aa[8];
250 
251     b->len=a->len;
252 }
253 
254 #endif
255 
256 #if MIRACL==8
257 
258 #define WORDS 34
259 #define NPW 2
260 #define ROMSZ 136
261 
262 /* For Pentanomial x^271+x^207+x^175+x^111+1 */
263 
264 #ifdef MR_AVR
265 __attribute__((__progmem__))
266 #endif
267 static const mr_small rom[]={
268 0x71,0x62,0xD6,0x98,0x14,0x40,0x1B,0x59,0xF6,0xF2,0x57,0x53,0x4E,0xC,0x6F,0xA1,0x10,0xE5,0x96,0x26,0x91,0xEF,0x6A,0xD7,0x1D,0x8D,0x77,0x8C,0x25,0x41,0xC0,0x75,0xB1,0x10,
269 0xBF,0x26,0x9C,0x5B,0x38,0x7F,0xDC,0x80,0xFD,0xB9,0xF3,0xBA,0xA7,0xC2,0x65,0x2B,0xB9,0xD8,0xF8,0x20,0x96,0xC1,0x84,0x6A,0x97,0x70,0x6E,0xE1,0x56,0xB8,0xD,0x6D,0x2,0x7C,
270 0x51,0x21,0xEE,0xD0,0x8F,0x42,0xDF,0x4E,0x38,0x31,0x1,0xD6,0xE6,0x9,0x45,0x8A,0xE7,0x68,0xC4,0xF7,0x66,0xBE,0x5F,0xBB,0x5E,0x32,0x52,0x16,0xF9,0xA,0x74,0xA2,0x67,0x2C,
271 0x26,0x30,0x83,0x3E,0x9A,0x86,0x9B,0x32,0x8,0x6,0xF8,0xD5,0xC7,0x6E,0x71,0xB3,0x59,0xAE,0x3,0x2B,0x89,0x5C,0xE3,0x3E,0x34,0xB1,0xAB,0x49,0x74,0x3E,0xF9,0x5A,0xFB,0x48
272 };
273 
fincr2(big a,big c)274 void fincr2(big a,big c)
275 {
276     mr_small *aa,*cc;
277     aa=a->w; cc=c->w;
278 
279     cc[0]^=aa[0];
280     cc[1]^=aa[1];
281     cc[2]^=aa[2];
282     cc[3]^=aa[3];
283     cc[4]^=aa[4];
284     cc[5]^=aa[5];
285     cc[6]^=aa[6];
286     cc[7]^=aa[7];
287     cc[8]^=aa[8];
288     cc[9]^=aa[9];
289     cc[10]^=aa[10];
290     cc[11]^=aa[11];
291     cc[12]^=aa[12];
292     cc[13]^=aa[13];
293     cc[14]^=aa[14];
294     cc[15]^=aa[15];
295     cc[16]^=aa[16];
296     cc[17]^=aa[17];
297     cc[18]^=aa[18];
298     cc[19]^=aa[19];
299     cc[20]^=aa[20];
300     cc[21]^=aa[21];
301     cc[22]^=aa[22];
302     cc[23]^=aa[23];
303     cc[24]^=aa[24];
304     cc[25]^=aa[25];
305     cc[26]^=aa[26];
306     cc[27]^=aa[27];
307     cc[28]^=aa[28];
308     cc[29]^=aa[29];
309     cc[30]^=aa[30];
310     cc[31]^=aa[31];
311     cc[32]^=aa[32];
312     cc[33]^=aa[33];
313 
314     c->len=WORDS;
315     if (cc[33]==0) mr_lzero(c);
316 }
317 
fadd2(big a,big b,big c)318 void fadd2(big a,big b,big c)
319 {
320     mr_small *aa,*bb,*cc;
321     aa=a->w; bb=b->w; cc=c->w;
322 
323     cc[0]=aa[0]^bb[0];
324     cc[1]=aa[1]^bb[1];
325     cc[2]=aa[2]^bb[2];
326     cc[3]=aa[3]^bb[3];
327     cc[4]=aa[4]^bb[4];
328     cc[5]=aa[5]^bb[5];
329     cc[6]=aa[6]^bb[6];
330     cc[7]=aa[7]^bb[7];
331     cc[8]=aa[8]^bb[8];
332     cc[9]=aa[9]^bb[9];
333     cc[10]=aa[10]^bb[10];
334     cc[11]=aa[11]^bb[11];
335     cc[12]=aa[12]^bb[12];
336     cc[13]=aa[13]^bb[13];
337     cc[14]=aa[14]^bb[14];
338     cc[15]=aa[15]^bb[15];
339     cc[16]=aa[16]^bb[16];
340     cc[17]=aa[17]^bb[17];
341     cc[18]=aa[18]^bb[18];
342     cc[19]=aa[19]^bb[19];
343     cc[20]=aa[20]^bb[20];
344     cc[21]=aa[21]^bb[21];
345     cc[22]=aa[22]^bb[22];
346     cc[23]=aa[23]^bb[23];
347     cc[24]=aa[24]^bb[24];
348     cc[25]=aa[25]^bb[25];
349     cc[26]=aa[26]^bb[26];
350     cc[27]=aa[27]^bb[27];
351     cc[28]=aa[28]^bb[28];
352     cc[29]=aa[29]^bb[29];
353     cc[30]=aa[30]^bb[30];
354     cc[31]=aa[31]^bb[31];
355     cc[32]=aa[32]^bb[32];
356     cc[33]=aa[33]^bb[33];
357 
358     c->len=WORDS;
359     if (cc[33]==0) mr_lzero(c);
360 }
361 
362 /* fast inlined copy code - replaces copy(.) */
363 
fcopy2(big a,big b)364 void fcopy2(big a,big b)
365 {
366     mr_small *aa,*bb;
367     aa=a->w; bb=b->w;
368 
369     bb[0]=aa[0];
370     bb[1]=aa[1];
371     bb[2]=aa[2];
372     bb[3]=aa[3];
373     bb[4]=aa[4];
374     bb[5]=aa[5];
375     bb[6]=aa[6];
376     bb[7]=aa[7];
377     bb[8]=aa[8];
378     bb[9]=aa[9];
379     bb[10]=aa[10];
380     bb[11]=aa[11];
381     bb[12]=aa[12];
382     bb[13]=aa[13];
383     bb[14]=aa[14];
384     bb[15]=aa[15];
385     bb[16]=aa[16];
386     bb[17]=aa[17];
387     bb[18]=aa[18];
388     bb[19]=aa[19];
389     bb[20]=aa[20];
390     bb[21]=aa[21];
391     bb[22]=aa[22];
392     bb[23]=aa[23];
393     bb[24]=aa[24];
394     bb[25]=aa[25];
395     bb[26]=aa[26];
396     bb[27]=aa[27];
397     bb[28]=aa[28];
398     bb[29]=aa[29];
399     bb[30]=aa[30];
400     bb[31]=aa[31];
401     bb[32]=aa[32];
402     bb[33]=aa[33];
403 
404     b->len=a->len;
405 }
406 
407 #endif
408 
409 /* Use internal workspace variables w1-w13 - must be careful doing this! - see comment below */
410 
mul(_MIPD_ big * a,big * b,big * r)411 void mul(_MIPD_ big *a,big *b,big *r)
412 {
413     /* Special multiplier for GF(2^{4m}) values of the form (x,y,y+1,0) */
414 
415     fcopy2(a[1],mr_mip->w2);
416     fcopy2(b[1],mr_mip->w3);
417     fadd2(a[1],a[0],mr_mip->w8);    /* e=w+p */
418     fadd2(b[1],b[0],mr_mip->w9);    /* s=t+q */
419 
420     /* only 3 modmults.. */
421 
422     modmult2(_MIPP_ mr_mip->w9,mr_mip->w8,mr_mip->w9);      /* z=(w+p)*(t+q) */
423     modmult2(_MIPP_ mr_mip->w3,mr_mip->w2,mr_mip->w4);      /* tw=t*w */
424     modmult2(_MIPP_ a[0],b[0],mr_mip->w8);            /* pq=p*q */
425     fincr2(mr_mip->w4,mr_mip->w9);                    /* z+=tw  */
426     fincr2(mr_mip->w8,mr_mip->w9);                    /* z+=pq  */
427     fincr2(mr_mip->w3,mr_mip->w2);                    /* w+=t   */
428     fadd2(mr_mip->w2,mr_mip->w4,mr_mip->w3);          /* t=w+tw */
429     incr2(mr_mip->w3,1,mr_mip->w3);                   /* t=w+tw+1  */
430 
431     fadd2(mr_mip->w9,a[0],mr_mip->w12);            /* x=z+p     */
432     fincr2(b[0],mr_mip->w12);                      /* x=z+p+q   */
433 
434     fadd2(mr_mip->w8,mr_mip->w3,r[0]);                /* r[0]=pq+t */
435     fadd2(mr_mip->w9,mr_mip->w3,r[1]);                /* r[1]=z+t  */
436     fadd2(mr_mip->w12,mr_mip->w4,r[2]);               /* r[2]=z+p+q+tw */
437     fcopy2(mr_mip->w2,r[3]);                       /* r[3]=w    */
438 }
439 
440 /* squaring GF(2^{4m}) values */
441 
square4(_MIPD_ big * a,big * c)442 void square4(_MIPD_ big *a,big *c)
443 {
444     if (a!=c)
445     {
446         fcopy2(a[0],c[0]);
447         fcopy2(a[1],c[1]);
448         fcopy2(a[2],c[2]);
449         fcopy2(a[3],c[3]);
450     }
451 
452     modsquare2(_MIPP_ c[3],c[3]);
453     fcopy2(c[2],mr_mip->w1);
454     modsquare2(_MIPP_ mr_mip->w1,mr_mip->w1);
455     fcopy2(c[1],c[2]);
456     modsquare2(_MIPP_ c[2],c[2]);
457     modsquare2(_MIPP_ c[0],c[0]);
458     fincr2(c[3],c[2]);
459     fincr2(mr_mip->w1,c[0]);
460     fcopy2(mr_mip->w1,c[1]);
461 
462     return;
463 }
464 
465 /* multiplying general GF(2^{4m}) values */
466 /* Uses karatsuba - 9 modmults - very time critical */
467 /* Use internal workspace variables w1-w13 - must be careful doing this! */
468 /* The thing is to ensure that none of the invoked miracl internal routines are using the same variables */
469 /* So first check the miracl source code.... I did... Its OK ... */
470 
mult4(_MIPD_ big * a,big * b,big * c)471 void mult4(_MIPD_ big *a,big *b,big *c)
472 {
473     fadd2(a[1],a[3],mr_mip->w3);
474     fadd2(a[0],a[2],mr_mip->w4);
475     fadd2(b[1],b[3],mr_mip->w8);
476     fadd2(b[0],b[2],mr_mip->w9);
477 
478     modmult2(_MIPP_ mr_mip->w8,mr_mip->w3,mr_mip->w10);
479     modmult2(_MIPP_ mr_mip->w9,mr_mip->w4,mr_mip->w11);
480     modmult2(_MIPP_ a[1],b[1],mr_mip->w2);
481     modmult2(_MIPP_ a[2],b[2],mr_mip->w1);
482 
483     fadd2(mr_mip->w2,mr_mip->w1,mr_mip->w13);
484 
485     fadd2(a[1],a[0],c[1]);
486     fadd2(b[0],b[1],mr_mip->w12);
487     modmult2(_MIPP_ c[1],mr_mip->w12,c[1]);
488     modmult2(_MIPP_ a[0],b[0],c[0]);
489     fincr2(c[0],c[1]);
490     fincr2(mr_mip->w2,c[1]);
491 
492     fcopy2(a[2],c[2]);
493     fadd2(a[2],a[3],mr_mip->w12);
494     fadd2(b[2],b[3],mr_mip->w2);
495     modmult2(_MIPP_ mr_mip->w12,mr_mip->w2,mr_mip->w12);
496 
497     fincr2(mr_mip->w12,mr_mip->w1);
498     modmult2(_MIPP_ a[3],b[3],mr_mip->w2);
499     fincr2(mr_mip->w2,mr_mip->w1);
500 
501     fadd2(mr_mip->w9,mr_mip->w8,mr_mip->w12);
502     fcopy2(mr_mip->w12,c[3]);
503     fadd2(mr_mip->w4,mr_mip->w3,mr_mip->w12);
504     modmult2(_MIPP_ c[3],mr_mip->w12,c[3]);
505 
506     fincr2(mr_mip->w2,mr_mip->w1);
507 
508     fincr2(mr_mip->w10,c[3]);
509     fincr2(mr_mip->w11,c[3]);
510     fincr2(mr_mip->w1,c[3]);
511     fincr2(c[1],c[3]);
512 
513     fadd2(mr_mip->w13,c[0],c[2]);
514     fincr2(mr_mip->w11,c[2]);
515     fincr2(mr_mip->w1,c[2]);
516 
517     fincr2(mr_mip->w1,c[1]);
518     fincr2(mr_mip->w10,mr_mip->w13);
519     fincr2(mr_mip->w13,c[1]);
520 
521     fincr2(mr_mip->w2,mr_mip->w13);
522     fincr2(mr_mip->w13,c[0]);
523 
524     return;
525 }
526 
527 /* Frobenius action - GF(2^{4m}) ^ 2^m */
528 
powq(_MIPD_ big * x)529 void powq(_MIPD_ big *x)
530 {
531     fcopy2(x[1],mr_mip->w1);
532     fincr2(x[1],x[0]);
533     fcopy2(x[2],x[1]);
534     fincr2(x[3],x[1]);
535     fcopy2(mr_mip->w1,x[2]);
536 }
537 
538 /* return degree of GF(2^{4m}) polynomial */
539 
degree(big * c)540 int degree(big *c)
541 {
542     if (size(c[3])!=0) return 3;
543     if (size(c[2])!=0) return 2;
544     if (size(c[1])!=0) return 1;
545     return 0;
546 }
547 
548 /* Raise a GF(2^{4m}) to a power */
549 
power4(_MIPD_ big * a,big k,big * u)550 void power4(_MIPD_ big *a,big k,big *u)
551 {
552     int i,j,m,nb,n,nbw,nzs;
553     big u2[4],t[4][4];
554 #ifndef MR_STATIC
555     char *mem=(char *)memalloc(_MIPP_ 20);
556 #else
557     char mem[MR_BIG_RESERVE(20)];
558     memset(mem,0,MR_BIG_RESERVE(20));
559 #endif
560     m=0;
561     for (i=0;i<4;i++)
562         for (j=0;j<4;j++) t[i][j]=mirvar_mem(_MIPP_ mem,m++);
563 
564     u2[0]=mirvar_mem(_MIPP_ mem,16);
565     u2[1]=mirvar_mem(_MIPP_ mem,17);
566     u2[2]=mirvar_mem(_MIPP_ mem,18);
567     u2[3]=mirvar_mem(_MIPP_ mem,19);
568 
569     if(size(k)==0)
570     {
571         convert(_MIPP_ 1,u[0]);
572         zero(u[1]);
573         zero(u[2]);
574         zero(u[3]);
575 #ifndef MR_STATIC
576         memkill(_MIPP_ mem,20);
577 #else
578         memset(mem,0,MR_BIG_RESERVE(20));
579 #endif
580         return;
581     }
582 
583     for (i=0;i<4;i++)
584         fcopy2(a[i],u[i]);
585 
586     if (size(k)==1)
587     {
588 #ifndef MR_STATIC
589         memkill(_MIPP_ mem,20);
590 #else
591         memset(mem,0,MR_BIG_RESERVE(20));
592 #endif
593         return;
594     }
595 
596     square4(_MIPP_ u,u2);
597     for (i=0;i<4;i++)
598         fcopy2(u[i],t[0][i]);
599 
600     for(i=1;i<4;i++)
601         mult4(_MIPP_ u2,t[i-1],t[i]);
602 
603     nb=logb2(_MIPP_ k);
604     if(nb>1) for(i=nb-2;i>=0;)
605     {
606         n=mr_window(_MIPP_ k,i,&nbw,&nzs,3); /* small 3 bit window to save RAM */
607         for(j=0;j<nbw;j++)
608             square4(_MIPP_ u,u);
609         if(n>0)	mult4(_MIPP_ u,t[n/2],u);
610         i-=nbw;
611         if(nzs!=0)
612         {
613             for (j=0;j<nzs;j++) square4(_MIPP_ u,u);
614             i-=nzs;
615         }
616     }
617 
618 #ifndef MR_STATIC
619     memkill(_MIPP_ mem,20);
620 #else
621     memset(mem,0,MR_BIG_RESERVE(20));
622 #endif
623 }
624 
625 /* Inversion of GF(2^{4m}) - fast but not time critical */
626 
invert(_MIPD_ big * x)627 void invert(_MIPD_ big *x)
628 {
629     int degF,degG,degB,degC,d,i,j;
630     big alpha,beta,gamma,BB[4],FF[5],CC[4],GG[5],t;
631     big *BP=BB,*CP=CC,*FP=FF,*GP=GG,*TP;
632 #ifndef MR_STATIC
633     char *mem=(char *)memalloc(_MIPP_ 22);
634 #else
635     char mem[MR_BIG_RESERVE(22)];       /* 972 bytes from the stack */
636     memset(mem,0,MR_BIG_RESERVE(22));
637 #endif
638 
639     alpha=mirvar_mem(_MIPP_ mem,0);
640     beta=mirvar_mem(_MIPP_ mem,1);
641     gamma=mirvar_mem(_MIPP_ mem,2);
642     t=mirvar_mem(_MIPP_ mem,3);
643 
644     for (i=0;i<4;i++)
645         BB[i]=mirvar_mem(_MIPP_ mem,4+i);
646     for (i=0;i<5;i++)
647         FF[i]=mirvar_mem(_MIPP_ mem,8+i);
648     for (i=0;i<4;i++)
649         CC[i]=mirvar_mem(_MIPP_ mem,13+i);
650     for (i=0;i<5;i++)
651         GG[i]=mirvar_mem(_MIPP_ mem,17+i);
652 
653     convert(_MIPP_ 1,CP[0]);
654     convert(_MIPP_ 1,FP[0]);
655     convert(_MIPP_ 1,FP[1]);
656     convert(_MIPP_ 1,FP[4]);   /* F = x^4+x+1 - irreducible polynomial for GF(2^{4m}) */
657 
658     degF=4; degG=degree(x); degC=0; degB=-1;
659 
660     if (degG==0)
661     {
662         inverse2(_MIPP_ x[0],x[0]);
663         return;
664     }
665 
666     for (i=0;i<4;i++)
667     {
668         fcopy2(x[i],GP[i]);
669         zero(x[i]);
670     }
671 
672     while (degF!=0)
673     {
674         if (degF<degG)
675         { /* swap */
676             TP=FP; FP=GP; GP=TP;  d=degF; degF=degG; degG=d;
677             TP=BP; BP=CP; CP=TP;  d=degB; degB=degC; degC=d;
678         }
679         j=degF-degG;
680         modsquare2(_MIPP_ GP[degG],alpha);
681         modmult2(_MIPP_ FP[degF],GP[degG],beta);
682         modmult2(_MIPP_ GP[degG],FP[degF-1],t);
683         modmult2(_MIPP_ FP[degF],GP[degG-1],gamma);
684         fincr2(t,gamma);
685         zero(t);
686 
687         for (i=0;i<=degF;i++ )
688         {
689             modmult2(_MIPP_ FP[i],alpha,FP[i]);
690             if (i>=j-1)
691             {
692                 modmult2(_MIPP_ gamma,GP[i-j+1],t);
693                 fincr2(t,FP[i]);
694             }
695             if (i>=j)
696             {
697                 modmult2(_MIPP_ beta,GP[i-j],t);
698                 fincr2(t,FP[i]);
699             }
700         }
701         for (i=0;i<=degB || i<=(degC+j);i++)
702         {
703             modmult2(_MIPP_ BP[i],alpha,BP[i]);
704             if (i>=j-1)
705             {
706                 modmult2(_MIPP_ gamma,CP[i-j+1],t);
707                 fincr2(t,BP[i]);
708             }
709             if (i>=j)
710             {
711                 modmult2(_MIPP_ beta,CP[i-j],t);
712                 fincr2(t,BP[i]);
713             }
714         }
715 
716         while (degF>=0 && size(FP[degF])==0) degF--;
717         if (degF==degG)
718         {
719             fcopy2(FP[degF],alpha);
720             for (i=0;i<=degF;i++)
721             {
722                 modmult2(_MIPP_ FP[i],GP[degF],FP[i]);
723                 modmult2(_MIPP_ alpha,GP[i],t);
724                 fincr2(t,FP[i]);
725             }
726             for (i=0;i<=4-degF;i++)
727             {
728                 modmult2(_MIPP_ BP[i],GP[degF],BP[i]);
729                 modmult2(_MIPP_ CP[i],alpha,t);
730                 fincr2(t,BP[i]);
731 
732             }
733             while (degF>=0 && size(FP[degF])==0) degF--;
734         }
735         degB=3;
736         while (degB>=0 && size(BP[degB])==0) degB--;
737     }
738 
739     inverse2(_MIPP_ FP[0],alpha);
740 
741     for (i=0;i<=degB;i++)
742         modmult2(_MIPP_ alpha,BP[i],x[i]);
743 #ifndef MR_STATIC
744     memkill(_MIPP_ mem,22);
745 #else
746     memset(mem,0,MR_BIG_RESERVE(22));
747 #endif
748 }
749 
750 /* Tate Pairing - calculated from eta_T pairing */
751 
eta_T(_MIPD_ epoint * P,epoint * Q,big * f,big * res)752 void eta_T(_MIPD_ epoint *P,epoint *Q,big *f,big *res)
753 {
754     big xp,yp,xq,yq,t,d;
755     big miller[4],v[4],u[4];
756     int i,m=M;
757     int promptr;
758 #ifndef MR_STATIC
759     char *mem=(char *)memalloc(_MIPP_ 18);
760 #else
761     char mem[MR_BIG_RESERVE(18)];            /* 972 bytes from stack */
762     memset(mem,0,MR_BIG_RESERVE(18));
763 #endif
764 
765     xp=mirvar_mem(_MIPP_ mem,0);
766     yp=mirvar_mem(_MIPP_ mem,1);
767     xq=mirvar_mem(_MIPP_ mem,2);
768     yq=mirvar_mem(_MIPP_ mem,3);
769     t=mirvar_mem(_MIPP_ mem,4);
770     d=mirvar_mem(_MIPP_ mem,5);
771     for (i=0;i<4;i++) miller[i]=mirvar_mem(_MIPP_ mem,6+i);
772     for (i=0;i<4;i++) v[i]=mirvar_mem(_MIPP_ mem,10+i);
773     for (i=0;i<4;i++) u[i]=mirvar_mem(_MIPP_ mem,14+i);
774 
775     fcopy2(P->X,xp);
776     fcopy2(P->Y,yp);
777     fcopy2(Q->X,xq);
778     fcopy2(Q->Y,yq);
779 
780     incr2(xp,1,t);       /* t=xp+1 */
781 
782     fadd2(xp,xq,d);
783     incr2(d,1,d);        /* xp+xq+1 */
784     modmult2(_MIPP_ d,t,d); /* t*(xp+xq+1) */
785     fincr2(yp,d);
786     fincr2(yq,d);
787     incr2(d,B,d);
788     incr2(d,1,f[0]);     /* f[0]=t*(xp+xq+1)+yp+yq+B+1 */
789 
790     convert(_MIPP_ 1,miller[0]);
791 
792     fadd2(t,xq,f[2]);       /* f[2]=t+xq   */
793     incr2(f[2],1,f[1]);     /* f[1]=t+xq+1 */
794     zero(f[3]);
795     promptr=0;
796 
797     for (i=0;i<(m-1)/2;i+=2)
798     {
799         fcopy2(xp,t);
800         sqroot2(_MIPP_ xp,xp);
801         sqroot2(_MIPP_ yp,yp);
802         fadd2(xp,xq,d);
803         modmult2(_MIPP_ d,t,d);
804         fincr2(yp,d);
805         fincr2(yq,d);
806         fadd2(d,xp,v[0]);       /* v[0]=t*(xp+xq)+yp+yq+xp */
807         fadd2(t,xq,v[2]);       /* v[2]=t+xq   */
808         incr2(v[2],1,v[1]);     /* v[1]=t+xq+1 */
809         modsquare2(_MIPP_ xq,xq);  /* xq*=xq */
810         modsquare2(_MIPP_ yq,yq);  /* yp*=yp */
811 
812         fcopy2(xp,t);           /* same again - unlooped times 2 */
813 
814         sqroot2(_MIPP_ xp,xp);
815         sqroot2(_MIPP_ yp,yp);
816 
817         fadd2(xp,xq,d);
818         modmult2(_MIPP_ d,t,d);
819         fincr2(yp,d);
820         fincr2(yq,d);
821         fadd2(d,xp,u[0]);
822         fadd2(t,xq,u[2]);
823         incr2(u[2],1,u[1]);
824         modsquare2(_MIPP_ xq,xq);
825         modsquare2(_MIPP_ yq,yq);
826 
827         mul(_MIPP_ u,v,u);         /* fast mul */
828         mult4(_MIPP_ miller,u,miller);
829     }
830 
831     mult4(_MIPP_ miller,f,miller);
832 
833     for (i=0;i<4;i++)
834     {
835         fcopy2(miller[i],u[i]);
836         fcopy2(miller[i],v[i]);
837         fcopy2(miller[i],f[i]);
838     }
839 
840     /* final exponentiation */
841 
842     for (i=0;i<(m+1)/2;i++) square4(_MIPP_ u,u);  /* u*=u */
843     powq(_MIPP_ u);
844     powq(_MIPP_ f);
845     for(i=0;i<4;i++) fcopy2(f[i],v[i]);
846     powq(_MIPP_ f);
847     for(i=0;i<4;i++) fcopy2(f[i],res[i]);
848     powq(_MIPP_ f);
849     mult4(_MIPP_ f,u,f);
850     mult4(_MIPP_ f,miller,f);
851     mult4(_MIPP_ res,v,res);
852     powq(_MIPP_ u);
853     powq(_MIPP_ u);
854     mult4(_MIPP_ res,u,res);
855 
856     /* doing inversion here could kill the stack... */
857 
858 #ifndef MR_STATIC
859     memkill(_MIPP_ mem,18);
860 #else
861     memset(mem,0,MR_BIG_RESERVE(18));
862 #endif
863 }
864 
865 /* ... so do inversion here to avoid stack overflow */
866 
tate(_MIPD_ epoint * P,epoint * Q,big * res)867 void tate(_MIPD_ epoint *P,epoint *Q,big *res)
868 {
869     int i;
870     big f[4];
871 #ifndef MR_STATIC
872     char *mem=(char *)memalloc(_MIPP_ 4);
873 #else
874     char mem[MR_BIG_RESERVE(4)];            /* 160 bytes from stack */
875     memset(mem,0,MR_BIG_RESERVE(4));
876 #endif
877     for (i=0;i<4;i++) f[i]=mirvar_mem(_MIPP_ mem,i);
878 
879     eta_T(_MIPP_ P,Q,f,res);
880 
881     invert(_MIPP_ f);
882     mult4(_MIPP_ res,f,res);
883 
884 #ifndef MR_STATIC
885     memkill(_MIPP_ mem,4);
886 #else
887     memset(mem,0,MR_BIG_RESERVE(4));
888 #endif
889 }
890 
891 /* Max stack requirement < 4K */
892 
main()893 int main()
894 {
895     big a2,a6,bx,r;
896     big res[4];
897     epoint *P,*Q;
898 
899     int i,romptr;
900     miracl instance;                           /* sizeof(miracl)= 2000 bytes from the stack */
901 #ifndef MR_STATIC
902 #ifdef MR_GENERIC_MT
903     miracl *mr_mip=mirsys(WORDS*NPW,16);
904 #else
905     miracl *mr_mip=mirsys(WORDS*NPW,16);
906 #endif
907     char *mem=(char *)memalloc(_MIPP_ 8);
908     char *mem1=(char *)ecp_memalloc(_MIPP_ 2);
909 #else
910 #ifdef MR_GENERIC_MT
911     miracl *mr_mip=mirsys(&instance,MR_STATIC*NPW,16); /* size of bigs is fixed */
912 #else
913     miracl *mr_mip=mirsys(&instance,MR_STATIC*NPW,16);
914 #endif
915     char mem[MR_BIG_RESERVE(8)];               /* reserve space on the stack for 8 bigs */
916     char mem1[MR_ECP_RESERVE(2)];              /* reserve space on stack for 2 curve points */
917     memset(mem,0,MR_BIG_RESERVE(8));           /* clear this memory */
918     memset(mem1,0,MR_ECP_RESERVE(2));          /* ~668 bytes in all  */
919 #endif
920 
921     /* Initialise bigs */
922 
923     a2=mirvar_mem(_MIPP_ mem,0);
924     a6=mirvar_mem(_MIPP_ mem,1);
925     bx=mirvar_mem(_MIPP_ mem,2);
926     for (i=0;i<4;i++)
927         res[i]=mirvar_mem(_MIPP_ mem,3+i);
928     r=mirvar_mem(_MIPP_ mem,7);
929 
930     /* printf("ROM size= %d\n",sizeof(rom)+sizeof(prom)); */
931 #ifndef MR_NO_STANDARD_IO
932 #ifdef MR_STATIC
933     printf("n Bigs require n*%d+%d bytes\n",MR_SIZE,MR_SL);
934     printf("n Points require n*%d+%d bytes\n",MR_ESIZE,MR_SL);
935     printf("sizeof(miracl)= %d\n",sizeof(miracl));
936 #endif
937 #endif
938     /* Initialise Elliptic curve points */
939 
940     P=epoint_init_mem(_MIPP_ mem1,0);
941     Q=epoint_init_mem(_MIPP_ mem1,1);
942 
943     /* Initialise supersingular curve */
944 
945     convert(_MIPP_ 1,a2);
946     convert(_MIPP_ B,a6);
947 
948     /* The -M tells MIRACL that this is a supersingular curve */
949 
950     if (!ecurve2_init(_MIPP_ -M,T,U,V,a2,a6,FALSE,MR_PROJECTIVE))
951     {
952 #ifndef MR_NO_STANDARD_IO
953         printf("Problem with the curve\n");
954 #endif
955         return 0;
956     }
957 
958     /* Get P and Q from ROM */
959     /* These should have been multiplied by the cofactor 487805 = 5*97561 */
960     /* 487805 is a cofactor of the group order 2^271+2^136+1 */
961 
962     romptr=0;
963     init_point_from_rom(P,WORDS,rom,ROMSZ,&romptr);
964     init_point_from_rom(Q,WORDS,rom,ROMSZ,&romptr);
965 
966 #ifndef MR_NO_STANDARD_IO
967     printf( "P= \n");
968     otnum(_MIPP_ P->X,stdout);
969     otnum(_MIPP_ P->Y,stdout);
970     printf( "Q= \n");
971     otnum(_MIPP_ Q->X,stdout);
972     otnum(_MIPP_ Q->Y,stdout);
973 #endif
974 
975     bigbits(_MIPP_ 160,r);
976 
977     /* Simple bilinearity test */
978 
979     tate(_MIPP_ P,Q,res);
980 
981     /* this could break the 4k stack, 2060+668+2996 >4K    */
982     /* so we cannot afford much precomputation in power4   */
983 
984     power4(_MIPP_ res,r,res);   /* res=res^{sr} */
985 
986 #ifndef MR_NO_STANDARD_IO
987     printf( "\ne(P,Q)^r= \n");
988     for (i=0;i<4;i++)
989     {
990         otnum(_MIPP_ res[i],stdout);
991         zero(res[i]);
992     }
993 #endif
994 
995     ecurve2_mult(_MIPP_ r,Q,Q);   /* Q=rQ */
996 
997     epoint2_norm(_MIPP_ Q);
998 
999     tate(_MIPP_ P,Q,res);         /* Now invert is taken out of Tate, and the stack should be OK */
1000 
1001 #ifndef MR_NO_STANDARD_IO
1002     printf( "\ne(P,rQ)= \n");
1003     for (i=0;i<4;i++)
1004         otnum(_MIPP_ res[i],stdout);
1005 #endif
1006 
1007     /* all done */
1008 
1009 #ifndef MR_STATIC
1010     memkill(_MIPP_ mem,8);
1011     ecp_memkill(_MIPP_ mem1,2);
1012 #else
1013     memset(mem,0,MR_BIG_RESERVE(8));        /* clear this stack memory */
1014     memset(mem1,0,MR_ECP_RESERVE(2));
1015 #endif
1016 
1017     mirexit(_MIPPO_ );  /* clears workspace memory */
1018     return 0;
1019 }
1020 
1021 
1022