/*
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 Weierstrass elliptic curve functions over FP2 */
/* SU=m, m is Stack Usage */

#include "ecp2_FP256BN.h"

int ECP2_FP256BN_isinf(ECP2_FP256BN *P)
{
    if (P->inf) return 1;
    P->inf=FP2_FP256BN_iszilch(&(P->x)) & FP2_FP256BN_iszilch(&(P->z));
    return P->inf;
}

/* Set P=Q */
/* SU= 16 */
void ECP2_FP256BN_copy(ECP2_FP256BN *P,ECP2_FP256BN *Q)
{
    P->inf=Q->inf;
    FP2_FP256BN_copy(&(P->x),&(Q->x));
    FP2_FP256BN_copy(&(P->y),&(Q->y));
    FP2_FP256BN_copy(&(P->z),&(Q->z));
}

/* set P to Infinity */
/* SU= 8 */
void ECP2_FP256BN_inf(ECP2_FP256BN *P)
{
    P->inf=1;
    FP2_FP256BN_zero(&(P->x));
    FP2_FP256BN_one(&(P->y));
    FP2_FP256BN_zero(&(P->z));
}

/* Conditional move Q to P dependant on d */
static void ECP2_FP256BN_cmove(ECP2_FP256BN *P,ECP2_FP256BN *Q,int d)
{
    FP2_FP256BN_cmove(&(P->x),&(Q->x),d);
    FP2_FP256BN_cmove(&(P->y),&(Q->y),d);
    FP2_FP256BN_cmove(&(P->z),&(Q->z),d);
    d=~(d-1);
    P->inf^=(P->inf^Q->inf)&d;
}

/* return 1 if b==c, no branching */
static int teq(sign32 b,sign32 c)
{
    sign32 x=b^c;
    x-=1;  // if x=0, x now -1
    return (int)((x>>31)&1);
}

/* Constant time select from pre-computed table */
static void ECP2_FP256BN_select(ECP2_FP256BN *P,ECP2_FP256BN W[],sign32 b)
{
    ECP2_FP256BN MP;
    sign32 m=b>>31;
    sign32 babs=(b^m)-m;

    babs=(babs-1)/2;

    ECP2_FP256BN_cmove(P,&W[0],teq(babs,0));  // conditional move
    ECP2_FP256BN_cmove(P,&W[1],teq(babs,1));
    ECP2_FP256BN_cmove(P,&W[2],teq(babs,2));
    ECP2_FP256BN_cmove(P,&W[3],teq(babs,3));
    ECP2_FP256BN_cmove(P,&W[4],teq(babs,4));
    ECP2_FP256BN_cmove(P,&W[5],teq(babs,5));
    ECP2_FP256BN_cmove(P,&W[6],teq(babs,6));
    ECP2_FP256BN_cmove(P,&W[7],teq(babs,7));

    ECP2_FP256BN_copy(&MP,P);
    ECP2_FP256BN_neg(&MP);  // minus P
    ECP2_FP256BN_cmove(P,&MP,(int)(m&1));
}

/* return 1 if P==Q, else 0 */
/* SU= 312 */
int ECP2_FP256BN_equals(ECP2_FP256BN *P,ECP2_FP256BN *Q)
{
    FP2_FP256BN a,b;
    if (ECP2_FP256BN_isinf(P) && ECP2_FP256BN_isinf(Q)) return 1;
    if (ECP2_FP256BN_isinf(P) || ECP2_FP256BN_isinf(Q)) return 0;

    FP2_FP256BN_mul(&a,&(P->x),&(Q->z));
    FP2_FP256BN_mul(&b,&(Q->x),&(P->z));
    if (!FP2_FP256BN_equals(&a,&b)) return 0;

    FP2_FP256BN_mul(&a,&(P->y),&(Q->z));
    FP2_FP256BN_mul(&b,&(Q->y),&(P->z));
    if (!FP2_FP256BN_equals(&a,&b)) return 0;
    return 1;
}

/* Make P affine (so z=1) */
/* SU= 232 */
void ECP2_FP256BN_affine(ECP2_FP256BN *P)
{
    FP2_FP256BN one,iz;
    if (ECP2_FP256BN_isinf(P)) return;

    FP2_FP256BN_one(&one);
    if (FP2_FP256BN_isunity(&(P->z)))
    {
        FP2_FP256BN_reduce(&(P->x));
        FP2_FP256BN_reduce(&(P->y));
        return;
    }

    FP2_FP256BN_inv(&iz,&(P->z));
    FP2_FP256BN_mul(&(P->x),&(P->x),&iz);
    FP2_FP256BN_mul(&(P->y),&(P->y),&iz);

    FP2_FP256BN_reduce(&(P->x));
    FP2_FP256BN_reduce(&(P->y));
    FP2_FP256BN_copy(&(P->z),&one);
}

