version3/cpp/fp24.cpp

/*
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 Fp^12 functions */
/* SU=m, m is Stack Usage (no lazy )*/
/* FP24 elements are of the form a+i.b+i^2.c */

#include "fp24_YYY.h"

using namespace XXX;

/* 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 FP24_select(YYY::FP24 *f,YYY::FP24 g[],sign32 b)
{
    YYY::FP24 invf;
    sign32 m=b>>31;
    sign32 babs=(b^m)-m;

    babs=(babs-1)/2;

    FP24_cmove(f,&g[0],teq(babs,0));  // conditional move
    FP24_cmove(f,&g[1],teq(babs,1));
    FP24_cmove(f,&g[2],teq(babs,2));
    FP24_cmove(f,&g[3],teq(babs,3));
    FP24_cmove(f,&g[4],teq(babs,4));
    FP24_cmove(f,&g[5],teq(babs,5));
    FP24_cmove(f,&g[6],teq(babs,6));
    FP24_cmove(f,&g[7],teq(babs,7));

    FP24_copy(&invf,f);
    FP24_conj(&invf,&invf);  // 1/f
    FP24_cmove(f,&invf,(int)(m&1));
}

/* test x==0 ? */
/* SU= 8 */
int YYY::FP24_iszilch(FP24 *x)
{
    if (FP8_iszilch(&(x->a)) && FP8_iszilch(&(x->b)) && FP8_iszilch(&(x->c))) return 1;
    return 0;
}

/* test x==1 ? */
/* SU= 8 */
int YYY::FP24_isunity(FP24 *x)
{
    if (FP8_isunity(&(x->a)) && FP8_iszilch(&(x->b)) && FP8_iszilch(&(x->c))) return 1;
    return 0;
}

/* FP24 copy w=x */
/* SU= 16 */
void YYY::FP24_copy(FP24 *w,FP24 *x)
{
    if (x==w) return;
    FP8_copy(&(w->a),&(x->a));
    FP8_copy(&(w->b),&(x->b));
    FP8_copy(&(w->c),&(x->c));
}

/* FP24 w=1 */
/* SU= 8 */
void YYY::FP24_one(FP24 *w)
{
    FP8_one(&(w->a));
    FP8_zero(&(w->b));
    FP8_zero(&(w->c));
}

/* return 1 if x==y, else 0 */
/* SU= 16 */
int YYY::FP24_equals(FP24 *x,FP24 *y)
{
    if (FP8_equals(&(x->a),&(y->a)) && FP8_equals(&(x->b),&(y->b)) && FP8_equals(&(x->b),&(y->b)))
        return 1;
    return 0;
}

/* Set w=conj(x) */
/* SU= 8 */
void YYY::FP24_conj(FP24 *w,FP24 *x)
{
    FP24_copy(w,x);
    FP8_conj(&(w->a),&(w->a));
    FP8_nconj(&(w->b),&(w->b));
    FP8_conj(&(w->c),&(w->c));
}

/* Create FP24 from FP8 */
/* SU= 8 */
void YYY::FP24_from_FP8(FP24 *w,FP8 *a)
{
    FP8_copy(&(w->a),a);
    FP8_zero(&(w->b));
    FP8_zero(&(w->c));
}

/* Create FP24 from 3 FP8's */
/* SU= 16 */
void YYY::FP24_from_FP8s(FP24 *w,FP8 *a,FP8 *b,FP8 *c)
{
    FP8_copy(&(w->a),a);
    FP8_copy(&(w->b),b);
    FP8_copy(&(w->c),c);
}

/* Granger-Scott Unitary Squaring. This does not benefit from lazy reduction */
/* SU= 600 */
void YYY::FP24_usqr(FP24 *w,FP24 *x)
{
    FP8 A,B,C,D;

    FP8_copy(&A,&(x->a));

    FP8_sqr(&(w->a),&(x->a));
    FP8_add(&D,&(w->a),&(w->a));
    FP8_add(&(w->a),&D,&(w->a));

    FP8_norm(&(w->a));
    FP8_nconj(&A,&A);

    FP8_add(&A,&A,&A);
    FP8_add(&(w->a),&(w->a),&A);
    FP8_sqr(&B,&(x->c));
    FP8_times_i(&B);

    FP8_add(&D,&B,&B);
    FP8_add(&B,&B,&D);
    FP8_norm(&B);

    FP8_sqr(&C,&(x->b));

    FP8_add(&D,&C,&C);
    FP8_add(&C,&C,&D);

    FP8_norm(&C);
    FP8_conj(&(w->b),&(x->b));
    FP8_add(&(w->b),&(w->b),&(w->b));
    FP8_nconj(&(w->c),&(x->c));

    FP8_add(&(w->c),&(w->c),&(w->c));
    FP8_add(&(w->b),&B,&(w->b));
    FP8_add(&(w->c),&C,&(w->c));

//	FP24_norm(w);
    FP24_reduce(w);	    /* reduce here as in pow function repeated squarings would trigger multiple reductions */
}

