1 
2 /***************************************************************************
3                                                                            *
4 Copyright 2012 CertiVox IOM Ltd.                                           *
5                                                                            *
6 This file is part of CertiVox MIRACL Crypto SDK.                           *
7                                                                            *
8 The CertiVox MIRACL Crypto SDK provides developers with an                 *
9 extensive and efficient set of cryptographic functions.                    *
10 For further information about its features and functionalities please      *
11 refer to http://www.certivox.com                                           *
12                                                                            *
13 * The CertiVox MIRACL Crypto SDK is free software: you can                 *
14   redistribute it and/or modify it under the terms of the                  *
15   GNU Affero General Public License as published by the                    *
16   Free Software Foundation, either version 3 of the License,               *
17   or (at your option) any later version.                                   *
18                                                                            *
19 * The CertiVox MIRACL Crypto SDK is distributed in the hope                *
20   that it will be useful, but WITHOUT ANY WARRANTY; without even the       *
21   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
22   See the GNU Affero General Public License for more details.              *
23                                                                            *
24 * You should have received a copy of the GNU Affero General Public         *
25   License along with CertiVox MIRACL Crypto SDK.                           *
26   If not, see <http://www.gnu.org/licenses/>.                              *
27                                                                            *
28 You can be released from the requirements of the license by purchasing     *
29 a commercial license. Buying such a license is mandatory as soon as you    *
30 develop commercial activities involving the CertiVox MIRACL Crypto SDK     *
31 without disclosing the source code of your own applications, or shipping   *
32 the CertiVox MIRACL Crypto SDK with a closed source product.               *
33                                                                            *
34 ***************************************************************************/
35 /*
36  *
37  * mnt_pair.cpp
38  *
39  * MNT curve, ate pairing embedding degree 6, ideal for security level AES-80
40  *
41  *
42  *  Irreducible binomial MUST be of the form x^6+2. This excludes many of the curves
43  *  found using the mnt utility!
44  *  NOTE: This version uses a "compositum". That is the ZZn6 class is a cubic tower over ZZn2, but can
45  *  also be considered as a quadratic tower over ZZn3. The routine shuffle converts from one form to the other.
46  *  The former is fastest for ZZn6 arithmetic, the latter form is required for handling the second parameter
47  *  to the pairing, which is on the quadratic twist E(Fp3)
48  *
49  * Provides high level interface to pairing functions
50  *
51  * GT=pairing(G2,G1)
52  *
53  * This is calculated on a Pairing Friendly Curve (PFC), which must first be defined.
54  *
55  * G1 is a point over the base field, and G2 is a point over an extension field of degree 3
56  * GT is a finite field point over the 6-th extension, where 6 is the embedding degree.
57  *
58  */
59 
60 #define MR_PAIRING_MNT
61 #include "pairing_3.h"
62 
63 // AES_SECURITY=80 bit curve
64 // MNT curve parameters, x,A,B
65 // Thanks to Drew Sutherland for providing the MNT curve
66 // irreducible poly is x^6+2
67 static char param[]="-D285DA0CFEF02F06F812";
68 static char curveB[]="77479D33943B5B1F590B54258B72F316B3261D45";
69 
read_only_error(void)70 void read_only_error(void)
71 {
72 	cout << "Attempt to write to read-only object" << endl;
73 	exit(0);
74 }
75 
set_frobenius_constant(ZZn2 & X)76 void set_frobenius_constant(ZZn2 &X)
77 {
78     Big p=get_modulus();
79     switch (get_mip()->pmod8)
80     {
81     case 5:
82          X.set((Big)0,(Big)1); // = (sqrt(-2)^(p-1)/2
83          break;
84     case 3:                     // = (1+sqrt(-1))^(p-1)/2
85          X.set((Big)1,(Big)1);
86          break;
87    case 7:
88          X.set((Big)2,(Big)1); // = (2+sqrt(-1))^(p-1)/2
89     default: break;
90     }
91     X=pow(X,(p-1)/3);
92 }
93 
94 // Using SHA as basic hash algorithm
95 //
96 // Hash function
97 //
98 
99 #define HASH_LEN 20
100 
H1(char * string)101 Big H1(char *string)
102 { // Hash a zero-terminated string to a number < modulus
103     Big h,p;
104     char s[HASH_LEN];
105     int i,j;
106     sha sh;
107 
108     shs_init(&sh);
109 
110     for (i=0;;i++)
111     {
112         if (string[i]==0) break;
113         shs_process(&sh,string[i]);
114     }
115     shs_hash(&sh,s);
116     p=get_modulus();
117     h=1; j=0; i=1;
118     forever
119     {
120         h*=256;
121         if (j==HASH_LEN)  {h+=i++; j=0;}
122         else         h+=s[j++];
123         if (h>=p) break;
124     }
125     h%=p;
126     return h;
127 }
128 
start_hash(void)129 void PFC::start_hash(void)
130 {
131 	shs_init(&SH);
132 }
133 
finish_hash_to_group(void)134 Big PFC::finish_hash_to_group(void)
135 {
136 	Big hash;
137 	char s[HASH_LEN];
138     shs_hash(&SH,s);
139     hash=from_binary(HASH_LEN,s);
140 	return hash%(*ord);
141 }
142 
add_to_hash(const GT & x)143 void PFC::add_to_hash(const GT& x)
144 {
145 	ZZn6 u=x.g;
146 	ZZn2 v;
147 	ZZn l,h;
148 	Big a,xx[2];
149 	int i,j,m;
150 
151 	u.get(v);
152 	v.get(l,h);
153 	xx[0]=l; xx[1]=h;
154 
155     for (i=0;i<2;i++)
156     {
157         a=xx[i];
158         while (a>0)
159         {
160             m=a%256;
161             shs_process(&SH,m);
162             a/=256;
163         }
164     }
165 
166 }
167 
add_to_hash(const G2 & x)168 void PFC::add_to_hash(const G2& x)
169 {
170 	ZZn3 X,Y;
171 	ECn3 v=x.g;
172 	Big a;
173 	ZZn xx[6];
174 
175 	int i,m;
176 
177 	v.get(X,Y);
178 	X.get(xx[0],xx[1],xx[2]);
179 	Y.get(xx[3],xx[4],xx[5]);
180 	for (i=0;i<6;i++)
181     {
182         a=(Big)xx[i];
183         while (a>0)
184         {
185             m=a%256;
186             shs_process(&SH,m);
187             a/=256;
188         }
189     }
190 }
191 
add_to_hash(const G1 & x)192 void PFC::add_to_hash(const G1& x)
193 {
194 	Big a,X,Y;
195 	int i,m;
196 	x.g.get(X,Y);
197 	a=X;
198     while (a>0)
199     {
200         m=a%256;
201         shs_process(&SH,m);
202         a/=256;
203     }
204 	a=Y;
205     while (a>0)
206     {
207         m=a%256;
208         shs_process(&SH,m);
209         a/=256;
210     }
211 }
212 
add_to_hash(const Big & x)213 void PFC::add_to_hash(const Big& x)
214 {
215 	int m;
216 	Big a=x;
217     while (a>0)
218     {
219         m=a%256;
220         shs_process(&SH,m);
221         a/=256;
222     }
223 }
224 
H2(ZZn6 y)225 Big H2(ZZn6 y)
226 { // Hash and compress an Fp6 to a big number
227     sha sh;
228     ZZn u,v,w;
229 	ZZn2 x;
230     Big a,h,xx[2];
231     char s[HASH_LEN];
232     int i,j,m;
233 
234     shs_init(&sh);
235 	y.get(x);
236     x.get(u,v);
237     xx[0]=u; xx[1]=v;
238 
239     for (i=0;i<2;i++)
240     {
241         a=xx[i];
242         while (a>0)
243         {
244             m=a%256;
245             shs_process(&sh,m);
246             a/=256;
247         }
248     }
249     shs_hash(&sh,s);
250     h=from_binary(HASH_LEN,s);
251     return h;
252 }
253 
254 #ifndef MR_AFFINE_ONLY
255 
force(ZZn & x,ZZn & y,ZZn & z,ECn & A)256 void force(ZZn& x,ZZn& y,ZZn& z,ECn& A)
257 {  // A=(x,y,z)
258     copy(getbig(x),A.get_point()->X);
259     copy(getbig(y),A.get_point()->Y);
260     copy(getbig(z),A.get_point()->Z);
261     A.get_point()->marker=MR_EPOINT_GENERAL;
262 }
263 
extract(ECn & A,ZZn & x,ZZn & y,ZZn & z)264 void extract(ECn &A, ZZn& x,ZZn& y,ZZn& z)
265 { // (x,y,z) <- A
266     big t;
267     x=(A.get_point())->X;
268     y=(A.get_point())->Y;
269     t=(A.get_point())->Z;
270     if (A.get_status()!=MR_EPOINT_GENERAL) z=1;
271     else                                   z=t;
272 }
273 
274 #endif
275 
force(ZZn & x,ZZn & y,ECn & A)276 void force(ZZn& x,ZZn& y,ECn& A)
277 { // A=(x,y)
278     copy(getbig(x),A.get_point()->X);
279     copy(getbig(y),A.get_point()->Y);
280     A.get_point()->marker=MR_EPOINT_NORMALIZED;
281 }
282 
extract(ECn & A,ZZn & x,ZZn & y)283 void extract(ECn& A,ZZn& x,ZZn& y)
284 { // (x,y) <- A
285     x=(A.get_point())->X;
286     y=(A.get_point())->Y;
287 }
288 
289 
shuffle(const ZZn3 & first,const ZZn3 & second)290 ZZn6 shuffle(const ZZn3 &first, const ZZn3 &second)
291 { // shuffle from a pair ZZn3's to three ZZn2's, as required by ZZn6
292 	ZZn6 w;
293 	ZZn x0,x1,x2,x3,x4,x5;
294 	ZZn2 t0,t1,t2;
295 	first.get(x0,x2,x4);
296 	second.get(x1,x3,x5);
297 	t0.set(x0,x3);
298 	t1.set(x1,x4);
299 	t2.set(x2,x5);
300 	w.set(t0,t1,t2);
301 	return w;
302 }
303 
unshuffle(ZZn6 & S,ZZn3 & first,ZZn3 & second)304 void unshuffle(ZZn6 &S,ZZn3 &first,ZZn3 &second)
305 { // unshuffle a ZZn6 into two ZZn3's
306 	ZZn x0,x1,x2,x3,x4,x5;
307 	ZZn2 t0,t1,t2;
308 	S.get(t0,t1,t2);
309 	t0.get(x0,x3);
310 	t1.get(x1,x4);
311 	t2.get(x2,x5);
312 	first.set(x0,x2,x4);
313 	second.set(x1,x3,x5);
314 }
315 
316 // Calculate q*P. P(X,Y) -> P(X^p,Y^p))
317 
q_power_frobenius(ECn3 & S,ZZn2 & X)318 void q_power_frobenius(ECn3 &S,ZZn2& X)
319 {
320 	ZZn6 X1,X2,Y1,Y2;
321 	ZZn3 Sx,Sy,T;
322 
323 	int qnr=get_mip()->cnr;
324 
325 	S.get(Sx,Sy);
326 
327 	// untwist
328     Sx=Sx/qnr;
329     Sy=tx(Sy);
330     Sy=Sy/(qnr*qnr);
331 
332 	X1=shuffle(Sx,(ZZn3)0); Y1=shuffle((ZZn3)0,Sy);
333 	X1.powq(X); Y1.powq(X);
334 	unshuffle(X1,Sx,T); unshuffle(Y1,T,Sy);
335 
336 	// twist
337 	Sx=qnr*Sx;
338 	Sy=txd(Sy*qnr*qnr);
339 	S.set(Sx,Sy);
340 }
341 
342 //
343 // Line from A to destination C. Let A=(x,y)
344 // Line Y-slope.X-c=0, through A, so intercept c=y-slope.x
345 // Line Y-slope.X-y+slope.x = (Y-y)-slope.(X-x) = 0
346 // Now evaluate at Q -> return (Qy-y)-slope.(Qx-x)
347 //
348 
line(ECn3 & A,ECn3 & C,ECn3 & B,int type,ZZn3 & slope,ZZn3 & ex1,ZZn3 & ex2,ZZn & Px,ZZn & Py)349 ZZn6 line(ECn3& A,ECn3& C,ECn3& B,int type,ZZn3& slope,ZZn3& ex1,ZZn3& ex2,ZZn& Px,ZZn& Py)
350 {
351     ZZn6 w;
352 	ZZn3 d;
353     ZZn3 x,y;
354 #ifdef MR_ECN3_PROJECTIVE
355 	ZZn3 z,z3,t;
356 	C.getZ(z3);
357 	d.set1(Py);
358 
359 	if (type==MR_ADD)
360 	{ // exploit that B is in affine
361 		ZZn3 x2,y2;
362 		B.get(x2,y2);
363 		y2*=z3; d*=z3;
364 		w=shuffle(y2-slope*(Px+x2),d);
365 	}
366 	if (type==MR_DOUBLE)
367 	{ // use extra information from point doubling
368 		A.get(x,y,z);
369 		w=shuffle(ex1-slope*(Px*ex2+x),d*z3*ex2);
370 	}
371 #else
372 	A.get(x,y);
373     d.set1(Py);
374 	w=shuffle(y-slope*(Px+x),d);
375 #endif
376     return w;
377 }
378 
379 //
380 // Add A=A+B  (or A=A+A)
381 // Return line function value
382 //
383 
g(ECn3 & A,ECn3 & B,ZZn & Px,ZZn & Py)384 ZZn6 g(ECn3& A,ECn3& B,ZZn& Px,ZZn& Py)
385 {
386     BOOL type;
387     ZZn3 lam,ex1,ex2;
388     ECn3 Q=A;
389 
390 // Evaluate line from A to A+B
391     type=A.add(B,lam,&ex1,&ex2);
392 
393     return line(Q,A,B,type,lam,ex1,ex2,Px,Py);
394 }
395 
396 // if multiples of G2 can be precalculated, its a lot faster!
397 
gp(ZZn3 * ptable,int & j,ZZn & Px,ZZn & Py)398 ZZn6 gp(ZZn3* ptable,int &j,ZZn& Px,ZZn& Py)
399 {
400 	ZZn6 w;
401 	ZZn3 d;
402 	d.set1(Py);
403 	w=shuffle(ptable[j]*Px+ptable[j+1],d);
404 	j+=2;
405 	return w;
406 }
407 
408 //
409 // Spill precomputation on pairing to byte array
410 //
411 
spill(G2 & w,char * & bytes)412 int PFC::spill(G2& w,char *& bytes)
413 {
414 	int i,j,len,m;
415 	int bytes_per_big=(MIRACL/8)*(get_mip()->nib-1);
416 
417 	ZZn a,b,c;
418 	Big X=*x;
419 	if (w.ptable==NULL) return 0;
420 
421 	m=2*(bits(X)-2+ham(X));
422 	len=m*3*bytes_per_big;
423 
424 	bytes=new char[len];
425 	for (i=j=0;i<m;i++)
426 	{
427 		w.ptable[i].get(a,b,c);
428 		to_binary((Big)a,bytes_per_big,&bytes[j],TRUE);
429 		j+=bytes_per_big;
430 		to_binary((Big)b,bytes_per_big,&bytes[j],TRUE);
431 		j+=bytes_per_big;
432 		to_binary((Big)c,bytes_per_big,&bytes[j],TRUE);
433 		j+=bytes_per_big;
434 	}
435 
436 	delete [] w.ptable;
437 	w.ptable=NULL;
438 	return len;
439 }
440 
441 //
442 // Restore precomputation on pairing to byte array
443 //
444 
restore(char * bytes,G2 & w)445 void PFC::restore(char * bytes,G2& w)
446 {
447 	int i,j,len,m;
448 	int bytes_per_big=(MIRACL/8)*(get_mip()->nib-1);
449 
450 	ZZn a,b,c;
451 	Big X=*x;
452 	if (w.ptable!=NULL) return;
453 
454 	m=2*(bits(X)-2+ham(X));
455 	len=m*3*bytes_per_big;
456 
457 	w.ptable=new ZZn3[m];
458 	for (i=j=0;i<m;i++)
459 	{
460 		a=from_binary(bytes_per_big,&bytes[j]);
461 		j+=bytes_per_big;
462 		b=from_binary(bytes_per_big,&bytes[j]);
463 		j+=bytes_per_big;
464 		c=from_binary(bytes_per_big,&bytes[j]);
465 		j+=bytes_per_big;
466 		w.ptable[i].set(a,b,c);
467 	}
468 	for (i=0;i<len;i++) bytes[i]=0;
469 
470 	delete [] bytes;
471 }
472 
473 // precompute G2 table for pairing
474 
precomp_for_pairing(G2 & w)475 int PFC::precomp_for_pairing(G2& w)
476 {
477 	int i,j,nb,type,len;
478 	ECn3 A,Q,B;
479 	ZZn3 lam,x1,y1;
480 	Big X=*x;
481 
482 	A=w.g;
483 	A.norm();
484 	B=A;
485 	nb=bits(X);
486 	j=0;
487 	len=2*(nb-2+ham(X));
488 	w.ptable=new ZZn3[len];
489 	get_mip()->coord=MR_AFFINE;  // switch to affine
490     for (i=nb-2;i>=0;i--)
491     {
492 		Q=A;
493 // Evaluate line from A to A+B
494 		A.add(A,lam,NULL,NULL);
495 		Q.get(x1,y1);
496 		w.ptable[j++]=-lam; w.ptable[j++]=y1-lam*x1;
497 
498 		if (bit(X,i)==1)
499 		{
500 			Q=A;
501 			type=A.add(B,lam,NULL,NULL);
502 			Q.get(x1,y1);
503 			w.ptable[j++]=-lam; w.ptable[j++]=y1-lam*x1;
504 		}
505     }
506 	get_mip()->coord=MR_PROJECTIVE;
507 	return len;
508 }
509 
multi_miller(int n,G2 ** QQ,G1 ** PP)510 GT PFC::multi_miller(int n,G2** QQ,G1** PP)
511 {
512 	GT z;
513     ZZn *Px,*Py;
514 	int i,j,*k,nb;
515     ECn3 *Q,*A;
516 	ECn P;
517     ZZn6 res;
518 	Big X=*x;
519 
520 	Px=new ZZn[n];
521 	Py=new ZZn[n];
522 	Q=new ECn3[n];
523 	A=new ECn3[n];
524 	k=new int[n];
525 
526     nb=bits(X);
527 	res=1;
528 
529 	for (j=0;j<n;j++)
530 	{
531 		k[j]=0;
532 		P=PP[j]->g; normalise(P); Q[j]=QQ[j]->g;
533 		extract(P,Px[j],Py[j]);
534 		Px[j]+=Px[j];
535 		Py[j]+=Py[j];
536 	}
537 
538 	for (j=0;j<n;j++)
539 	{
540 #ifdef MR_ECN3_PROJECTIVE
541 		Q[j].norm();
542 #endif
543 		A[j]=Q[j];
544 	}
545 
546 	for (i=nb-2;i>=0;i--)
547 	{
548 		res*=res;
549 		for (j=0;j<n;j++)
550 		{
551 			if (QQ[j]->ptable==NULL)
552 				res*=g(A[j],A[j],Px[j],Py[j]);
553 			else
554 				res*=gp(QQ[j]->ptable,k[j],Px[j],Py[j]);
555 		}
556 		if (bit(X,i)==1)
557 			for (j=0;j<n;j++)
558 			{
559 				if (QQ[j]->ptable==NULL)
560 					res*=g(A[j],Q[j],Px[j],Py[j]);
561 				else
562 					res*=gp(QQ[j]->ptable,k[j],Px[j],Py[j]);
563 			}
564 		if (res.iszero()) return 0;
565 	}
566 
567 	delete [] k;
568 	delete [] A;
569 	delete [] Q;
570 	delete [] Py;
571 	delete [] Px;
572 
573 	z.g=res;
574 	return z;
575 }
576 
577 //
578 // R-ate Pairing G2 x G1 -> GT
579 //
580 // P is a point of order q in G1. Q(x,y) is a point of order q in G2.
581 // Note that Q is a point on the sextic twist of the curve over Fp^2, P(x,y) is a point on the
582 // curve over the base field Fp
583 //
584 
miller_loop(const G2 & QQ,const G1 & PP)585 GT PFC::miller_loop(const G2& QQ,const G1& PP)
586 {
587 	GT z;
588     int i,j,n,nb,nbw,nzs;
589     ECn3 A,Q;
590 	ECn P;
591 	ZZn Px,Py;
592 	BOOL precomp;
593     ZZn6 res;
594 	Big X=*x;
595 
596 	P=PP.g; Q=QQ.g;
597 #ifdef MR_ECN3_PROJECTIVE
598 	Q.norm();
599 #endif
600 	precomp=FALSE;
601 	if (QQ.ptable!=NULL) precomp=TRUE;
602 
603 	normalise(P);
604 	extract(P,Px,Py);
605 
606     Px+=Px;  // because x^6+2 is irreducible.. simplifies line function calculation
607     Py+=Py;
608 
609     res=1;
610     A=Q;    // reset A
611     nb=bits(X);
612 	res.mark_as_miller();
613 	j=0;
614 
615     for (i=nb-2;i>=0;i--)
616     {
617 		res*=res;
618 		if (precomp) res*=gp(QQ.ptable,j,Px,Py);
619 		else         res*=g(A,A,Px,Py);
620 
621 		if (bit(X,i)==1)
622 		{
623 			if (precomp) res*=gp(QQ.ptable,j,Px,Py);
624 			else         res*=g(A,Q,Px,Py);
625 		}
626     }
627 
628 	z.g=res;
629 	return z;
630 }
631 
final_exp(const GT & z)632 GT PFC::final_exp(const GT& z)
633 {
634 	GT y;
635 	ZZn6 w,res;
636 	Big X=*x;
637 
638 	res=z.g;
639 
640     w=res;
641     w.powq(*frob);
642     res*=w;                        // ^(p+1)
643 
644     w=res;
645     w.powq(*frob); w.powq(*frob); w.powq(*frob);
646     res=w/res;                     // ^(p^3-1)
647 
648 // exploit the clever "trick" for a half-length exponentiation!
649 
650     res.mark_as_unitary();
651 
652     w=res;
653     res.powq(*frob);  // res*=res;  // res=pow(res,CF);
654 
655     if (X<0) res/=powu(w,-X);
656     else res*=powu(w,X);
657 
658     y.g=res;
659 
660 	return y;
661 }
662 
PFC(int s,csprng * rng)663 PFC::PFC(int s, csprng *rng)
664 {
665 	int mod_bits,words;
666 	if (s!=80)
667 	{
668 		cout << "No suitable curve available" << endl;
669 		exit(0);
670 	}
671 	mod_bits=2*s;
672 
673 	if (mod_bits%MIRACL==0)
674 		words=(mod_bits/MIRACL);
675 	else
676 		words=(mod_bits/MIRACL)+1;
677 
678 #ifdef MR_SIMPLE_BASE
679 	miracl *mip=mirsys((MIRACL/4)*words,16);
680 #else
681 	miracl *mip=mirsys(words,0);
682 	mip->IOBASE=16;
683 #endif
684 
685 	B=new Big;
686 	x=new Big;
687 	mod=new Big;
688 	ord=new Big;
689 	cof=new Big;
690 	npoints=new Big;
691 	trace=new Big;
692 	frob=new ZZn2;
693 
694 	*B=curveB;
695 	S=s;
696 	*x=param;
697 	Big X=*x;
698 
699 	*mod=X*X+1;
700 	*npoints=X*X-X+1;
701 	*trace=X+1;
702 	*cof=X*X+X+1;
703 	*ord=*npoints;
704 	ecurve(-3,*B,*mod,MR_PROJECTIVE);
705 	set_frobenius_constant(*frob);
706 	Big sru=pow((ZZn)-2,(*mod-1)/6);   // x^6+2 is irreducible
707     set_zzn3(-2,sru);
708 	mip->TWIST=MR_QUADRATIC;   // twisted curve E'(ZZn3)
709 
710 	RNG = rng;
711 }
712 
~PFC()713 PFC::~PFC()
714 {
715 	delete B;
716 	delete x;
717 	delete mod;
718 	delete ord;
719 	delete cof;
720 	delete npoints;
721 	delete trace;
722 	delete frob;
723 	mirexit();
724 }
725 
mult(const G1 & w,const Big & k)726 G1 PFC::mult(const G1& w,const Big& k)
727 {
728 	G1 z;
729 	if (w.mtable!=NULL)
730 	{ // we have precomputed values
731 		Big e=k;
732 		if (k<0) e=-e;
733 
734 		int i,j,t=w.mtbits; //MR_ROUNDUP(2*S,WINDOW_SIZE);
735 		j=recode(e,t,WINDOW_SIZE,t-1);
736 		z.g=w.mtable[j];
737 		for (i=t-2;i>=0;i--)
738 		{
739 			j=recode(e,t,WINDOW_SIZE,i);
740 			z.g+=z.g;
741 			if (j>0) z.g+=w.mtable[j];
742 		}
743 		if (k<0) z.g=-z.g;
744 	}
745 	else
746 	{
747 		z.g=w.g;
748 		z.g*=k;
749 	}
750 	return z;
751 }
752 
753 // GLV + Galbraith-Scott
754 
mult(const G2 & w,const Big & k)755 G2 PFC::mult(const G2& w,const Big& k)
756 {
757 	G2 z;
758 	Big X=*x;
759 	if (w.mtable!=NULL)
760 	{ // we have precomputed values
761 		Big e=k;
762 		if (k<0) e=-e;
763 
764 		int i,j,t=w.mtbits; //MR_ROUNDUP(2*S,WINDOW_SIZE);
765 		j=recode(e,t,WINDOW_SIZE,t-1);
766 		z.g=w.mtable[j];
767 		for (i=t-2;i>=0;i--)
768 		{
769 			j=recode(e,t,WINDOW_SIZE,i);
770 			z.g+=z.g;
771 			if (j>0) z.g+=w.mtable[j];
772 		}
773 		if (k<0) z.g=-z.g;
774 	}
775 	else
776 	{
777 		ECn3 v=w.g;
778 		q_power_frobenius(v,*frob);
779 		z.g=mul(v,k/X,w.g,k%X);
780 	}
781 	return z;
782 }
783 
784 // GLV method + Galbraith-Scott idea
785 
power(const GT & w,const Big & k)786 GT PFC::power(const GT& w,const Big& k)
787 {
788 	GT z;
789 	Big X=*x;
790 	if (w.etable!=NULL)
791 	{ // precomputation is available
792 		Big e=k;
793 		if (k<0) e=-e;
794 
795 		int i,j,t=w.etbits; //MR_ROUNDUP(2*S,WINDOW_SIZE);
796 		j=recode(e,t,WINDOW_SIZE,t-1);
797 		z.g=w.etable[j];
798 		for (i=t-2;i>=0;i--)
799 		{
800 			j=recode(e,t,WINDOW_SIZE,i);
801 			z.g*=z.g;
802 			if (j>0) z.g*=w.etable[j];
803 		}
804 		if (k<0) z.g=inverse(z.g);
805 	}
806 	else
807 	{
808 		ZZn6 y=w.g;
809 		y.powq(*frob);
810 		z.g=powu(y,k/X,w.g,k%X);
811 	}
812 	return z;
813 }
814 
815 // Use Scott et al. idea - http://eprint.iacr.org/2008/530.pdf
816 // Map to point of correct order
817 
map(ECn3 & S,Big x,ZZn2 & X)818 void map(ECn3 &S,Big x, ZZn2& X)
819 { // S=Phi(2xP)+phi^2(2xP)
820 	ZZn6 X1,X2,Y1,Y2;
821 	ZZn3 Sx,Sy,T;
822 	ECn3 S2;
823 	int qnr=get_mip()->cnr;
824 
825 	S*=x; S+=S; // hard work done here
826 
827 	S.get(Sx,Sy);
828 
829 	// untwist
830     Sx=Sx/qnr;
831     Sy=tx(Sy);
832     Sy=Sy/(qnr*qnr);
833 
834 	X1=shuffle(Sx,(ZZn3)0); Y1=shuffle((ZZn3)0,Sy);
835 	X1.powq(X); Y1.powq(X);
836 	X2=X1; Y2=Y1;
837 	X2.powq(X); Y2.powq(X);
838 	unshuffle(X1,Sx,T); unshuffle(Y1,T,Sy);
839 
840 	// twist
841 	Sx=qnr*Sx;
842 	Sy=txd(Sy*qnr*qnr);
843 	S.set(Sx,Sy);
844 	unshuffle(X2,Sx,T); unshuffle(Y2,T,Sy);
845 
846 	//twist (again, like we did last summer...)
847 	Sx=qnr*Sx;
848 	Sy=txd(Sy*qnr*qnr);
849 	S2.set(Sx,Sy);
850 	S+=S2;
851 }
852 
853 // random group element
854 
random(Big & w)855 void PFC::random(Big& w)
856 {
857 	if (RNG==NULL) w=rand(*ord);
858 	else w=strong_rand(RNG,*ord);
859 }
860 
861 // random AES key
862 
rankey(Big & k)863 void PFC::rankey(Big& k)
864 {
865 		if (RNG==NULL) k=rand(S,2);
866 		else k=strong_rand(RNG,S,2);
867 }
868 
hash_and_map(G2 & w,char * ID)869 void PFC::hash_and_map(G2& w,char *ID)
870 {
871     int i;
872     ZZn3 XX;
873 	Big X=*x;
874 
875     Big x0=H1(ID);
876     forever
877     {
878         x0+=1;
879         XX.set2((ZZn)x0);
880         if (!w.g.set(XX)) continue;
881 
882         break;
883     }
884 	map(w.g,X,*frob);
885 }
886 
random(G2 & w)887 void PFC::random(G2& w)
888 {
889     int i;
890     ZZn3 XX;
891 	Big X=*x;
892 	Big x0;
893 
894     if (RNG==NULL) x0=rand(*mod);
895 	else x0=strong_rand(RNG,*mod);
896     forever
897     {
898         x0+=1;
899         XX.set2((ZZn)x0);
900         if (!w.g.set(XX)) continue;
901 
902         break;
903     }
904 	map(w.g,X,*frob);
905 }
906 
hash_and_map(G1 & w,char * ID)907 void PFC::hash_and_map(G1& w,char *ID)
908 {
909     Big x0=H1(ID);
910     while (!w.g.set(x0,x0)) x0+=1;
911 }
912 
random(G1 & w)913 void PFC::random(G1& w)
914 {
915 	Big x0;
916 	if (RNG==NULL) x0=rand(*mod);
917 	else x0=strong_rand(RNG,*mod);
918 
919 	while (!w.g.set(x0,x0)) x0+=1;
920 }
921 
hash_to_aes_key(const GT & w)922 Big PFC::hash_to_aes_key(const GT& w)
923 {
924 	Big m=pow((Big)2,S);
925 	return H2(w.g)%m;
926 }
927 
hash_to_group(char * ID)928 Big PFC::hash_to_group(char *ID)
929 {
930 	Big m=H1(ID);
931 	return m%(*ord);
932 }
933 
operator *(const GT & x,const GT & y)934 GT operator*(const GT& x,const GT& y)
935 {
936 	GT z=x;
937 	z.g*=y.g;
938 	return z;
939 }
940 
operator /(const GT & x,const GT & y)941 GT operator/(const GT& x,const GT& y)
942 {
943 	GT z=x;
944 	z.g/=y.g;
945 	return z;
946 }
947 
948 //
949 // spill precomputation on GT to byte array
950 //
951 
spill(char * & bytes)952 int GT::spill(char *& bytes)
953 {
954 	int i,j,n=(1<<WINDOW_SIZE);
955 	int bytes_per_big=(MIRACL/8)*(get_mip()->nib-1);
956 	int len=n*6*bytes_per_big;
957 	ZZn2 a,b,c;
958 	Big x,y;
959 
960 	if (etable==NULL) return 0;
961 
962 	bytes=new char[len];
963 	for (i=j=0;i<n;i++)
964 	{
965 		etable[i].get(a,b,c);
966 		a.get(x,y);
967 		to_binary(x,bytes_per_big,&bytes[j],TRUE);
968 		j+=bytes_per_big;
969 		to_binary(y,bytes_per_big,&bytes[j],TRUE);
970 		j+=bytes_per_big;
971 		b.get(x,y);
972 		to_binary(x,bytes_per_big,&bytes[j],TRUE);
973 		j+=bytes_per_big;
974 		to_binary(y,bytes_per_big,&bytes[j],TRUE);
975 		j+=bytes_per_big;
976 		c.get(x,y);
977 		to_binary(x,bytes_per_big,&bytes[j],TRUE);
978 		j+=bytes_per_big;
979 		to_binary(y,bytes_per_big,&bytes[j],TRUE);
980 		j+=bytes_per_big;
981 	}
982 	delete [] etable;
983 	etable=NULL;
984 	return len;
985 }
986 
987 //
988 // restore precomputation for GT from byte array
989 //
990 
restore(char * bytes)991 void GT::restore(char *bytes)
992 {
993 	int i,j,n=(1<<WINDOW_SIZE);
994 	int bytes_per_big=(MIRACL/8)*(get_mip()->nib-1);
995 	int len=n*6*bytes_per_big;
996 	ZZn2 a,b,c;
997 	Big x,y;
998 	if (etable!=NULL) return;
999 
1000 	etable=new ZZn6[1<<WINDOW_SIZE];
1001 	for (i=j=0;i<n;i++)
1002 	{
1003 		x=from_binary(bytes_per_big,&bytes[j]);
1004 		j+=bytes_per_big;
1005 		y=from_binary(bytes_per_big,&bytes[j]);
1006 		j+=bytes_per_big;
1007 		a.set(x,y);
1008 		x=from_binary(bytes_per_big,&bytes[j]);
1009 		j+=bytes_per_big;
1010 		y=from_binary(bytes_per_big,&bytes[j]);
1011 		j+=bytes_per_big;
1012 		b.set(x,y);
1013 		x=from_binary(bytes_per_big,&bytes[j]);
1014 		j+=bytes_per_big;
1015 		y=from_binary(bytes_per_big,&bytes[j]);
1016 		j+=bytes_per_big;
1017 		c.set(x,y);
1018 		etable[i].set(a,b,c);
1019 	}
1020 	delete [] bytes;
1021 }
1022 
1023 
operator +(const G1 & x,const G1 & y)1024 G1 operator+(const G1& x,const G1& y)
1025 {
1026 	G1 z=x;
1027 	z.g+=y.g;
1028 	return z;
1029 }
1030 
operator -(const G1 & x)1031 G1 operator-(const G1& x)
1032 {
1033 	G1 z=x;
1034 	z.g=-z.g;
1035 	return z;
1036 }
1037 
1038 //
1039 // spill precomputation on G1 to byte array
1040 //
1041 
spill(char * & bytes)1042 int G1::spill(char *& bytes)
1043 {
1044 	int i,j,n=(1<<WINDOW_SIZE);
1045 	int bytes_per_big=(MIRACL/8)*(get_mip()->nib-1);
1046 	int len=n*2*bytes_per_big;
1047 	Big x,y;
1048 
1049 	if (mtable==NULL) return 0;
1050 
1051 	bytes=new char[len];
1052 	for (i=j=0;i<n;i++)
1053 	{
1054 		mtable[i].get(x,y);
1055 		to_binary(x,bytes_per_big,&bytes[j],TRUE);
1056 		j+=bytes_per_big;
1057 		to_binary(y,bytes_per_big,&bytes[j],TRUE);
1058 		j+=bytes_per_big;
1059 	}
1060 	delete [] mtable;
1061 	mtable=NULL;
1062 	return len;
1063 }
1064 
1065 //
1066 // restore precomputation for G1 from byte array
1067 //
1068 
restore(char * bytes)1069 void G1::restore(char *bytes)
1070 {
1071 	int i,j,n=(1<<WINDOW_SIZE);
1072 	int bytes_per_big=(MIRACL/8)*(get_mip()->nib-1);
1073 	int len=n*2*bytes_per_big;
1074 	Big x,y;
1075 	if (mtable!=NULL) return;
1076 
1077 	mtable=new ECn[1<<WINDOW_SIZE];
1078 	for (i=j=0;i<n;i++)
1079 	{
1080 		x=from_binary(bytes_per_big,&bytes[j]);
1081 		j+=bytes_per_big;
1082 		y=from_binary(bytes_per_big,&bytes[j]);
1083 		j+=bytes_per_big;
1084 		mtable[i].set(x,y);
1085 	}
1086 	delete [] bytes;
1087 }
1088 
operator +(const G2 & x,const G2 & y)1089 G2 operator+(const G2& x,const G2& y)
1090 {
1091 	G2 z=x;
1092 	z.g+=y.g;
1093 	return z;
1094 }
1095 
operator -(const G2 & x)1096 G2 operator-(const G2& x)
1097 {
1098 	G2 z=x;
1099 	z.g=-z.g;
1100 	return z;
1101 }
1102 
1103 //
1104 // spill precomputation on G2 to byte array
1105 //
1106 
spill(char * & bytes)1107 int G2::spill(char *& bytes)
1108 {
1109 	int i,j,n=(1<<WINDOW_SIZE);
1110 	int bytes_per_big=(MIRACL/8)*(get_mip()->nib-1);
1111 	int len=n*6*bytes_per_big;
1112 	ZZn3 x,y;
1113 	ZZn a,b,c;
1114 
1115 	if (mtable==NULL) return 0;
1116 
1117 	bytes=new char[len];
1118 	for (i=j=0;i<n;i++)
1119 	{
1120 		mtable[i].get(x,y);
1121 		x.get(a,b,c);
1122 		to_binary((Big)a,bytes_per_big,&bytes[j],TRUE);
1123 		j+=bytes_per_big;
1124 		to_binary((Big)b,bytes_per_big,&bytes[j],TRUE);
1125 		j+=bytes_per_big;
1126 		to_binary((Big)c,bytes_per_big,&bytes[j],TRUE);
1127 		j+=bytes_per_big;
1128 		y.get(a,b,c);
1129 		to_binary((Big)a,bytes_per_big,&bytes[j],TRUE);
1130 		j+=bytes_per_big;
1131 		to_binary((Big)b,bytes_per_big,&bytes[j],TRUE);
1132 		j+=bytes_per_big;
1133 		to_binary((Big)c,bytes_per_big,&bytes[j],TRUE);
1134 		j+=bytes_per_big;
1135 	}
1136 	delete [] mtable;
1137 	mtable=NULL;
1138 	return len;
1139 }
1140 
1141 //
1142 // restore precomputation for G2 from byte array
1143 //
1144 
restore(char * bytes)1145 void G2::restore(char *bytes)
1146 {
1147 	int i,j,n=(1<<WINDOW_SIZE);
1148 	int bytes_per_big=(MIRACL/8)*(get_mip()->nib-1);
1149 	int len=n*6*bytes_per_big;
1150 	ZZn3 x,y;
1151 	ZZn a,b,c;
1152 	if (mtable!=NULL) return;
1153 
1154 	mtable=new ECn3[1<<WINDOW_SIZE];
1155 	for (i=j=0;i<n;i++)
1156 	{
1157 		a=from_binary(bytes_per_big,&bytes[j]);
1158 		j+=bytes_per_big;
1159 		b=from_binary(bytes_per_big,&bytes[j]);
1160 		j+=bytes_per_big;
1161 		c=from_binary(bytes_per_big,&bytes[j]);
1162 		j+=bytes_per_big;
1163 		x.set(a,b,c);
1164 		a=from_binary(bytes_per_big,&bytes[j]);
1165 		j+=bytes_per_big;
1166 		b=from_binary(bytes_per_big,&bytes[j]);
1167 		j+=bytes_per_big;
1168 		c=from_binary(bytes_per_big,&bytes[j]);
1169 		j+=bytes_per_big;
1170 		y.set(a,b,c);
1171 		mtable[i].set(x,y);
1172 	}
1173 	delete [] bytes;
1174 }
1175 
member(const GT & z)1176 BOOL PFC::member(const GT& z)
1177 {
1178 	ZZn6 r=z.g;
1179 	ZZn6 w=z.g;
1180 	Big X=*x;
1181 	if (!r.is_unitary()) return FALSE;
1182 	if (r*conj(r)!=(ZZn6)1) return FALSE; // not unitary
1183 	w.powq(*frob);
1184 	if (X<0) r=powu(inverse(r),-X);
1185 	else     r=powu(r,X);
1186 	if (r==w) return TRUE;
1187 	return FALSE;
1188 }
1189 
pairing(const G2 & x,const G1 & y)1190 GT PFC::pairing(const G2& x,const G1& y)
1191 {
1192 	GT z;
1193 	z=miller_loop(x,y);
1194 	z=final_exp(z);
1195 	return z;
1196 }
1197 
multi_pairing(int n,G2 ** y,G1 ** x)1198 GT PFC::multi_pairing(int n,G2 **y,G1 **x)
1199 {
1200 	GT z;
1201 	z=multi_miller(n,y,x);
1202 	z=final_exp(z);
1203 	return z;
1204 
1205 }
1206 
precomp_for_mult(G1 & w,BOOL small)1207 int PFC::precomp_for_mult(G1& w,BOOL small)
1208 {
1209 	ECn v=w.g;
1210 	int i,j,k,bp,is,t;
1211 	if (small) t=MR_ROUNDUP(2*S,WINDOW_SIZE);
1212 	else       t=MR_ROUNDUP(bits(*ord),WINDOW_SIZE);
1213 	w.mtable=new ECn[1<<WINDOW_SIZE];
1214 	w.mtable[1]=v;
1215 	w.mtbits=t;
1216 	for (j=0;j<t;j++)
1217         v+=v;
1218     k=1;
1219     for (i=2;i<(1<<WINDOW_SIZE);i++)
1220     {
1221         if (i==(1<<k))
1222         {
1223             k++;
1224 			normalise(v);
1225 			w.mtable[i]=v;
1226             for (j=0;j<t;j++)
1227 				v+=v;
1228             continue;
1229         }
1230         bp=1;
1231         for (j=0;j<k;j++)
1232         {
1233             if (i&bp)
1234 			{
1235 				is=1<<j;
1236 				w.mtable[i]+=w.mtable[is];
1237 			}
1238             bp<<=1;
1239         }
1240         normalise(w.mtable[i]);
1241     }
1242 	return (1<<WINDOW_SIZE);
1243 }
1244 
precomp_for_mult(G2 & w,BOOL small)1245 int PFC::precomp_for_mult(G2& w,BOOL small)
1246 {
1247 	ECn3 v;
1248 
1249 	ZZn3 x,y;
1250 	int i,j,k,bp,is,t;
1251 	if (small) t=MR_ROUNDUP(2*S,WINDOW_SIZE);
1252 	else       t=MR_ROUNDUP(bits(*ord),WINDOW_SIZE);
1253 	w.g.norm();
1254 	v=w.g;
1255 	w.mtable=new ECn3[1<<WINDOW_SIZE];
1256 	v.norm();
1257 	w.mtable[1]=v;
1258 	w.mtbits=t;
1259 	for (j=0;j<t;j++)
1260         v+=v;
1261     k=1;
1262 
1263     for (i=2;i<(1<<WINDOW_SIZE);i++)
1264     {
1265         if (i==(1<<k))
1266         {
1267             k++;
1268 			v.norm();
1269 			w.mtable[i]=v;
1270             for (j=0;j<t;j++)
1271 				v+=v;
1272             continue;
1273         }
1274         bp=1;
1275         for (j=0;j<k;j++)
1276         {
1277             if (i&bp)
1278 			{
1279 				is=1<<j;
1280 				w.mtable[i]+=w.mtable[is];
1281 			}
1282             bp<<=1;
1283         }
1284 		w.mtable[i].norm();
1285     }
1286 	return (1<<WINDOW_SIZE);
1287 }
1288 
precomp_for_power(GT & w,BOOL small)1289 int PFC::precomp_for_power(GT& w,BOOL small)
1290 {
1291 	ZZn6 v=w.g;
1292 	int i,j,k,bp,is,t;
1293 	if (small) t=MR_ROUNDUP(2*S,WINDOW_SIZE);
1294 	else       t=MR_ROUNDUP(bits(*ord),WINDOW_SIZE);
1295 	w.etable=new ZZn6[1<<WINDOW_SIZE];
1296 	w.etable[0]=1;
1297 	w.etable[1]=v;
1298 	w.etbits=t;
1299 	for (j=0;j<t;j++)
1300         v*=v;
1301     k=1;
1302 
1303     for (i=2;i<(1<<WINDOW_SIZE);i++)
1304     {
1305         if (i==(1<<k))
1306         {
1307             k++;
1308 			w.etable[i]=v;
1309             for (j=0;j<t;j++)
1310 				v*=v;
1311             continue;
1312         }
1313         bp=1;
1314 		w.etable[i]=1;
1315         for (j=0;j<k;j++)
1316         {
1317             if (i&bp)
1318 			{
1319 				is=1<<j;
1320 				w.etable[i]*=w.etable[is];
1321 			}
1322             bp<<=1;
1323         }
1324     }
1325 	return (1<<WINDOW_SIZE);
1326 }
1327