optimized.sol

pragma solidity ^0.4.0;

contract ECops {
    uint256 constant n = 0x30644E72E131A029B85045B68181585D97816A916871CA8D3C208C16D87CFD47;
    uint256 constant a = 0;
    uint256 constant b = 3;
    
    // Returns the inverse in the field of modulo n 
    function inverse(uint256 num) public pure returns(uint256 invNum) {
        uint256 t = 0;
        uint256 newT = 1;
        uint256 r = n;
        uint256 newR = num;
        uint256 q;
        while (newR != 0) {
            (q, r) = divmod(r, newR, n);
            (t, newT) = (newT, addmod(t, (n - mulmod(q, newT, n)), n));
            (r, newR) = (newR, r - q * newR );
        }

        invNum = t;
    }
    
    // Double an elliptic curve point
    function twiceProj(uint256 x0, uint256 y0, uint256 z0) public pure
        returns(uint256 x1, uint256 y1, uint256 z1)
    {
        uint256 t;
        uint256 u;
        uint256 v;
        uint256 w;

        if(x0 == 0 && y0 == 0) {
            return (0,0,1);
        }

        u = y0 * z0 << 1;

        v = u * x0 * y0 << 1;

        x0 = mulmod(x0, x0, n);
        t = mulmod(x0, 3, n);
        z0 = mulmod(z0, z0, n);
        z0 = mulmod(z0, a, n);
        t = addmod(t, z0, n);

        w = mulmod(t, t, n);
        x0 = v << 1;
        w = addmod(w, n-x0, n);

        x0 = addmod(v, n-w, n);
        x0 = mulmod(t, x0, n);
        y0 = mulmod(y0, u, n);
        y0 = mulmod(y0, y0, n);
        y0 = y0 << 1;
        y1 = addmod(x0, n-y0, n);

        x1 = mulmod(u, w, n);

        z1 = mulmod(u, u, n);
        z1 = mulmod(z1, u, n);
    }

    // Add elliptic curve points
    function addProj(uint256 x0, uint256 y0, uint256 z0,
                     uint256 x1, uint256 y1, uint256 z1) public pure
        returns(uint256 x2, uint256 y2, uint256 z2)
    {
        uint256 u0;
        uint256 u1;
        uint256 v = mulmod(z0, z1, n);

        if (x0 == 0 && y0 == 0) {
            return (x1, y1, z1);
        } 
        else if (x1 == 0 && y1 == 0) {
            return (x0, y0, z0);
        }
        
        u0 = mulmod(x0, z1, n);
        u1 = mulmod(x1, z0, n);

        if (u0 == u1) {
            if (mulmod(y0, z1, n) == mulmod(y1, z0, n)) {
                return twiceProj(x0, y0, z0);
            }
            else {
                return (0, 0, 1);
            }
        }

        (x2, y2, z2) = addProj2(v, u0, u1, mulmod(y1, z0, n), mulmod(y0, z1, n));
    }

    // Multiple an elliptic curve point in a 2 power base (i.e., (2^exp)*P))
    function multiplyPowerBase2(uint256 x0, uint256 y0, uint exp) public pure
        returns(uint256 x1, uint256 y1)
    {
        uint256 base2X = x0;
        uint256 base2Y = y0;
        uint256 base2Z = 1;

        for(uint i = 0; i < exp; i++) {
            (base2X, base2Y, base2Z) = twiceProj(base2X, base2Y, base2Z);
        }
        return toAffinePoint(base2X, base2Y, base2Z);
    }

    // Multiply an elliptic curve point in a scalar
    function multiplyScalar(uint256 x0, uint256 y0,
                            uint scalar) public pure
        returns(uint256 x1, uint256 y1)
    {
        uint256 base2X = x0;
        uint256 base2Y = y0;
        uint256 base2Z = 1;
        uint256 z1 = 1;

        if(scalar == 0) {
            return (0, 0);
        } else if(scalar == 1) {
            return (x0, y0);
        } else if(scalar == 2) {
            return twice(x0, y0);
        }

        if(scalar%2 == 0) {
            (x1, y1) = (0, 0);
        } else {
            (x1, y1) = (x0, y0);
        }

        scalar = scalar >> 1;

        while(scalar > 0) {
            (base2X, base2Y, base2Z) = twiceProj(base2X, base2Y, base2Z);

            if(scalar%2 == 1) {
                (x1, y1, z1) = addProj(base2X, base2Y, base2Z, x1, y1, z1);
            }

            scalar = scalar >> 1;
        }

        return toAffinePoint(x1, y1, z1);
    }
}