/* FP24 squaring w=x^2 */
/* SU= 600 */
void YYY::FP24_sqr(FP24 *w,FP24 *x)
{
    /* Use Chung-Hasan SQR2 method from http://cacr.uwaterloo.ca/techreports/2006/cacr2006-24.pdf */

    FP8 A,B,C,D;

    FP8_sqr(&A,&(x->a));
    FP8_mul(&B,&(x->b),&(x->c));
    FP8_add(&B,&B,&B);
FP8_norm(&B);
    FP8_sqr(&C,&(x->c));

    FP8_mul(&D,&(x->a),&(x->b));
    FP8_add(&D,&D,&D);

    FP8_add(&(w->c),&(x->a),&(x->c));
    FP8_add(&(w->c),&(x->b),&(w->c));
FP8_norm(&(w->c));	

    FP8_sqr(&(w->c),&(w->c));

    FP8_copy(&(w->a),&A);
    FP8_add(&A,&A,&B);

    FP8_norm(&A);

    FP8_add(&A,&A,&C);
    FP8_add(&A,&A,&D);

    FP8_norm(&A);

    FP8_neg(&A,&A);
    FP8_times_i(&B);
    FP8_times_i(&C);

    FP8_add(&(w->a),&(w->a),&B);
    FP8_add(&(w->b),&C,&D);
    FP8_add(&(w->c),&(w->c),&A);

    FP24_norm(w);
}

/* FP24 full multiplication w=w*y */


/* SU= 896 */
/* FP24 full multiplication w=w*y */
void YYY::FP24_mul(FP24 *w,FP24 *y)
{
    FP8 z0,z1,z2,z3,t0,t1;

    FP8_mul(&z0,&(w->a),&(y->a));
    FP8_mul(&z2,&(w->b),&(y->b));  //

    FP8_add(&t0,&(w->a),&(w->b));
    FP8_add(&t1,&(y->a),&(y->b));  //

FP8_norm(&t0);
FP8_norm(&t1);

    FP8_mul(&z1,&t0,&t1);
    FP8_add(&t0,&(w->b),&(w->c));
    FP8_add(&t1,&(y->b),&(y->c));  //

FP8_norm(&t0);
FP8_norm(&t1);

    FP8_mul(&z3,&t0,&t1);

    FP8_neg(&t0,&z0);
    FP8_neg(&t1,&z2);

    FP8_add(&z1,&z1,&t0);   // z1=z1-z0
//    FP8_norm(&z1);
    FP8_add(&(w->b),&z1,&t1);
// z1=z1-z2
    FP8_add(&z3,&z3,&t1);        // z3=z3-z2
    FP8_add(&z2,&z2,&t0);        // z2=z2-z0

    FP8_add(&t0,&(w->a),&(w->c));
    FP8_add(&t1,&(y->a),&(y->c));

FP8_norm(&t0);
FP8_norm(&t1);

    FP8_mul(&t0,&t1,&t0);
    FP8_add(&z2,&z2,&t0);

    FP8_mul(&t0,&(w->c),&(y->c));
    FP8_neg(&t1,&t0);

    FP8_add(&(w->c),&z2,&t1);
    FP8_add(&z3,&z3,&t1);
    FP8_times_i(&t0);
    FP8_add(&(w->b),&(w->b),&t0);
FP8_norm(&z3);
    FP8_times_i(&z3);
    FP8_add(&(w->a),&z0,&z3);

    FP24_norm(w);
}

