/*
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.
*/

/* MiotCL Fp^12 functions */
/* FP12 elements are of the form a+i.b+i^2.c */

package XXX

//import "fmt"

type FP48 struct {
	a *FP16
	b *FP16
	c *FP16
	stype int
}

/* Constructors */
func NewFP48fp16(d *FP16) *FP48 {
	F := new(FP48)
	F.a = NewFP16copy(d)
	F.b = NewFP16int(0)
	F.c = NewFP16int(0)
	F.stype=FP_SPARSER
	return F
}

func NewFP48int(d int) *FP48 {
	F := new(FP48)
	F.a = NewFP16int(d)
	F.b = NewFP16int(0)
	F.c = NewFP16int(0)
	if d==1 {
		F.stype=FP_ONE
	} else {
		F.stype=FP_SPARSER
	}
	return F
}

func NewFP48fp16s(d *FP16, e *FP16, f *FP16) *FP48 {
	F := new(FP48)
	F.a = NewFP16copy(d)
	F.b = NewFP16copy(e)
	F.c = NewFP16copy(f)
	F.stype=FP_DENSE
	return F
}

func NewFP48copy(x *FP48) *FP48 {
	F := new(FP48)
	F.a = NewFP16copy(x.a)
	F.b = NewFP16copy(x.b)
	F.c = NewFP16copy(x.c)
	F.stype=x.stype
	return F
}

/* reduce all components of this mod Modulus */
func (F *FP48) reduce() {
	F.a.reduce()
	F.b.reduce()
	F.c.reduce()
}

/* normalise all components of this */
func (F *FP48) norm() {
	F.a.norm()
	F.b.norm()
	F.c.norm()
}

/* test x==0 ? */
func (F *FP48) iszilch() bool {
	return (F.a.iszilch() && F.b.iszilch() && F.c.iszilch())
}

/* Conditional move */
func (F *FP48) cmove(g *FP48, d int) {
	F.a.cmove(g.a, d)
	F.b.cmove(g.b, d)
	F.c.cmove(g.c, d)
	d=^(d-1)
	F.stype^=(F.stype^g.stype)&d
}

/* Constant time select from pre-computed table */
func (F *FP48) selector(g []*FP48, b int32) {

	m := b >> 31
	babs := (b ^ m) - m

	babs = (babs - 1) / 2

	F.cmove(g[0], teq(babs, 0)) // conditional move
	F.cmove(g[1], teq(babs, 1))
	F.cmove(g[2], teq(babs, 2))
	F.cmove(g[3], teq(babs, 3))
	F.cmove(g[4], teq(babs, 4))
	F.cmove(g[5], teq(babs, 5))
	F.cmove(g[6], teq(babs, 6))
	F.cmove(g[7], teq(babs, 7))

	invF := NewFP48copy(F)
	invF.conj()
	F.cmove(invF, int(m&1))
}

/* test x==1 ? */
func (F *FP48) Isunity() bool {
	one := NewFP16int(1)
	return (F.a.Equals(one) && F.b.iszilch() && F.c.iszilch())
}

/* return 1 if x==y, else 0 */
func (F *FP48) Equals(x *FP48) bool {
	return (F.a.Equals(x.a) && F.b.Equals(x.b) && F.c.Equals(x.c))
}

/* extract a from this */
func (F *FP48) geta() *FP16 {
	return F.a
}

/* extract b */
func (F *FP48) getb() *FP16 {
	return F.b
}

/* extract c */
func (F *FP48) getc() *FP16 {
	return F.c
}

/* copy this=x */
func (F *FP48) Copy(x *FP48) {
	F.a.copy(x.a)
	F.b.copy(x.b)
	F.c.copy(x.c)
	F.stype=x.stype
}

/* set this=1 */
func (F *FP48) one() {
	F.stype=FP_ONE
	F.a.one()
	F.b.zero()
	F.c.zero()
}

/* this=conj(this) */
func (F *FP48) conj() {
	F.a.conj()
	F.b.nconj()
	F.c.conj()
}