/* extract x, y from point P */
/* SU= 16 */
int ECP2_FP256BN_get(FP2_FP256BN *x,FP2_FP256BN *y,ECP2_FP256BN *P)
{
    if (P->inf) return -1;
    ECP2_FP256BN_affine(P);
    FP2_FP256BN_copy(y,&(P->y));
    FP2_FP256BN_copy(x,&(P->x));
    return 0;
}

/* SU= 152 */
/* Output point P */
void ECP2_FP256BN_output(ECP2_FP256BN *P)
{
    FP2_FP256BN x,y;
    if (P->inf)
    {
        printf("Infinity\n");
        return;
    }
    ECP2_FP256BN_get(&x,&y,P);
    printf("(");
    FP2_FP256BN_output(&x);
    printf(",");
    FP2_FP256BN_output(&y);
    printf(")\n");
}

/* SU= 232 */
void ECP2_FP256BN_outputxyz(ECP2_FP256BN *P)
{
    ECP2_FP256BN Q;
    if (P->inf)
    {
        printf("Infinity\n");
        return;
    }
    ECP2_FP256BN_copy(&Q,P);
    printf("(");
    FP2_FP256BN_output(&(Q.x));
    printf(",");
    FP2_FP256BN_output(&(Q.y));
    printf(",");
    FP2_FP256BN_output(&(Q.z));
    printf(")\n");
}

/* SU= 168 */
/* Convert Q to octet string */
void ECP2_FP256BN_toOctet(octet *W,ECP2_FP256BN *Q)
{
    BIG_256_56 b;
    FP2_FP256BN qx,qy;
    ECP2_FP256BN_get(&qx,&qy,Q);

    FP_FP256BN_redc(b,&(qx.a));
    BIG_256_56_toBytes(&(W->val[0]),b);
    FP_FP256BN_redc(b,&(qx.b));
    BIG_256_56_toBytes(&(W->val[MODBYTES_256_56]),b);
    FP_FP256BN_redc(b,&(qy.a));
    BIG_256_56_toBytes(&(W->val[2*MODBYTES_256_56]),b);
    FP_FP256BN_redc(b,&(qy.b));
    BIG_256_56_toBytes(&(W->val[3*MODBYTES_256_56]),b);

    W->len=4*MODBYTES_256_56;

}

/* SU= 176 */
/* restore Q from octet string */
int ECP2_FP256BN_fromOctet(ECP2_FP256BN *Q,octet *W)
{
    BIG_256_56 b;
    FP2_FP256BN qx,qy;
    BIG_256_56_fromBytes(b,&(W->val[0]));
    FP_FP256BN_nres(&(qx.a),b);
    BIG_256_56_fromBytes(b,&(W->val[MODBYTES_256_56]));
    FP_FP256BN_nres(&(qx.b),b);
    BIG_256_56_fromBytes(b,&(W->val[2*MODBYTES_256_56]));
    FP_FP256BN_nres(&(qy.a),b);
    BIG_256_56_fromBytes(b,&(W->val[3*MODBYTES_256_56]));
    FP_FP256BN_nres(&(qy.b),b);

    if (ECP2_FP256BN_set(Q,&qx,&qy)) return 1;
    return 0;
}

/* SU= 128 */
/* Calculate RHS of twisted curve equation x^3+B/i or x^3+Bi*/
void ECP2_FP256BN_rhs(FP2_FP256BN *rhs,FP2_FP256BN *x)
{
    /* calculate RHS of elliptic curve equation */
    FP2_FP256BN t;
    BIG_256_56 b;
    FP2_FP256BN_sqr(&t,x);

    FP2_FP256BN_mul(rhs,&t,x);

    /* Assuming CURVE_A=0 */

    BIG_256_56_rcopy(b,CURVE_B_FP256BN);

    FP2_FP256BN_from_BIG(&t,b);

#if SEXTIC_TWIST_FP256BN == D_TYPE
    FP2_FP256BN_div_ip(&t);   /* IMPORTANT - here we use the correct SEXTIC twist of the curve */
#endif

#if SEXTIC_TWIST_FP256BN == M_TYPE
    FP2_FP256BN_norm(&t);
    FP2_FP256BN_mul_ip(&t);   /* IMPORTANT - here we use the correct SEXTIC twist of the curve */
    FP2_FP256BN_norm(&t);

#endif


    FP2_FP256BN_add(rhs,&t,rhs);
    FP2_FP256BN_reduce(rhs);
}