/* FP24 multiplication w=w*y */
/* SU= 744 */
/* catering for special case that arises from special form of ATE pairing line function */
void YYY::FP24_smul(FP24 *w,FP24 *y,int type)
{
    FP8 z0,z1,z2,z3,t0,t1;

	if (type==D_TYPE)
	{ // y->c is 0

		FP8_copy(&z3,&(w->b));
		FP8_mul(&z0,&(w->a),&(y->a));

		FP8_pmul(&z2,&(w->b),&(y->b).a);
		FP8_add(&(w->b),&(w->a),&(w->b));
		FP8_copy(&t1,&(y->a));
		FP4_add(&t1.a,&t1.a,&(y->b).a);

		FP8_norm(&t1);
		FP8_norm(&(w->b));

		FP8_mul(&(w->b),&(w->b),&t1);
		FP8_add(&z3,&z3,&(w->c));
		FP8_norm(&z3);
		FP8_pmul(&z3,&z3,&(y->b).a);
		FP8_neg(&t0,&z0);
		FP8_neg(&t1,&z2);

		FP8_add(&(w->b),&(w->b),&t0);   // z1=z1-z0
//    FP8_norm(&(w->b));
		FP8_add(&(w->b),&(w->b),&t1);   // z1=z1-z2

		FP8_add(&z3,&z3,&t1);        // z3=z3-z2
		FP8_add(&z2,&z2,&t0);        // z2=z2-z0

		FP8_add(&t0,&(w->a),&(w->c));

		FP8_norm(&t0);
		FP8_norm(&z3);

		FP8_mul(&t0,&(y->a),&t0);
		FP8_add(&(w->c),&z2,&t0);

		FP8_times_i(&z3);
		FP8_add(&(w->a),&z0,&z3);
	}

	if (type==M_TYPE)
	{ // y->b is zero
		FP8_mul(&z0,&(w->a),&(y->a));
		FP8_add(&t0,&(w->a),&(w->b));
		FP8_norm(&t0);

		FP8_mul(&z1,&t0,&(y->a));
		FP8_add(&t0,&(w->b),&(w->c));
		FP8_norm(&t0);

		FP8_pmul(&z3,&t0,&(y->c).b);
		FP8_times_i(&z3);

		FP8_neg(&t0,&z0);
		FP8_add(&z1,&z1,&t0);   // z1=z1-z0

		FP8_copy(&(w->b),&z1);

		FP8_copy(&z2,&t0);

		FP8_add(&t0,&(w->a),&(w->c));
		FP8_add(&t1,&(y->a),&(y->c));

		FP8_norm(&t0);
		FP8_norm(&t1);

		FP8_mul(&t0,&t1,&t0);
		FP8_add(&z2,&z2,&t0);

		FP8_pmul(&t0,&(w->c),&(y->c).b);
		FP8_times_i(&t0);
		FP8_neg(&t1,&t0);
		FP8_times_i(&t0);

		FP8_add(&(w->c),&z2,&t1);
		FP8_add(&z3,&z3,&t1);

		FP8_add(&(w->b),&(w->b),&t0);
		FP8_norm(&z3);
		FP8_times_i(&z3);
		FP8_add(&(w->a),&z0,&z3);
	}
    FP24_norm(w);
}

/* Set w=1/x */
/* SU= 600 */
void YYY::FP24_inv(FP24 *w,FP24 *x)
{
    FP8 f0,f1,f2,f3;
//    FP24_norm(x);

    FP8_sqr(&f0,&(x->a));
    FP8_mul(&f1,&(x->b),&(x->c));
    FP8_times_i(&f1);
    FP8_sub(&f0,&f0,&f1);  /* y.a */
	FP8_norm(&f0); 		

    FP8_sqr(&f1,&(x->c));
    FP8_times_i(&f1);
    FP8_mul(&f2,&(x->a),&(x->b));
    FP8_sub(&f1,&f1,&f2);  /* y.b */
	FP8_norm(&f1); 

    FP8_sqr(&f2,&(x->b));
    FP8_mul(&f3,&(x->a),&(x->c));
    FP8_sub(&f2,&f2,&f3);  /* y.c */
	FP8_norm(&f2); 

    FP8_mul(&f3,&(x->b),&f2);
    FP8_times_i(&f3);
    FP8_mul(&(w->a),&f0,&(x->a));
    FP8_add(&f3,&(w->a),&f3);
    FP8_mul(&(w->c),&f1,&(x->c));
    FP8_times_i(&(w->c));



    FP8_add(&f3,&(w->c),&f3);
	FP8_norm(&f3);
	
    FP8_inv(&f3,&f3);
    FP8_mul(&(w->a),&f0,&f3);
    FP8_mul(&(w->b),&f1,&f3);
    FP8_mul(&(w->c),&f2,&f3);

}

/* constant time powering by small integer of max length bts */

void YYY::FP24_pinpow(FP24 *r,int e,int bts)
{
    int i,b;
    FP24 R[2];

    FP24_one(&R[0]);
    FP24_copy(&R[1],r);

    for (i=bts-1; i>=0; i--)
    {
        b=(e>>i)&1;
        FP24_mul(&R[1-b],&R[b]);
        FP24_usqr(&R[b],&R[b]);
    }
    FP24_copy(r,&R[0]);
}

/* Compressed powering of unitary elements y=x^(e mod r) */

