// File: BlockSynthesis.sol



//SPDX-License-Identifier: Unlicense
pragma solidity ^0.8.17;

/*
    Describes a Bitcoin block with a single transaction: the concatenation of
    genTx0, extraNonce1, extraNonce2, and genTx1.
*/
struct SingleTxBitcoinBlock {
    bytes genTx0;
    bytes4 extraNonce1;
    bytes extraNonce2;
    bytes genTx1;

    bytes4 nonce;
    bytes4 bits;
    bytes4 nTime;
    bytes32 previousBlockHash;
    bytes4 version;
}

contract BlockSynthesis {
    function reverseUint256(uint256 input) internal pure returns (uint256 v) {
        v = input;

        // swap bytes
        v = ((v & 0xFF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00) >> 8) |
            ((v & 0x00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF00FF) << 8);

        // swap 2-byte long pairs
        v = ((v & 0xFFFF0000FFFF0000FFFF0000FFFF0000FFFF0000FFFF0000FFFF0000FFFF0000) >> 16) |
            ((v & 0x0000FFFF0000FFFF0000FFFF0000FFFF0000FFFF0000FFFF0000FFFF0000FFFF) << 16);

        // swap 4-byte long pairs
        v = ((v & 0xFFFFFFFF00000000FFFFFFFF00000000FFFFFFFF00000000FFFFFFFF00000000) >> 32) |
            ((v & 0x00000000FFFFFFFF00000000FFFFFFFF00000000FFFFFFFF00000000FFFFFFFF) << 32);

        // swap 8-byte long pairs
        v = ((v & 0xFFFFFFFFFFFFFFFF0000000000000000FFFFFFFFFFFFFFFF0000000000000000) >> 64) |
            ((v & 0x0000000000000000FFFFFFFFFFFFFFFF0000000000000000FFFFFFFFFFFFFFFF) << 64);

        // swap 16-byte long pairs
        v = (v >> 128) | (v << 128);
    }

    function reverseBytes(bytes memory input) internal pure returns (bytes memory) {
        bytes memory newBytes = new bytes(input.length);
        for (uint i = 0; i < input.length; i++) {
            newBytes[input.length - i - 1] = input[i];
        }
        return newBytes;
    }

    function sha256d(bytes memory data) public pure returns (bytes32) {
        return sha256(bytes.concat(sha256(data)));
    }

    function createCoinbaseTx(
        bytes memory genTx0,
        bytes4 extraNonce1,
        bytes memory extraNonce2,
        bytes memory genTx1
    ) public pure returns (bytes memory) {
        return bytes.concat(genTx0, extraNonce1, extraNonce2, genTx1);
    }

    function coinbaseHash(bytes memory coinbaseTx) public pure returns (bytes32) {
        return bytes32(reverseUint256(uint256(sha256d(coinbaseTx))));
    }

    function createBlockHeader(
        bytes4 nonce,
        bytes4 bits,
        bytes4 nTime,
        bytes32 merkleRoot,
        bytes32 previousBlockHash,
        bytes4 version
    ) public pure returns (bytes memory) {
        bytes memory headerPreReverse = bytes.concat(nonce, bits, nTime, merkleRoot, previousBlockHash, version);
        return reverseBytes(headerPreReverse);
    }

    function createSingleTxHeader(SingleTxBitcoinBlock calldata  blockData) public pure returns (bytes memory) {
        bytes memory coinbaseTx = createCoinbaseTx(blockData.genTx0, blockData.extraNonce1, blockData.extraNonce2, blockData.genTx1);
        bytes32 merkleRoot = coinbaseHash(coinbaseTx);
        return createBlockHeader(blockData.nonce, blockData.bits, blockData.nTime,
            merkleRoot, blockData.previousBlockHash, blockData.version);
    }

    function blockHash(bytes memory blockHeader) public pure returns (bytes32) {
        return bytes32(reverseUint256(uint256(sha256d(blockHeader))));
    }

    function blockDifficulty(bytes32 hash) public pure returns (uint256) {
        return 0x00000000ffffffffffffffffffffffffffffffffffffffffffffffffffffffff / uint256(hash);
    }
}

// File: CKVerifier.sol



//SPDX-License-Identifier: Unlicense
pragma solidity ^0.8.17;

import "./EllipticCurve.sol";
import "./BlockSynthesis.sol";