/* Set P=(x,y). Return 1 if (x,y) is on the curve, else return 0*/
/* SU= 232 */
int ECP2_FP256BN_set(ECP2_FP256BN *P,FP2_FP256BN *x,FP2_FP256BN *y)
{
    FP2_FP256BN rhs,y2;

    FP2_FP256BN_sqr(&y2,y);
    ECP2_FP256BN_rhs(&rhs,x);

    if (!FP2_FP256BN_equals(&y2,&rhs))
    {
        P->inf=1;
        return 0;
    }

    P->inf=0;
    FP2_FP256BN_copy(&(P->x),x);
    FP2_FP256BN_copy(&(P->y),y);

    FP2_FP256BN_one(&(P->z));
    return 1;
}

/* Set P=(x,y). Return 1 if (x,.) is on the curve, else return 0 */
/* SU= 232 */
int ECP2_FP256BN_setx(ECP2_FP256BN *P,FP2_FP256BN *x)
{
    FP2_FP256BN y;
    ECP2_FP256BN_rhs(&y,x);

    if (!FP2_FP256BN_sqrt(&y,&y))
    {
        P->inf=1;
        return 0;
    }

    P->inf=0;
    FP2_FP256BN_copy(&(P->x),x);
    FP2_FP256BN_copy(&(P->y),&y);
    FP2_FP256BN_one(&(P->z));
    return 1;
}

/* Set P=-P */
/* SU= 8 */
void ECP2_FP256BN_neg(ECP2_FP256BN *P)
{
//	if (ECP2_FP256BN_isinf(P)) return;
    FP2_FP256BN_norm(&(P->y));
    FP2_FP256BN_neg(&(P->y),&(P->y));
    FP2_FP256BN_norm(&(P->y));
}

/* R+=R */
/* return -1 for Infinity, 0 for addition, 1 for doubling */
/* SU= 448 */
int ECP2_FP256BN_dbl(ECP2_FP256BN *P)
{
    FP2_FP256BN t0,t1,t2,iy,x3,y3;
    if (P->inf) return -1;

    FP2_FP256BN_copy(&iy,&(P->y));		//FP2 iy=new FP2(y);
#if SEXTIC_TWIST_FP256BN==D_TYPE
    FP2_FP256BN_mul_ip(&iy);			//iy.mul_ip();
    FP2_FP256BN_norm(&iy);				//iy.norm();
#endif
    //FP2_FP256BN_copy(&t0,&(P->y));		//FP2 t0=new FP2(y);                  //***** Change
    FP2_FP256BN_sqr(&t0,&(P->y));			//t0.sqr();
#if SEXTIC_TWIST_FP256BN==D_TYPE
    FP2_FP256BN_mul_ip(&t0);			//t0.mul_ip();
#endif
    //FP2_FP256BN_copy(&t1,&iy);				//FP2 t1=new FP2(iy);
    FP2_FP256BN_mul(&t1,&iy,&(P->z));	//t1.mul(z);
    //FP2_FP256BN_copy(&t2,&(P->z));		//FP2 t2=new FP2(z);
    FP2_FP256BN_sqr(&t2,&(P->z));				//t2.sqr();

    //FP2_FP256BN_copy(&(P->z),&t0);		//z.copy(t0);
    FP2_FP256BN_add(&(P->z),&t0,&t0);	//z.add(t0);
    FP2_FP256BN_norm(&(P->z));				//z.norm();
    FP2_FP256BN_add(&(P->z),&(P->z),&(P->z));	//z.add(z);
    FP2_FP256BN_add(&(P->z),&(P->z),&(P->z));	//z.add(z);
    FP2_FP256BN_norm(&(P->z));			//z.norm();

    FP2_FP256BN_imul(&t2,&t2,3*CURVE_B_I_FP256BN);	//t2.imul(3*ROM.CURVE_B_I);
#if SEXTIC_TWIST_FP256BN==M_TYPE
    FP2_FP256BN_mul_ip(&t2);
    FP2_FP256BN_norm(&t2);
#endif

    //FP2_FP256BN_copy(&x3,&t2);			//FP2 x3=new FP2(t2);
    FP2_FP256BN_mul(&x3,&t2,&(P->z));	//x3.mul(z);

    //FP2_FP256BN_copy(&y3,&t0);			//FP2 y3=new FP2(t0);

    FP2_FP256BN_add(&y3,&t0,&t2);		//y3.add(t2);
    FP2_FP256BN_norm(&y3);				//y3.norm();
    FP2_FP256BN_mul(&(P->z),&(P->z),&t1);	//z.mul(t1);

    //FP2_FP256BN_copy(&t1,&t2);			//t1.copy(t2);
    FP2_FP256BN_add(&t1,&t2,&t2);		//t1.add(t2);
    FP2_FP256BN_add(&t2,&t2,&t1);		//t2.add(t1);
    FP2_FP256BN_norm(&t2);				//t2.norm();
    FP2_FP256BN_sub(&t0,&t0,&t2);		//t0.sub(t2);
    FP2_FP256BN_norm(&t0);				//t0.norm();                           //y^2-9bz^2
    FP2_FP256BN_mul(&y3,&y3,&t0);		//y3.mul(t0);
    FP2_FP256BN_add(&(P->y),&y3,&x3);		//y3.add(x3);                          //(y^2+3z*2)(y^2-9z^2)+3b.z^2.8y^2
    //FP2_FP256BN_copy(&t1,&(P->x));		//t1.copy(x);
    FP2_FP256BN_mul(&t1,&(P->x),&iy);		//t1.mul(iy);						//
    //FP2_FP256BN_copy(&(P->x),&t0);		//x.copy(t0);
    FP2_FP256BN_norm(&t0);			//x.norm();
    FP2_FP256BN_mul(&(P->x),&t0,&t1);	//x.mul(t1);
    FP2_FP256BN_add(&(P->x),&(P->x),&(P->x));	//x.add(x);       //(y^2-9bz^2)xy2

    FP2_FP256BN_norm(&(P->x));			//x.norm();
    //FP2_FP256BN_copy(&(P->y),&y3);		//y.copy(y3);
    FP2_FP256BN_norm(&(P->y));			//y.norm();

    return 1;
}

