| /* |
| Licensed to the Apache Software Foundation (ASF) under one |
| or more contributor license agreements. See the NOTICE file |
| distributed with this work for additional information |
| regarding copyright ownership. The ASF licenses this file |
| to you under the Apache License, Version 2.0 (the |
| "License"); you may not use this file except in compliance |
| with the License. You may obtain a copy of the License at |
| |
| http://www.apache.org/licenses/LICENSE-2.0 |
| |
| Unless required by applicable law or agreed to in writing, |
| software distributed under the License is distributed on an |
| "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY |
| KIND, either express or implied. See the License for the |
| specific language governing permissions and limitations |
| under the License. |
| */ |
| |
| /* AMCL basic functions for Large Finite Field support */ |
| |
| #include "amcl.h" |
| |
| /* Arazi and Qi inversion mod 256 */ |
| static int invmod256(int a) |
| { |
| int U,t1,t2,b,c; |
| t1=0; |
| c=(a>>1)&1; |
| t1+=c; |
| t1&=1; |
| t1=2-t1; |
| t1<<=1; |
| U=t1+1; |
| |
| // i=2 |
| b=a&3; |
| t1=U*b; |
| t1>>=2; |
| c=(a>>2)&3; |
| t2=(U*c)&3; |
| t1+=t2; |
| t1*=U; |
| t1&=3; |
| t1=4-t1; |
| t1<<=2; |
| U+=t1; |
| |
| // i=4 |
| b=a&15; |
| t1=U*b; |
| t1>>=4; |
| c=(a>>4)&15; |
| t2=(U*c)&15; |
| t1+=t2; |
| t1*=U; |
| t1&=15; |
| t1=16-t1; |
| t1<<=4; |
| U+=t1; |
| |
| return U; |
| } |
| |
| /* a=1/a mod 2^BIGBITS. This is very fast! */ |
| void BIG_invmod2m(BIG a) |
| { |
| int i; |
| BIG U,t1,b,c; |
| BIG_zero(U); |
| BIG_inc(U,invmod256(BIG_lastbits(a,8))); |
| for (i=8; i<BIGBITS; i<<=1) |
| { |
| BIG_copy(b,a); |
| BIG_mod2m(b,i); // bottom i bits of a |
| |
| BIG_smul(t1,U,b); |
| BIG_shr(t1,i); // top i bits of U*b |
| |
| BIG_copy(c,a); |
| BIG_shr(c,i); |
| BIG_mod2m(c,i); // top i bits of a |
| |
| BIG_smul(b,U,c); |
| BIG_mod2m(b,i); // bottom i bits of U*c |
| |
| BIG_add(t1,t1,b); |
| BIG_smul(b,t1,U); |
| BIG_copy(t1,b); // (t1+b)*U |
| BIG_mod2m(t1,i); // bottom i bits of (t1+b)*U |
| |
| BIG_one(b); |
| BIG_shl(b,i); |
| BIG_sub(t1,b,t1); |
| BIG_norm(t1); |
| |
| BIG_shl(t1,i); |
| |
| BIG_add(U,U,t1); |
| } |
| BIG_copy(a,U); |
| BIG_norm(a); |
| BIG_mod2m(a,BIGBITS); |
| } |
| |
| /* |
| void FF_rcopy(BIG x[],const BIG y[],int n) |
| { |
| int i; |
| for (i=0;i<n;i++) |
| BIG_rcopy(x[i],y[i]); |
| } |
| */ |
| |
| /* x=y */ |
| void FF_copy(BIG x[],BIG y[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_copy(x[i],y[i]); |
| } |
| |
| /* x=y<<n */ |
| static void FF_dsucopy(BIG x[],BIG y[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| { |
| BIG_copy(x[n+i],y[i]); |
| BIG_zero(x[i]); |
| } |
| } |
| |
| /* x=y */ |
| static void FF_dscopy(BIG x[],BIG y[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| { |
| BIG_copy(x[i],y[i]); |
| BIG_zero(x[n+i]); |
| } |
| } |
| |
| /* x=y>>n */ |
| static void FF_sducopy(BIG x[],BIG y[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_copy(x[i],y[n+i]); |
| } |
| |
| /* set to zero */ |
| void FF_zero(BIG x[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_zero(x[i]); |
| } |
| |
| /* test equals 0 */ |
| int FF_iszilch(BIG x[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| if (!BIG_iszilch(x[i])) return 0; |
| return 1; |
| } |
| |
| /* shift right by BIGBITS-bit words */ |
| static void FF_shrw(BIG a[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| { |
| BIG_copy(a[i],a[i+n]); |
| BIG_zero(a[i+n]); |
| } |
| } |
| |
| /* shift left by BIGBITS-bit words */ |
| static void FF_shlw(BIG a[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| { |
| BIG_copy(a[i+n],a[i]); |
| BIG_zero(a[i]); |
| } |
| } |
| |
| /* extract last bit */ |
| int FF_parity(BIG x[]) |
| { |
| return BIG_parity(x[0]); |
| } |
| |
| /* extract last m bits */ |
| int FF_lastbits(BIG x[],int m) |
| { |
| return BIG_lastbits(x[0],m); |
| } |
| |
| /* x=1 */ |
| void FF_one(BIG x[],int n) |
| { |
| int i; |
| BIG_one(x[0]); |
| for (i=1; i<n; i++) |
| BIG_zero(x[i]); |
| } |
| |
| /* x=m, where m is 32-bit int */ |
| void FF_init(BIG x[],sign32 m,int n) |
| { |
| int i; |
| BIG_zero(x[0]); |
| #if CHUNK<64 |
| x[0][0]=(chunk)(m&BMASK); |
| x[0][1]=(chunk)(m>>BASEBITS); |
| #else |
| x[0][0]=(chunk)m; |
| #endif |
| for (i=1; i<n; i++) |
| BIG_zero(x[i]); |
| } |
| |
| /* compare x and y - must be normalised */ |
| int FF_comp(BIG x[],BIG y[],int n) |
| { |
| int i,j; |
| for (i=n-1; i>=0; i--) |
| { |
| j=BIG_comp(x[i],y[i]); |
| if (j!=0) return j; |
| } |
| return 0; |
| } |
| |
| /* recursive add */ |
| static void FF_radd(BIG z[],int zp,BIG x[],int xp,BIG y[],int yp,int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_add(z[zp+i],x[xp+i],y[yp+i]); |
| } |
| |
| /* recursive inc */ |
| static void FF_rinc(BIG z[],int zp,BIG y[],int yp,int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_add(z[zp+i],z[zp+i],y[yp+i]); |
| } |
| |
| /* recursive sub */ |
| /* |
| static void FF_rsub(BIG z[],int zp,BIG x[],int xp,BIG y[],int yp,int n) |
| { |
| int i; |
| for (i=0;i<n;i++) |
| BIG_sub(z[zp+i],x[xp+i],y[yp+i]); |
| } |
| */ |
| |
| /* recursive dec */ |
| static void FF_rdec(BIG z[],int zp,BIG y[],int yp,int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_sub(z[zp+i],z[zp+i],y[yp+i]); |
| } |
| |
| /* simple add */ |
| void FF_add(BIG z[],BIG x[],BIG y[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_add(z[i],x[i],y[i]); |
| } |
| |
| /* simple sub */ |
| void FF_sub(BIG z[],BIG x[],BIG y[],int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_sub(z[i],x[i],y[i]); |
| } |
| |
| /* increment/decrement by a small integer */ |
| void FF_inc(BIG x[],int m,int n) |
| { |
| BIG_inc(x[0],m); |
| FF_norm(x,n); |
| } |
| |
| void FF_dec(BIG x[],int m,int n) |
| { |
| BIG_dec(x[0],m); |
| FF_norm(x,n); |
| } |
| |
| /* normalise - but hold any overflow in top part unless n<0 */ |
| static void FF_rnorm(BIG z[],int zp,int n) |
| { |
| int i,trunc=0; |
| chunk carry; |
| if (n<0) |
| { |
| /* -v n signals to do truncation */ |
| n=-n; |
| trunc=1; |
| } |
| for (i=0; i<n-1; i++) |
| { |
| carry=BIG_norm(z[zp+i]); |
| |
| z[zp+i][NLEN-1]^=carry<<P_TBITS; /* remove it */ |
| z[zp+i+1][0]+=carry; |
| } |
| carry=BIG_norm(z[zp+n-1]); |
| if (trunc) z[zp+n-1][NLEN-1]^=carry<<P_TBITS; |
| } |
| |
| void FF_norm(BIG z[],int n) |
| { |
| FF_rnorm(z,0,n); |
| } |
| |
| /* shift left by one bit */ |
| void FF_shl(BIG x[],int n) |
| { |
| int i; |
| int carry,delay_carry=0; |
| for (i=0; i<n-1; i++) |
| { |
| carry=BIG_fshl(x[i],1); |
| x[i][0]|=delay_carry; |
| x[i][NLEN-1]^=(chunk)carry<<P_TBITS; |
| delay_carry=carry; |
| } |
| BIG_fshl(x[n-1],1); |
| x[n-1][0]|=delay_carry; |
| } |
| |
| /* shift right by one bit */ |
| void FF_shr(BIG x[],int n) |
| { |
| int i; |
| int carry; |
| for (i=n-1; i>0; i--) |
| { |
| carry=BIG_fshr(x[i],1); |
| x[i-1][NLEN-1]|=(chunk)carry<<P_TBITS; |
| } |
| BIG_fshr(x[0],1); |
| } |
| |
| void FF_output(BIG x[],int n) |
| { |
| int i; |
| FF_norm(x,n); |
| for (i=n-1; i>=0; i--) |
| { |
| BIG_output(x[i]); |
| printf(" "); |
| } |
| } |
| |
| void FF_rawoutput(BIG x[],int n) |
| { |
| int i; |
| for (i=n-1; i>=0; i--) |
| { |
| BIG_rawoutput(x[i]); |
| printf(" "); |
| } |
| } |
| |
| /* Convert FFs to/from octet strings */ |
| void FF_toOctet(octet *w,BIG x[],int n) |
| { |
| int i; |
| w->len=n*MODBYTES; |
| for (i=0; i<n; i++) |
| { |
| BIG_toBytes(&(w->val[(n-i-1)*MODBYTES]),x[i]); |
| } |
| } |
| |
| void FF_fromOctet(BIG x[],octet *w,int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| { |
| BIG_fromBytes(x[i],&(w->val[(n-i-1)*MODBYTES])); |
| } |
| } |
| |
| /* in-place swapping using xor - side channel resistant */ |
| static void FF_cswap(BIG a[],BIG b[],int d,int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| BIG_cswap(a[i],b[i],d); |
| return; |
| } |
| |
| /* z=x*y, t is workspace */ |
| static void FF_karmul(BIG z[],int zp,BIG x[],int xp,BIG y[],int yp,BIG t[],int tp,int n) |
| { |
| int nd2; |
| if (n==1) |
| { |
| BIG_mul(t[tp],x[xp],y[yp]); |
| BIG_split(z[zp+1],z[zp],t[tp],BIGBITS); |
| return; |
| } |
| |
| nd2=n/2; |
| FF_radd(z,zp,x,xp,x,xp+nd2,nd2); |
| FF_rnorm(z,zp,nd2); /* needs this if recursion level too deep */ |
| |
| FF_radd(z,zp+nd2,y,yp,y,yp+nd2,nd2); |
| FF_rnorm(z,zp+nd2,nd2); |
| FF_karmul(t,tp,z,zp,z,zp+nd2,t,tp+n,nd2); |
| FF_karmul(z,zp,x,xp,y,yp,t,tp+n,nd2); |
| FF_karmul(z,zp+n,x,xp+nd2,y,yp+nd2,t,tp+n,nd2); |
| FF_rdec(t,tp,z,zp,n); |
| FF_rdec(t,tp,z,zp+n,n); |
| FF_rinc(z,zp+nd2,t,tp,n); |
| FF_rnorm(z,zp,2*n); |
| } |
| |
| static void FF_karsqr(BIG z[],int zp,BIG x[],int xp,BIG t[],int tp,int n) |
| { |
| int nd2; |
| if (n==1) |
| { |
| BIG_sqr(t[tp],x[xp]); |
| BIG_split(z[zp+1],z[zp],t[tp],BIGBITS); |
| return; |
| } |
| nd2=n/2; |
| FF_karsqr(z,zp,x,xp,t,tp+n,nd2); |
| FF_karsqr(z,zp+n,x,xp+nd2,t,tp+n,nd2); |
| FF_karmul(t,tp,x,xp,x,xp+nd2,t,tp+n,nd2); |
| FF_rinc(z,zp+nd2,t,tp,n); |
| FF_rinc(z,zp+nd2,t,tp,n); |
| |
| FF_rnorm(z,zp+nd2,n); /* was FF_rnorm(z,zp,2*n) */ |
| } |
| |
| static void FF_karmul_lower(BIG z[],int zp,BIG x[],int xp,BIG y[],int yp,BIG t[],int tp,int n) |
| { |
| /* Calculates Least Significant bottom half of x*y */ |
| int nd2; |
| if (n==1) |
| { |
| /* only calculate bottom half of product */ |
| BIG_smul(z[zp],x[xp],y[yp]); |
| return; |
| } |
| nd2=n/2; |
| FF_karmul(z,zp,x,xp,y,yp,t,tp+n,nd2); |
| FF_karmul_lower(t,tp,x,xp+nd2,y,yp,t,tp+n,nd2); |
| FF_rinc(z,zp+nd2,t,tp,nd2); |
| FF_karmul_lower(t,tp,x,xp,y,yp+nd2,t,tp+n,nd2); |
| FF_rinc(z,zp+nd2,t,tp,nd2); |
| FF_rnorm(z,zp+nd2,-nd2); /* truncate it */ |
| } |
| |
| static void FF_karmul_upper(BIG z[],BIG x[],BIG y[],BIG t[],int n) |
| { |
| /* Calculates Most Significant upper half of x*y, given lower part */ |
| int nd2; |
| |
| nd2=n/2; |
| FF_radd(z,n,x,0,x,nd2,nd2); |
| FF_radd(z,n+nd2,y,0,y,nd2,nd2); |
| FF_rnorm(z,n,nd2); |
| FF_rnorm(z,n+nd2,nd2); |
| |
| FF_karmul(t,0,z,n+nd2,z,n,t,n,nd2); /* t = (a0+a1)(b0+b1) */ |
| FF_karmul(z,n,x,nd2,y,nd2,t,n,nd2); /* z[n]= a1*b1 */ |
| /* z[0-nd2]=l(a0b0) z[nd2-n]= h(a0b0)+l(t)-l(a0b0)-l(a1b1) */ |
| FF_rdec(t,0,z,n,n); /* t=t-a1b1 */ |
| FF_rinc(z,nd2,z,0,nd2); /* z[nd2-n]+=l(a0b0) = h(a0b0)+l(t)-l(a1b1) */ |
| FF_rdec(z,nd2,t,0,nd2); /* z[nd2-n]=h(a0b0)+l(t)-l(a1b1)-l(t-a1b1)=h(a0b0) */ |
| FF_rnorm(z,0,-n); /* a0b0 now in z - truncate it */ |
| FF_rdec(t,0,z,0,n); /* (a0+a1)(b0+b1) - a0b0 */ |
| FF_rinc(z,nd2,t,0,n); |
| |
| FF_rnorm(z,nd2,n); |
| } |
| |
| /* z=x*y */ |
| void FF_mul(BIG z[],BIG x[],BIG y[],int n) |
| { |
| #ifndef C99 |
| BIG t[2*FFLEN]; |
| #else |
| BIG t[2*n]; |
| #endif |
| // FF_norm(x,n); /* change here */ |
| // FF_norm(y,n); /* change here */ |
| FF_karmul(z,0,x,0,y,0,t,0,n); |
| } |
| |
| /* return low part of product */ |
| static void FF_lmul(BIG z[],BIG x[],BIG y[],int n) |
| { |
| #ifndef C99 |
| BIG t[2*FFLEN]; |
| #else |
| BIG t[2*n]; |
| #endif |
| // FF_norm(x,n); /* change here */ |
| // FF_norm(y,n); /* change here */ |
| FF_karmul_lower(z,0,x,0,y,0,t,0,n); |
| } |
| |
| /* Set b=b mod c */ |
| void FF_mod(BIG b[],BIG c[],int n) |
| { |
| int k=0; |
| |
| FF_norm(b,n); |
| if (FF_comp(b,c,n)<0) |
| return; |
| do |
| { |
| FF_shl(c,n); |
| k++; |
| } |
| while (FF_comp(b,c,n)>=0); |
| |
| while (k>0) |
| { |
| FF_shr(c,n); |
| if (FF_comp(b,c,n)>=0) |
| { |
| FF_sub(b,b,c,n); |
| FF_norm(b,n); |
| } |
| k--; |
| } |
| } |
| |
| /* z=x^2 */ |
| void FF_sqr(BIG z[],BIG x[],int n) |
| { |
| #ifndef C99 |
| BIG t[2*FFLEN]; |
| #else |
| BIG t[2*n]; |
| #endif |
| // FF_norm(x,n); /* change here */ |
| FF_karsqr(z,0,x,0,t,0,n); |
| } |
| |
| /* r=t mod modulus, N is modulus, ND is Montgomery Constant */ |
| static void FF_reduce(BIG r[],BIG T[],BIG N[],BIG ND[],int n) |
| { |
| /* fast karatsuba Montgomery reduction */ |
| #ifndef C99 |
| BIG t[2*FFLEN]; |
| BIG m[FFLEN]; |
| #else |
| BIG t[2*n]; |
| BIG m[n]; |
| #endif |
| FF_sducopy(r,T,n); /* keep top half of T */ |
| //FF_norm(T,n); /* change here */ |
| FF_karmul_lower(m,0,T,0,ND,0,t,0,n); /* m=T.(1/N) mod R */ |
| |
| //FF_norm(N,n); /* change here */ |
| FF_karmul_upper(T,N,m,t,n); /* T=mN */ |
| FF_sducopy(m,T,n); |
| |
| FF_add(r,r,N,n); |
| FF_sub(r,r,m,n); |
| FF_norm(r,n); |
| } |
| |
| |
| /* Set r=a mod b */ |
| /* a is of length - 2*n */ |
| /* r,b is of length - n */ |
| void FF_dmod(BIG r[],BIG a[],BIG b[],int n) |
| { |
| int k; |
| #ifndef C99 |
| BIG m[2*FFLEN]; |
| BIG x[2*FFLEN]; |
| #else |
| BIG m[2*n]; |
| BIG x[2*n]; |
| #endif |
| FF_copy(x,a,2*n); |
| FF_norm(x,2*n); |
| FF_dsucopy(m,b,n); |
| k=BIGBITS*n; |
| |
| while (FF_comp(x,m,2*n)>=0) |
| { |
| FF_sub(x,x,m,2*n); |
| FF_norm(x,2*n); |
| } |
| |
| while (k>0) |
| { |
| FF_shr(m,2*n); |
| |
| if (FF_comp(x,m,2*n)>=0) |
| { |
| FF_sub(x,x,m,2*n); |
| FF_norm(x,2*n); |
| } |
| |
| k--; |
| } |
| FF_copy(r,x,n); |
| FF_mod(r,b,n); |
| } |
| |
| /* Set r=1/a mod p. Binary method - a<p on entry */ |
| |
| void FF_invmodp(BIG r[],BIG a[],BIG p[],int n) |
| { |
| #ifndef C99 |
| BIG u[FFLEN],v[FFLEN],x1[FFLEN],x2[FFLEN],t[FFLEN],one[FFLEN]; |
| #else |
| BIG u[n],v[n],x1[n],x2[n],t[n],one[n]; |
| #endif |
| FF_copy(u,a,n); |
| FF_copy(v,p,n); |
| FF_one(one,n); |
| FF_copy(x1,one,n); |
| FF_zero(x2,n); |
| |
| // reduce n in here as well! |
| while (FF_comp(u,one,n)!=0 && FF_comp(v,one,n)!=0) |
| { |
| while (FF_parity(u)==0) |
| { |
| FF_shr(u,n); |
| if (FF_parity(x1)!=0) |
| { |
| FF_add(x1,p,x1,n); |
| FF_norm(x1,n); |
| } |
| FF_shr(x1,n); |
| } |
| while (FF_parity(v)==0) |
| { |
| FF_shr(v,n); |
| if (FF_parity(x2)!=0) |
| { |
| FF_add(x2,p,x2,n); |
| FF_norm(x2,n); |
| } |
| FF_shr(x2,n); |
| } |
| if (FF_comp(u,v,n)>=0) |
| { |
| |
| FF_sub(u,u,v,n); |
| FF_norm(u,n); |
| if (FF_comp(x1,x2,n)>=0) FF_sub(x1,x1,x2,n); |
| else |
| { |
| FF_sub(t,p,x2,n); |
| FF_add(x1,x1,t,n); |
| } |
| FF_norm(x1,n); |
| } |
| else |
| { |
| FF_sub(v,v,u,n); |
| FF_norm(v,n); |
| if (FF_comp(x2,x1,n)>=0) FF_sub(x2,x2,x1,n); |
| else |
| { |
| FF_sub(t,p,x1,n); |
| FF_add(x2,x2,t,n); |
| } |
| FF_norm(x2,n); |
| } |
| } |
| if (FF_comp(u,one,n)==0) |
| FF_copy(r,x1,n); |
| else |
| FF_copy(r,x2,n); |
| } |
| |
| /* nesidue mod m */ |
| static void FF_nres(BIG a[],BIG m[],int n) |
| { |
| #ifndef C99 |
| BIG d[2*FFLEN]; |
| #else |
| BIG d[2*n]; |
| #endif |
| |
| if (n==1) |
| { |
| BIG_dscopy(d[0],a[0]); |
| BIG_dshl(d[0],NLEN*BASEBITS); |
| BIG_dmod(a[0],d[0],m[0]); |
| } |
| else |
| { |
| FF_dsucopy(d,a,n); |
| FF_dmod(a,d,m,n); |
| } |
| } |
| |
| static void FF_redc(BIG a[],BIG m[],BIG ND[],int n) |
| { |
| #ifndef C99 |
| BIG d[2*FFLEN]; |
| #else |
| BIG d[2*n]; |
| #endif |
| if (n==1) |
| { |
| BIG_dzero(d[0]); |
| BIG_dscopy(d[0],a[0]); |
| BIG_monty(a[0],m[0],((chunk)1<<BASEBITS)-ND[0][0],d[0]); |
| } |
| else |
| { |
| FF_mod(a,m,n); |
| FF_dscopy(d,a,n); |
| FF_reduce(a,d,m,ND,n); |
| FF_mod(a,m,n); |
| } |
| } |
| |
| /* U=1/a mod 2^m - Arazi & Qi */ |
| static void FF_invmod2m(BIG U[],BIG a[],int n) |
| { |
| int i; |
| #ifndef C99 |
| BIG t1[FFLEN],b[FFLEN],c[FFLEN]; |
| #else |
| BIG t1[2*n],b[n],c[n]; |
| #endif |
| |
| FF_zero(U,n); |
| FF_zero(b,n); |
| FF_zero(c,n); |
| FF_zero(t1,2*n); |
| |
| BIG_copy(U[0],a[0]); |
| BIG_invmod2m(U[0]); |
| for (i=1; i<n; i<<=1) |
| { |
| FF_copy(b,a,i); |
| FF_mul(t1,U,b,i); |
| FF_shrw(t1,i); // top half to bottom half, top half=0 |
| |
| FF_copy(c,a,2*i); |
| FF_shrw(c,i); // top half of c |
| FF_lmul(b,U,c,i); // should set top half of b=0 |
| FF_add(t1,t1,b,i); |
| FF_norm(t1,2*i); |
| FF_lmul(b,t1,U,i); |
| FF_copy(t1,b,i); |
| FF_one(b,i); |
| FF_shlw(b,i); |
| FF_sub(t1,b,t1,2*i); |
| FF_norm(t1,2*i); |
| FF_shlw(t1,i); |
| FF_add(U,U,t1,2*i); |
| } |
| |
| FF_norm(U,n); |
| } |
| |
| void FF_random(BIG x[],csprng *rng,int n) |
| { |
| int i; |
| for (i=0; i<n; i++) |
| { |
| BIG_random(x[i],rng); |
| } |
| /* make sure top bit is 1 */ |
| while (BIG_nbits(x[n-1])<MODBYTES*8) BIG_random(x[n-1],rng); |
| } |
| |
| /* generate random x mod p */ |
| void FF_randomnum(BIG x[],BIG p[],csprng *rng,int n) |
| { |
| int i; |
| #ifndef C99 |
| BIG d[2*FFLEN]; |
| #else |
| BIG d[2*n]; |
| #endif |
| for (i=0; i<2*n; i++) |
| { |
| BIG_random(d[i],rng); |
| } |
| FF_dmod(x,d,p,n); |
| } |
| |
| static void FF_modmul(BIG z[],BIG x[],BIG y[],BIG p[],BIG ND[],int n) |
| { |
| #ifndef C99 |
| BIG d[2*FFLEN]; |
| #else |
| BIG d[2*n]; |
| #endif |
| chunk ex=P_EXCESS(x[n-1]); |
| chunk ey=P_EXCESS(y[n-1]); |
| #ifdef dchunk |
| if ((dchunk)(ex+1)*(ey+1)>(dchunk)P_FEXCESS) |
| #else |
| if ((ex+1)>P_FEXCESS/(ey+1)) |
| #endif |
| { |
| #ifdef DEBUG_REDUCE |
| printf("Product too large - reducing it %d %d\n",ex,ey); |
| #endif |
| FF_mod(x,p,n); |
| } |
| |
| if (n==1) |
| { |
| BIG_mul(d[0],x[0],y[0]); |
| BIG_monty(z[0],p[0],((chunk)1<<BASEBITS)-ND[0][0],d[0]); |
| } |
| else |
| { |
| FF_mul(d,x,y,n); |
| FF_reduce(z,d,p,ND,n); |
| } |
| } |
| |
| static void FF_modsqr(BIG z[],BIG x[],BIG p[],BIG ND[],int n) |
| { |
| #ifndef C99 |
| BIG d[2*FFLEN]; |
| #else |
| BIG d[2*n]; |
| #endif |
| chunk ex=P_EXCESS(x[n-1]); |
| #ifdef dchunk |
| if ((dchunk)(ex+1)*(ex+1)>(dchunk)P_FEXCESS) |
| #else |
| if ((ex+1)>P_FEXCESS/(ex+1)) |
| #endif |
| { |
| #ifdef DEBUG_REDUCE |
| printf("Product too large - reducing it %d\n",ex); |
| #endif |
| FF_mod(x,p,n); |
| } |
| |
| if (n==1) |
| { |
| BIG_sqr(d[0],x[0]); |
| BIG_monty(z[0],p[0],((chunk)1<<BASEBITS)-ND[0][0],d[0]); |
| } |
| else |
| { |
| FF_sqr(d,x,n); |
| FF_reduce(z,d,p,ND,n); |
| } |
| } |
| |
| /* r=x^e mod p using side-channel resistant Montgomery Ladder, for large e */ |
| void FF_skpow(BIG r[],BIG x[],BIG e[],BIG p[],int n) |
| { |
| int i,b; |
| #ifndef C99 |
| BIG R0[FFLEN],R1[FFLEN],ND[FFLEN]; |
| #else |
| BIG R0[n],R1[n],ND[n]; |
| #endif |
| FF_invmod2m(ND,p,n); |
| |
| FF_one(R0,n); |
| FF_copy(R1,x,n); |
| FF_nres(R0,p,n); |
| FF_nres(R1,p,n); |
| |
| for (i=8*MODBYTES*n-1; i>=0; i--) |
| { |
| b=BIG_bit(e[i/BIGBITS],i%BIGBITS); |
| FF_modmul(r,R0,R1,p,ND,n); |
| |
| FF_cswap(R0,R1,b,n); |
| FF_modsqr(R0,R0,p,ND,n); |
| |
| FF_copy(R1,r,n); |
| FF_cswap(R0,R1,b,n); |
| } |
| FF_copy(r,R0,n); |
| FF_redc(r,p,ND,n); |
| } |
| |
| /* r=x^e mod p using side-channel resistant Montgomery Ladder, for short e */ |
| void FF_skspow(BIG r[],BIG x[],BIG e,BIG p[],int n) |
| { |
| int i,b; |
| #ifndef C99 |
| BIG R0[FFLEN],R1[FFLEN],ND[FFLEN]; |
| #else |
| BIG R0[n],R1[n],ND[n]; |
| #endif |
| FF_invmod2m(ND,p,n); |
| FF_one(R0,n); |
| FF_copy(R1,x,n); |
| FF_nres(R0,p,n); |
| FF_nres(R1,p,n); |
| for (i=8*MODBYTES-1; i>=0; i--) |
| { |
| b=BIG_bit(e,i); |
| FF_modmul(r,R0,R1,p,ND,n); |
| FF_cswap(R0,R1,b,n); |
| FF_modsqr(R0,R0,p,ND,n); |
| FF_copy(R1,r,n); |
| FF_cswap(R0,R1,b,n); |
| } |
| FF_copy(r,R0,n); |
| FF_redc(r,p,ND,n); |
| } |
| |
| /* raise to an integer power - right-to-left method */ |
| void FF_power(BIG r[],BIG x[],int e,BIG p[],int n) |
| { |
| int f=1; |
| #ifndef C99 |
| BIG w[FFLEN],ND[FFLEN]; |
| #else |
| BIG w[n],ND[n]; |
| #endif |
| FF_invmod2m(ND,p,n); |
| |
| FF_copy(w,x,n); |
| FF_nres(w,p,n); |
| |
| if (e==2) |
| { |
| FF_modsqr(r,w,p,ND,n); |
| } |
| else for (;;) |
| { |
| if (e%2==1) |
| { |
| if (f) FF_copy(r,w,n); |
| else FF_modmul(r,r,w,p,ND,n); |
| f=0; |
| } |
| e>>=1; |
| if (e==0) break; |
| FF_modsqr(w,w,p,ND,n); |
| } |
| |
| FF_redc(r,p,ND,n); |
| } |
| |
| /* r=x^e mod p, faster but not side channel resistant */ |
| void FF_pow(BIG r[],BIG x[],BIG e[],BIG p[],int n) |
| { |
| int i,b; |
| #ifndef C99 |
| BIG w[FFLEN],ND[FFLEN]; |
| #else |
| BIG w[n],ND[n]; |
| #endif |
| FF_invmod2m(ND,p,n); |
| |
| FF_copy(w,x,n); |
| FF_one(r,n); |
| FF_nres(r,p,n); |
| FF_nres(w,p,n); |
| |
| for (i=8*MODBYTES*n-1; i>=0; i--) |
| { |
| FF_modsqr(r,r,p,ND,n); |
| b=BIG_bit(e[i/BIGBITS],i%BIGBITS); |
| if (b==1) FF_modmul(r,r,w,p,ND,n); |
| } |
| FF_redc(r,p,ND,n); |
| } |
| |
| /* double exponentiation r=x^e.y^f mod p */ |
| void FF_pow2(BIG r[],BIG x[],BIG e,BIG y[],BIG f,BIG p[],int n) |
| { |
| int i,eb,fb; |
| #ifndef C99 |
| BIG xn[FFLEN],yn[FFLEN],xy[FFLEN],ND[FFLEN]; |
| #else |
| BIG xn[n],yn[n],xy[n],ND[n]; |
| #endif |
| |
| FF_invmod2m(ND,p,n); |
| |
| FF_copy(xn,x,n); |
| FF_copy(yn,y,n); |
| FF_nres(xn,p,n); |
| FF_nres(yn,p,n); |
| FF_modmul(xy,xn,yn,p,ND,n); |
| FF_one(r,n); |
| FF_nres(r,p,n); |
| |
| for (i=8*MODBYTES-1; i>=0; i--) |
| { |
| eb=BIG_bit(e,i); |
| fb=BIG_bit(f,i); |
| FF_modsqr(r,r,p,ND,n); |
| if (eb==1) |
| { |
| if (fb==1) FF_modmul(r,r,xy,p,ND,n); |
| else FF_modmul(r,r,xn,p,ND,n); |
| } |
| else |
| { |
| if (fb==1) FF_modmul(r,r,yn,p,ND,n); |
| } |
| } |
| FF_redc(r,p,ND,n); |
| } |
| |
| static sign32 igcd(sign32 x,sign32 y) |
| { |
| /* integer GCD, returns GCD of x and y */ |
| sign32 r; |
| if (y==0) return x; |
| while ((r=x%y)!=0) |
| x=y,y=r; |
| return y; |
| } |
| |
| /* quick and dirty check for common factor with s */ |
| int FF_cfactor(BIG w[],sign32 s,int n) |
| { |
| int r; |
| sign32 g; |
| #ifndef C99 |
| BIG x[FFLEN],y[FFLEN]; |
| #else |
| BIG x[n],y[n]; |
| #endif |
| FF_init(y,s,n); |
| FF_copy(x,w,n); |
| FF_norm(x,n); |
| |
| // if (FF_parity(x)==0) return 1; |
| do |
| { |
| FF_sub(x,x,y,n); |
| FF_norm(x,n); |
| while (!FF_iszilch(x,n) && FF_parity(x)==0) FF_shr(x,n); |
| } |
| while (FF_comp(x,y,n)>0); |
| #if CHUNK<32 |
| g=x[0][0]+((sign32)(x[0][1])<<BASEBITS); |
| #else |
| g=(sign32)x[0][0]; |
| #endif |
| r=igcd(s,g); |
| if (r>1) return 1; |
| return 0; |
| } |
| |
| /* Miller-Rabin test for primality. Slow. */ |
| int FF_prime(BIG p[],csprng *rng,int n) |
| { |
| int i,j,loop,s=0; |
| #ifndef C99 |
| BIG d[FFLEN],x[FFLEN],unity[FFLEN],nm1[FFLEN]; |
| #else |
| BIG d[n],x[n],unity[n],nm1[n]; |
| #endif |
| sign32 sf=4849845;/* 3*5*.. *19 */ |
| |
| FF_norm(p,n); |
| |
| if (FF_cfactor(p,sf,n)) return 0; |
| |
| FF_one(unity,n); |
| FF_sub(nm1,p,unity,n); |
| FF_norm(nm1,n); |
| FF_copy(d,nm1,n); |
| while (FF_parity(d)==0) |
| { |
| FF_shr(d,n); |
| s++; |
| } |
| if (s==0) return 0; |
| |
| for (i=0; i<10; i++) |
| { |
| FF_randomnum(x,p,rng,n); |
| FF_pow(x,x,d,p,n); |
| if (FF_comp(x,unity,n)==0 || FF_comp(x,nm1,n)==0) continue; |
| loop=0; |
| for (j=1; j<s; j++) |
| { |
| FF_power(x,x,2,p,n); |
| if (FF_comp(x,unity,n)==0) return 0; |
| if (FF_comp(x,nm1,n)==0 ) |
| { |
| loop=1; |
| break; |
| } |
| } |
| if (loop) continue; |
| return 0; |
| } |
| |
| return 1; |
| } |
| |
| /* |
| BIG P[4]= {{0x1670957,0x1568CD3C,0x2595E5,0xEED4F38,0x1FC9A971,0x14EF7E62,0xA503883,0x9E1E05E,0xBF59E3},{0x1844C908,0x1B44A798,0x3A0B1E7,0xD1B5B4E,0x1836046F,0x87E94F9,0x1D34C537,0xF7183B0,0x46D07},{0x17813331,0x19E28A90,0x1473A4D6,0x1CACD01F,0x1EEA8838,0xAF2AE29,0x1F85292A,0x1632585E,0xD945E5},{0x919F5EF,0x1567B39F,0x19F6AD11,0x16CE47CF,0x9B36EB1,0x35B7D3,0x483B28C,0xCBEFA27,0xB5FC21}}; |
| |
| int main() |
| { |
| int i; |
| BIG p[4],e[4],x[4],r[4]; |
| csprng rng; |
| char raw[100]; |
| for (i=0;i<100;i++) raw[i]=i; |
| RAND_seed(&rng,100,raw); |
| |
| |
| FF_init(x,3,4); |
| |
| FF_copy(p,P,4); |
| FF_copy(e,p,4); |
| FF_dec(e,1,4); |
| FF_norm(e,4); |
| |
| |
| |
| printf("p= ");FF_output(p,4); printf("\n"); |
| if (FF_prime(p,&rng,4)) printf("p is a prime\n"); |
| printf("e= ");FF_output(e,4); printf("\n"); |
| |
| FF_skpow(r,x,e,p,4); |
| printf("r= ");FF_output(r,4); printf("\n"); |
| } |
| |
| */ |