Vyper: VIP: Square root function

[jacqueswww]

What's your issue about?

We don't have sqrt function, this is useful to have for bonding curves.

How can it be fixed?

Embed this code:
```python

def test_sqrt(get_contract):

code = """

@public
def sqrt256(x: uint256) -> uint256:
z: uint256 = (x + 1) / 2
y: uint256 = x

for i in range(256):
    if (z > y):
        break
    y = z
    z = (x / z + z) / 2
return y

@public
def sqrt_decimal(x: decimal) -> decimal:
assert x >= 0.0
if x == 0.0:
return 0.0
z: decimal = (x + 1.0) / 2.0
y: decimal = x

for i in range(256):
    if (z == y):
        break
    y = z
    z = (x / z + z) / 2.0
return y

"""

c = get_contract(code)

assert c.sqrt256(2) == 1
# math.sqrt(2000) == 44.721359549995796
assert c.sqrt_decimal(2000) == Decimal('44.7213595499')

[fubuloubu]

Should this be a VIP? It's a small function, but we can discuss how to write it and use it.

[charles-cooper]

Neat. I think a more efficient implementation would be in terms of (MOD)EXP and a vyper-defined log2 function. I believe log2 can be implemented in O(log2(256)) steps instead of O(256); plus, log2 is more generally useful. Another thing might be to let the user control the number of iterations of the approximation.

[fubuloubu]

:+1: on log2. Would so love to have that. I was working on an integer version a bit ago.

[charles-cooper]

I think log2 should be something like

@public
def log2(x: uint256) -> uint256 :
  ret = 0
  for i in range(0, 8) :
    p = 2**(8 - i)
    mask = 2**p - 1
    cond = iszero(x & mask) # encode branchless 'if'
    x >>= cond # who needs JUMP
    ret += cond * p
  return ret

[jacqueswww]

Sure but that's log2, I just want a sqrt for now ;) I have not found a more efficient version than above. Or am I missing something?

[charles-cooper]

@jacqueswww sqrt(x) can be represented as exp_b(log_b(x)/2)

$ python <<EOF
import math
print(math.exp(math.log(2000)/2))
print(math.sqrt(2000))
EOF

44.72135955
44.72135955

[fubuloubu]

@charlee-cooper would that be the most gas-efficient algo for calculating square root though?

[jacqueswww]

exp gas is defined as: (exp == 0) ? 10 : (10 + 10 * (1 + log256(exp))) hmmmmm . I think that could be small enough

[jacqueswww]

We need a nice gas comparison too, also we need to use the same tool to compile & compare calls to solidity.

[charles-cooper]

Also no matter the implementation it would be nice to also be able to return the trailing bits. Another use case for bigmath!

[charles-cooper]

@fubuloubu I think it's the most gas efficient. Each loop iteration would have to be 32x less efficient in order for it to be slower. Fwiw, the fast inverse square root (https://en.m.wikipedia.org/wiki/Fast_inverse_square_root) is fast because it uses a log approximation + 1 Newton step instead of a full Newton search. While we don't have access to the same trick they use for log, instead we can just calculate it.

[jacqueswww]

@charles-cooper sounds good, I would suggest we whip up a log2 PoC (just as a vyper functions?) and we compare?

[charles-cooper]

@jacqueswww @fubuloubu here's a PoC that I got to compile. The problem is a bunch of bits get truncated because it's all integer log/exp. I will think more about how to get those bits ;)

from decimal import Decimal
def test_sqrt(get_contract):
    code = """
@public
def sqrt_decimal(x: decimal) -> decimal:
    # sqrt(x) = exp( log(x * 10e10) / 2 ) / 10e5
    DECIMAL_DEN: decimal = convert(10**10, decimal)
    SQRT_DEN: decimal = convert(10**5, decimal)
    assert x >= 0.0
    if x == 0.0:
        return 0.0
        #return 0

    # calculate log2
    log2_intermediate: uint256 = 0
    x_i: uint256 = convert(x * DECIMAL_DEN, uint256) # convert to int for bit stuff

    for i in range(8):
        p: int128 = convert(shift(1, 7 - i), int128)
        mask: uint256 = shift(1, p) - 1
        y_i: uint256 = shift(x_i, -p)
        cond: uint256 = bitwise_and(y_i, mask)
        if cond > 0 : # real code would have no jump
          x_i = y_i
          log2_intermediate += convert(p, uint256)

    exp2: uint256 = 2 ** (log2_intermediate / 2)
    return convert(exp2, decimal) / SQRT_DEN
    """

    c = get_contract(code)

    assert c.sqrt_decimal(Decimal('7.70636789')) == Decimal('2.7760345') # sqrt

[charles-cooper]

Maybe the catch is you can only get 8 bits of fidelity? :thinking:

[fubuloubu]

"fidelity" is an awkward thing to explain to users. truncate the result or provide the same fidelity that the output type expects is better IMO.

[jacqueswww]

Meeting discussion:

• Sqrt will be used for log2 implementation.
• Use Babylonian method
• Only Decimal version for the time being. Integer can always be added later.

Decision yet to be made:

• None, let's do this!

[jacqueswww]

@charles-cooper @fubuloubu Check https://github.com/andytudhope/Recollections/issues/4#issuecomment-472098452.

The second algorithm is very similar to what @charles-cooper was referring to when we had the call (adjusting precision). I could definitely see having sqrt_hfp function builtin to be useful?

[fubuloubu]

I feel like sqrt should only take a decimal and produce a decimal, so we guarantee 10 decimal places of precision.

[fubuloubu]

Having a "special" version where sqrt(2) produces a _larger_ number than the operand is highly non-intuitive I think.

[jacqueswww]

I completely agree, it should be normalised to using a uint256 and outputting a uint256.
Something like:

def sqrt_hfp(value: uint256, decimal_places: value: uint256, floating_accuracy: int128) -> uint256

Perhaps ?

[fubuloubu]

Noooooo

def sqrt(value: decimal) -> decimal

[jacqueswww]

Haha to elaborate further: I fully agree on the default on being def sqrt(value: decimal) -> decimal and will implement it accordingly. But I don't see a reason not to have a second type - if it's required by contract devs using ERC20 tokens (which is a big part of contract out there, and will remain so).

[fubuloubu]

Why not convert? It'll be more explicit about this.

[jacqueswww]

@fubuloubu because you'll loose the the other 8 digits of precision? Let's try out the one method first, and then see from there.

[fubuloubu]

Oh, I think I see what you're saying. I think a better solution for that would be a different kind of conversion that takes the integer value and interprets it as having X decimal places for this kind of functionality. (The current convert truncates the integer when converting back)

I think this would be better in an extension to the current decimal functionality by allowing to pick an arbitrary amount of decimal places (for example, most tokens are 18 decimal places not 10)

[charles-cooper]

Here's an idea for the API, a user could specify the number of bits of precision (i.e. from the left) with a positive number of places, and the number of 'places' from the right with a negative number. E.g. 3 decimal places past 0 could be represented by places=-7.

A couple more late-night thoughts, the user might want to control the rounding mode. I'm not sure how complicated this would be to implement/how much extra gas it would cost. Also it might be easier implementation-wise to specify the number of bits rather than the number of decimal places (and also easier for the user to reason about when the input is uint256).

def sqrt_hfp(value: uint256, decimal_places: value: uint256, floating_accuracy: int128) -> uint256

Oh, I think I see what you're saying. I think a better solution for that would be a different kind of conversion that takes the integer value and interprets it as having X decimal places for this kind of functionality. (The current convert truncates the integer when converting back)

We should start having a conversation about a floating point standard ;). It keeps popping up oh so innocuously and it would be better to have a standard (and not necessarily opinionated) framework for talking about these kinds of issues and the tradeoffs between different approaches.

[fubuloubu]

We should start having a conversation about a floating point standard

I didn't saying anything about floating point numbers lol. I definitely think that we should not have floating point numbers in the EVM due to their ambiguity (you can represent a given number many different ways), lossy-ness, and unique behaviors (IEEE NaN and Inf, anyone?).

I do think the fixed-point standardization of the ABI and mathematical functions is extremely underdeveloped. Fixed-point math is important for financial applications as lossy math fundamentally leads to financial harm. We do need to give people tools to understand the effects that different algorithms have on the precision of the results, but hopefully there is only so much math people are doing that 10 decimal pts of precision is enough for most cases.

There are far more corner cases and undefined behavior in using floats, and as soon as we give access to these data types a lot more people would use them without understanding the problems with their use. Everything will be all fine and good until the first person who runs into a significant inconsistency error caused by floating point math which artificially inflates the supply of some token, or reduces it through "theft" of people's holdings.

Fixed point is a much safer approach for all of our users.

[pipermerriam]

My intuition is that any constant precision here is going to be inflexible enough to be a problem.

I'm not intimate with vyper's decimal facilities, but implementing something like the context API from the python decimal library might be worth considering. This could be as simple as having an optional context object that any of functions that operate on decimals take, and later implementing some of the nicer contextmanager style APIs.

https://docs.python.org/3.5/library/decimal.html#decimal.Context

[fubuloubu]

Should a precision failure cause a revert?

[pipermerriam]

Should a precision failure cause a revert?

There probably isn't a global right answer to that, hence the concept of decimal contexts to allow this to be controlled.

Admittedly there is a counter argument here which is that APIs that allow lots of configuration can also be surprising or lead to errors if you don't understand the configuration rules.

Assuming we like the idea of contexts, then we potentially expose a way for a contract to define it's default context as well as something like the decimal.localcontext() API which allows for contextmanager style APIs for isolated context rules like revert if there is a precision error.

[charles-cooper]

@fubuloubu I'm not suggesting that the EVM should implement floating point operations or even that we should adopt any particular specification, but rather that we should have a framework where we can talk about these kinds of issues, the tradeoffs and ways to work with and reason about those tradeoffs. (And there are tradeoffs, e.g. fixed point arithmetic isn't entirely lossless either, which is why #1086 was proposed). It is much better IMO to be able to recognize and say, 'this problem is a specific case of this general class of problem, and here are the standard ways in which we approach this problem', rather than having an ad hoc conversation and approach every time a specific case comes up. But that's just my 2c.

[fubuloubu]

Yup, I agree. Maybe we need that framework before this PR can progress?

[charles-cooper]

Maybe we need that framework before this PR can progress?

I think the decimal case is fairly clear-cut, but yeah, maybe for the variable precision version we do.

[fubuloubu]

I agree. The PR for this should implement a decimal-only function for now

[davesque]

@charles-cooper @fubuloubu Can we clarify what we mean with the term "framework"? I'm a bit fuzzy on what's being suggested there. Is that just what @pipermerriam was talking about with having a precision "context" API a la python's standard decimal library? Or were you thinking about developing some specific process for proposing and working out the details of these particular kinds of features? If so, that feels like too much scope creep to me.

[fubuloubu]

Very much like scope creep :smile: (for either suggestion)

It could be important later on, if we see demand for more advanced mathematical functionality.

I'd say any of these suggestions (other than adding a decimal-only sqrt function) are probably way out of scope for our v1 release.