/* Set P+=Q */
/* SU= 400 */
int ECP2_FP256BN_add(ECP2_FP256BN *P,ECP2_FP256BN *Q)
{
    FP2_FP256BN t0,t1,t2,t3,t4,x3,y3,z3;
    int b3=3*CURVE_B_I_FP256BN;
    if (Q->inf) return 0;
    if (P->inf)
    {
        ECP2_FP256BN_copy(P,Q);
        return 0;
    }

    //FP2_FP256BN_copy(&t0,&(P->x));		//FP2 t0=new FP2(x);
    FP2_FP256BN_mul(&t0,&(P->x),&(Q->x));	//t0.mul(Q.x);         // x.Q.x
    //FP2_FP256BN_copy(&t1,&(P->y));		//FP2 t1=new FP2(y);
    FP2_FP256BN_mul(&t1,&(P->y),&(Q->y));	//t1.mul(Q.y);		 // y.Q.y

    //FP2_FP256BN_copy(&t2,&(P->z));		//FP2 t2=new FP2(z);
    FP2_FP256BN_mul(&t2,&(P->z),&(Q->z));	//t2.mul(Q.z);
    //FP2_FP256BN_copy(&t3,&(P->x));		//FP2 t3=new FP2(x);
    FP2_FP256BN_add(&t3,&(P->x),&(P->y));	//t3.add(y);
    FP2_FP256BN_norm(&t3);				//t3.norm();          //t3=X1+Y1
    //FP2_FP256BN_copy(&t4,&(Q->x));		//FP2 t4=new FP2(Q.x);
    FP2_FP256BN_add(&t4,&(Q->x),&(Q->y));	//t4.add(Q.y);
    FP2_FP256BN_norm(&t4);				//t4.norm();			//t4=X2+Y2
    FP2_FP256BN_mul(&t3,&t3,&t4);		//t3.mul(t4);						//t3=(X1+Y1)(X2+Y2)
    //FP2_FP256BN_copy(&t4,&t0);			//t4.copy(t0);
    FP2_FP256BN_add(&t4,&t0,&t1);		//t4.add(t1);		//t4=X1.X2+Y1.Y2

    FP2_FP256BN_sub(&t3,&t3,&t4);		//t3.sub(t4);
    FP2_FP256BN_norm(&t3);				//t3.norm();
#if SEXTIC_TWIST_FP256BN==D_TYPE
    FP2_FP256BN_mul_ip(&t3);			//t3.mul_ip();
    FP2_FP256BN_norm(&t3);				//t3.norm();         //t3=(X1+Y1)(X2+Y2)-(X1.X2+Y1.Y2) = X1.Y2+X2.Y1
#endif
    //FP2_FP256BN_copy(&t4,&(P->y));		//t4.copy(y);
    FP2_FP256BN_add(&t4,&(P->y),&(P->z));	//t4.add(z);
    FP2_FP256BN_norm(&t4);				//t4.norm();			//t4=Y1+Z1
    //FP2_FP256BN_copy(&x3,&(Q->y));		//FP2 x3=new FP2(Q.y);
    FP2_FP256BN_add(&x3,&(Q->y),&(Q->z));	//x3.add(Q.z);
    FP2_FP256BN_norm(&x3);				//x3.norm();			//x3=Y2+Z2

    FP2_FP256BN_mul(&t4,&t4,&x3);		//t4.mul(x3);						//t4=(Y1+Z1)(Y2+Z2)
    //FP2_FP256BN_copy(&x3,&t1);			//x3.copy(t1);					//
    FP2_FP256BN_add(&x3,&t1,&t2);		//x3.add(t2);						//X3=Y1.Y2+Z1.Z2

    FP2_FP256BN_sub(&t4,&t4,&x3);		//t4.sub(x3);
    FP2_FP256BN_norm(&t4);				//t4.norm();
#if SEXTIC_TWIST_FP256BN==D_TYPE
    FP2_FP256BN_mul_ip(&t4);			//t4.mul_ip();
    FP2_FP256BN_norm(&t4);				//t4.norm();          //t4=(Y1+Z1)(Y2+Z2) - (Y1.Y2+Z1.Z2) = Y1.Z2+Y2.Z1
#endif
    //FP2_FP256BN_copy(&x3,&(P->x));		//x3.copy(x);
    FP2_FP256BN_add(&x3,&(P->x),&(P->z));	//x3.add(z);
    FP2_FP256BN_norm(&x3);				//x3.norm();	// x3=X1+Z1
    //FP2_FP256BN_copy(&y3,&(Q->x));		//FP2 y3=new FP2(Q.x);
    FP2_FP256BN_add(&y3,&(Q->x),&(Q->z));	//y3.add(Q.z);
    FP2_FP256BN_norm(&y3);				//y3.norm();				// y3=X2+Z2
    FP2_FP256BN_mul(&x3,&x3,&y3);		//x3.mul(y3);							// x3=(X1+Z1)(X2+Z2)
    //FP2_FP256BN_copy(&y3,&t0);			//y3.copy(t0);
    FP2_FP256BN_add(&y3,&t0,&t2);		//y3.add(t2);							// y3=X1.X2+Z1+Z2
    FP2_FP256BN_sub(&y3,&x3,&y3);		//y3.rsub(x3);
    FP2_FP256BN_norm(&y3);				//y3.norm();				// y3=(X1+Z1)(X2+Z2) - (X1.X2+Z1.Z2) = X1.Z2+X2.Z1
#if SEXTIC_TWIST_FP256BN==D_TYPE
    FP2_FP256BN_mul_ip(&t0);			//t0.mul_ip();
    FP2_FP256BN_norm(&t0);				//t0.norm(); // x.Q.x
    FP2_FP256BN_mul_ip(&t1);			//t1.mul_ip();
    FP2_FP256BN_norm(&t1);				//t1.norm(); // y.Q.y
#endif
    //FP2_FP256BN_copy(&x3,&t0);			//x3.copy(t0);
    FP2_FP256BN_add(&x3,&t0,&t0);		//x3.add(t0);
    FP2_FP256BN_add(&t0,&t0,&x3);		//t0.add(x3);
    FP2_FP256BN_norm(&t0);				//t0.norm();
    FP2_FP256BN_imul(&t2,&t2,b3);		//t2.imul(b);
#if SEXTIC_TWIST_FP256BN==M_TYPE
    FP2_FP256BN_mul_ip(&t2);
#endif
    //FP2_FP256BN_copy(&z3,&t1);			//FP2 z3=new FP2(t1);
    FP2_FP256BN_add(&z3,&t1,&t2);		//z3.add(t2);
    FP2_FP256BN_norm(&z3);				//z3.norm();
    FP2_FP256BN_sub(&t1,&t1,&t2);		//t1.sub(t2);
    FP2_FP256BN_norm(&t1);				//t1.norm();
    FP2_FP256BN_imul(&y3,&y3,b3);		//y3.imul(b);
#if SEXTIC_TWIST_FP256BN==M_TYPE
    FP2_FP256BN_mul_ip(&y3);
    FP2_FP256BN_norm(&y3);
#endif
    //FP2_FP256BN_copy(&x3,&y3);			//x3.copy(y3);
    FP2_FP256BN_mul(&x3,&y3,&t4);		//x3.mul(t4);
    //FP2_FP256BN_copy(&t2,&t3);			//t2.copy(t3);
    FP2_FP256BN_mul(&t2,&t3,&t1);		//t2.mul(t1);
    FP2_FP256BN_sub(&(P->x),&t2,&x3);		//x3.rsub(t2);
    FP2_FP256BN_mul(&y3,&y3,&t0);		//y3.mul(t0);
    FP2_FP256BN_mul(&t1,&t1,&z3);		//t1.mul(z3);
    FP2_FP256BN_add(&(P->y),&y3,&t1);		//y3.add(t1);
    FP2_FP256BN_mul(&t0,&t0,&t3);		//t0.mul(t3);
    FP2_FP256BN_mul(&z3,&z3,&t4);		//z3.mul(t4);
    FP2_FP256BN_add(&(P->z),&z3,&t0);		//z3.add(t0);

    //FP2_FP256BN_copy(&(P->x),&x3);		//x.copy(x3);
    FP2_FP256BN_norm(&(P->x));			//x.norm();
    //FP2_FP256BN_copy(&(P->y),&y3);		//y.copy(y3);
    FP2_FP256BN_norm(&(P->y));			//y.norm();
    //FP2_FP256BN_copy(&(P->z),&z3);		//z.copy(z3);
    FP2_FP256BN_norm(&(P->z));			//z.norm();

    return 0;
}