/* Granger-Scott Unitary Squaring */
func (F *FP48) usqr() {
	A := NewFP16copy(F.a)
	B := NewFP16copy(F.c)
	C := NewFP16copy(F.b)
	D := NewFP16int(0)

	F.a.sqr()
	D.copy(F.a)
	D.add(F.a)
	F.a.add(D)

	F.a.norm()
	A.nconj()

	A.add(A)
	F.a.add(A)
	B.sqr()
	B.times_i()

	D.copy(B)
	D.add(B)
	B.add(D)
	B.norm()

	C.sqr()
	D.copy(C)
	D.add(C)
	C.add(D)
	C.norm()

	F.b.conj()
	F.b.add(F.b)
	F.c.nconj()

	F.c.add(F.c)
	F.b.add(B)
	F.c.add(C)
	F.reduce()
	F.stype=FP_DENSE
}

/* Chung-Hasan SQR2 method from http://cacr.uwaterloo.ca/techreports/2006/cacr2006-24.pdf */
func (F *FP48) sqr() {
	if F.stype==FP_ONE {
		return
	}
	A := NewFP16copy(F.a)
	B := NewFP16copy(F.b)
	C := NewFP16copy(F.c)
	D := NewFP16copy(F.a)

	A.sqr()
	B.mul(F.c)
	B.add(B)
	B.norm()
	C.sqr()
	D.mul(F.b)
	D.add(D)

	F.c.add(F.a)
	F.c.add(F.b)
	F.c.norm()
	F.c.sqr()

	F.a.copy(A)

	A.add(B)
	A.norm()
	A.add(C)
	A.add(D)
	A.norm()

	A.neg()
	B.times_i()
	C.times_i()

	F.a.add(B)

	F.b.copy(C)
	F.b.add(D)
	F.c.add(A)
	if F.stype==FP_SPARSER {
		F.stype=FP_SPARSE
	} else {
		F.stype=FP_DENSE
	}
	F.norm()
}

/* FP48 full multiplication this=this*y */
func (F *FP48) Mul(y *FP48) {
	z0 := NewFP16copy(F.a)
	z1 := NewFP16int(0)
	z2 := NewFP16copy(F.b)
	z3 := NewFP16int(0)
	t0 := NewFP16copy(F.a)
	t1 := NewFP16copy(y.a)

	z0.mul(y.a)
	z2.mul(y.b)

	t0.add(F.b)
	t0.norm()
	t1.add(y.b)
	t1.norm()

	z1.copy(t0)
	z1.mul(t1)
	t0.copy(F.b)
	t0.add(F.c)
	t0.norm()

	t1.copy(y.b)
	t1.add(y.c)
	t1.norm()
	z3.copy(t0)
	z3.mul(t1)

	t0.copy(z0)
	t0.neg()
	t1.copy(z2)
	t1.neg()

	z1.add(t0)
	//z1.norm();
	F.b.copy(z1)
	F.b.add(t1)

	z3.add(t1)
	z2.add(t0)

	t0.copy(F.a)
	t0.add(F.c)
	t0.norm()
	t1.copy(y.a)
	t1.add(y.c)
	t1.norm()
	t0.mul(t1)
	z2.add(t0)

	t0.copy(F.c)
	t0.mul(y.c)
	t1.copy(t0)
	t1.neg()

	F.c.copy(z2)
	F.c.add(t1)
	z3.add(t1)
	t0.times_i()
	F.b.add(t0)
	z3.norm()
	z3.times_i()
	F.a.copy(z0)
	F.a.add(z3)
	F.stype=FP_DENSE
	F.norm()
}

