bn256pairing.scrypt

import "bn256.scrypt";


struct LineFuncRes {
    FQ2 a;
    FQ2 b;
    FQ2 c;
    TwistPoint rOut;
}


library BN256Pairing {

    static function lineFuncAdd(TwistPoint r, TwistPoint p, CurvePoint q, FQ2 r2) : LineFuncRes {
        // See the mixed addition algorithm from "Faster Computation of the
        // Tate Pairing", http://arxiv.org/pdf/0904.0854v3.pdf

        FQ2 B = BN256.mulFQ2(p.x, r.t);
        FQ2 D = BN256.addFQ2(p.y, r.z);
        D = BN256.squareFQ2(D);
        D = BN256.subFQ2(D, r2);
        D = BN256.subFQ2(D, r.t);
        D = BN256.mulFQ2(D, r.t);

        FQ2 H = BN256.subFQ2(B, r.x);
        FQ2 I = BN256.squareFQ2(H);

        FQ2 E = BN256.addFQ2(I, I);
        E = BN256.addFQ2(E, E);

        FQ2 J = BN256.mulFQ2(H, E);

        FQ2 L1 = BN256.subFQ2(D, r.y);
        L1 = BN256.subFQ2(L1, r.y);

        FQ2 V = BN256.mulFQ2(r.x, E);

        FQ2 rOutX = BN256.squareFQ2(L1);
        rOutX = BN256.subFQ2(rOutX, J);
        rOutX = BN256.subFQ2(rOutX, V);
        rOutX = BN256.subFQ2(rOutX, V);

        FQ2 rOutZ = BN256.addFQ2(r.z, H);
        rOutZ = BN256.squareFQ2(rOutZ);
        rOutZ = BN256.subFQ2(rOutZ, r.t);
        rOutZ = BN256.subFQ2(rOutZ, I);

        FQ2 t = BN256.subFQ2(V, rOutX);
        t = BN256.mulFQ2(t, L1);
        FQ2 t2 = BN256.mulFQ2(r.y, J);
        t2 = BN256.addFQ2(t2, t2);
        FQ2 rOutY = BN256.subFQ2(t, t2);

        FQ2 rOutT = BN256.squareFQ2(rOutZ);

        t = BN256.addFQ2(p.y, rOutZ);
        t = BN256.squareFQ2(t);
        t = BN256.subFQ2(t, r2);
        t = BN256.subFQ2(t, rOutT);

        t2 = BN256.mulFQ2(L1, p.x);
        t2 = BN256.addFQ2(t2, t2);
        FQ2 a = BN256.subFQ2(t2, t);

        FQ2 c = BN256.mulScalarFQ2(rOutZ, q.y);
        c = BN256.addFQ2(c, c);

        FQ2 b = BN256.subFQ2(BN256.FQ2Zero, L1);
        b = BN256.mulScalarFQ2(b, q.x);
        b = BN256.addFQ2(b, b);

        TwistPoint rOut = {
            rOutX, rOutY, rOutZ, rOutT
        };

        return {a, b, c, rOut};
    }

    static function lineFuncDouble(TwistPoint r, CurvePoint q) : LineFuncRes {
        // See the doubling algorithm for a=0 from "Faster Computation of the
        // Tate Pairing", http://arxiv.org/pdf/0904.0854v3.pdf

        FQ2 A = BN256.squareFQ2(r.x);
        FQ2 B = BN256.squareFQ2(r.y);
        FQ2 C = BN256.squareFQ2(B);

        FQ2 D = BN256.addFQ2(r.x, B);
        D = BN256.squareFQ2(D);
        D = BN256.subFQ2(D, A);
        D = BN256.subFQ2(D, C);
        D = BN256.addFQ2(D, D);

        FQ2 E = BN256.addFQ2(A, A);
        E = BN256.addFQ2(E, A);

        FQ2 G = BN256.squareFQ2(E);

        FQ2 rOutX = BN256.subFQ2(G, D);
        rOutX = BN256.subFQ2(rOutX, D);

        FQ2 rOutZ = BN256.addFQ2(r.y, r.z);
        rOutZ = BN256.squareFQ2(rOutZ);
        rOutZ = BN256.subFQ2(rOutZ, B);
        rOutZ = BN256.subFQ2(rOutZ, r.t);

        FQ2 rOutY = BN256.subFQ2(D, rOutX);
        rOutY = BN256.mulFQ2(rOutY, E);
        FQ2 t = BN256.addFQ2(C, C);
        t = BN256.addFQ2(t, t);
        t = BN256.addFQ2(t, t);
        rOutY = BN256.subFQ2(rOutY, t);

        FQ2 rOutT = BN256.squareFQ2(rOutZ);

        t = BN256.mulFQ2(E, r.t);
        t = BN256.addFQ2(t, t);
        FQ2 b = BN256.subFQ2(BN256.FQ2Zero, t);
        b = BN256.mulScalarFQ2(b, q.x);

        FQ2 a = BN256.addFQ2(r.x, E);
        a = BN256.squareFQ2(a);
        a = BN256.subFQ2(a, A);
        a = BN256.subFQ2(a, G);
        t = BN256.addFQ2(B, B);
        t = BN256.addFQ2(t, t);
        a = BN256.subFQ2(a, t);

        FQ2 c = BN256.mulFQ2(rOutZ, r.t);
        c = BN256.addFQ2(c, c);
        c = BN256.mulScalarFQ2(c, q.y);

        TwistPoint rOut = {
            rOutX, rOutY, rOutZ, rOutT
        };

        return {a, b, c, rOut};
    }

    static function mulLine(FQ12 ret, FQ2 a, FQ2 b, FQ2 c) : FQ12 {
        FQ6 a2 = {BN256.FQ2Zero, a, b};
        a2 = BN256.mulFQ6(a2, ret.x);
        FQ6 t3 = BN256.mulScalarFQ6(ret.y, c);

        FQ2 t = BN256.addFQ2(b, c);
        FQ6 t2 = {BN256.FQ2Zero, a, t};

        FQ6 resX = BN256.addFQ6(ret.x, ret.y);
        FQ6 resY = t3;

        resX = BN256.mulFQ6(resX, t2);
        resX = BN256.subFQ6(resX, a2);
        resX = BN256.subFQ6(resX, resY);
        a2 = BN256.mulTauFQ6(a2);
        resY = BN256.addFQ6(resY, a2);

        return {resX, resY};
    }

    static function miller(TwistPoint q, CurvePoint p) : FQ12 {
        FQ12 ret = BN256.FQ12One;

        TwistPoint aAffine = BN256.makeAffineTwistPoint(q);
        CurvePoint bAffine = BN256.makeAffineCurvePoint(p);

        TwistPoint minusA = BN256.negTwistPoint(aAffine);

        TwistPoint r = aAffine;

        FQ2 r2 = BN256.squareFQ2(aAffine.y);

        // sixuPlus2NAF is 6u+2 in non-adjacent form.
        // Unrolled loop to get rid of in-loop branching. Reference impl.:
        // https://github.com/ethereum/go-ethereum/blob/bd6879ac518431174a490ba42f7e6e822dcb3ee1/crypto/bn256/google/optate.go#L213
        // var sixuPlus2NAF = {}int8{0, 0, 0, 1, 0, 1, 0, -1, 0, 0, 1, -1, 0, 0, 1, 0,
        //                           0, 1, 1, 0, -1, 0, 0, 1, 0, -1, 0, 0, 0, 0, 1, 1,
        //                           1, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 1,
        //                           1, 0, 0, -1, 0, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0, 1, 1}
        
        //---- 1
        LineFuncRes lfr = lineFuncDouble(r, bAffine);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- -1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, minusA, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- -1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, minusA, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- -1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, minusA, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- -1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, minusA, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- -1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, minusA, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- -1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, minusA, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- -1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, minusA, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- -1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, minusA, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 1
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        lfr = lineFuncAdd(r, aAffine, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;
        //---- 0
        lfr = lineFuncDouble(r, bAffine);
        ret = BN256.squareFQ12(ret);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        // In order to calculate Q1 we have to convert q from the sextic twist
        // to the full GF(p^12) group, apply the Frobenius there, and convert
        // back.
        //
        // The twist isomorphism is (x', y') -> (xω², yω³). If we consider just
        // x for a moment, then after applying the Frobenius, we have x̄ω^(2p)
        // where x̄ is the conjugate of x. If we are going to apply the inverse
        // isomorphism we need a value with a single coefficient of ω² so we
        // rewrite this as x̄ω^(2p-2)ω². ξ⁶ = ω and, due to the construction of
        // p, 2p-2 is a multiple of six. Therefore we can rewrite as
        // x̄ξ^((p-1)/3)ω² and applying the inverse isomorphism eliminates the
        // ω².
        //
        // A similar argument can be made for the y value.

        FQ2 q1x = BN256.conjugateFQ2(aAffine.x);
        q1x = BN256.mulFQ2(q1x, BN256.xiToPMinus1Over3);
        FQ2 q1y = BN256.conjugateFQ2(aAffine.y);
        q1y = BN256.mulFQ2(q1y, BN256.xiToPMinus1Over2);

        TwistPoint q1 = {
            q1x, q1y, {0, 1}, {0, 1}
        };

        // For Q2 we are applying the p² Frobenius. The two conjugations cancel
        // out and we are left only with the factors from the isomorphism. In
        // the case of x, we end up with a pure number which is why
        // xiToPSquaredMinus1Over3 is ∈ GF(p). With y we get a factor of -1. We
        // ignore this to end up with -Q2.

        FQ2 minusQ2x = BN256.mulScalarFQ2(aAffine.x, BN256.xiToPSquaredMinus1Over3);
        TwistPoint minusQ2 = {
            minusQ2x, aAffine.y, {0, 1}, {0, 1}
        };

        r2 = BN256.squareFQ2(q1.y);
        lfr = lineFuncAdd(r, q1, bAffine, r2);
        ret = mulLine(ret, lfr.a, lfr.b, lfr.c);
        r = lfr.rOut;

        r2 = BN256.squareFQ2(minusQ2.y);
        lfr = lineFuncAdd(r, minusQ2, bAffine, r2);
        return mulLine(ret, lfr.a, lfr.b, lfr.c);
    }

    static function finalExponentiation(FQ12 in) : FQ12 {
        FQ12 t1 = {
            BN256.negFQ6(in.x),
            in.y
        }; 

        FQ12 inv = BN256.inverseFQ12(in);
        t1 = BN256.mulFQ12(t1, inv);

        FQ12 t2 = BN256.frobeniusP2FQ12(t1);
        t1 = BN256.mulFQ12(t1, t2);

        FQ12 fp = BN256.frobeniusFQ12(t1);
        FQ12 fp2 = BN256.frobeniusP2FQ12(t1);
        FQ12 fp3 = BN256.frobeniusFQ12(fp2);

        FQ12 fu = BN256.expFQ12_u(t1);
        FQ12 fu2 = BN256.expFQ12_u(fu);
        FQ12 fu3 = BN256.expFQ12_u(fu2);

        FQ12 y3 = BN256.frobeniusFQ12(fu);
        FQ12 fu2p = BN256.frobeniusFQ12(fu2);
        FQ12 fu3p = BN256.frobeniusFQ12(fu3);
        FQ12 y2 = BN256.frobeniusP2FQ12(fu2);

        FQ12 y0 = BN256.mulFQ12(fp, fp2);
        y0 = BN256.mulFQ12(y0, fp3);

        FQ12 y1 = BN256.conjugateFQ12(t1);
        FQ12 y5 = BN256.conjugateFQ12(fu2);
        y3 = BN256.conjugateFQ12(y3);
        FQ12 y4 = BN256.mulFQ12(fu ,fu2p);
        y4 = BN256.conjugateFQ12(y4);

        FQ12 y6 = BN256.mulFQ12(fu3, fu3p);
        y6 = BN256.conjugateFQ12(y6);

        FQ12 t0 = BN256.squareFQ12(y6);
        t0 = BN256.mulFQ12(t0, y4);
        t0 = BN256.mulFQ12(t0, y5);
        t1 = BN256.mulFQ12(y3, y5);
        t1 = BN256.mulFQ12(t1, t0);
        t0 = BN256.mulFQ12(t0, y2);
        t1 = BN256.squareFQ12(t1);
        t1 = BN256.mulFQ12(t1, t0);
        t1 = BN256.squareFQ12(t1);
        t0 = BN256.mulFQ12(t1, y1);
        t1 = BN256.mulFQ12(t1, y0);
        t0 = BN256.squareFQ12(t0);
        t0 = BN256.mulFQ12(t0, t1);

        return t0;
    }

    static function pairInternal(CurvePoint g1, TwistPoint g2) : FQ12 {
        FQ12 e = miller(g2, g1);
        FQ12 ret = finalExponentiation(e);

        if (BN256.isInfTwistPoint(g2) || BN256.isInfCurvePoint(g1)) {
            ret = BN256.FQ12One;
        }

        return ret;
    }

    static function pair(G1Point g1, G2Point g2) : FQ12 {
        return pairInternal(
                BN256.createCurvePoint(g1), 
                BN256.createTwistPoint(g2)
            );
    }

    static function pairCheckP4PrecalcInternal(
            CurvePoint a0, TwistPoint b0,
            FQ12 millerBetaAlpha,
            CurvePoint a2, TwistPoint b2,
            CurvePoint a3, TwistPoint b3
            ) : bool {
        FQ12 acc = BN256.FQ12One;

        a0 = BN256.makeAffineCurvePoint(a0);
        a2 = BN256.makeAffineCurvePoint(a2);
        a3 = BN256.makeAffineCurvePoint(a3);

        if (!BN256.isInfCurvePoint(a0) && !BN256.isInfTwistPoint(b0)) {
            acc = BN256.mulFQ12(acc, miller(b0, a0));
        }
        acc = BN256.mulFQ12(acc, millerBetaAlpha);
        if (!BN256.isInfCurvePoint(a2) && !BN256.isInfTwistPoint(b2)) {
            acc = BN256.mulFQ12(acc, miller(b2, a2));
        }
        if (!BN256.isInfCurvePoint(a3) && !BN256.isInfTwistPoint(b3)) {
            acc = BN256.mulFQ12(acc, miller(b3, a3));
        }

        acc = finalExponentiation(acc);

        return acc == BN256.FQ12One;
    }

    // Check four pairs.
    // A * B + inputs * (-gamma) + C * (-delta) == alpha * beta
    // where alpha * beta is precalculated.
    static function pairCheckP4Precalc(
            G1Point a0, G2Point b0,
            FQ12 millerBetaAlpha,
            G1Point a2, G2Point b2,
            G1Point a3, G2Point b3) : bool {
        return pairCheckP4PrecalcInternal(
                BN256.createCurvePoint(a0), BN256.createTwistPoint(b0),
                millerBetaAlpha,
                BN256.createCurvePoint(a2), BN256.createTwistPoint(b2),
                BN256.createCurvePoint(a3), BN256.createTwistPoint(b3)
            );
    }

}