/* Set P-=Q */
/* SU= 16 */
void ECP2_FP256BN_sub(ECP2_FP256BN *P,ECP2_FP256BN *Q)
{
    ECP2_FP256BN_neg(Q);
    ECP2_FP256BN_add(P,Q);
    ECP2_FP256BN_neg(Q);
}

/* P*=e */
/* SU= 280 */
void ECP2_FP256BN_mul(ECP2_FP256BN *P,BIG_256_56 e)
{
    /* fixed size windows */
    int i,nb,s,ns;
    BIG_256_56 mt,t;
    ECP2_FP256BN Q,W[8],C;
    sign8 w[1+(NLEN_256_56*BASEBITS_256_56+3)/4];

    if (ECP2_FP256BN_isinf(P)) return;
    ECP2_FP256BN_affine(P);


    /* precompute table */

    ECP2_FP256BN_copy(&Q,P);
    ECP2_FP256BN_dbl(&Q);
    ECP2_FP256BN_copy(&W[0],P);

    for (i=1; i<8; i++)
    {
        ECP2_FP256BN_copy(&W[i],&W[i-1]);
        ECP2_FP256BN_add(&W[i],&Q);
    }

    /* make exponent odd - add 2P if even, P if odd */
    BIG_256_56_copy(t,e);
    s=BIG_256_56_parity(t);
    BIG_256_56_inc(t,1);
    BIG_256_56_norm(t);
    ns=BIG_256_56_parity(t);
    BIG_256_56_copy(mt,t);
    BIG_256_56_inc(mt,1);
    BIG_256_56_norm(mt);
    BIG_256_56_cmove(t,mt,s);
    ECP2_FP256BN_cmove(&Q,P,ns);
    ECP2_FP256BN_copy(&C,&Q);

    nb=1+(BIG_256_56_nbits(t)+3)/4;

    /* convert exponent to signed 4-bit window */
    for (i=0; i<nb; i++)
    {
        w[i]=BIG_256_56_lastbits(t,5)-16;
        BIG_256_56_dec(t,w[i]);
        BIG_256_56_norm(t);
        BIG_256_56_fshr(t,4);
    }
    w[nb]=BIG_256_56_lastbits(t,5);

    ECP2_FP256BN_copy(P,&W[(w[nb]-1)/2]);
    for (i=nb-1; i>=0; i--)
    {
        ECP2_FP256BN_select(&Q,W,w[i]);
        ECP2_FP256BN_dbl(P);
        ECP2_FP256BN_dbl(P);
        ECP2_FP256BN_dbl(P);
        ECP2_FP256BN_dbl(P);
        ECP2_FP256BN_add(P,&Q);
    }
    ECP2_FP256BN_sub(P,&C); /* apply correction */
    ECP2_FP256BN_affine(P);
}