/* FP48 full multiplication w=w*y */
/* Supports sparse multiplicands */
/* Usually w is denser than y */
func (F *FP48) ssmul(y *FP48) {
	if F.stype==FP_ONE {
		F.Copy(y)
		return
	}
	if y.stype==FP_ONE {
		return
	}
	if y.stype>=FP_SPARSE {
		z0:=NewFP16copy(F.a)
		z1:=NewFP16int(0)
		z2:=NewFP16int(0)
		z3:=NewFP16int(0)
		z0.mul(y.a)

		if SEXTIC_TWIST==M_TYPE {
			if y.stype==FP_SPARSE || F.stype==FP_SPARSE {
				z2.getb().copy(F.b.getb())
				z2.getb().mul(y.b.getb())
				z2.geta().zero()
				if y.stype!=FP_SPARSE {
					z2.geta().copy(F.b.getb())
					z2.geta().mul(y.b.geta())
				}
				if F.stype!=FP_SPARSE {
					z2.geta().copy(F.b.geta())
					z2.geta().mul(y.b.getb())
				}
				z2.times_i()
			} else {
				z2.copy(F.b)
				z2.mul(y.b)
			}
		} else {
			z2.copy(F.b)
			z2.mul(y.b)
		}
		t0:=NewFP16copy(F.a)
		t1:=NewFP16copy(y.a)
		t0.add(F.b); t0.norm();
		t1.add(y.b); t1.norm()

		z1.copy(t0); z1.mul(t1)
		t0.copy(F.b); t0.add(F.c); t0.norm()
		t1.copy(y.b); t1.add(y.c); t1.norm()

		z3.copy(t0); z3.mul(t1)

		t0.copy(z0); t0.neg()
		t1.copy(z2); t1.neg()

		z1.add(t0)
		F.b.copy(z1); F.b.add(t1)

		z3.add(t1)
		z2.add(t0)

		t0.copy(F.a); t0.add(F.c); t0.norm()
		t1.copy(y.a); t1.add(y.c); t1.norm()
	
		t0.mul(t1)
		z2.add(t0)

		if SEXTIC_TWIST==D_TYPE {
			if y.stype==FP_SPARSE || F.stype==FP_SPARSE {
				t0.geta().copy(F.c.geta())
				t0.geta().mul(y.c.geta())
				t0.getb().zero()
				if y.stype!=FP_SPARSE {
					t0.getb().copy(F.c.geta())
					t0.getb().mul(y.c.getb())
				}
				if F.stype!=FP_SPARSE {
					t0.getb().copy(F.c.getb())
					t0.getb().mul(y.c.geta())
				}
			} else {
				t0.copy(F.c)
				t0.mul(y.c)
			}
		} else {
			t0.copy(F.c)
			t0.mul(y.c)
		}
		t1.copy(t0); t1.neg()

		F.c.copy(z2); F.c.add(t1)
		z3.add(t1)
		t0.times_i()
		F.b.add(t0)
		z3.norm()
		z3.times_i()
		F.a.copy(z0); F.a.add(z3);
	} else {
		if F.stype==FP_SPARSER {
			F.smul(y)
			return
		}
		if SEXTIC_TWIST==D_TYPE { // dense by sparser - 13m 
			z0:=NewFP16copy(F.a)
			z2:=NewFP16copy(F.b)
			z3:=NewFP16copy(F.b)
			t0:=NewFP16int(0)
			t1:=NewFP16copy(y.a)
			z0.mul(y.a)
			z2.pmul(y.b.real())
			F.b.add(F.a)
			t1.real().add(y.b.real())

			t1.norm()
			F.b.norm()
			F.b.mul(t1)
			z3.add(F.c)
			z3.norm()
			z3.pmul(y.b.real())

			t0.copy(z0); t0.neg()
			t1.copy(z2); t1.neg()

			F.b.add(t0)

			F.b.add(t1)
			z3.add(t1)
			z2.add(t0)

			t0.copy(F.a); t0.add(F.c); t0.norm()
			z3.norm()
			t0.mul(y.a)
			F.c.copy(z2); F.c.add(t0)

			z3.times_i()
			F.a.copy(z0); F.a.add(z3)
		}
		if SEXTIC_TWIST==M_TYPE {
			z0:=NewFP16copy(F.a)
			z1:=NewFP16int(0)
			z2:=NewFP16int(0)
			z3:=NewFP16int(0)
			t0:=NewFP16copy(F.a)
			t1:=NewFP16int(0)
		
			z0.mul(y.a)
			t0.add(F.b); t0.norm()

			z1.copy(t0); z1.mul(y.a)
			t0.copy(F.b); t0.add(F.c)
			t0.norm()

			z3.copy(t0)
			z3.pmul(y.c.getb())
			z3.times_i()

			t0.copy(z0); t0.neg()
			z1.add(t0)
			F.b.copy(z1)
			z2.copy(t0)

			t0.copy(F.a); t0.add(F.c); t0.norm()
			t1.copy(y.a); t1.add(y.c); t1.norm()

			t0.mul(t1)
			z2.add(t0)
			t0.copy(F.c)
			
			t0.pmul(y.c.getb())
			t0.times_i()
			t1.copy(t0); t1.neg()

			F.c.copy(z2); F.c.add(t1)
			z3.add(t1)
			t0.times_i()
			F.b.add(t0)
			z3.norm()
			z3.times_i()
			F.a.copy(z0); F.a.add(z3)
		}	
	}
	F.stype=FP_DENSE
	F.norm()
}