void YYY::FP24_compow(FP8 *c,FP24 *x,BIG e,BIG r)
{
    FP24 g1,g2;
	FP8 cp,cpm1,cpm2;
    FP2 f;
	BIG q,a,b,m;

    BIG_rcopy(a,Fra);
    BIG_rcopy(b,Frb);
    FP2_from_BIGs(&f,a,b);

    BIG_rcopy(q,Modulus);

    FP24_copy(&g1,x);
	FP24_copy(&g2,x);

    BIG_copy(m,q);
    BIG_mod(m,r);

    BIG_copy(a,e);
    BIG_mod(a,m);

    BIG_copy(b,e);
    BIG_sdiv(b,m);

    FP24_trace(c,&g1);

	if (BIG_iszilch(b))
	{
		FP8_xtr_pow(c,c,e);
		return;
	}

    FP24_frob(&g2,&f,1);
    FP24_trace(&cp,&g2);
    FP24_conj(&g1,&g1);
    FP24_mul(&g2,&g1);
    FP24_trace(&cpm1,&g2);
    FP24_mul(&g2,&g1);

    FP24_trace(&cpm2,&g2);

    FP8_xtr_pow2(c,&cp,c,&cpm1,&cpm2,a,b);

}

/* Note this is simple square and multiply, so not side-channel safe */

void YYY::FP24_pow(FP24 *r,FP24 *a,BIG b)
{
    FP24 w,sf;
    BIG b1,b3;
    int i,nb,bt;
	BIG_copy(b1,b);
    BIG_norm(b1);
	BIG_pmul(b3,b1,3);
	BIG_norm(b3);

	FP24_copy(&sf,a);
	FP24_norm(&sf);
    FP24_copy(&w,&sf);


	nb=BIG_nbits(b3);
	for (i=nb-2;i>=1;i--)
	{
		FP24_usqr(&w,&w);
		bt=BIG_bit(b3,i)-BIG_bit(b1,i);
		if (bt==1)
			FP24_mul(&w,&sf);
		if (bt==-1)
		{
			FP24_conj(&sf,&sf);
			FP24_mul(&w,&sf);
			FP24_conj(&sf,&sf);
		}
	}

	FP24_copy(r,&w);
	FP24_reduce(r);
}


/* SU= 528 */
/* set r=a^b */
/* Note this is simple square and multiply, so not side-channel safe 

void YYY::FP24_ppow(FP24 *r,FP24 *a,BIG b)
{
    FP24 w;
    BIG z,zilch;
    int bt;
    BIG_zero(zilch);
    BIG_norm(b);
    BIG_copy(z,b);
    FP24_copy(&w,a);
    FP24_one(r);

    while(1)
    {
        bt=BIG_parity(z);
        BIG_shr(z,1);
        if (bt)
		{
			//printf("In mul\n");
            FP24_mul(r,&w);
			//printf("Out of mul\n");
		}
        if (BIG_comp(z,zilch)==0) break;
		//printf("In sqr\n");
        FP24_sqr(&w,&w);
		//printf("Out of sqr\n");
    }

    FP24_reduce(r);
}  */


/* p=q0^u0.q1^u1.q2^u2.q3^u3... */
/* Side channel attack secure */
// Bos & Costello https://eprint.iacr.org/2013/458.pdf
// Faz-Hernandez & Longa & Sanchez  https://eprint.iacr.org/2013/158.pdf

