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

//
//  fp48.swift
//
//  Created by Michael Scott on 07/07/2015.
//  Copyright (c) 2015 Michael Scott. All rights reserved.
//

/* AMCL Fp^48 functions */
/* FP48 elements are of the form a+i.b+i^2.c */

public struct FP48
{
    private var a:FP16
    private var b:FP16
    private var c:FP16
    
    /* reduce all components of this mod Modulus */
    mutating func reduce()
    {
        a.reduce()
        b.reduce()
        c.reduce()
    }
    /* normalise all components of this */
    mutating func norm()
    {
        a.norm();
        b.norm();
        c.norm();
    }
    /* Constructors */
    init(_ d:FP16)
    {
        a=FP16(d)
        b=FP16(0)
        c=FP16(0)
    }
    
    init(_ d:Int)
    {
        a=FP16(d)
        b=FP16(0)
        c=FP16(0)
    }

    init(_ d:FP16,_ e:FP16,_ f:FP16)
    {
        a=FP16(d)
        b=FP16(e)
        c=FP16(f)
    }

    init(_ x:FP48)
    {
        a=FP16(x.a)
        b=FP16(x.b)
        c=FP16(x.c)
    }
    /* test x==0 ? */
    func iszilch() -> Bool
    {
        //reduce();
        return a.iszilch() && b.iszilch() && c.iszilch()
    }

    mutating func cmove(_ g:FP48,_ d:Int)
    {
        a.cmove(g.a,d)
        b.cmove(g.b,d)
        c.cmove(g.c,d)
    }

    /* return 1 if b==c, no branching */
    private static func teq(_ b:Int32,_ c:Int32) -> Int
    {
        var x=b^c
        x-=1  // if x=0, x now -1
        return Int((x>>31)&1)
    }
    /* Constant time select from pre-computed table */
    mutating func select(_ g:[FP48],_ b:Int32)
    {

        let m=b>>31
        var babs=(b^m)-m
        
        babs=(babs-1)/2
    
        cmove(g[0],FP48.teq(babs,0)) // conditional move
        cmove(g[1],FP48.teq(babs,1))
        cmove(g[2],FP48.teq(babs,2))
        cmove(g[3],FP48.teq(babs,3))
        cmove(g[4],FP48.teq(babs,4))
        cmove(g[5],FP48.teq(babs,5))
        cmove(g[6],FP48.teq(babs,6))
        cmove(g[7],FP48.teq(babs,7))
    
        var invf=FP48(self)
        invf.conj()
        cmove(invf,Int(m&1))
    }

    /* test x==1 ? */
    public func isunity() -> Bool
    {
        let one=FP16(1)
        return a.equals(one) && b.iszilch() && c.iszilch()
    }
    /* return 1 if x==y, else 0 */
    public func equals(_ x:FP48) -> Bool
    {
        return a.equals(x.a) && b.equals(x.b) && c.equals(x.c)
    }

    /* extract a from self */
    func geta() -> FP16
    {
        return a
    }
    /* extract b */
    func getb()  -> FP16
    {
        return b
    }
    /* extract c */
    func getc() -> FP16
    {
        return c
    }
    /* copy self=x */
    public mutating func copy(_ x:FP48)
    {
        a.copy(x.a)
        b.copy(x.b)
        c.copy(x.c)
    }
    /* set self=1 */
    mutating func one()
    {
        a.one()
        b.zero()
        c.zero()
    }
    /* self=conj(self) */
    mutating func conj()
    {
        a.conj()
        b.nconj()
        c.conj()
    }

    /* Granger-Scott Unitary Squaring */
    mutating func usqr()
    {
        var A=FP16(a)
        var B=FP16(c)
        var C=FP16(b)
        var D=FP16(0)
    
        a.sqr()
        D.copy(a); D.add(a)
        a.add(D)
    
        a.norm()
        A.nconj()
    
        A.add(A)
        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()
    
        b.conj()
        b.add(b)
        c.nconj()
    
        c.add(c)
        b.add(B)
        c.add(C)
        reduce()
    
    }

    /* Chung-Hasan SQR2 method from http://cacr.uwaterloo.ca/techreports/2006/cacr2006-24.pdf */
    mutating func sqr()
    {
        var A=FP16(a)
        var B=FP16(b)
        var C=FP16(c)
        var D=FP16(a)
    
        A.sqr()
        B.mul(c)
        B.add(B); B.norm()
        C.sqr()
        D.mul(b)
        D.add(D)

        c.add(a)
        c.add(b); c.norm()
        c.sqr()
    
        a.copy(A)
    
        A.add(B)
        A.norm()
        A.add(C)
        A.add(D)
        A.norm()
    
        A.neg()
        B.times_i()
        C.times_i()
    
        a.add(B)

        b.copy(C); b.add(D)
        c.add(A)
    
        norm()
    }