/* Special case of multiplication arises from special form of ATE pairing line function */
func (F *FP48) smul(y *FP48) {
	if SEXTIC_TWIST==D_TYPE {	
		w1:=NewFP8copy(F.a.geta())
		w2:=NewFP8copy(F.a.getb())
		w3:=NewFP8copy(F.b.geta())

		w1.mul(y.a.geta())
		w2.mul(y.a.getb())
		w3.mul(y.b.geta())

		ta:=NewFP8copy(F.a.geta())
		tb:=NewFP8copy(y.a.geta())
		ta.add(F.a.getb()); ta.norm()
		tb.add(y.a.getb()); tb.norm()
		tc:=NewFP8copy(ta)
		tc.mul(tb);
		t:=NewFP8copy(w1)
		t.add(w2)
		t.neg()
		tc.add(t)

		ta.copy(F.a.geta()); ta.add(F.b.geta()); ta.norm()
		tb.copy(y.a.geta()); tb.add(y.b.geta()); tb.norm()
		td:=NewFP8copy(ta)
		td.mul(tb)
		t.copy(w1)
		t.add(w3)
		t.neg()
		td.add(t)

		ta.copy(F.a.getb()); ta.add(F.b.geta()); ta.norm()
		tb.copy(y.a.getb()); tb.add(y.b.geta()); tb.norm()
		te:=NewFP8copy(ta)
		te.mul(tb)
		t.copy(w2)
		t.add(w3)
		t.neg()
		te.add(t)

		w2.times_i()
		w1.add(w2)

		F.a.geta().copy(w1); F.a.getb().copy(tc)
		F.b.geta().copy(td); F.b.getb().copy(te)
		F.c.geta().copy(w3); F.c.getb().zero()

		F.a.norm()
		F.b.norm()
	} else {
		w1:=NewFP8copy(F.a.geta())
		w2:=NewFP8copy(F.a.getb())
		w3:=NewFP8copy(F.c.getb())

		w1.mul(y.a.geta())
		w2.mul(y.a.getb())
		w3.mul(y.c.getb())

		ta:=NewFP8copy(F.a.geta())
		tb:=NewFP8copy(y.a.geta())
		ta.add(F.a.getb()); ta.norm()
		tb.add(y.a.getb()); tb.norm()
		tc:=NewFP8copy(ta)
		tc.mul(tb);
		t:=NewFP8copy(w1)
		t.add(w2)
		t.neg()
		tc.add(t)

		ta.copy(F.a.geta()); ta.add(F.c.getb()); ta.norm()
		tb.copy(y.a.geta()); tb.add(y.c.getb()); tb.norm()
		td:=NewFP8copy(ta)
		td.mul(tb)
		t.copy(w1)
		t.add(w3)
		t.neg()
		td.add(t)

		ta.copy(F.a.getb()); ta.add(F.c.getb()); ta.norm()
		tb.copy(y.a.getb()); tb.add(y.c.getb()); tb.norm()
		te:=NewFP8copy(ta)
		te.mul(tb)
		t.copy(w2)
		t.add(w3)
		t.neg()
		te.add(t)

		w2.times_i()
		w1.add(w2)
		F.a.geta().copy(w1); F.a.getb().copy(tc)

		w3.times_i()
		w3.norm()
		F.b.geta().zero(); F.b.getb().copy(w3)

		te.norm()
		te.times_i()
		F.c.geta().copy(te)
		F.c.getb().copy(td)

		F.a.norm()
		F.c.norm()

	}
	F.stype=FP_SPARSE;
}