/* Calculates q.P using Frobenius constant X */
/* SU= 96 */
void ECP2_FP256BN_frob(ECP2_FP256BN *P,FP2_FP256BN *X)
{
    FP2_FP256BN X2;
    if (P->inf) return;
//printf("X= "); FP2_FP256BN_output(X); printf("\n");
    FP2_FP256BN_sqr(&X2,X);
//printf("X2= "); FP2_FP256BN_output(&X2); printf("\n");
    FP2_FP256BN_conj(&(P->x),&(P->x));
    FP2_FP256BN_conj(&(P->y),&(P->y));
    FP2_FP256BN_conj(&(P->z),&(P->z));
    FP2_FP256BN_reduce(&(P->z));

    FP2_FP256BN_mul(&(P->x),&X2,&(P->x));
    FP2_FP256BN_mul(&(P->y),&X2,&(P->y));
    FP2_FP256BN_mul(&(P->y),X,&(P->y));


//printf("Px= "); FP2_FP256BN_output(&(P->x)); printf("\n");
//printf("Py= "); FP2_FP256BN_output(&(P->y)); printf("\n");
//printf("Pz= "); FP2_FP256BN_output(&(P->z)); printf("\n");
}

void ECP2_FP256BN_mul4(ECP2_FP256BN *P,ECP2_FP256BN Q[4],BIG_256_56 u[4])
{
    int i,j,a[4],nb;
    ECP2_FP256BN W[8],T,C;
    BIG_256_56 mt,t[4];
    sign8 w[NLEN_256_56*BASEBITS_256_56+1];

    for (i=0; i<4; i++)
    {
        BIG_256_56_copy(t[i],u[i]);
        ECP2_FP256BN_affine(&Q[i]);
    }

    /* precompute table */

    ECP2_FP256BN_copy(&W[0],&Q[0]);
    ECP2_FP256BN_sub(&W[0],&Q[1]);  /* P-Q */
    ECP2_FP256BN_copy(&W[1],&W[0]);
    ECP2_FP256BN_copy(&W[2],&W[0]);
    ECP2_FP256BN_copy(&W[3],&W[0]);
    ECP2_FP256BN_copy(&W[4],&Q[0]);
    ECP2_FP256BN_add(&W[4],&Q[1]);  /* P+Q */
    ECP2_FP256BN_copy(&W[5],&W[4]);
    ECP2_FP256BN_copy(&W[6],&W[4]);
    ECP2_FP256BN_copy(&W[7],&W[4]);

    ECP2_FP256BN_copy(&T,&Q[2]);
    ECP2_FP256BN_sub(&T,&Q[3]);       /* R-S */
    ECP2_FP256BN_sub(&W[1],&T);
    ECP2_FP256BN_add(&W[2],&T);
    ECP2_FP256BN_sub(&W[5],&T);
    ECP2_FP256BN_add(&W[6],&T);
    ECP2_FP256BN_copy(&T,&Q[2]);
    ECP2_FP256BN_add(&T,&Q[3]);      /* R+S */
    ECP2_FP256BN_sub(&W[0],&T);
    ECP2_FP256BN_add(&W[3],&T);
    ECP2_FP256BN_sub(&W[4],&T);
    ECP2_FP256BN_add(&W[7],&T);

    /* if multiplier is even add 1 to multiplier, and add P to correction */
    ECP2_FP256BN_inf(&C);

    BIG_256_56_zero(mt);
    for (i=0; i<4; i++)
    {
        if (BIG_256_56_parity(t[i])==0)
        {
            BIG_256_56_inc(t[i],1);
            BIG_256_56_norm(t[i]);
            ECP2_FP256BN_add(&C,&Q[i]);
        }
        BIG_256_56_add(mt,mt,t[i]);
        BIG_256_56_norm(mt);
    }

    nb=1+BIG_256_56_nbits(mt);

    /* convert exponent to signed 1-bit window */
    for (j=0; j<nb; j++)
    {
        for (i=0; i<4; i++)
        {
            a[i]=BIG_256_56_lastbits(t[i],2)-2;
            BIG_256_56_dec(t[i],a[i]);
            BIG_256_56_norm(t[i]);
            BIG_256_56_fshr(t[i],1);
        }
        w[j]=8*a[0]+4*a[1]+2*a[2]+a[3];
    }
    w[nb]=8*BIG_256_56_lastbits(t[0],2)+4*BIG_256_56_lastbits(t[1],2)+2*BIG_256_56_lastbits(t[2],2)+BIG_256_56_lastbits(t[3],2);

    ECP2_FP256BN_copy(P,&W[(w[nb]-1)/2]);
    for (i=nb-1; i>=0; i--)
    {
        ECP2_FP256BN_select(&T,W,w[i]);
        ECP2_FP256BN_dbl(P);
        ECP2_FP256BN_add(P,&T);
    }
    ECP2_FP256BN_sub(P,&C); /* apply correction */

    ECP2_FP256BN_affine(P);
}