struct SingleTxBitcoinBlockSplit {
    bytes genTx0;
    bytes4 extraNonce1;
    bytes extraNonce2;

    // response is the first 32 bytes of genTx1
    bytes32 response;
    // remainingTx is the remaining bytes of genTx1
    bytes remainingTx;

    bytes4 nonce;
    bytes4 bits;
    bytes4 nTime;
    bytes32 previousBlockHash;
    bytes4 version;
}

/**
 * @dev A contract that attests to proofs of complete knowledge
 */
interface ICKVerifier {
    /**
     * @dev Returns true if the address has shown a proof of complete knowledge
     * to this verifier.
     */
    function isCKVerified(address addr) external returns (bool);
}

contract CKVerifier is BlockSynthesis, ICKVerifier {

    uint public constant GX = 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798;
    uint public constant GY = 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8;
    uint public constant AA = 0;
    uint public constant BB = 7;
    uint public constant PP = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F;
    uint public difficulty;
    uint public timeThreshold;
    uint public currJob = 0;
    uint public nonceRange;
    uint public numberOfRounds;

    mapping (uint => uint[]) public commitmentsX; // job id to commitment
    mapping (uint => uint[]) public commitmentsY; // job id to commitment
    mapping (uint => uint) public randomnessInputs; // job id to trusted random input
    mapping (uint => uint) public pubKeysX;
    mapping (uint => uint) public pubKeysY;
    mapping (uint => uint) public startTimes;
    mapping (address => bool) public isCKVerified;
    
    event NewJob(uint indexed jobId, uint indexed pubKeyX, uint indexed pubKeyY);
    event ChallengeInitialized(uint indexed jobId, uint indexed randomness, uint startTime);
    event CKVerified(uint indexed jobId, address indexed addr);

    constructor(uint _d, uint _tau, uint nRounds, uint _nonceRange) {
        difficulty = _d;
        timeThreshold = _tau;
        numberOfRounds = nRounds;
        // Length of extraNonce2, in bytes
        nonceRange = _nonceRange;
    }
    
    function registerJob(uint[] calldata _aX, uint[] calldata _aY, uint _pkX, uint _pkY) public returns(uint) {
        require(_aX.length == numberOfRounds);
        require(_aY.length == numberOfRounds);
        currJob += 1;
        commitmentsX[currJob] = _aX;
        commitmentsY[currJob] = _aY;
        pubKeysX[currJob] = _pkX;
        pubKeysY[currJob] = _pkY;
        emit NewJob(currJob, _pkX, _pkY);
        return currJob;
    }

    function initChallenge(uint jobId) public returns(uint) {
        require(jobId <= currJob, "Job not registered");
        require(startTimes[jobId] == 0, "Challenge already initialized");
        randomnessInputs[jobId] = uint(keccak256(abi.encodePacked(jobId, block.difficulty)));
        uint startTime = block.timestamp;
        startTimes[jobId] = startTime;
        emit ChallengeInitialized(jobId, randomnessInputs[jobId], startTime);
        return startTime;
    }

    function powAccept(bytes32 dataHash) private view returns(bool) {
        if (difficulty == 0) return true;
        return blockDifficulty(dataHash) >= difficulty;
    }

    function zkAccept(uint aX, uint aY, uint challenge, uint response, uint pkX, uint pkY) private pure returns(bool) {
        (uint gsX, uint gsY) = EllipticCurve.ecMul(response, GX, GY, AA, PP);
        (uint tempX, uint tempY) = EllipticCurve.ecMul(challenge, pkX, pkY, AA, PP);
        (uint rhsX, uint rhsY) = EllipticCurve.ecAdd(aX, aY, tempX, tempY, AA, PP);
        return (gsX == rhsX) && (gsY == rhsY);
    }
    
    function deriveAddress(uint jobId) private view returns (address) {
        uint pkx = pubKeysX[jobId];
        uint pky = pubKeysY[jobId];
        return address(uint160(uint256(keccak256(bytes.concat(bytes32(pkx), bytes32(pky))))));
    }

    function verify(uint jobId, SingleTxBitcoinBlock[] calldata blocks) public returns (bool) {
        require(startTimes[jobId] > 0, "Challenge not initiated");
        uint randomInput = randomnessInputs[jobId];
        bool accepted = wouldVerify(jobId, blocks, randomInput); 
        if (accepted) {
            address addr = deriveAddress(jobId);
            isCKVerified[addr] = true;
            emit CKVerified(jobId, addr);
        }
        // Maybe emit Verify log for ease of use only
        return accepted;
    }
	
    // Bitcoin block struct containing everything
    // accept BitcoinBlock[] - numberOfRounds in length or more
    // [genTx0, genTx1, extraNonce1, extraNonce2] -> merkle hash, previousBlockHash, nonce, bits, nTime, version
    // Instead of genTx1: have response, and remaining tx as arguments
    //    	genTx1 = bytes.concat(response, remainingTx)
    
    function wouldVerify(uint jobId, SingleTxBitcoinBlock[] calldata blocks, uint randomInput) public view returns (bool) {
        require(blocks.length == numberOfRounds);
        if (block.timestamp - startTimes[jobId] > timeThreshold) {
            return false;
        }
        // response = gentx1 (first 32 bytes)
        // bitcoin block hash
        for (uint i = 0; i < blocks.length; i++) {
	        bytes32 dataHash = blockHash(createSingleTxHeader(blocks[i]));
	        if (!powAccept(dataHash)) {
	            return false;
	        }
	        if (blocks[i].extraNonce2.length > nonceRange) {
	            return false;
	        }
        }
        // challenge can be determined by all arguments except the response (just hash all arguments together, including randomness - concat at end)
        for (uint i = 0; i < blocks.length; i++) {
            SingleTxBitcoinBlock calldata currentBlock = blocks[i];
            bytes32 challenge = sha256(bytes.concat(currentBlock.version, currentBlock.previousBlockHash,
                currentBlock.genTx0, currentBlock.extraNonce1, currentBlock.genTx1[32:],
                currentBlock.nTime, currentBlock.bits, bytes32(randomInput + i)));

	        bytes32 response = bytes32(currentBlock.genTx1[:32]);
	        
	        uint aXi = commitmentsX[jobId][i];
	        uint aYi = commitmentsY[jobId][i];
	        uint pkX = pubKeysX[jobId];
	        uint pkY = pubKeysY[jobId];
	        if (!zkAccept(aXi, aYi, uint(challenge), uint(response), pkX, pkY)) {
	        	return false;
	        }
        }
        return true;
    }
}