/* this=1/this */
func (F *FP48) Inverse() {
	f0 := NewFP16copy(F.a)
	f1 := NewFP16copy(F.b)
	f2 := NewFP16copy(F.a)
	f3 := NewFP16int(0)

	//F.norm()
	f0.sqr()
	f1.mul(F.c)
	f1.times_i()
	f0.sub(f1)
	f0.norm()

	f1.copy(F.c)
	f1.sqr()
	f1.times_i()
	f2.mul(F.b)
	f1.sub(f2)
	f1.norm()

	f2.copy(F.b)
	f2.sqr()
	f3.copy(F.a)
	f3.mul(F.c)
	f2.sub(f3)
	f2.norm()

	f3.copy(F.b)
	f3.mul(f2)
	f3.times_i()
	F.a.mul(f0)
	f3.add(F.a)
	F.c.mul(f1)
	F.c.times_i()

	f3.add(F.c)
	f3.norm()
	f3.inverse()
	F.a.copy(f0)
	F.a.mul(f3)
	F.b.copy(f1)
	F.b.mul(f3)
	F.c.copy(f2)
	F.c.mul(f3)
	F.stype=FP_DENSE
}

/* this=this^p using Frobenius */
func (F *FP48) frob(f *FP2, n int) {
	f2 := NewFP2copy(f)
	f3 := NewFP2copy(f)

	f2.sqr()
	f3.mul(f2)

	f3.mul_ip()
	f3.norm()
	f3.mul_ip()
	f3.norm()

	for i := 0; i < n; i++ {
		F.a.frob(f3)
		F.b.frob(f3)
		F.c.frob(f3)

		F.b.qmul(f)
		F.b.times_i4()
		F.b.times_i2()
		F.c.qmul(f2)
		F.c.times_i4()
		F.c.times_i4()
		F.c.times_i4()
	}
	F.stype=FP_DENSE
}

/* trace function */
func (F *FP48) trace() *FP16 {
	t := NewFP16int(0)
	t.copy(F.a)
	t.imul(3)
	t.reduce()
	return t
}

