Solady Cube Root

Can You Prove the Code is Correct?

/// @dev Returns the cube root of `x`, rounded down.
function cbrt(uint256 x) internal pure returns (uint256 z) {
  /// @solidity memory-safe-assembly
  assembly {
      let r := shl(7, lt(0xffffffffffffffffffffffffffffffff, x))
      r := or(r, shl(6, lt(0xffffffffffffffff, shr(r, x))))
      r := or(r, shl(5, lt(0xffffffff, shr(r, x))))
      r := or(r, shl(4, lt(0xffff, shr(r, x))))
      r := or(r, shl(3, lt(0xff, shr(r, x))))

      z := div(shl(div(r, 3), shl(lt(0xf, shr(r, x)), 0xf)), xor(7, mod(r, 3)))

      z := div(add(add(div(x, mul(z, z)), z), z), 3)
      z := div(add(add(div(x, mul(z, z)), z), z), 3)
      z := div(add(add(div(x, mul(z, z)), z), z), 3)
      z := div(add(add(div(x, mul(z, z)), z), z), 3)
      z := div(add(add(div(x, mul(z, z)), z), z), 3)
      z := div(add(add(div(x, mul(z, z)), z), z), 3)
      z := div(add(add(div(x, mul(z, z)), z), z), 3)

      z := sub(z, lt(div(x, mul(z, z)), z))
  }
}

Questions A

1. What's the purpose of the first part of the code?
2. If one line of z := div(add(add(div(x, mul(z, z)), z), z), 3) was removed, can you find an counterexample?
3. If the last line z := sub(z, lt(div(x, mul(z, z)), z)) was removed, can you find an counterexample?
4. The EVM can only handle integers. If the variables are floats, what conclusion can be reached about the precision of the estimation?

Questions B

/// @dev Returns the square root of `x`, denominated in `WAD`, rounded down.
function sqrtWad(uint256 x) internal pure returns (uint256 z) {
  unchecked {
      if (x <= type(uint256).max / 10 ** 18) return sqrt(x * 10 ** 18);
      z = (1 + sqrt(x)) * 10 ** 9;
      z = (fullMulDivUnchecked(x, 10 ** 18, z) + z) >> 1;
  }
  /// @solidity memory-safe-assembly
  assembly {
      z := sub(z, gt(999999999999999999, sub(mulmod(z, z, x), 1)))
  }
}

1. What's the purpose of the last line? Why does it work?

/// @dev Returns the cube root of `x`, denominated in `WAD`, rounded down.
function cbrtWad(uint256 x) internal pure returns (uint256 z) {
  unchecked {
      if (x <= type(uint256).max / 10 ** 36) return cbrt(x * 10 ** 36);
      z = (1 + cbrt(x)) * 10 ** 12;
      z = (fullMulDivUnchecked(x, 10 ** 36, z * z) + z + z) / 3;
      x = fullMulDivUnchecked(x, 10 ** 36, z * z);
  }
  /// @solidity memory-safe-assembly
  assembly {
      z := sub(z, lt(x, z))
  }
}

2. fullMulDivUnchecked is expensive, can you use the same trick to replace the second fullMulDivUnchecked?

Answer - Optimized cbrtWad

https://github.com/Vectorized/solady/blob/e6d2a1c848921ca90a9ddac78ab24ffa5c7fb54f/src/utils/FixedPointMathLib.sol#L740-L764

Answer - The Proof

https://github.com/xuwinnie/reports/blob/main/solady-cbrt-proof-xuwinnie-202407.pdf

Questions C

1. Did you find any mistakes in the proof?
2. The optimized cbrtWad currently handles four cases. Can you reduce it to two cases?
3. z := sub(z, gt(lt(t, shr(1, p)), iszero(t))) Can you replace p/2 with a magical constant, just like sqrtWad did?