bn256.scrypt

type FQ = int;

struct FQ2 {
    FQ x;
    FQ y;
}

struct FQ6 {
    FQ2 x;
    FQ2 y;
    FQ2 z;
}

struct CurvePoint {
    FQ x;
    FQ y;
    FQ z;
    FQ t;
}

struct TwistPoint {
    FQ2 x;
    FQ2 y;
    FQ2 z;
    FQ2 t;
}

// These two are just to make it easier for users to interface with the code
// by not having them to deal with z and t coords.
struct G1Point {
    FQ x;
    FQ y;
}

struct G2Point {
    FQ2 x;
    FQ2 y;
}

// FQ12 implements the field of size p¹² as a quadratic extension of FQ6
// where ω²=τ.
struct FQ12 {
    FQ6 x;
    FQ6 y;
}

library BN256 {

    // Curve bits:
    static const int CURVE_BITS = 256; 
    static const int CURVE_BITS_P8 = 264; // +8 bits

    // Key int size:
    static const int S = 33;    // 32 bytes plus sign byte
    static const bytes mask = reverseBytes(num2bin(1, 33), 33);
    static const bytes zero = reverseBytes(num2bin(0, 33), 33);

    // Generator of G1:
    static const CurvePoint G1 = {1, 2, 1, 1};

    // Generator of G2:
    static const TwistPoint G2 = {
        {
            0x198e9393920d483a7260bfb731fb5d25f1aa493335a9e71297e485b7aef312c2,
            0x1800deef121f1e76426a00665e5c4479674322d4f75edadd46debd5cd992f6ed
        },
        {
            0x090689d0585ff075ec9e99ad690c3395bc4b313370b38ef355acdadcd122975b,
            0x12c85ea5db8c6deb4aab71808dcb408fe3d1e7690c43d37b4ce6cc0166fa7daa
        },
        {0, 1},
        {0, 1}
    };

    static const FQ2 FQ2Zero = {0, 0};
    static const FQ2 FQ2One = {0, 1};

    static const FQ6 FQ6Zero = {
        FQ2Zero, FQ2Zero, FQ2Zero
    };

    static const FQ6 FQ6One = {
        FQ2Zero, FQ2Zero, FQ2One
    };

    static const FQ12 FQ12Zero = {FQ6Zero, FQ6Zero};
    static const FQ12 FQ12One = {FQ6Zero, FQ6One};

    // Curve field modulus:
    static int P = 0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47;

    // xiToPMinus1Over6 is ξ^((p-1)/6) where ξ = i+9.
    static const FQ2 xiToPMinus1Over6 = {
        0x246996f3b4fae7e6a6327cfe12150b8e747992778eeec7e5ca5cf05f80f362ac,
        0x1284b71c2865a7dfe8b99fdd76e68b605c521e08292f2176d60b35dadcc9e470
    };

    // xiTo2PMinus2Over3 is ξ^((2p-2)/3) where ξ = i+9.
    static const FQ2 xiTo2PMinus2Over3 = {
        0x2c145edbe7fd8aee9f3a80b03b0b1c923685d2ea1bdec763c13b4711cd2b8126,
        0x05b54f5e64eea80180f3c0b75a181e84d33365f7be94ec72848a1f55921ea762
    };

    // xiToPMinus1Over2 is ξ^((p-1)/2) where ξ = i+9.
    static const FQ2 xiToPMinus1Over2 = {
        0x07c03cbcac41049a0704b5a7ec796f2b21807dc98fa25bd282d37f632623b0e3,
        0x063cf305489af5dcdc5ec698b6e2f9b9dbaae0eda9c95998dc54014671a0135a
    };

    // xiToPMinus1Over3 is ξ^((p-1)/3) where ξ = i+9.
    static const FQ2 xiToPMinus1Over3 = {
        0x16c9e55061ebae204ba4cc8bd75a079432ae2a1d0b7c9dce1665d51c640fcba2,
        0x2fb347984f7911f74c0bec3cf559b143b78cc310c2c3330c99e39557176f553d
    };

    // xiTo2PSquaredMinus2Over3 is ξ^((2p²-2)/3) where ξ = i+9 (a cubic root of unity, mod p).
    static const FQ xiTo2PSquaredMinus2Over3 = 0x59e26bcea0d48bacd4f263f1acdb5c4f5763473177fffffe;
    
    // xiToPSquaredMinus1Over3 is ξ^((p²-1)/3) where ξ = i+9.
    static const FQ xiToPSquaredMinus1Over3 = 0x30644e72e131a0295e6dd9e7e0acccb0c28f069fbb966e3de4bd44e5607cfd48;

    // xiToPSquaredMinus1Over6 is ξ^((1p²-1)/6) where ξ = i+9 (a cubic root of -1, mod p).
    static const FQ xiToPSquaredMinus1Over6 = 0x30644e72e131a0295e6dd9e7e0acccb0c28f069fbb966e3de4bd44e5607cfd49;

    static function modReduce(int k, int modulus) : int {
        int res = k % modulus;
        return (res < 0) ? res + modulus : res;
    }

    static function mulFQ2(FQ2 a, FQ2 b) : FQ2 {
        int tx = a.x * b.y;
        int t =  b.x * a.y;
        int tx_2 = tx + t;

        int ty = a.y * b.y;
        int t_2 = a.x * b.x;
        int ty_2 = ty - t_2;
        
        return {modReduce(tx_2, P), modReduce(ty_2, P)};
    }

    static function mulXiFQ2(FQ2 a) : FQ2 {
        // (xi+y)(i+3) = (9x+y)i+(9y-x)
        FQ tx = (a.x << 3) + a.x + a.y;
        FQ ty = (a.y << 3) + a.y - a.x;
        return {modReduce(tx, P), modReduce(ty, P)};
    }

    static function mulScalarFQ2(FQ2 a, int scalar) : FQ2 {
        return {
            modReduce(a.x * scalar, P),
            modReduce(a.y * scalar, P)
        };
    }

    static function addFQ2(FQ2 a, FQ2 b) : FQ2 {
        return {
            modReduce(a.x + b.x, P), 
            modReduce(a.y + b.y, P)
        };
    }

    static function subFQ2(FQ2 a, FQ2 b) : FQ2 {
        return {
            modReduce(a.x - b.x, P),
            modReduce(a.y - b.y, P)
        };
    }

    static function negFQ2(FQ2 a) : FQ2 {
        return {
            modReduce(a.x * -1, P), 
            modReduce(a.y * -1, P)
        };
    }

    static function conjugateFQ2(FQ2 a) : FQ2 {
        return {
            modReduce(a.x * -1, P), 
            modReduce(a.y, P)
        };
    }

    static function doubleFQ2(FQ2 a) : FQ2 {
        return {
            modReduce(a.x * 2, P),
            modReduce(a.y * 2, P)
        };
    }

    static function squareFQ2(FQ2 a) : FQ2 {
        int tx = a.y - a.x;
        int ty = a.x + a.y;
        int ty_2 = ty * tx;

        int tx_2 = (a.x * a.y) * 2;

        return {modReduce(tx_2, P), modReduce(ty_2, P)};
    }

    static function modInverseBranchlessP(int x) : int {
        asm {
47fd7cd8168c203c8dca7168916a81975d588181b64550b829a031e1724e6430 OP_SWAP OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_MUL OP_OVER OP_MOD OP_DUP OP_TOALTSTACK OP_DUP OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_FROMALTSTACK OP_MUL OP_OVER OP_MOD OP_NIP
        }
    }

    static function inverseFQ2(FQ2 a) : FQ2 {
        int t2 = a.y * a.y; 
        int t1 = (a.x * a.x) + t2;

        int inv = modInverseBranchlessP(t1);

        int axNeg = a.x * -1;

        return {
            modReduce(axNeg * inv, P),
            modReduce(a.y * inv, P)
        };
    }

    static function mulFQ6(FQ6 a, FQ6 b) : FQ6 {
        // "Multiplication and Squaring on Pairing-Friendly Fields"
        // Section 4, Karatsuba method.
        // http://eprint.iacr.org/2006/471.pdf

        FQ2 v0 = mulFQ2(a.z, b.z);
        FQ2 v1 = mulFQ2(a.y, b.y);
        FQ2 v2 = mulFQ2(a.x, b.x);

        FQ2 t0 = addFQ2(a.x, a.y);
        FQ2 t1 = addFQ2(b.x, b.y);
        FQ2 tz = mulFQ2(t0, t1);

        tz = subFQ2(tz, v1);
        tz = subFQ2(tz, v2);
        tz = mulXiFQ2(tz);
        tz = addFQ2(tz, v0);

        t0 = addFQ2(a.y, a.z);
        t1 = addFQ2(b.y, b.z);
        
        FQ2 ty = mulFQ2(t0, t1);
        ty = subFQ2(ty, v0);
        ty = subFQ2(ty, v1);
        t0 = mulXiFQ2(v2);
        ty = addFQ2(ty, t0);

        t0 = addFQ2(a.x, a.z);
        t1 = addFQ2(b.x, b.z);
        FQ2 tx = mulFQ2(t0, t1);
        tx = subFQ2(tx, v0);
        tx = addFQ2(tx, v1);
        tx = subFQ2(tx, v2);

        return {tx, ty, tz};
    }

    static function doubleFQ6(FQ6 a) : FQ6 {
        return {
            doubleFQ2(a.x),
            doubleFQ2(a.y),
            doubleFQ2(a.z)
        };
    }

    static function mulScalarFQ6(FQ6 a, FQ2 scalar) : FQ6 {
        return {
            mulFQ2(a.x, scalar),
            mulFQ2(a.y, scalar),
            mulFQ2(a.z, scalar)
        };
    }

    static function addFQ6(FQ6 a, FQ6 b) : FQ6 {
        return {
            addFQ2(a.x, b.x),
            addFQ2(a.y, b.y),
            addFQ2(a.z, b.z)
        };
    }

    static function subFQ6(FQ6 a, FQ6 b) : FQ6 {
        return {
            subFQ2(a.x, b.x),
            subFQ2(a.y, b.y),
            subFQ2(a.z, b.z)
        };
    }

    static function negFQ6(FQ6 a) : FQ6 {
        return {
            negFQ2(a.x),
            negFQ2(a.y),
            negFQ2(a.z)
        };
    }

    static function squareFQ6(FQ6 a) : FQ6 {
        FQ2 v0 = squareFQ2(a.z);
        FQ2 v1 = squareFQ2(a.y);
        FQ2 v2 = squareFQ2(a.x);

        FQ2 c0 = addFQ2(a.x, a.y);
        c0 = squareFQ2(c0);
        c0 = subFQ2(c0, v1);
        c0 = subFQ2(c0, v2);
        c0 = mulXiFQ2(c0);
        c0 = addFQ2(c0, v0);

        FQ2 c1 = addFQ2(a.y, a.z);
        c1 = squareFQ2(c1);
        c1 = subFQ2(c1, v0);
        c1 = subFQ2(c1, v1);
        FQ2 xiV2 = mulXiFQ2(v2);
        c1 = addFQ2(c1, xiV2);

        FQ2 c2 = addFQ2(a.x, a.z);
        c2 = squareFQ2(c2);
        c2 = subFQ2(c2, v0);
        c2 = addFQ2(c2, v1);
        c2 = subFQ2(c2, v2);

        return {c2, c1, c0};
    }

    static function mulTauFQ6(FQ6 a) : FQ6 {
        // MulTau computes τ·(aτ²+bτ+c) = bτ²+cτ+aξ
        return {
            a.y,
            a.z,
            mulXiFQ2(a.x)
        };
    }

    static function inverseFQ6(FQ6 a) : FQ6 {
        // See "Implementing cryptographic pairings", M. Scott, section 3.2.
        // ftp://136.206.11.249/pub/crypto/pairings.pdf

        // Here we can give a short explanation of how it works: let j be a cubic root of
        // unity in GF(p²) so that 1+j+j²=0.
        // Then (xτ² + yτ + z)(xj²τ² + yjτ + z)(xjτ² + yj²τ + z)
        // = (xτ² + yτ + z)(Cτ²+Bτ+A)
        // = (x³ξ²+y³ξ+z³-3ξxyz) = F is an element of the base field (the norm).
        //
        // On the other hand (xj²τ² + yjτ + z)(xjτ² + yj²τ + z)
        // = τ²(y²-ξxz) + τ(ξx²-yz) + (z²-ξxy)
        //
        // So that's why A = (z²-ξxy), B = (ξx²-yz), C = (y²-ξxz)

        FQ2 A = squareFQ2(a.z);
        FQ2 t1 = mulFQ2(a.x, a.y);
        t1 = mulXiFQ2(t1);
        A = subFQ2(A, t1);
        
        FQ2 B = squareFQ2(a.x);
        B = mulXiFQ2(B);
        t1 = mulFQ2(a.y, a.z);
        B = subFQ2(B, t1);

        FQ2 C = squareFQ2(a.y);
        t1 = mulFQ2(a.x, a.z);
        C = subFQ2(C, t1);
        
        FQ2 F = mulFQ2(C, a.y);
        F = mulXiFQ2(F);
        t1 = mulFQ2(A, a.z);
        F = addFQ2(F, t1);
        t1 = mulFQ2(B, a.x);
        t1 = mulXiFQ2(t1);
        F = addFQ2(F, t1);

        F = inverseFQ2(F);

        return {
            mulFQ2(C, F),
            mulFQ2(B, F),
            mulFQ2(A, F)
        };
    }

    static function mulScalarFQ12(FQ12 a, FQ6 b) : FQ12 {
        return {
            mulFQ6(a.x, b),
            mulFQ6(a.y, b)
        };
    }

    static function inverseFQ12(FQ12 a) : FQ12 {
        // See "Implementing cryptographic pairings", M. Scott, section 3.2.
        // ftp://136.206.11.249/pub/crypto/pairings.pdf

        FQ6 t1 = squareFQ6(a.x);
        FQ6 t2 = squareFQ6(a.y);
        FQ6 t1_2 = mulTauFQ6(t1);
        FQ6 t1_3 = subFQ6(t2, t1_2);
        FQ6 t2_2 = inverseFQ6(t1_3);

        FQ12 e = {
            negFQ6(a.x),
            a.y
        };
        
        return mulScalarFQ12(e, t2_2);
    }

    static function mulFQ12(FQ12 a, FQ12 b) : FQ12 {
        FQ6 tx = mulFQ6(a.x, b.y);
        FQ6 t = mulFQ6(b.x, a.y);
        FQ6 tx2 = addFQ6(tx, t);

        FQ6 ty = mulFQ6(a.y, b.y);
        FQ6 t2 = mulFQ6(a.x, b.x);
        FQ6 t3 = mulTauFQ6(t2);
        
        return {tx2, addFQ6(ty, t3)};
    }

    static function frobeniusFQ6(FQ6 a) : FQ6 {
        return {
            mulFQ2(conjugateFQ2(a.x), xiTo2PMinus2Over3),
            mulFQ2(conjugateFQ2(a.y), xiToPMinus1Over3),
            conjugateFQ2(a.z)
        };
    }

    static function frobeniusP2FQ6(FQ6 a) : FQ6 {
        // FrobeniusP2 computes (xτ²+yτ+z)^(p²) = xτ^(2p²) + yτ^(p²) + z
        return {
            // τ^(2p²) = τ²τ^(2p²-2) = τ²ξ^((2p²-2)/3)
            mulScalarFQ2(a.x, xiTo2PSquaredMinus2Over3),
            // τ^(p²) = ττ^(p²-1) = τξ^((p²-1)/3)
            mulScalarFQ2(a.y, xiToPSquaredMinus1Over3),
            a.z
        };
    }

    static function mulGFP(FQ6 a, int b) : FQ6 {
        return {
            mulScalarFQ2(a.x, b),
            mulScalarFQ2(a.y, b),
            mulScalarFQ2(a.z, b)
        };
    }

    static function conjugateFQ12(FQ12 a) : FQ12 {
        return {
            negFQ6(a.x),
            a.y
        };
    }

    static function frobeniusFQ12(FQ12 a) : FQ12 {
        // Frobenius computes (xω+y)^p = x^p ω·ξ^((p-1)/6) + y^p
        return {
            mulScalarFQ6(frobeniusFQ6(a.x), xiToPMinus1Over6),
            frobeniusFQ6(a.y)
        };
    }

    static function frobeniusP2FQ12(FQ12 a) : FQ12 {
        // FrobeniusP2 computes (xω+y)^p² = x^p² ω·ξ^((p²-1)/6) + y^p²
        return {
            mulGFP(frobeniusP2FQ6(a.x), xiToPSquaredMinus1Over6),
            frobeniusP2FQ6(a.y)
        };
    }

    static function squareFQ12(FQ12 a) : FQ12 {
        // Complex squaring algorithm
        FQ6 v0 = mulFQ6(a.x, a.y);

        FQ6 t = mulTauFQ6(a.x);
        FQ6 t2 = addFQ6(a.y, t);
        FQ6 ty = addFQ6(a.x, a.y);
        FQ6 ty2 = mulFQ6(ty, t2);
        FQ6 ty3 = subFQ6(ty2, v0);
        FQ6 t3 = mulTauFQ6(v0);

        FQ6 ty4 = subFQ6(ty3, t3);

        return {
            doubleFQ6(v0),
            ty4
        };
    }

    static function expFQ12_u(FQ12 a) : FQ12 {
        // u is the BN parameter that determines the prime.
        // u = 4965661367192848881;

        FQ12 sum = FQ12One;

        // Unrolled loop. Reference impl.:
        // https://github.com/ethereum/go-ethereum/blob/bd6879ac518431174a490ba42f7e6e822dcb3ee1/crypto/bn256/google/gfp12.go#L138
        FQ12 sum0 = squareFQ12(sum);
        FQ12 sum1 = mulFQ12(sum0, a);
        FQ12 sum2 = squareFQ12(sum1);
        FQ12 sum3 = squareFQ12(sum2);
        FQ12 sum4 = squareFQ12(sum3);
        FQ12 sum5 = squareFQ12(sum4);
        FQ12 sum6 = mulFQ12(sum5, a);
        FQ12 sum7 = squareFQ12(sum6);
        FQ12 sum8 = squareFQ12(sum7);
        FQ12 sum9 = squareFQ12(sum8);
        FQ12 sum10 = mulFQ12(sum9, a);
        FQ12 sum11 = squareFQ12(sum10);
        FQ12 sum12 = mulFQ12(sum11, a);
        FQ12 sum13 = squareFQ12(sum12);
        FQ12 sum14 = mulFQ12(sum13, a);
        FQ12 sum15 = squareFQ12(sum14);
        FQ12 sum16 = squareFQ12(sum15);
        FQ12 sum17 = mulFQ12(sum16, a);
        FQ12 sum18 = squareFQ12(sum17);
        FQ12 sum19 = squareFQ12(sum18);
        FQ12 sum20 = squareFQ12(sum19);
        FQ12 sum21 = mulFQ12(sum20, a);
        FQ12 sum22 = squareFQ12(sum21);
        FQ12 sum23 = mulFQ12(sum22, a);
        FQ12 sum24 = squareFQ12(sum23);
        FQ12 sum25 = squareFQ12(sum24);
        FQ12 sum26 = squareFQ12(sum25);
        FQ12 sum27 = mulFQ12(sum26, a);
        FQ12 sum28 = squareFQ12(sum27);
        FQ12 sum29 = squareFQ12(sum28);
        FQ12 sum30 = squareFQ12(sum29);
        FQ12 sum31 = mulFQ12(sum30, a);
        FQ12 sum32 = squareFQ12(sum31);
        FQ12 sum33 = squareFQ12(sum32);
        FQ12 sum34 = mulFQ12(sum33, a);
        FQ12 sum35 = squareFQ12(sum34);
        FQ12 sum36 = squareFQ12(sum35);
        FQ12 sum37 = mulFQ12(sum36, a);
        FQ12 sum38 = squareFQ12(sum37);
        FQ12 sum39 = mulFQ12(sum38, a);
        FQ12 sum40 = squareFQ12(sum39);
        FQ12 sum41 = squareFQ12(sum40);
        FQ12 sum42 = mulFQ12(sum41, a);
        FQ12 sum43 = squareFQ12(sum42);
        FQ12 sum44 = squareFQ12(sum43);
        FQ12 sum45 = squareFQ12(sum44);
        FQ12 sum46 = squareFQ12(sum45);
        FQ12 sum47 = mulFQ12(sum46, a);
        FQ12 sum48 = squareFQ12(sum47);
        FQ12 sum49 = squareFQ12(sum48);
        FQ12 sum50 = squareFQ12(sum49);
        FQ12 sum51 = mulFQ12(sum50, a);
        FQ12 sum52 = squareFQ12(sum51);
        FQ12 sum53 = squareFQ12(sum52);
        FQ12 sum54 = mulFQ12(sum53, a);
        FQ12 sum55 = squareFQ12(sum54);
        FQ12 sum56 = squareFQ12(sum55);
        FQ12 sum57 = squareFQ12(sum56);
        FQ12 sum58 = mulFQ12(sum57, a);
        FQ12 sum59 = squareFQ12(sum58);
        FQ12 sum60 = mulFQ12(sum59, a);
        FQ12 sum61 = squareFQ12(sum60);
        FQ12 sum62 = squareFQ12(sum61);
        FQ12 sum63 = mulFQ12(sum62, a);
        FQ12 sum64 = squareFQ12(sum63);
        FQ12 sum65 = squareFQ12(sum64);
        FQ12 sum66 = squareFQ12(sum65);
        FQ12 sum67 = mulFQ12(sum66, a);
        FQ12 sum68 = squareFQ12(sum67);
        FQ12 sum69 = squareFQ12(sum68);
        FQ12 sum70 = squareFQ12(sum69);
        FQ12 sum71 = squareFQ12(sum70);
        FQ12 sum72 = squareFQ12(sum71);
        FQ12 sum73 = mulFQ12(sum72, a);
        FQ12 sum74 = squareFQ12(sum73);
        FQ12 sum75 = squareFQ12(sum74);
        FQ12 sum76 = squareFQ12(sum75);
        FQ12 sum77 = mulFQ12(sum76, a);
        FQ12 sum78 = squareFQ12(sum77);
        FQ12 sum79 = mulFQ12(sum78, a);
        FQ12 sum80 = squareFQ12(sum79);
        FQ12 sum81 = mulFQ12(sum80, a);
        FQ12 sum82 = squareFQ12(sum81);
        FQ12 sum83 = mulFQ12(sum82, a);
        FQ12 sum84 = squareFQ12(sum83);
        FQ12 sum85 = mulFQ12(sum84, a);
        FQ12 sum86 = squareFQ12(sum85);
        FQ12 sum87 = squareFQ12(sum86);
        FQ12 sum88 = squareFQ12(sum87);
        FQ12 sum89 = squareFQ12(sum88);
        FQ12 sum90 = mulFQ12(sum89, a);
        
        return sum90;

    }

    static function expFQ12(FQ12 a, int power) : FQ12 {
        FQ12 sum = FQ12One;
        FQ12 t = FQ12One;

        bytes mb = reverseBytes(num2bin(power, S), S);
        bool firstOne = false;

        loop (CURVE_BITS_P8) : i {
            if (firstOne) {
                t = squareFQ12(sum);
            }

            if ((mb & (mask << ((CURVE_BITS_P8 - 1) - i))) != zero) {
                firstOne = true;
                sum = mulFQ12(t, a);
            } else {
                sum = t;
            }
        }

        return sum;
    }

    // ----------------------------------------------------

    static function doubleG1Point(G1Point a) : G1Point {
        CurvePoint res = doubleCurvePoint(
                createCurvePoint(a)
        );
        
        return getG1Point(res);
    }

    static function doubleCurvePoint(CurvePoint a) : CurvePoint {
        // See http://hyperelliptic.org/EFD/g1p/auto-code/shortw/jacobian-0/doubling/dbl-2009-l.op3
        CurvePoint res = {0, 0, 0, 0};

        int A = modReduce(a.x * a.x, P);
        int B = modReduce(a.y * a.y, P);
        int C = modReduce(B * B, P);

        int t = a.x + B;
        int t2 = modReduce(t * t, P);
        t = t2 - A;
        t2 = t - C;

        int d = t2 * 2;
        t = A * 2;
        int e = t + A;
        int f = modReduce(e * e, P);

        t = d * 2;
        res.x = f - t;

        t = C * 2;
        t2 = t * 2;
        t = t2 * 2;
        res.y = d - res.x;
        t2 = modReduce(e * res.y, P);
        res.y = t2 - t;

        //int prod = res.y * a.z;
        //if (a.t != 0) {
        //    prod = a.y * a.z;
        //}
        //res.z = modReduce(prod, P) * 2;

        int prod = a.y * a.z;
        res.z = modReduce(prod, P) * 2;

        return res;
    }

    static function addG1Points(G1Point a, G1Point b) : G1Point {
        CurvePoint res = addCurvePoints(
                createCurvePoint(a),
                createCurvePoint(b)
        );
        
        return getG1Point(res);
    }

    static function addCurvePoints(CurvePoint a, CurvePoint b) : CurvePoint {
        // See http://hyperelliptic.org/EFD/g1p/auto-code/shortw/jacobian-0/addition/add-2007-bl.op3
        CurvePoint res = {0, 0, 0, 0};
        
        if (a.z == 0) {
            res = b;
        } else if (b.z == 0) {
            res = a;
        } else {
            // Normalize the points by replacing a = [x1:y1:z1] and b = [x2:y2:z2]
            // by [u1:s1:z1·z2] and [u2:s2:z1·z2]
            // where u1 = x1·z2², s1 = y1·z2³ and u1 = x2·z1², s2 = y2·z1³
            
            int z12 = modReduce(a.z * a.z, P);
            int z22 = modReduce(b.z * b.z, P);
            
            int u1 = modReduce(a.x * z22, P);
            int u2 = modReduce(b.x * z12, P);

            int t = modReduce(b.z * z22, P);
            int s1 = modReduce(a.y * t, P);

            t = modReduce(a.z * z12, P);
            int s2 = modReduce(b.y * t, P);

            // Compute x = (2h)²(s²-u1-u2)
            // where s = (s2-s1)/(u2-u1) is the slope of the line through
            // (u1,s1) and (u2,s2). The extra factor 2h = 2(u2-u1) comes from the value of z below.
            // This is also:
            // 4(s2-s1)² - 4h²(u1+u2) = 4(s2-s1)² - 4h³ - 4h²(2u1)
            //                        = r² - j - 2v
            // with the notations below.

            int h = u2 - u1;
            bool xEqual = h == 0;

            t = h * 2;
            // i = 4h²
            int i = modReduce(t * t, P);
            // j = 4h³
            int j = modReduce(h * i, P);

            t = s2 - s1;
            bool yEqual = t == 0;

            if (xEqual && yEqual) {
                res = doubleCurvePoint(a);
            } else {
                int r = t + t;
                int v = modReduce(u1 * i, P);

                // t4 = 4(s2-s1)²
                int t4 = modReduce(r * r, P);
                int t6 = t4 - j;
                t = v * 2;

                res.x = t6 - t;

                // Set y = -(2h)³(s1 + s*(x/4h²-u1))
                // This is also
                // y = - 2·s1·j - (s2-s1)(2x - 2i·u1) = r(v-x) - 2·s1·j
                t = v - res.x;
                t4 = modReduce(s1 * j, P);
                t6 = t4 * 2;
                t4 = modReduce(r * t, P);
                res.y = t4 - t6;
                
                // Set z = 2(u2-u1)·z1·z2 = 2h·z1·z2
                t = a.z + b.z;
                t4 = modReduce(t * t, P);
                t = t4 - z12;
                t4 = t - z22;
                res.z = modReduce(t4 * h, P);
            }
        }

        return res;
    }

    static function mulG1Point(G1Point a, int m) : G1Point {
        CurvePoint res = mulCurvePoint(
                createCurvePoint(a),
                m
        );
        
        return getG1Point(res);
    }
        
    static function mulCurvePoint(CurvePoint a, int m) : CurvePoint {
        CurvePoint res = {0, 1, 0, 0};

        if (m != 0) {
            // Double and add method.
            // Lowest bit to highest.
            CurvePoint t =   {0, 0, 0, 0};
            CurvePoint sum = {0, 0, 0, 0};

            bytes mb = reverseBytes(num2bin(m, S), S);
            bool firstOne = false;

            loop (CURVE_BITS_P8) : i {
                if (firstOne) {
                    t = doubleCurvePoint(sum);
                }

                if ((mb & (mask << ((CURVE_BITS_P8 - 1) - i))) != zero) {
                    firstOne = true;
                    sum = addCurvePoints(t, a);
                } else {
                    sum = t;
                }
            }

            res = sum;
        }

        return res;
    }

    static function makeAffineCurvePoint(CurvePoint a) : CurvePoint {
        // MakeAffine converts a to affine form. If c is ∞, then it sets
        // a to 0 : 1 : 0.

        CurvePoint res = a;
        if (a.z != 1) {
            if (a.z == 0) {
                res = {0, 1, 0, 0};
            } else {
                FQ zInv = modInverseBranchlessP(a.z);
                FQ t = modReduce(a.y * zInv, P);
                FQ zInv2 = modReduce(zInv * zInv, P);
                FQ ay = modReduce(t * zInv2, P);
                FQ ax = modReduce(a.x * zInv2, P);
                
                res = {ax, ay, 1, 1};
            }
        }

        return res;
    }

    static function negCurvePoint(CurvePoint a) : CurvePoint {
        return {
            a.x,
            -a.y,
            a.z,
            0
        };
    }

    static function isInfCurvePoint(CurvePoint a) : bool {
        return a.z == 0;
    }

    static function createCurvePoint(G1Point ccp) : CurvePoint {
        return ccp == {0, 0} ? {0, 1, 0, 0} : {ccp.x, ccp.y, 1, 1};
    }

    static function getG1Point(CurvePoint cp) : G1Point {
        CurvePoint acp = makeAffineCurvePoint(cp);
        return acp == {0, 1, 0, 0} ? {0, 0} : {acp.x, acp.y};
    }

    // ----------------------------------------------------

    static function doubleG2Point(G2Point a) : G2Point {
        TwistPoint res = doubleTwistPoint(
                createTwistPoint(a)
        );
        
        return getG2Point(res);
    }

    static function doubleTwistPoint(TwistPoint a) : TwistPoint {
        // See http://hyperelliptic.org/EFD/g1p/auto-code/shortw/jacobian-0/doubling/dbl-2009-l.op3
        TwistPoint res = {FQ2Zero, FQ2Zero, FQ2Zero, FQ2Zero};

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

        FQ2 t = addFQ2(a.x, B);
        FQ2 t2 = squareFQ2(t);
        t = subFQ2(t2, A);
        t2 = subFQ2(t, C);

        FQ2 d = mulScalarFQ2(t2, 2);
        t = mulScalarFQ2(A, 2);
        FQ2 e = addFQ2(t, A);
        FQ2 f = squareFQ2(e);

        t = mulScalarFQ2(d, 2);
        res.x = subFQ2(f, t);

        t = mulScalarFQ2(C, 2);
        t2 = mulScalarFQ2(t, 2);
        t = mulScalarFQ2(t2, 2);
        res.y = subFQ2(d, res.x);
        t2 = mulFQ2(e, res.y);
        res.y = subFQ2(t2, t);

        res.z = mulScalarFQ2(mulFQ2(a.y, a.z), 2);

        return res;
    }

    static function addG2Points(G2Point a, G2Point b) : G2Point {
        TwistPoint res = addTwistPoints(
                createTwistPoint(a),
                createTwistPoint(b)
        );
        
        return getG2Point(res);
    }

    static function addTwistPoints(TwistPoint a, TwistPoint b) : TwistPoint {
        TwistPoint res = {FQ2Zero, FQ2Zero, FQ2Zero, a.t};
        
        if (a.z == FQ2Zero) {
            res = b;
        } else if (b.z == FQ2Zero) {
            res = a;
        } else {
            // See http://hyperelliptic.org/EFD/g1p/auto-code/shortw/jacobian-0/addition/add-2007-bl.op3

            // Normalize the points by replacing a = [x1:y1:z1] and b = [x2:y2:z2]
            // by [u1:s1:z1·z2] and [u2:s2:z1·z2]
            // where u1 = x1·z2², s1 = y1·z2³ and u1 = x2·z1², s2 = y2·z1³
            
            FQ2 z12 = squareFQ2(a.z);
            FQ2 z22 = squareFQ2(b.z);
            
            FQ2 u1 = mulFQ2(a.x, z22);
            FQ2 u2 = mulFQ2(b.x, z12);

            FQ2 t = mulFQ2(b.z, z22);
            FQ2 s1 = mulFQ2(a.y, t);

            t = mulFQ2(a.z, z12);
            FQ2 s2 = mulFQ2(b.y, t);

            // Compute x = (2h)²(s²-u1-u2)
            // where s = (s2-s1)/(u2-u1) is the slope of the line through
            // (u1,s1) and (u2,s2). The extra factor 2h = 2(u2-u1) comes from the value of z below.
            // This is also:
            // 4(s2-s1)² - 4h²(u1+u2) = 4(s2-s1)² - 4h³ - 4h²(2u1)
            //                        = r² - j - 2v
            // with the notations below.

            FQ2 h = subFQ2(u2, u1);
            bool xEqual = h == FQ2Zero;

            t = mulScalarFQ2(h, 2);
            // i = 4h²
            FQ2 i = squareFQ2(t);
            // j = 4h³
            FQ2 j = mulFQ2(h, i);

            t = subFQ2(s2, s1);
            bool yEqual = t == FQ2Zero;
            if (xEqual && yEqual) {
                res = doubleTwistPoint(a);
            } else {
                FQ2 r = mulScalarFQ2(t, 2);
                FQ2 v = mulFQ2(u1, i);

                // t4 = 4(s2-s1)²
                FQ2 t4 = squareFQ2(r);
                FQ2 t6 = subFQ2(t4, j);
                t = mulScalarFQ2(v, 2);

                res.x = subFQ2(t6, t);

                // Set y = -(2h)³(s1 + s*(x/4h²-u1))
                // This is also
                // y = - 2·s1·j - (s2-s1)(2x - 2i·u1) = r(v-x) - 2·s1·j
                t = subFQ2(v, res.x);
                t4 = mulFQ2(s1, j);
                t6 = mulScalarFQ2(t4, 2);
                t4 = mulFQ2(r, t);
                res.y = subFQ2(t4, t6);
                
                // Set z = 2(u2-u1)·z1·z2 = 2h·z1·z2
                t = addFQ2(a.z, b.z);
                t4 = squareFQ2(t);
                t = subFQ2(t4, z12);
                t4 = subFQ2(t, z22);
                res.z = mulFQ2(t4, h);
            }
        }

        return res;
    }

    static function makeAffineTwistPoint(TwistPoint a) : TwistPoint {
        TwistPoint res = a; 
        if (a.z != {0, 1}) {
            if (a.z == FQ2Zero) {
                res = {
                    FQ2Zero,
                    FQ2One,
                    FQ2Zero,
                    FQ2Zero
                };
            } else {
                FQ2 zInv = inverseFQ2(a.z);
                FQ2 t = mulFQ2(a.y, zInv);
                FQ2 zInv2 = squareFQ2(zInv);
                res.y = mulFQ2(t, zInv2);
                t = mulFQ2(a.x, zInv2);
                res.x = t;
                res.z = {0, 1};
                res.t = {0, 1};
            }
        }

        return res;
    }

    static function negTwistPoint(TwistPoint a) : TwistPoint {
        return {
            a.x,
            subFQ2(FQ2Zero, a.y),
            a.z,
            FQ2Zero
        };
    }

    static function isInfTwistPoint(TwistPoint a) : bool {
        return a.z == FQ2Zero;
    }

    static function createTwistPoint(G2Point ctp) : TwistPoint {
        // Output x and y coords from ZoKrates need to be swaped to work correctly.
        FQ2 x = {ctp.x.y, ctp.x.x};
        FQ2 y = {ctp.y.y, ctp.y.x};
        return ctp == {FQ2Zero, FQ2Zero} ? 
                      {FQ2Zero, FQ2One, FQ2Zero, FQ2Zero} : 
                      {x, y, FQ2One, FQ2One};
    }

    static function getG2Point(TwistPoint tp) : G2Point {
        TwistPoint atp = makeAffineTwistPoint(tp);
        return atp == {FQ2Zero, FQ2One, FQ2Zero, FQ2Zero} ?
                      {FQ2Zero, FQ2Zero} :
                      {atp.x, atp.y};
    }

}