/* convert from byte array to FP48 */
func FP48_fromBytes(w []byte) *FP48 {
	var t [int(MODBYTES)]byte
	MB := int(MODBYTES)

	for i := 0; i < MB; i++ {
		t[i] = w[i]
	}
	a := FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+MB]
	}
	b := FromBytes(t[:])
	c := NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+2*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+3*MB]
	}
	b = FromBytes(t[:])
	d := NewFP2bigs(a, b)

	ea := NewFP4fp2s(c, d)

	for i := 0; i < MB; i++ {
		t[i] = w[i+4*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+5*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+6*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+7*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	eb := NewFP4fp2s(c, d)

	e := NewFP8fp4s(ea, eb)

	for i := 0; i < MB; i++ {
		t[i] = w[i+8*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+9*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+10*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+11*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	ea = NewFP4fp2s(c, d)

	for i := 0; i < MB; i++ {
		t[i] = w[i+12*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+13*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+14*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+15*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	eb = NewFP4fp2s(c, d)

	f := NewFP8fp4s(ea, eb)

	g := NewFP16fp8s(e, f)

	for i := 0; i < MB; i++ {
		t[i] = w[i+16*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+17*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+18*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+19*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	ea = NewFP4fp2s(c, d)

	for i := 0; i < MB; i++ {
		t[i] = w[i+20*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+21*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+22*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+23*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	eb = NewFP4fp2s(c, d)

	e = NewFP8fp4s(ea, eb)

	for i := 0; i < MB; i++ {
		t[i] = w[i+24*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+25*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+26*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+27*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	ea = NewFP4fp2s(c, d)

	for i := 0; i < MB; i++ {
		t[i] = w[i+28*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+29*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+30*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+31*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	eb = NewFP4fp2s(c, d)

	f = NewFP8fp4s(ea, eb)

	h := NewFP16fp8s(e, f)

	for i := 0; i < MB; i++ {
		t[i] = w[i+32*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+33*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+34*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+35*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	ea = NewFP4fp2s(c, d)

	for i := 0; i < MB; i++ {
		t[i] = w[i+36*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+37*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+38*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+39*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	eb = NewFP4fp2s(c, d)

	e = NewFP8fp4s(ea, eb)

	for i := 0; i < MB; i++ {
		t[i] = w[i+40*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+41*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+42*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+43*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	ea = NewFP4fp2s(c, d)

	for i := 0; i < MB; i++ {
		t[i] = w[i+44*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+45*MB]
	}
	b = FromBytes(t[:])
	c = NewFP2bigs(a, b)

	for i := 0; i < MB; i++ {
		t[i] = w[i+46*MB]
	}
	a = FromBytes(t[:])
	for i := 0; i < MB; i++ {
		t[i] = w[i+47*MB]
	}
	b = FromBytes(t[:])
	d = NewFP2bigs(a, b)

	eb = NewFP4fp2s(c, d)

	f = NewFP8fp4s(ea, eb)

	i := NewFP16fp8s(e, f)

	return NewFP48fp16s(g, h, i)
}

/* convert this to byte array */
func (F *FP48) ToBytes(w []byte) {
	var t [int(MODBYTES)]byte
	MB := int(MODBYTES)
	F.a.a.a.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i] = t[i]
	}
	F.a.a.a.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+MB] = t[i]
	}
	F.a.a.a.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+2*MB] = t[i]
	}
	F.a.a.a.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+3*MB] = t[i]
	}
	F.a.a.b.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+4*MB] = t[i]
	}
	F.a.a.b.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+5*MB] = t[i]
	}
	F.a.a.b.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+6*MB] = t[i]
	}
	F.a.a.b.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+7*MB] = t[i]
	}

	F.a.b.a.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+8*MB] = t[i]
	}
	F.a.b.a.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+9*MB] = t[i]
	}
	F.a.b.a.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+10*MB] = t[i]
	}
	F.a.b.a.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+11*MB] = t[i]
	}
	F.a.b.b.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+12*MB] = t[i]
	}
	F.a.b.b.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+13*MB] = t[i]
	}
	F.a.b.b.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+14*MB] = t[i]
	}
	F.a.b.b.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+15*MB] = t[i]
	}

	F.b.a.a.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+16*MB] = t[i]
	}
	F.b.a.a.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+17*MB] = t[i]
	}
	F.b.a.a.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+18*MB] = t[i]
	}
	F.b.a.a.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+19*MB] = t[i]
	}
	F.b.a.b.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+20*MB] = t[i]
	}
	F.b.a.b.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+21*MB] = t[i]
	}
	F.b.a.b.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+22*MB] = t[i]
	}
	F.b.a.b.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+23*MB] = t[i]
	}

	F.b.b.a.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+24*MB] = t[i]
	}
	F.b.b.a.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+25*MB] = t[i]
	}
	F.b.b.a.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+26*MB] = t[i]
	}
	F.b.b.a.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+27*MB] = t[i]
	}
	F.b.b.b.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+28*MB] = t[i]
	}
	F.b.b.b.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+29*MB] = t[i]
	}
	F.b.b.b.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+30*MB] = t[i]
	}
	F.b.b.b.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+31*MB] = t[i]
	}

	F.c.a.a.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+32*MB] = t[i]
	}
	F.c.a.a.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+33*MB] = t[i]
	}
	F.c.a.a.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+34*MB] = t[i]
	}
	F.c.a.a.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+35*MB] = t[i]
	}
	F.c.a.b.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+36*MB] = t[i]
	}
	F.c.a.b.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+37*MB] = t[i]
	}
	F.c.a.b.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+38*MB] = t[i]
	}
	F.c.a.b.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+39*MB] = t[i]
	}

	F.c.b.a.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+40*MB] = t[i]
	}
	F.c.b.a.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+41*MB] = t[i]
	}
	F.c.b.a.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+42*MB] = t[i]
	}
	F.c.b.a.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+43*MB] = t[i]
	}
	F.c.b.b.geta().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+44*MB] = t[i]
	}
	F.c.b.b.geta().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+45*MB] = t[i]
	}
	F.c.b.b.getb().GetA().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+46*MB] = t[i]
	}
	F.c.b.b.getb().GetB().ToBytes(t[:])
	for i := 0; i < MB; i++ {
		w[i+47*MB] = t[i]
	}
}

/* convert to hex string */
func (F *FP48) ToString() string {
	return ("[" + F.a.toString() + "," + F.b.toString() + "," + F.c.toString() + "]")
}