void YYY::FP24_pow8(FP24 *p,FP24 *q,BIG u[8])
{
    int i,j,k,nb,pb1,pb2,bt;
	FP24 g1[8],g2[8],r;
	BIG t[8],mt;
    sign8 w1[NLEN_XXX*BASEBITS_XXX+1];
    sign8 s1[NLEN_XXX*BASEBITS_XXX+1];
    sign8 w2[NLEN_XXX*BASEBITS_XXX+1];
    sign8 s2[NLEN_XXX*BASEBITS_XXX+1];
    FP fx,fy;
	FP2 X;

    FP_rcopy(&fx,Fra);
    FP_rcopy(&fy,Frb);
    FP2_from_FPs(&X,&fx,&fy);

    for (i=0; i<8; i++)
        BIG_copy(t[i],u[i]);

// Precomputed table
    FP24_copy(&g1[0],&q[0]); // q[0]
    FP24_copy(&g1[1],&g1[0]);
	FP24_mul(&g1[1],&q[1]);	// q[0].q[1]
    FP24_copy(&g1[2],&g1[0]);
	FP24_mul(&g1[2],&q[2]);	// q[0].q[2]
	FP24_copy(&g1[3],&g1[1]);
	FP24_mul(&g1[3],&q[2]);	// q[0].q[1].q[2]
	FP24_copy(&g1[4],&g1[0]);
	FP24_mul(&g1[4],&q[3]);  // q[0].q[3]
	FP24_copy(&g1[5],&g1[1]);
	FP24_mul(&g1[5],&q[3]);	// q[0].q[1].q[3]
	FP24_copy(&g1[6],&g1[2]);
	FP24_mul(&g1[6],&q[3]);	// q[0].q[2].q[3]
	FP24_copy(&g1[7],&g1[3]);
	FP24_mul(&g1[7],&q[3]);	// q[0].q[1].q[2].q[3]

// Use Frobenius

	for (i=0;i<8;i++)
	{
		FP24_copy(&g2[i],&g1[i]);
		FP24_frob(&g2[i],&X,4);
	}

// Make it odd
	pb1=1-BIG_parity(t[0]);
	BIG_inc(t[0],pb1);
	BIG_norm(t[0]);

	pb2=1-BIG_parity(t[4]);
	BIG_inc(t[4],pb2);
	BIG_norm(t[4]);

// Number of bits
    BIG_zero(mt);
    for (i=0; i<8; i++)
    {
        BIG_or(mt,mt,t[i]);
    }
    nb=1+BIG_nbits(mt);

// Sign pivot 
	s1[nb-1]=1;
	s2[nb-1]=1;
	for (i=0;i<nb-1;i++)
	{
        BIG_fshr(t[0],1);
		s1[i]=2*BIG_parity(t[0])-1;
        BIG_fshr(t[4],1);
		s2[i]=2*BIG_parity(t[4])-1;
	}

// Recoded exponents
    for (i=0; i<nb; i++)
    {
		w1[i]=0;
		k=1;
		for (j=1; j<4; j++)
		{
			bt=s1[i]*BIG_parity(t[j]);
			BIG_fshr(t[j],1);

			BIG_dec(t[j],(bt>>1));
			BIG_norm(t[j]);
			w1[i]+=bt*k;
			k*=2;
        }

		w2[i]=0;
		k=1;
		for (j=5; j<8; j++)
		{
			bt=s2[i]*BIG_parity(t[j]);
			BIG_fshr(t[j],1);

			BIG_dec(t[j],(bt>>1));
			BIG_norm(t[j]);
			w2[i]+=bt*k;
			k*=2;
        }
    }	

// Main loop
	FP24_select(p,g1,2*w1[nb-1]+1);
	FP24_select(&r,g2,2*w2[nb-1]+1);
	FP24_mul(p,&r);
    for (i=nb-2; i>=0; i--)
    {
		FP24_usqr(p,p);
        FP24_select(&r,g1,2*w1[i]+s1[i]);
        FP24_mul(p,&r);
        FP24_select(&r,g2,2*w2[i]+s2[i]);
        FP24_mul(p,&r);
    }

// apply correction
	FP24_conj(&r,&q[0]);   
	FP24_mul(&r,p);
	FP24_cmove(p,&r,pb1);
	FP24_conj(&r,&q[4]);   
	FP24_mul(&r,p);
	FP24_cmove(p,&r,pb2);

	FP24_reduce(p);
}