// File: EllipticCurve.sol



// SPDX-License-Identifier: MIT

pragma solidity >=0.5.3 <0.9.0;


/**
 * @title Elliptic Curve Library
 * @dev Library providing arithmetic operations over elliptic curves.
 * This library does not check whether the inserted points belong to the curve
 * `isOnCurve` function should be used by the library user to check the aforementioned statement.
 * @author Witnet Foundation
 */
library EllipticCurve {

  // Pre-computed constant for 2 ** 255
  uint256 constant private U255_MAX_PLUS_1 = 57896044618658097711785492504343953926634992332820282019728792003956564819968;

  /// @dev Modular euclidean inverse of a number (mod p).
  /// @param _x The number
  /// @param _pp The modulus
  /// @return q such that x*q = 1 (mod _pp)
  function invMod(uint256 _x, uint256 _pp) internal pure returns (uint256) {
    require(_x != 0 && _x != _pp && _pp != 0, "Invalid number");
    uint256 q = 0;
    uint256 newT = 1;
    uint256 r = _pp;
    uint256 t;
    while (_x != 0) {
      t = r / _x;
      (q, newT) = (newT, addmod(q, (_pp - mulmod(t, newT, _pp)), _pp));
      (r, _x) = (_x, r - t * _x);
    }

    return q;
  }

  /// @dev Modular exponentiation, b^e % _pp.
  /// Source: https://github.com/androlo/standard-contracts/blob/master/contracts/src/crypto/ECCMath.sol
  /// @param _base base
  /// @param _exp exponent
  /// @param _pp modulus
  /// @return r such that r = b**e (mod _pp)
  function expMod(uint256 _base, uint256 _exp, uint256 _pp) internal pure returns (uint256) {
    require(_pp!=0, "Modulus is zero");

    if (_base == 0)
      return 0;
    if (_exp == 0)
      return 1;

    uint256 r = 1;
    uint256 bit = U255_MAX_PLUS_1;
    assembly {
      for { } gt(bit, 0) { }{
        r := mulmod(mulmod(r, r, _pp), exp(_base, iszero(iszero(and(_exp, bit)))), _pp)
        r := mulmod(mulmod(r, r, _pp), exp(_base, iszero(iszero(and(_exp, div(bit, 2))))), _pp)
        r := mulmod(mulmod(r, r, _pp), exp(_base, iszero(iszero(and(_exp, div(bit, 4))))), _pp)
        r := mulmod(mulmod(r, r, _pp), exp(_base, iszero(iszero(and(_exp, div(bit, 8))))), _pp)
        bit := div(bit, 16)
      }
    }

    return r;
  }

  /// @dev Converts a point (x, y, z) expressed in Jacobian coordinates to affine coordinates (x', y', 1).
  /// @param _x coordinate x
  /// @param _y coordinate y
  /// @param _z coordinate z
  /// @param _pp the modulus
  /// @return (x', y') affine coordinates
  function toAffine(
    uint256 _x,
    uint256 _y,
    uint256 _z,
    uint256 _pp)
  internal pure returns (uint256, uint256)
  {
    uint256 zInv = invMod(_z, _pp);
    uint256 zInv2 = mulmod(zInv, zInv, _pp);
    uint256 x2 = mulmod(_x, zInv2, _pp);
    uint256 y2 = mulmod(_y, mulmod(zInv, zInv2, _pp), _pp);

    return (x2, y2);
  }

  /// @dev Derives the y coordinate from a compressed-format point x [[SEC-1]](https://www.secg.org/SEC1-Ver-1.0.pdf).
  /// @param _prefix parity byte (0x02 even, 0x03 odd)
  /// @param _x coordinate x
  /// @param _aa constant of curve
  /// @param _bb constant of curve
  /// @param _pp the modulus
  /// @return y coordinate y
  function deriveY(
    uint8 _prefix,
    uint256 _x,
    uint256 _aa,
    uint256 _bb,
    uint256 _pp)
  internal pure returns (uint256)
  {
    require(_prefix == 0x02 || _prefix == 0x03, "Invalid compressed EC point prefix");

    // x^3 + ax + b
    uint256 y2 = addmod(mulmod(_x, mulmod(_x, _x, _pp), _pp), addmod(mulmod(_x, _aa, _pp), _bb, _pp), _pp);
    y2 = expMod(y2, (_pp + 1) / 4, _pp);
    // uint256 cmp = yBit ^ y_ & 1;
    uint256 y = (y2 + _prefix) % 2 == 0 ? y2 : _pp - y2;

    return y;
  }

  /// @dev Check whether point (x,y) is on curve defined by a, b, and _pp.
  /// @param _x coordinate x of P1
  /// @param _y coordinate y of P1
  /// @param _aa constant of curve
  /// @param _bb constant of curve
  /// @param _pp the modulus
  /// @return true if x,y in the curve, false else
  function isOnCurve(
    uint _x,
    uint _y,
    uint _aa,
    uint _bb,
    uint _pp)
  internal pure returns (bool)
  {
    if (0 == _x || _x >= _pp || 0 == _y || _y >= _pp) {
      return false;
    }
    // y^2
    uint lhs = mulmod(_y, _y, _pp);
    // x^3
    uint rhs = mulmod(mulmod(_x, _x, _pp), _x, _pp);
    if (_aa != 0) {
      // x^3 + a*x
      rhs = addmod(rhs, mulmod(_x, _aa, _pp), _pp);
    }
    if (_bb != 0) {
      // x^3 + a*x + b
      rhs = addmod(rhs, _bb, _pp);
    }

    return lhs == rhs;
  }

  /// @dev Calculate inverse (x, -y) of point (x, y).
  /// @param _x coordinate x of P1
  /// @param _y coordinate y of P1
  /// @param _pp the modulus
  /// @return (x, -y)
  function ecInv(
    uint256 _x,
    uint256 _y,
    uint256 _pp)
  internal pure returns (uint256, uint256)
  {
    return (_x, (_pp - _y) % _pp);
  }

  /// @dev Add two points (x1, y1) and (x2, y2) in affine coordinates.
  /// @param _x1 coordinate x of P1
  /// @param _y1 coordinate y of P1
  /// @param _x2 coordinate x of P2
  /// @param _y2 coordinate y of P2
  /// @param _aa constant of the curve
  /// @param _pp the modulus
  /// @return (qx, qy) = P1+P2 in affine coordinates
  function ecAdd(
    uint256 _x1,
    uint256 _y1,
    uint256 _x2,
    uint256 _y2,
    uint256 _aa,
    uint256 _pp)
    internal pure returns(uint256, uint256)
  {
    uint x = 0;
    uint y = 0;
    uint z = 0;

    // Double if x1==x2 else add
    if (_x1==_x2) {
      // y1 = -y2 mod p
      if (addmod(_y1, _y2, _pp) == 0) {
        return(0, 0);
      } else {
        // P1 = P2
        (x, y, z) = jacDouble(
          _x1,
          _y1,
          1,
          _aa,
          _pp);
      }
    } else {
      (x, y, z) = jacAdd(
        _x1,
        _y1,
        1,
        _x2,
        _y2,
        1,
        _pp);
    }
    // Get back to affine
    return toAffine(
      x,
      y,
      z,
      _pp);
  }

  /// @dev Substract two points (x1, y1) and (x2, y2) in affine coordinates.
  /// @param _x1 coordinate x of P1
  /// @param _y1 coordinate y of P1
  /// @param _x2 coordinate x of P2
  /// @param _y2 coordinate y of P2
  /// @param _aa constant of the curve
  /// @param _pp the modulus
  /// @return (qx, qy) = P1-P2 in affine coordinates
  function ecSub(
    uint256 _x1,
    uint256 _y1,
    uint256 _x2,
    uint256 _y2,
    uint256 _aa,
    uint256 _pp)
  internal pure returns(uint256, uint256)
  {
    // invert square
    (uint256 x, uint256 y) = ecInv(_x2, _y2, _pp);
    // P1-square
    return ecAdd(
      _x1,
      _y1,
      x,
      y,
      _aa,
      _pp);
  }

  /// @dev Multiply point (x1, y1, z1) times d in affine coordinates.
  /// @param _k scalar to multiply
  /// @param _x coordinate x of P1
  /// @param _y coordinate y of P1
  /// @param _aa constant of the curve
  /// @param _pp the modulus
  /// @return (qx, qy) = d*P in affine coordinates
  function ecMul(
    uint256 _k,
    uint256 _x,
    uint256 _y,
    uint256 _aa,
    uint256 _pp)
  internal pure returns(uint256, uint256)
  {
    // Jacobian multiplication
    (uint256 x1, uint256 y1, uint256 z1) = jacMul(
      _k,
      _x,
      _y,
      1,
      _aa,
      _pp);
    // Get back to affine
    return toAffine(
      x1,
      y1,
      z1,
      _pp);
  }

  /// @dev Adds two points (x1, y1, z1) and (x2 y2, z2).
  /// @param _x1 coordinate x of P1
  /// @param _y1 coordinate y of P1
  /// @param _z1 coordinate z of P1
  /// @param _x2 coordinate x of square
  /// @param _y2 coordinate y of square
  /// @param _z2 coordinate z of square
  /// @param _pp the modulus
  /// @return (qx, qy, qz) P1+square in Jacobian
  function jacAdd(
    uint256 _x1,
    uint256 _y1,
    uint256 _z1,
    uint256 _x2,
    uint256 _y2,
    uint256 _z2,
    uint256 _pp)
  internal pure returns (uint256, uint256, uint256)
  {
    if (_x1==0 && _y1==0)
      return (_x2, _y2, _z2);
    if (_x2==0 && _y2==0)
      return (_x1, _y1, _z1);

    // We follow the equations described in https://pdfs.semanticscholar.org/5c64/29952e08025a9649c2b0ba32518e9a7fb5c2.pdf Section 5
    uint[4] memory zs; // z1^2, z1^3, z2^2, z2^3
    zs[0] = mulmod(_z1, _z1, _pp);
    zs[1] = mulmod(_z1, zs[0], _pp);
    zs[2] = mulmod(_z2, _z2, _pp);
    zs[3] = mulmod(_z2, zs[2], _pp);

    // u1, s1, u2, s2
    zs = [
      mulmod(_x1, zs[2], _pp),
      mulmod(_y1, zs[3], _pp),
      mulmod(_x2, zs[0], _pp),
      mulmod(_y2, zs[1], _pp)
    ];

    // In case of zs[0] == zs[2] && zs[1] == zs[3], double function should be used
    require(zs[0] != zs[2] || zs[1] != zs[3], "Use jacDouble function instead");

    uint[4] memory hr;
    //h
    hr[0] = addmod(zs[2], _pp - zs[0], _pp);
    //r
    hr[1] = addmod(zs[3], _pp - zs[1], _pp);
    //h^2
    hr[2] = mulmod(hr[0], hr[0], _pp);
    // h^3
    hr[3] = mulmod(hr[2], hr[0], _pp);
    // qx = -h^3  -2u1h^2+r^2
    uint256 qx = addmod(mulmod(hr[1], hr[1], _pp), _pp - hr[3], _pp);
    qx = addmod(qx, _pp - mulmod(2, mulmod(zs[0], hr[2], _pp), _pp), _pp);
    // qy = -s1*z1*h^3+r(u1*h^2 -x^3)
    uint256 qy = mulmod(hr[1], addmod(mulmod(zs[0], hr[2], _pp), _pp - qx, _pp), _pp);
    qy = addmod(qy, _pp - mulmod(zs[1], hr[3], _pp), _pp);
    // qz = h*z1*z2
    uint256 qz = mulmod(hr[0], mulmod(_z1, _z2, _pp), _pp);
    return(qx, qy, qz);
  }

  /// @dev Doubles a points (x, y, z).
  /// @param _x coordinate x of P1
  /// @param _y coordinate y of P1
  /// @param _z coordinate z of P1
  /// @param _aa the a scalar in the curve equation
  /// @param _pp the modulus
  /// @return (qx, qy, qz) 2P in Jacobian
  function jacDouble(
    uint256 _x,
    uint256 _y,
    uint256 _z,
    uint256 _aa,
    uint256 _pp)
  internal pure returns (uint256, uint256, uint256)
  {
    if (_z == 0)
      return (_x, _y, _z);

    // We follow the equations described in https://pdfs.semanticscholar.org/5c64/29952e08025a9649c2b0ba32518e9a7fb5c2.pdf Section 5
    // Note: there is a bug in the paper regarding the m parameter, M=3*(x1^2)+a*(z1^4)
    // x, y, z at this point represent the squares of _x, _y, _z
    uint256 x = mulmod(_x, _x, _pp); //x1^2
    uint256 y = mulmod(_y, _y, _pp); //y1^2
    uint256 z = mulmod(_z, _z, _pp); //z1^2

    // s
    uint s = mulmod(4, mulmod(_x, y, _pp), _pp);
    // m
    uint m = addmod(mulmod(3, x, _pp), mulmod(_aa, mulmod(z, z, _pp), _pp), _pp);

    // x, y, z at this point will be reassigned and rather represent qx, qy, qz from the paper
    // This allows to reduce the gas cost and stack footprint of the algorithm
    // qx
    x = addmod(mulmod(m, m, _pp), _pp - addmod(s, s, _pp), _pp);
    // qy = -8*y1^4 + M(S-T)
    y = addmod(mulmod(m, addmod(s, _pp - x, _pp), _pp), _pp - mulmod(8, mulmod(y, y, _pp), _pp), _pp);
    // qz = 2*y1*z1
    z = mulmod(2, mulmod(_y, _z, _pp), _pp);

    return (x, y, z);
  }

  /// @dev Multiply point (x, y, z) times d.
  /// @param _d scalar to multiply
  /// @param _x coordinate x of P1
  /// @param _y coordinate y of P1
  /// @param _z coordinate z of P1
  /// @param _aa constant of curve
  /// @param _pp the modulus
  /// @return (qx, qy, qz) d*P1 in Jacobian
  function jacMul(
    uint256 _d,
    uint256 _x,
    uint256 _y,
    uint256 _z,
    uint256 _aa,
    uint256 _pp)
  internal pure returns (uint256, uint256, uint256)
  {
    // Early return in case that `_d == 0`
    if (_d == 0) {
      return (_x, _y, _z);
    }

    uint256 remaining = _d;
    uint256 qx = 0;
    uint256 qy = 0;
    uint256 qz = 1;

    // Double and add algorithm
    while (remaining != 0) {
      if ((remaining & 1) != 0) {
        (qx, qy, qz) = jacAdd(
          qx,
          qy,
          qz,
          _x,
          _y,
          _z,
          _pp);
      }
      remaining = remaining / 2;
      (_x, _y, _z) = jacDouble(
        _x,
        _y,
        _z,
        _aa,
        _pp);
    }
    return (qx, qy, qz);
  }
}