/* this=this^e */
func (F *FP48) Pow(e *BIG) *FP48 {
	sf := NewFP48copy(F)
	sf.norm()
	e1 := NewBIGcopy(e)
	e1.norm()
	e3 := NewBIGcopy(e1)
	e3.pmul(3)
	e3.norm()

	w := NewFP48copy(sf)

	nb := e3.nbits()
	for i := nb - 2; i >= 1; i-- {
		w.usqr()
		bt := e3.bit(i) - e1.bit(i)
		if bt == 1 {
			w.Mul(sf)
		}
		if bt == -1 {
			sf.conj()
			w.Mul(sf)
			sf.conj()
		}
	}
	w.reduce()
	return w

}

/* constant time powering by small integer of max length bts */
func (F *FP48) pinpow(e int, bts int) {
	var R []*FP48
	R = append(R, NewFP48int(1))
	R = append(R, NewFP48copy(F))

	for i := bts - 1; i >= 0; i-- {
		b := (e >> uint(i)) & 1
		R[1-b].Mul(R[b])
		R[b].usqr()
	}
	F.Copy(R[0])
}

/* Fast compressed FP16 power of unitary FP48 */
func (F *FP48) Compow(e *BIG, r *BIG) *FP16 {
	q := NewBIGints(Modulus)
	f := NewFP2bigs(NewBIGints(Fra), NewBIGints(Frb))

	m := NewBIGcopy(q)
	m.Mod(r)

	a := NewBIGcopy(e)
	a.Mod(m)

	b := NewBIGcopy(e)
	b.div(m)

	g1 := NewFP48copy(F)
	c := g1.trace()

	if b.iszilch() {
		c = c.xtr_pow(e)
		return c
	}

	g2 := NewFP48copy(F)
	g2.frob(f, 1)
	cp := g2.trace()

	g1.conj()
	g2.Mul(g1)
	cpm1 := g2.trace()
	g2.Mul(g1)
	cpm2 := g2.trace()

	c = c.xtr_pow2(cp, cpm1, cpm2, a, b)
	return c
}

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