    /* FP48 full multiplication this=this*y */
    mutating func mul(_ y:FP48)
    {
        var z0=FP16(a)
        var z1=FP16(0)
        var z2=FP16(b)
        var z3=FP16(0)
        var t0=FP16(a)
        var t1=FP16(y.a)
    
        z0.mul(y.a)
        z2.mul(y.b)
    
        t0.add(b)
        t1.add(y.b)

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

        t0.norm(); t1.norm()

        z3.copy(t0); z3.mul(t1)
    
        t0.copy(z0); t0.neg()
        t1.copy(z2); t1.neg()
    
        z1.add(t0)
        b.copy(z1); b.add(t1)
    
        z3.add(t1)
        z2.add(t0)
    
        t0.copy(a); t0.add(c)
        t1.copy(y.a); t1.add(y.c)

        t0.norm(); t1.norm()

        t0.mul(t1)
        z2.add(t0)
    
        t0.copy(c); t0.mul(y.c)
        t1.copy(t0); t1.neg()
    
        c.copy(z2); c.add(t1)
        z3.add(t1)
        t0.times_i()
        b.add(t0)
    
        z3.norm()
        z3.times_i()
        a.copy(z0); a.add(z3)
    
        norm()
    }

    /* Special case of multiplication arises from special form of ATE pairing line function */
    mutating func smul(_ y:FP48,_ twist:Int)
    {
        if twist == ECP.D_TYPE {
            var z0=FP16(a)
            var z2=FP16(b)
            var z3=FP16(b)
            var t0=FP16(0)
            var t1=FP16(y.a)
    
            z0.mul(y.a)
            z2.pmul(y.b.real())
            b.add(a)
            t1.adds(y.b.real())
    
            b.norm(); t1.norm()

            b.mul(t1)
            z3.add(c); z3.norm()
            z3.pmul(y.b.real())
    
            t0.copy(z0); t0.neg()
            t1.copy(z2); t1.neg()
    
            b.add(t0)
    
            b.add(t1)
            z3.add(t1)
            z2.add(t0)
    
            t0.copy(a); t0.add(c)
            t0.norm(); z3.norm()
            t0.mul(y.a)
            c.copy(z2); c.add(t0)
    
            z3.times_i()
            a.copy(z0); a.add(z3)
        }
        if twist == ECP.M_TYPE {
            var z0=FP16(a)
            var z1=FP16(0)
            var z2=FP16(0)
            var z3=FP16(0)
            var t0=FP16(a)
            var t1=FP16(0)
        
            z0.mul(y.a)
            t0.add(b)
            t0.norm()

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

            z3.copy(t0); //z3.mul(y.c);
            z3.pmul(y.c.getb())
            z3.times_i()

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

            z1.add(t0)
            b.copy(z1);
            z2.copy(t0)

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

            t0.norm()
            t1.norm()
    
            t0.mul(t1)
            z2.add(t0)

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

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

            c.copy(z2); c.add(t1)
            z3.add(t1)
            t0.times_i()
            b.add(t0)
            z3.norm();
            z3.times_i()
            a.copy(z0); a.add(z3)      
        }
        norm()
    }
    /* self=1/self */
    mutating func inverse()
    {
        var f0=FP16(a)
        var f1=FP16(b)
        var f2=FP16(a)
        var f3=FP16(0)
    
        norm()
        f0.sqr()
        f1.mul(c)
        f1.times_i()
        f0.sub(f1); f0.norm()
    
        f1.copy(c); f1.sqr()
        f1.times_i()
        f2.mul(b)
        f1.sub(f2); f1.norm()
    
        f2.copy(b); f2.sqr()
        f3.copy(a); f3.mul(c)
        f2.sub(f3); f2.norm()
    
        f3.copy(b); f3.mul(f2)
        f3.times_i()
        a.mul(f0)
        f3.add(a)
        c.mul(f1)
        c.times_i()
    
        f3.add(c); f3.norm()
        f3.inverse()
        a.copy(f0); a.mul(f3)
        b.copy(f1); b.mul(f3)
        c.copy(f2); c.mul(f3)
    }