/*
void YYY::FP24_pow8(FP24 *p,FP24 *q,BIG u[8])
{
    int i,j,a[4],nb,m;
    FP24 g[8],f[8],c,s[2];
    BIG t[8],mt;
    sign8 w[NLEN_XXX*BASEBITS_XXX+1];
    sign8 z[NLEN_XXX*BASEBITS_XXX+1];
    FP fx,fy;
	FP2 X;

    FP_rcopy(&fx,Fra);
    FP_rcopy(&fy,Frb);
    FP2_from_FPs(&X,&fx,&fy);

    for (i=0; i<8; i++)
        BIG_copy(t[i],u[i]);

    FP24_copy(&g[0],&q[0]);
    FP24_conj(&s[0],&q[1]);
    FP24_mul(&g[0],&s[0]);  // P/Q 
    FP24_copy(&g[1],&g[0]);
    FP24_copy(&g[2],&g[0]);
    FP24_copy(&g[3],&g[0]);
    FP24_copy(&g[4],&q[0]);
    FP24_mul(&g[4],&q[1]);  // P*Q 
    FP24_copy(&g[5],&g[4]);
    FP24_copy(&g[6],&g[4]);
    FP24_copy(&g[7],&g[4]);

    FP24_copy(&s[1],&q[2]);
    FP24_conj(&s[0],&q[3]);
    FP24_mul(&s[1],&s[0]);       // R/S 
    FP24_conj(&s[0],&s[1]);
    FP24_mul(&g[1],&s[0]);
    FP24_mul(&g[2],&s[1]);
    FP24_mul(&g[5],&s[0]);
    FP24_mul(&g[6],&s[1]);
    FP24_copy(&s[1],&q[2]);
    FP24_mul(&s[1],&q[3]);      // R*S 
    FP24_conj(&s[0],&s[1]);
    FP24_mul(&g[0],&s[0]);
    FP24_mul(&g[3],&s[1]);
    FP24_mul(&g[4],&s[0]);
    FP24_mul(&g[7],&s[1]);

// Use Frobenius

	for (i=0;i<8;i++)
	{
		FP24_copy(&f[i],&g[i]);
		FP24_frob(&f[i],&X,4);
	}


    // if power is even add 1 to power, and add q to correction 
    FP24_one(&c);

    BIG_zero(mt);
    for (i=0; i<8; i++)
    {
        if (BIG_parity(t[i])==0)
        {
            BIG_inc(t[i],1);
            BIG_norm(t[i]);
            FP24_mul(&c,&q[i]);
        }
        BIG_add(mt,mt,t[i]);
        BIG_norm(mt);
    }

    FP24_conj(&c,&c);
    nb=1+BIG_nbits(mt);

    // convert exponents to signed 1-bit windows 
    for (j=0; j<nb; j++)
    {
        for (i=0; i<4; i++)
        {
            a[i]=BIG_lastbits(t[i],2)-2;
            BIG_dec(t[i],a[i]);
            BIG_norm(t[i]);
            BIG_fshr(t[i],1);
        }
        w[j]=8*a[0]+4*a[1]+2*a[2]+a[3];
    }
    w[nb]=8*BIG_lastbits(t[0],2)+4*BIG_lastbits(t[1],2)+2*BIG_lastbits(t[2],2)+BIG_lastbits(t[3],2);


    for (j=0; j<nb; j++)
    {
        for (i=0; i<4; i++)
        {
            a[i]=BIG_lastbits(t[i+4],2)-2;
            BIG_dec(t[i+4],a[i]);
            BIG_norm(t[i+4]);
            BIG_fshr(t[i+4],1);
        }
        z[j]=8*a[0]+4*a[1]+2*a[2]+a[3];
    }
    z[nb]=8*BIG_lastbits(t[4],2)+4*BIG_lastbits(t[5],2)+2*BIG_lastbits(t[6],2)+BIG_lastbits(t[7],2);


    FP24_copy(p,&g[(w[nb]-1)/2]);
    FP24_mul(p,&f[(z[nb]-1)/2]);
    for (i=nb-1; i>=0; i--)
    {
		FP24_usqr(p,p);

        m=w[i]>>7;
        j=(w[i]^m)-m;  // j=abs(w[i]) 
        j=(j-1)/2;
        FP24_copy(&s[0],&g[j]);
        FP24_conj(&s[1],&g[j]);
        FP24_mul(p,&s[m&1]);

        m=z[i]>>7;
        j=(z[i]^m)-m;  // j=abs(w[i]) 
        j=(j-1)/2;
        FP24_copy(&s[0],&f[j]);
        FP24_conj(&s[1],&f[j]);
        FP24_mul(p,&s[m&1]);

    }
    FP24_mul(p,&c); // apply correction 
    FP24_reduce(p);
}
*/

/* Set w=w^p using Frobenius */
/* SU= 160 */
void YYY::FP24_frob(FP24 *w,FP2 *f,int n)
{
	int i;
	FP4 X2,X4;
    FP2 f3,f2;				// f=(1+i)^(p-7)/12
    FP2_sqr(&f2,f);     // 
    FP2_mul(&f3,&f2,f); // f3=f^3=(1+i)^(p-7)/4

	FP2_mul_ip(&f3);    // f3 = (1+i).f3 = (1+i)^(p-3)/4 
	FP2_norm(&f3);

	for (i=0;i<n;i++)
	{
		FP8_frob(&(w->a),&f3);   // a=a^p
		FP8_frob(&(w->b),&f3);   // b=b^p
		FP8_frob(&(w->c),&f3);   // c=c^p
  
		FP8_qmul(&(w->b),&(w->b),f); FP8_times_i2(&(w->b));
		FP8_qmul(&(w->c),&(w->c),&f2); FP8_times_i2(&(w->c)); FP8_times_i2(&(w->c));
	}
}


/* SU= 8 */
/* normalise all components of w */
void YYY::FP24_norm(FP24 *w)
{
    FP8_norm(&(w->a));
    FP8_norm(&(w->b));
    FP8_norm(&(w->c));
}

/* SU= 8 */
/* reduce all components of w */
void YYY::FP24_reduce(FP24 *w)
{
    FP8_reduce(&(w->a));
    FP8_reduce(&(w->b));
    FP8_reduce(&(w->c));
}

/* trace function w=trace(x) */
/* SU= 8 */
void YYY::FP24_trace(FP8 *w,FP24 *x)
{
    FP8_imul(w,&(x->a),3);
    FP8_reduce(w);
}

/* SU= 8 */
/* Output w in hex */
void YYY::FP24_output(FP24 *w)
{
    printf("[");
    FP8_output(&(w->a));
    printf(",");
    FP8_output(&(w->b));
    printf(",");
    FP8_output(&(w->c));
    printf("]");
}

