library ModUtils {

    /**
     * @dev Wrap the modular exponent pre-compile introduced in Byzantium.
     * Returns base^exponent mod p.
     */
    function modExp(uint256 base, uint256 exponent, uint256 p)
        internal
        view returns(uint256 o)
    {
        /* solium-disable-next-line */
        assembly {
            // Args for the precompile: [<length_of_BASE> <length_of_EXPONENT>
            // <length_of_MODULUS> <BASE> <EXPONENT> <MODULUS>]
            let output := mload(0x40)
            let args := add(output, 0x20)
            mstore(args, 0x20)
            mstore(add(args, 0x20), 0x20)
            mstore(add(args, 0x40), 0x20)
            mstore(add(args, 0x60), base)
            mstore(add(args, 0x80), exponent)
            mstore(add(args, 0xa0), p)

            // 0x05 is the modular exponent contract address
            if iszero(staticcall(not(0), 0x05, args, 0xc0, output, 0x20)) {
                revert(0, 0)
            }
            o := mload(output)
        }
    }

    /**
     * @dev Calculates and returns the square root of a mod p,
     * or 0 if there is no such root.
     */
    function modSqrt(uint256 a, uint256 p)
        internal
        view returns(uint256)
    {

        if (legendre(a, p) != 1) {
            return 0;
        }

        if (a == 0) {
            return 0;
        }

        if (p == 2) {
            return p;
        }

        if (p % 4 == 3) {
            return modExp(a, (p + 1) / 4, p);
        }

        uint256 s = p - 1;
        uint256 e = 0;

        while (s % 2 == 0) {
            s = s / 2;
            e = e + 1;
        }

        // Note the smaller int- finding n with Legendre symbol or -1
        // should be quick
        uint256 n = 2;
        while (legendre(n, p) != -1) {
            n = n + 1;
        }

        uint256 x = modExp(a, (s + 1) / 2, p);
        uint256 b = modExp(a, s, p);
        uint256 g = modExp(n, s, p);
        uint256 r = e;
        uint256 gs = 0;
        uint256 m = 0;
        uint256 t = b;

        while (true) {
            t = b;
            m = 0;

            for (m = 0; m < r; m++) {
                if (t == 1) {
                    break;
                }
                t = modExp(t, 2, p);
            }

            if (m == 0) {
                return x;
            }

            gs = modExp(g, uint256(2) ** (r - m - 1), p);
            g = (gs * gs) % p;
            x = (x * gs) % p;
            b = (b * g) % p;
            r = m;
        }
    }

    /**
     * @dev Calculates the Legendre symbol of the given a mod p.
     * @return Returns 1 if a is a quadratic residue mod p, -1 if it is
     * a non-quadratic residue, and 0 if a is 0.
     */
    function legendre(uint256 a, uint256 p)
        internal
        view returns(int256)
    {
        uint256 raised = modExp(a, (p - 1) / uint256(2), p);

        if (raised == 0 || raised == 1) {
            return int256(raised);
        } else if (raised == p - 1) {
            return -1;
        }

        require(false, "Failed to calculate legendre.");
    }
}
function modSqrt(uint256 a, uint256 p)
    internal
    view returns(uint256)
{

    if (legendre(a, p) != 1) {
        return 0;
    }

    if (a == 0) {
        return 0;
    }

    if (p == 2) {
        return p;
    }

    if (p % 4 == 3) {
        return modExp(a, (p + 1) / 4, p);
    }

    uint256 s = p - 1;
    uint256 e = 0;

    while (s % 2 == 0) {
        s = s / 2;
        e = e + 1;
    }

    // Note the smaller int- finding n with Legendre symbol or -1
    // should be quick
    uint256 n = 2;
    while (legendre(n, p) != -1) {
        n = n + 1;
    }

    uint256 x = modExp(a, (s + 1) / 2, p);
    uint256 b = modExp(a, s, p);
    uint256 g = modExp(n, s, p);
    uint256 r = e;
    uint256 gs = 0;
    uint256 m = 0;
    uint256 t = b;

    while (true) {
        t = b;
        m = 0;

        for (m = 0; m < r; m++) {
            if (t == 1) {
                break;
            }
            t = modExp(t, 2, p);
        }

        if (m == 0) {
            return x;
        }

        gs = modExp(g, uint256(2) ** (r - m - 1), p);
        g = (gs * gs) % p;
        x = (x * gs) % p;
        b = (b * g) % p;
        r = m;
    }
}

/**
 * @dev Calculates and returns the square root of a mod p,
 * or 0 if there is no such root.
 */