    /* self=self^p using Frobenius */
    mutating func frob(_ f:FP2,_ n:Int)
    {

        var f2=FP2(f)
        var f3=FP2(f)
    
        f2.sqr()
        f3.mul(f2)

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

        for _ in 0 ..< n {
            a.frob(f3)
            b.frob(f3)
            c.frob(f3)

            b.qmul(f); b.times_i4(); b.times_i2()
            c.qmul(f2); c.times_i4(); c.times_i4(); c.times_i4()
        }        

    }

    /* trace function */
    func trace() -> FP16
    {
        var t=FP16(0)
        t.copy(a)
        t.imul(3)
        t.reduce()
        return t
    }

    /* convert from byte array to FP48 */
    static func fromBytes(_ w:[UInt8]) -> FP48
    {
        let RM=Int(BIG.MODBYTES)
        var t=[UInt8](repeating: 0,count: RM)
    
        for i in 0 ..< RM {t[i]=w[i]}
        var a=BIG.fromBytes(t)
        for i in 0 ..< RM {t[i]=w[i+RM]}
        var b=BIG.fromBytes(t)
        var c=FP2(a,b)
    
        for i in 0 ..< RM {t[i]=w[i+2*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+3*RM]}
        b.copy(BIG.fromBytes(t))
        var d=FP2(a,b)
    
        var ea=FP4(c,d)

        for i in 0 ..< RM {t[i]=w[i+4*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+5*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+6*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+7*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))
    
        var eb=FP4(c,d)

        var ea4=FP8(ea,eb);

        for i in 0 ..< RM {t[i]=w[i+8*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+9*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+10*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+11*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))
    
        ea.copy(FP4(c,d))

        for i in 0 ..< RM {t[i]=w[i+12*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+13*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+14*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+15*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))
    
        eb.copy(FP4(c,d))

        var eb4=FP8(ea,eb);

        let e=FP16(ea4,eb4)


        for i in 0 ..< RM {t[i]=w[i+16*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+17*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+18*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+19*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))

        ea.copy(FP4(c,d))


        for i in 0 ..< RM {t[i]=w[i+20*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+21*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+22*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+23*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))

        eb.copy(FP4(c,d))

        ea4.copy(FP8(ea,eb))

        for i in 0 ..< RM {t[i]=w[i+24*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+25*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+26*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+27*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))

        ea.copy(FP4(c,d))


        for i in 0 ..< RM {t[i]=w[i+28*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+29*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+30*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+31*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))

        eb.copy(FP4(c,d))

        eb4.copy(FP8(ea,eb))

        let f=FP16(ea4,eb4)
    
    
        for i in 0 ..< RM {t[i]=w[i+32*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+33*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+34*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+35*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))

        ea.copy(FP4(c,d))

    
        for i in 0 ..< RM {t[i]=w[i+36*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+37*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+38*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+39*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))

        eb.copy(FP4(c,d))

        ea4.copy(FP8(ea,eb))

        for i in 0 ..< RM {t[i]=w[i+40*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+41*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+42*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+43*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))

        ea.copy(FP4(c,d))

    
        for i in 0 ..< RM {t[i]=w[i+44*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+45*RM]}
        b.copy(BIG.fromBytes(t))
        c.copy(FP2(a,b))
    
        for i in 0 ..< RM {t[i]=w[i+46*RM]}
        a.copy(BIG.fromBytes(t))
        for i in 0 ..< RM {t[i]=w[i+47*RM]}
        b.copy(BIG.fromBytes(t))
        d.copy(FP2(a,b))

        eb.copy(FP4(c,d))

        eb4.copy(FP8(ea,eb))

        let g=FP16(ea4,eb4)
    
        return FP48(e,f,g)
    }