func pow16(q []*FP48, u []*BIG) *FP48 {
	var g1 []*FP48
	var g2 []*FP48
	var g3 []*FP48
	var g4 []*FP48
	var w1 [NLEN*int(BASEBITS) + 1]int8
	var s1 [NLEN*int(BASEBITS) + 1]int8
	var w2 [NLEN*int(BASEBITS) + 1]int8
	var s2 [NLEN*int(BASEBITS) + 1]int8
	var w3 [NLEN*int(BASEBITS) + 1]int8
	var s3 [NLEN*int(BASEBITS) + 1]int8
	var w4 [NLEN*int(BASEBITS) + 1]int8
	var s4 [NLEN*int(BASEBITS) + 1]int8
	var t []*BIG
	r := NewFP48int(0)
	p := NewFP48int(0)
	mt := NewBIGint(0)
	var bt int8
	var k int

	for i := 0; i < 16; i++ {
		t = append(t, NewBIGcopy(u[i]))
	}

	g1 = append(g1, NewFP48copy(q[0])) // q[0]
	g1 = append(g1, NewFP48copy(g1[0]))
	g1[1].Mul(q[1]) // q[0].q[1]
	g1 = append(g1, NewFP48copy(g1[0]))
	g1[2].Mul(q[2]) // q[0].q[2]
	g1 = append(g1, NewFP48copy(g1[1]))
	g1[3].Mul(q[2]) // q[0].q[1].q[2]
	g1 = append(g1, NewFP48copy(g1[0]))
	g1[4].Mul(q[3]) // q[0].q[3]
	g1 = append(g1, NewFP48copy(g1[1]))
	g1[5].Mul(q[3]) // q[0].q[1].q[3]
	g1 = append(g1, NewFP48copy(g1[2]))
	g1[6].Mul(q[3]) // q[0].q[2].q[3]
	g1 = append(g1, NewFP48copy(g1[3]))
	g1[7].Mul(q[3]) // q[0].q[1].q[2].q[3]

	Fra := NewBIGints(Fra)
	Frb := NewBIGints(Frb)
	X := NewFP2bigs(Fra, Frb)

	// Use Frobenius
	for i := 0; i < 8; i++ {
		g2 = append(g2, NewFP48copy(g1[i]))
		g2[i].frob(X, 4)
		g3 = append(g3, NewFP48copy(g2[i]))
		g3[i].frob(X, 4)
		g4 = append(g4, NewFP48copy(g3[i]))
		g4[i].frob(X, 4)
	}

	// Make them odd
	pb1 := 1 - t[0].parity()
	t[0].inc(pb1)
	//	t[0].norm();

	pb2 := 1 - t[4].parity()
	t[4].inc(pb2)
	//	t[4].norm();

	pb3 := 1 - t[8].parity()
	t[8].inc(pb3)
	//	t[8].norm();

	pb4 := 1 - t[12].parity()
	t[12].inc(pb4)
	//	t[12].norm();

	// Number of bits
	mt.zero()
	for i := 0; i < 16; i++ {
		t[i].norm()
		mt.or(t[i])
	}

	nb := 1 + mt.nbits()

	// Sign pivot
	s1[nb-1] = 1
	s2[nb-1] = 1
	s3[nb-1] = 1
	s4[nb-1] = 1
	for i := 0; i < nb-1; i++ {
		t[0].fshr(1)
		s1[i] = 2*int8(t[0].parity()) - 1
		t[4].fshr(1)
		s2[i] = 2*int8(t[4].parity()) - 1
		t[8].fshr(1)
		s3[i] = 2*int8(t[8].parity()) - 1
		t[12].fshr(1)
		s4[i] = 2*int8(t[12].parity()) - 1

	}

	// Recoded exponents
	for i := 0; i < nb; i++ {
		w1[i] = 0
		k = 1
		for j := 1; j < 4; j++ {
			bt = s1[i] * int8(t[j].parity())
			t[j].fshr(1)
			t[j].dec(int(bt) >> 1)
			t[j].norm()
			w1[i] += bt * int8(k)
			k *= 2
		}
		w2[i] = 0
		k = 1
		for j := 5; j < 8; j++ {
			bt = s2[i] * int8(t[j].parity())
			t[j].fshr(1)
			t[j].dec(int(bt) >> 1)
			t[j].norm()
			w2[i] += bt * int8(k)
			k *= 2
		}
		w3[i] = 0
		k = 1
		for j := 9; j < 12; j++ {
			bt = s3[i] * int8(t[j].parity())
			t[j].fshr(1)
			t[j].dec(int(bt) >> 1)
			t[j].norm()
			w3[i] += bt * int8(k)
			k *= 2
		}
		w4[i] = 0
		k = 1
		for j := 13; j < 16; j++ {
			bt = s4[i] * int8(t[j].parity())
			t[j].fshr(1)
			t[j].dec(int(bt) >> 1)
			t[j].norm()
			w4[i] += bt * int8(k)
			k *= 2
		}
	}

	// Main loop
	p.selector(g1, int32(2*w1[nb-1]+1))
	r.selector(g2, int32(2*w2[nb-1]+1))
	p.Mul(r)
	r.selector(g3, int32(2*w3[nb-1]+1))
	p.Mul(r)
	r.selector(g4, int32(2*w4[nb-1]+1))
	p.Mul(r)
	for i := nb - 2; i >= 0; i-- {
		p.usqr()
		r.selector(g1, int32(2*w1[i]+s1[i]))
		p.Mul(r)
		r.selector(g2, int32(2*w2[i]+s2[i]))
		p.Mul(r)
		r.selector(g3, int32(2*w3[i]+s3[i]))
		p.Mul(r)
		r.selector(g4, int32(2*w4[i]+s4[i]))
		p.Mul(r)
	}

	// apply correction
	r.Copy(q[0])
	r.conj()
	r.Mul(p)
	p.cmove(r, pb1)
	r.Copy(q[4])
	r.conj()
	r.Mul(p)
	p.cmove(r, pb2)
	r.Copy(q[8])
	r.conj()
	r.Mul(p)
	p.cmove(r, pb3)
	r.Copy(q[12])
	r.conj()
	r.Mul(p)
	p.cmove(r, pb4)

	p.reduce()
	return p
}