/* Map to hash value to point on G2 from random BIG */
void ECP2_FP256BN_mapit(ECP2_FP256BN *Q,octet *W)
{
    BIG_256_56 q,one,Fx,Fy,x,hv;
    FP2_FP256BN X;
#if (PAIRING_FRIENDLY_FP256BN == BN)
    ECP2_FP256BN T,K;
#elif (PAIRING_FRIENDLY_FP256BN == BLS)
    ECP2_FP256BN xQ, x2Q;
#endif
    BIG_256_56_fromBytes(hv,W->val);
    BIG_256_56_rcopy(q,Modulus_FP256BN);
    BIG_256_56_one(one);
    BIG_256_56_mod(hv,q);

    for (;;)
    {
        FP2_FP256BN_from_BIGs(&X,one,hv);
        if (ECP2_FP256BN_setx(Q,&X)) break;
        BIG_256_56_inc(hv,1);
    }

    BIG_256_56_rcopy(Fx,Fra_FP256BN);
    BIG_256_56_rcopy(Fy,Frb_FP256BN);
    FP2_FP256BN_from_BIGs(&X,Fx,Fy);

#if SEXTIC_TWIST_FP256BN==M_TYPE
    FP2_FP256BN_inv(&X,&X);
    FP2_FP256BN_norm(&X);
#endif

    BIG_256_56_rcopy(x,CURVE_Bnx_FP256BN);

#if (PAIRING_FRIENDLY_FP256BN == BN)

    /* Faster Hashing to G2 - Fuentes-Castaneda, Knapp and Rodriguez-Henriquez */
    /* Q -> xQ + F(3xQ) + F(F(xQ)) + F(F(F(Q))). */
    ECP2_FP256BN_copy(&T,Q);
    ECP2_FP256BN_mul(&T,x);
#if SIGN_OF_X_FP256BN==NEGATIVEX
    ECP2_FP256BN_neg(&T);   // our x is negative
#endif
    ECP2_FP256BN_copy(&K,&T);
    ECP2_FP256BN_dbl(&K);
    ECP2_FP256BN_add(&K,&T);
    //ECP2_FP256BN_affine(&K);

    ECP2_FP256BN_frob(&K,&X);
    ECP2_FP256BN_frob(Q,&X);
    ECP2_FP256BN_frob(Q,&X);
    ECP2_FP256BN_frob(Q,&X);
    ECP2_FP256BN_add(Q,&T);
    ECP2_FP256BN_add(Q,&K);
    ECP2_FP256BN_frob(&T,&X);
    ECP2_FP256BN_frob(&T,&X);
    ECP2_FP256BN_add(Q,&T);
    ECP2_FP256BN_affine(Q);

#elif (PAIRING_FRIENDLY_FP256BN == BLS)

    /* Efficient hash maps to G2 on BLS curves - Budroni, Pintore */
    /* Q -> x2Q -xQ -Q +F(xQ -Q) +F(F(2Q)) */

    ECP2_FP256BN_copy(&xQ,Q);
    ECP2_FP256BN_mul(&xQ,x);

    ECP2_FP256BN_copy(&x2Q,&xQ);
    ECP2_FP256BN_mul(&x2Q,x);

#if SIGN_OF_X_FP256BN==NEGATIVEX
    ECP2_FP256BN_neg(&xQ);
#endif

    ECP2_FP256BN_sub(&x2Q,&xQ);
    ECP2_FP256BN_sub(&x2Q,Q);

    ECP2_FP256BN_sub(&xQ,Q);
    ECP2_FP256BN_frob(&xQ,&X);

    ECP2_FP256BN_dbl(Q);
    ECP2_FP256BN_frob(Q,&X);
    ECP2_FP256BN_frob(Q,&X);

    ECP2_FP256BN_add(Q,&x2Q);
    ECP2_FP256BN_add(Q,&xQ);

    ECP2_FP256BN_affine(Q);

#endif
}