    /* convert this to byte array */
    func toBytes(_ w:inout [UInt8])
    {
        let RM=Int(BIG.MODBYTES)
        var t=[UInt8](repeating: 0,count: RM)

        a.geta().geta().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i]=t[i]}
        a.geta().geta().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+RM]=t[i]}
        a.geta().geta().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+2*RM]=t[i]}
        a.geta().geta().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+3*RM]=t[i]}
    
        a.geta().getb().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+4*RM]=t[i]}
        a.geta().getb().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+5*RM]=t[i]}
        a.geta().getb().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+6*RM]=t[i]}
        a.geta().getb().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+7*RM]=t[i]}
  
        a.getb().geta().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+8*RM]=t[i]}
        a.getb().geta().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+9*RM]=t[i]}
        a.getb().geta().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+10*RM]=t[i]}
        a.getb().geta().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+11*RM]=t[i]}
    
        a.getb().getb().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+12*RM]=t[i]}
        a.getb().getb().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+13*RM]=t[i]}
        a.getb().getb().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+14*RM]=t[i]}
        a.getb().getb().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+15*RM]=t[i]}


        b.geta().geta().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+16*RM]=t[i]}
        b.geta().geta().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+17*RM]=t[i]}
        b.geta().geta().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+18*RM]=t[i]}
        b.geta().geta().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+19*RM]=t[i]}
        b.geta().getb().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+20*RM]=t[i]}
        b.geta().getb().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+21*RM]=t[i]}
        b.geta().getb().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+22*RM]=t[i]}
        b.geta().getb().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+23*RM]=t[i]}


        b.getb().geta().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+24*RM]=t[i]}
        b.getb().geta().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+25*RM]=t[i]}
        b.getb().geta().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+26*RM]=t[i]}
        b.getb().geta().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+27*RM]=t[i]}
        b.getb().getb().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+28*RM]=t[i]}
        b.getb().getb().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+29*RM]=t[i]}
        b.getb().getb().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+30*RM]=t[i]}
        b.getb().getb().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+31*RM]=t[i]}


        c.geta().geta().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+32*RM]=t[i]}
        c.geta().geta().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+33*RM]=t[i]}
        c.geta().geta().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+34*RM]=t[i]}
        c.geta().geta().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+35*RM]=t[i]}
        c.geta().getb().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+36*RM]=t[i]}
        c.geta().getb().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+37*RM]=t[i]}
        c.geta().getb().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+38*RM]=t[i]}
        c.geta().getb().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+39*RM]=t[i]}        


        c.getb().geta().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+40*RM]=t[i]}
        c.getb().geta().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+41*RM]=t[i]}
        c.getb().geta().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+42*RM]=t[i]}
        c.getb().geta().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+43*RM]=t[i]}
        c.getb().getb().geta().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+44*RM]=t[i]}
        c.getb().getb().geta().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+45*RM]=t[i]}
        c.getb().getb().getb().getA().toBytes(&t)
        for i in 0 ..< RM {w[i+46*RM]=t[i]}
        c.getb().getb().getb().getB().toBytes(&t)
        for i in 0 ..< RM {w[i+47*RM]=t[i]}        

    }

    /* convert to hex string */
    public func toString() -> String
    {
        return ("["+a.toString()+","+b.toString()+","+c.toString()+"]")
    }
    
    /* self=self^e */
    /* Note this is simple square and multiply, so not side-channel safe */
    func pow(_ e:BIG) -> FP48
    {
        var sf=FP48(self)
        sf.norm()
        
        var e1=BIG(e)
        e1.norm()        
        var e3=BIG(e1)
        e3.pmul(3)
        e3.norm();

        var w=FP48(sf)
        let nb=e3.nbits()
 
        for i in (1...nb-2).reversed()
        {
            w.usqr()
            let bt=e3.bit(UInt(i))-e1.bit(UInt(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 */
    mutating func pinpow(_ e:Int32,_ bts:Int32)
    {
        var R=[FP48]()
        R.append(FP48(1))
        R.append(FP48(self))

        //for var i=bts-1;i>=0;i--
        for i in (0...bts-1).reversed()
        {
            let b=Int((e>>i)&1)
            R[1-b].mul(R[b])
            R[b].usqr()
        }
        copy(R[0])
    }

    public func compow(_ e :BIG,_ r :BIG) -> FP16
    {
        let f=FP2(BIG(ROM.Fra),BIG(ROM.Frb))
        let q=BIG(ROM.Modulus)       

        var g1=FP48(self)
        var g2=FP48(self)

        var m=BIG(q)
        m.mod(r)
    
        var a=BIG(e)
        a.mod(m)
    
        var b=BIG(e)
        b.div(m);
    
        var c=g1.trace()

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

        g2.frob(f,1)
        let cp=g2.trace()
        g1.conj()
        g2.mul(g1)
        let cpm1=g2.trace()
        g2.mul(g1)
        let cpm2=g2.trace()
    
        c=c.xtr_pow2(cp,cpm1,cpm2,a,b)

        return c
    }


    /* P=u0.Q0+u1*Q1+u2*Q2+u3*Q3.. */
    // Bos & Costello https://eprint.iacr.org/2013/458.pdf
    // Faz-Hernandez & Longa & Sanchez  https://eprint.iacr.org/2013/158.pdf
    // Side channel attack secure 

    static func pow16(_ q:[FP48],_ u:[BIG]) -> FP48
    {
        var g1=[FP48]()
        var g2=[FP48]()        
        var g3=[FP48]()
        var g4=[FP48]()        

        
        for _ in 0 ..< 8 {
            g1.append(FP48(0))
            g2.append(FP48(0))       
            g3.append(FP48(0))
            g4.append(FP48(0))                    
        }
        
        var r=FP48(0)
        var p=FP48(0)
        
        var t=[BIG]()
        for i in 0 ..< 16 {
            t.append(BIG(u[i]))
            t[i].norm()
        }

        var mt=BIG(0);
        var w1=[Int8](repeating: 0,count: BIG.NLEN*Int(BIG.BASEBITS)+1)           
        var s1=[Int8](repeating: 0,count: BIG.NLEN*Int(BIG.BASEBITS)+1)   
        var w2=[Int8](repeating: 0,count: BIG.NLEN*Int(BIG.BASEBITS)+1)           
        var s2=[Int8](repeating: 0,count: BIG.NLEN*Int(BIG.BASEBITS)+1)   
        var w3=[Int8](repeating: 0,count: BIG.NLEN*Int(BIG.BASEBITS)+1)           
        var s3=[Int8](repeating: 0,count: BIG.NLEN*Int(BIG.BASEBITS)+1)   
        var w4=[Int8](repeating: 0,count: BIG.NLEN*Int(BIG.BASEBITS)+1)           
        var s4=[Int8](repeating: 0,count: BIG.NLEN*Int(BIG.BASEBITS)+1)   


// precompute table 
        g1[0].copy(q[0])   // q[0]
        g1[1].copy(g1[0]); g1[1].mul(q[1])   // q[0].q[1]
        g1[2].copy(g1[0]); g1[2].mul(q[2])   // q[0].q[2]
        g1[3].copy(g1[1]); g1[3].mul(q[2])   // q[0].q[1].q[2]
        g1[4].copy(g1[0]); g1[4].mul(q[3])   // q[0].q[3]
        g1[5].copy(g1[1]); g1[5].mul(q[3])   // q[0].q[1].q[3]
        g1[6].copy(g1[2]); g1[6].mul(q[3])   // q[0].q[2].q[3]
        g1[7].copy(g1[3]); g1[7].mul(q[3])   // q[0].q[1].q[2].q[3]

// Use Frobenius
    let f=FP2(BIG(ROM.Fra),BIG(ROM.Frb))
    for i in 0 ..< 8  {
        g2[i].copy(g1[i]); g2[i].frob(f,4)       
        g3[i].copy(g2[i]); g3[i].frob(f,4)     
        g4[i].copy(g3[i]); g4[i].frob(f,4)                     
    }

// Make it odd
        let pb1=1-t[0].parity()
        t[0].inc(pb1)
        t[0].norm()  

        let pb2=1-t[4].parity()
        t[4].inc(pb2)
        t[4].norm()          

        let pb3=1-t[8].parity()
        t[8].inc(pb3)
        t[8].norm()  

        let pb4=1-t[12].parity()
        t[12].inc(pb4)
        t[12].norm()    

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

        let nb=1+mt.nbits()

// Sign pivot 

        s1[nb-1]=1
        s2[nb-1]=1
        s3[nb-1]=1
        s4[nb-1]=1        
        for i in 0 ..< nb-1 {
            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 exponent
        for i in 0 ..< nb {
            w1[i]=0
            var k=1
            for j in 1 ..< 4 {
                let 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*k
            }
            w2[i]=0
            k=1
            for j in 5 ..< 8 {
                let 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*k
            }            
            w3[i]=0
            k=1
            for j in 9 ..< 12 {
                let 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*k
            }            
            w4[i]=0
            k=1
            for j in 13 ..< 16 {
                let 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*k
            }                                    
        }   

// Main loop
        p.select(g1,Int32(2*w1[nb-1]+1))
        r.select(g2,Int32(2*w2[nb-1]+1))
        p.mul(r)        
        r.select(g3,Int32(2*w3[nb-1]+1))
        p.mul(r)    
        r.select(g4,Int32(2*w4[nb-1]+1))
        p.mul(r)                    
        for i in (0 ..< nb-1).reversed() {
            p.usqr()
            r.select(g1,Int32(2*w1[i]+s1[i]))
            p.mul(r)
            r.select(g2,Int32(2*w2[i]+s2[i]))
            p.mul(r)            
            r.select(g3,Int32(2*w3[i]+s3[i]))
            p.mul(r)
            r.select(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
    }

}