/* SU= 64 */
/* Convert g to octet string w */
void YYY::FP24_toOctet(octet *W,FP24 *g)
{
    BIG a;
    W->len=24*MODBYTES_XXX;

    FP_redc(a,&(g->a.a.a.a));
    BIG_toBytes(&(W->val[0]),a);
    FP_redc(a,&(g->a.a.a.b));
    BIG_toBytes(&(W->val[MODBYTES_XXX]),a);
    FP_redc(a,&(g->a.a.b.a));
    BIG_toBytes(&(W->val[2*MODBYTES_XXX]),a);
    FP_redc(a,&(g->a.a.b.b));
    BIG_toBytes(&(W->val[3*MODBYTES_XXX]),a);
    FP_redc(a,&(g->a.b.a.a));
    BIG_toBytes(&(W->val[4*MODBYTES_XXX]),a);
    FP_redc(a,&(g->a.b.a.b));
    BIG_toBytes(&(W->val[5*MODBYTES_XXX]),a);
    FP_redc(a,&(g->a.b.b.a));
    BIG_toBytes(&(W->val[6*MODBYTES_XXX]),a);
    FP_redc(a,&(g->a.b.b.b));
    BIG_toBytes(&(W->val[7*MODBYTES_XXX]),a);

    FP_redc(a,&(g->b.a.a.a));
    BIG_toBytes(&(W->val[8*MODBYTES_XXX]),a);
    FP_redc(a,&(g->b.a.a.b));
    BIG_toBytes(&(W->val[9*MODBYTES_XXX]),a);
    FP_redc(a,&(g->b.a.b.a));
    BIG_toBytes(&(W->val[10*MODBYTES_XXX]),a);
    FP_redc(a,&(g->b.a.b.b));
    BIG_toBytes(&(W->val[11*MODBYTES_XXX]),a);
    FP_redc(a,&(g->b.b.a.a));
    BIG_toBytes(&(W->val[12*MODBYTES_XXX]),a);
    FP_redc(a,&(g->b.b.a.b));
    BIG_toBytes(&(W->val[13*MODBYTES_XXX]),a);
    FP_redc(a,&(g->b.b.b.a));
    BIG_toBytes(&(W->val[14*MODBYTES_XXX]),a);
    FP_redc(a,&(g->b.b.b.b));
    BIG_toBytes(&(W->val[15*MODBYTES_XXX]),a);

    FP_redc(a,&(g->c.a.a.a));
    BIG_toBytes(&(W->val[16*MODBYTES_XXX]),a);
    FP_redc(a,&(g->c.a.a.b));
    BIG_toBytes(&(W->val[17*MODBYTES_XXX]),a);
    FP_redc(a,&(g->c.a.b.a));
    BIG_toBytes(&(W->val[18*MODBYTES_XXX]),a);
    FP_redc(a,&(g->c.a.b.b));
    BIG_toBytes(&(W->val[19*MODBYTES_XXX]),a);
    FP_redc(a,&(g->c.b.a.a));
    BIG_toBytes(&(W->val[20*MODBYTES_XXX]),a);
    FP_redc(a,&(g->c.b.a.b));
    BIG_toBytes(&(W->val[21*MODBYTES_XXX]),a);
    FP_redc(a,&(g->c.b.b.a));
    BIG_toBytes(&(W->val[22*MODBYTES_XXX]),a);
    FP_redc(a,&(g->c.b.b.b));
    BIG_toBytes(&(W->val[23*MODBYTES_XXX]),a);
}