/*

int main()
{
	int i;
	ECP2_FP256BN G,P;
	ECP2_FP256BN *W;
	FP2_FP256BN x,y,w,z,f;
	BIG_256_56 r,xa,xb,ya,yb;

	BIG_256_56_rcopy(xa,CURVE_Pxa_FP256BN);
	BIG_256_56_rcopy(xb,CURVE_Pxb_FP256BN);
	BIG_256_56_rcopy(ya,CURVE_Pya_FP256BN);
	BIG_256_56_rcopy(yb,CURVE_Pyb_FP256BN);

	FP2_FP256BN_from_BIGs(&x,xa,xb);
	FP2_FP256BN_from_BIGs(&y,ya,yb);
	ECP2_FP256BN_set(&G,&x,&y);
	if (G.inf) printf("Failed to set - point not on curve\n");
	else printf("set success\n");

	ECP2_FP256BN_output(&G);

//	BIG_256_56_copy(r,CURVE_Order_FP256BN);
	BIG_256_56_rcopy(r,Modulus_FP256BN);

	ECP2_FP256BN_copy(&P,&G);

	ECP2_FP256BN_mul(&P,r);

	ECP2_FP256BN_output(&P);

	FP2_FP256BN_gfc(&f,12);

	ECP2_FP256BN_frob(&G,&f);

	ECP2_FP256BN_output(&G);

	return 0;
}

*/