/*
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");
}

*/