/* SU= 24 */
/* Restore g from octet string w */
void YYY::FP24_fromOctet(FP24 *g,octet *W)
{
	BIG b;

    BIG_fromBytes(b,&W->val[0]);
    FP_nres(&(g->a.a.a.a),b);
    BIG_fromBytes(b,&W->val[MODBYTES_XXX]);
    FP_nres(&(g->a.a.a.b),b);
    BIG_fromBytes(b,&W->val[2*MODBYTES_XXX]);
    FP_nres(&(g->a.a.b.a),b);
    BIG_fromBytes(b,&W->val[3*MODBYTES_XXX]);
    FP_nres(&(g->a.a.b.b),b);
    BIG_fromBytes(b,&W->val[4*MODBYTES_XXX]);
    FP_nres(&(g->a.b.a.a),b);
    BIG_fromBytes(b,&W->val[5*MODBYTES_XXX]);
    FP_nres(&(g->a.b.a.b),b);
    BIG_fromBytes(b,&W->val[6*MODBYTES_XXX]);
    FP_nres(&(g->a.b.b.a),b);
    BIG_fromBytes(b,&W->val[7*MODBYTES_XXX]);
    FP_nres(&(g->a.b.b.b),b);

    BIG_fromBytes(b,&W->val[8*MODBYTES_XXX]);
    FP_nres(&(g->b.a.a.a),b);
    BIG_fromBytes(b,&W->val[9*MODBYTES_XXX]);
    FP_nres(&(g->b.a.a.b),b);
    BIG_fromBytes(b,&W->val[10*MODBYTES_XXX]);
    FP_nres(&(g->b.a.b.a),b);
    BIG_fromBytes(b,&W->val[11*MODBYTES_XXX]);
    FP_nres(&(g->b.a.b.b),b);
    BIG_fromBytes(b,&W->val[12*MODBYTES_XXX]);
    FP_nres(&(g->b.b.a.a),b);
    BIG_fromBytes(b,&W->val[13*MODBYTES_XXX]);
    FP_nres(&(g->b.b.a.b),b);
    BIG_fromBytes(b,&W->val[14*MODBYTES_XXX]);
    FP_nres(&(g->b.b.b.a),b);
    BIG_fromBytes(b,&W->val[15*MODBYTES_XXX]);
    FP_nres(&(g->b.b.b.b),b);

    BIG_fromBytes(b,&W->val[16*MODBYTES_XXX]);
    FP_nres(&(g->c.a.a.a),b);
    BIG_fromBytes(b,&W->val[17*MODBYTES_XXX]);
    FP_nres(&(g->c.a.a.b),b);
    BIG_fromBytes(b,&W->val[18*MODBYTES_XXX]);
    FP_nres(&(g->c.a.b.a),b);
    BIG_fromBytes(b,&W->val[19*MODBYTES_XXX]);
    FP_nres(&(g->c.a.b.b),b);
    BIG_fromBytes(b,&W->val[20*MODBYTES_XXX]);
    FP_nres(&(g->c.b.a.a),b);
    BIG_fromBytes(b,&W->val[21*MODBYTES_XXX]);
    FP_nres(&(g->c.b.a.b),b);
    BIG_fromBytes(b,&W->val[22*MODBYTES_XXX]);
    FP_nres(&(g->c.b.b.a),b);
    BIG_fromBytes(b,&W->val[23*MODBYTES_XXX]);
    FP_nres(&(g->c.b.b.b),b);
}

/* Move b to a if d=1 */
void YYY::FP24_cmove(FP24 *f,FP24 *g,int d)
{
    FP8_cmove(&(f->a),&(g->a),d);
    FP8_cmove(&(f->b),&(g->b),d);
    FP8_cmove(&(f->c),&(g->c),d);
}

/*
using namespace YYY;

int main() {
	int i;
	FP2 f,w0,w1,X;
	FP4 f0,f1;
	FP8 t0,t1,t2;
	FP24 w,t,lv;
	BIG a,b;
	BIG p;


	char raw[100];
	csprng RNG;                // Crypto Strong RNG 

	for (i=0; i<100; i++) raw[i]=i;

	BIG_rcopy(a,Fra);
    BIG_rcopy(b,Frb);
	FP2_from_BIGs(&X,a,b);



    RAND_seed(&RNG,100,raw);   // initialise strong RNG 

	BIG_rcopy(p,Modulus);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w0,a,b);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w1,a,b);

	FP4_from_FP2s(&f0,&w0,&w1);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w0,a,b);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w1,a,b);

	FP4_from_FP2s(&f1,&w0,&w1);
	FP8_from_FP4s(&t0,&f0,&f1);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w0,a,b);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w1,a,b);

	FP4_from_FP2s(&f0,&w0,&w1);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w0,a,b);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w1,a,b);

	FP4_from_FP2s(&f1,&w0,&w1);
	FP8_from_FP4s(&t1,&f0,&f1);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w0,a,b);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w1,a,b);

	FP4_from_FP2s(&f0,&w0,&w1);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w0,a,b);

	BIG_randomnum(a,p,&RNG);
	BIG_randomnum(b,p,&RNG);
	FP2_from_BIGs(&w1,a,b);

	FP4_from_FP2s(&f1,&w0,&w1);
	FP8_from_FP4s(&t2,&f0,&f1);

	FP24_from_FP8s(&w,&t0,&t1,&t2);


	FP24_copy(&t,&w);

	printf("w= ");
	FP24_output(&w);
	printf("\n");

	FP24_norm(&w);

	printf("w^p= ");
	FP24_frob(&w,&X);
	FP24_output(&w);
	printf("\n");	

//	printf("p.w= ");
//	FP24_ppow(&t,&t,p);
//	FP24_output(&t);
//	printf("\n");	

	printf("1/w= ");
	FP24_inv(&t,&w);
	FP24_output(&t);
	printf("\n");	

	printf("w= ");
	FP24_inv(&w,&t);
	FP24_output(&w);
	printf("\n");	

	return 0;
}

*/