Eth2.0-specs: [BLS] Modular square root issue

[mratsim]

Problem

The BLS test vectors in https://github.com/ethereum/eth2.0-tests/pull/12 and https://github.com/ethereum/eth2.0-tests/pull/13 reuse the Trinity modular_squarerootimplementation at https://github.com/ethereum/trinity/blob/a1b0f058e7bc8e385c8dac3164c49098967fd5bb/eth2/_utils/bls.py#L63-L77

def modular_squareroot(value: FQ2) -> FQP:
    """
    ``modular_squareroot(x)`` returns the value ``y`` such that ``y**2 % q == x``,
    and None if this is not possible. In cases where there are two solutions,
    the value with higher imaginary component is favored;
    if both solutions have equal imaginary component the value with higher real
    component is favored.
    """
    candidate_squareroot = value ** ((FQ2_order + 8) // 16)
    check = candidate_squareroot ** 2 / value
    if check in eighth_roots_of_unity[::2]:
        x1 = candidate_squareroot / eighth_roots_of_unity[eighth_roots_of_unity.index(check) // 2]
        x2 = FQ2([-x1.coeffs[0], -x1.coeffs[1]])  # x2 = -x1
        return x1 if (x1.coeffs[1], x1.coeffs[0]) > (x2.coeffs[1], x2.coeffs[0]) else x2
    return None

Those are a copy paste of the current spec https://github.com/ethereum/eth2.0-specs/blob/928f9772fa51ed2e251b99c4795097da976c4bde/specs/bls_signature.md#modular_squareroot

Fq2_order = q ** 2 - 1
eighth_roots_of_unity = [Fq2([1,1]) ** ((Fq2_order * k) // 8) for k in range(8)]

def modular_squareroot(value: int) -> int:
    candidate_squareroot = value ** ((Fq2_order + 8) // 16)
    check = candidate_squareroot ** 2 / value
    if check in eighth_roots_of_unity[::2]:
        x1 = candidate_squareroot / eighth_roots_of_unity[eighth_roots_of_unity.index(check) // 2]
        x2 = -x1
        return x1 if (x1.coeffs[1].n, x1.coeffs[0].n) > (x2.coeffs[1].n, x2.coeffs[0].n) else x2
    return None

However, our experiments using Milagro modular square root showed that the spec does not match with established crypto library.

The main issue is that this comparison (x1.coeffs[1], x1.coeffs[0]) > (x2.coeffs[1], x2.coeffs[0]) is always true.

Analysis

This line x2 = FQ2([-x1.coeffs[0], -x1.coeffs[1]]) # x2 = -x1 is only negating the coefficient.
Given that x1 is always positive x2 will never be returned.

Probably Py-ECC should check the input coefficient and compute the modular inverse if one coefficient is negative.

Reference code: https://github.com/ethereum/py_ecc/blob/e9c42c76f1cc4f535dedec01b4c7ab2a59e46119/py_ecc/optimized_bls12_381/optimized_field_elements.py#L240-L246

class FQ2(FQP):
    def __init__(self, coeffs):
        self.coeffs = tuple(coeffs)
        self.modulus_coeffs = (1, 0)
        self.mc_tuples = [(0, 1)]
        self.degree = 2
        self.__class__.degree = 2

Note that Py-ECC does the correct initialization for FQ

class FQ(object):
    def __init__(self, n):
        if isinstance(n, self.__class__):
            self.n = n.n
        else:
            self.n = n % field_modulus
        assert isinstance(self.n, int)

In the specs, it should be clear that the - operations in x2 = FQ2([-x1.coeffs[0], -x1.coeffs[1]]) are doing modular inverse.

Test case

from py_ecc.optimized_bls12_381 import FQ, FQ2

x = FQ(-1)
y = FQ2([-1,-1])

print(x) # 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559786
print(y) # (-1, -1)

cc @cheatfate, @vbuterin, @pipermerriam, @hwwhww, @ChihChengLiang

Implementation question

Also I tried to find the algorithm used in the specs in papers but I couldn't find it in this extensive review of 10 modular square root over extension fields algorithms

[vbuterin]

This line x2 = FQ2([-x1.coeffs[0], -x1.coeffs[1]]) # x2 = -x1 is only negating the coefficient.
Given that x1 is always positive x2 will never be returned.

We're talking about numbers mod p here, so all of them are positive. The intention is that if either (i) x1.coeffs[1].n > p/2 or (ii) x1.coeffs[1].n == 0 and x1.coeffs[0].n > p/2 then return x1, otherwise return x2.

Or are you saying that the implementation doesn't actually store the coefficients as FQ elements so there's a bug and -x1.coeffs[0] actually is below zero? If the latter then it should be a simple fix.

I don't think it should matter what algorithm is used, as long as it chooses between the two square roots according to a consistent method, and picking the greater of root and -root seems like the most natural thing to do. But changing it to be consistent with whatever established crypto libraries prefer seems fine to me (unless established crypto libraries don't have a consistent formula and require something that requires switching to their very specific algorithm for modular square roots to know which root to use, in which case we should stick to picking the greater and wrap the crypto libraries to confirm).

Edit: I checked:

>>> x = py_ecc.bls12_381.FQ2([5,7]) ** 1
>>> y = -x
>>> x.coeffs
(5, 7)
>>> y.coeffs
(4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559782, 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559780)

So the implementation does work correctly. I'll modify the spec to make it clearer that it's a modular additive inverse.

[vbuterin]

https://github.com/ethereum/eth2.0-specs/pull/509

[mratsim]

I rechecked with FQ inputs instead of plain integers, it does store each coeff as FQ elements:

>>> import py_ecc
>>> x = py_ecc.bls12_381.FQ2([5,7]) ** 1
>>> y = py_ecc.bls12_381.FQ2([-x.coeffs[0], -x.coeffs[1]])
>>> y.coeffs
(4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559782, 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559780)
>>> 

I guess next step is to check the FQ2 modular square root implementation of Relic/Zcash/Milagro.

[mratsim]

Further investigation shows that it's the optimized_bls12_381 of Py-ECC that has an issue:

import py_ecc
x = py_ecc.bls12_381.FQ2([5,7])
y = py_ecc.bls12_381.FQ2([-x.coeffs[0], -x.coeffs[1]])
print(y.coeffs) # (4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559782, 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559780)


u = py_ecc.optimized_bls12_381.FQ2([5,7])
v = py_ecc.optimized_bls12_381.FQ2([-u.coeffs[0], -u.coeffs[1]])
print(v.coeffs) # (-5, -7)

Concretely, Trinity implementation has an issue due to this, it will always returns True or None.

Test case

from py_ecc.optimized_bls12_381 import (
    FQ, FQ2,
    field_modulus as q
)

FQ2_order = q ** 2 - 1
eighth_roots_of_unity = [
    FQ2([1, 1]) ** ((FQ2_order * k) // 8)
    for k in range(8)
]

# Trinity implementation,
# Modified to print intermediate result and if x1 > x2
def modular_squareroot(value: FQ2) -> bool:
    """
    ``modular_squareroot(x)`` returns the value ``y`` such that ``y**2 % q == x``,
    and None if this is not possible. In cases where there are two solutions,
    the value with higher imaginary component is favored;
    if both solutions have equal imaginary component the value with higher real
    component is favored.
    """
    candidate_squareroot = value ** ((FQ2_order + 8) // 16)
    check = candidate_squareroot ** 2 / value
    if check in eighth_roots_of_unity[::2]:
        x1 = candidate_squareroot / eighth_roots_of_unity[eighth_roots_of_unity.index(check) // 2]
        x2 = FQ2([-x1.coeffs[0], -x1.coeffs[1]])  # x2 = -x1
        print(f"x2: {x2}")
        # return x1 if (x1.coeffs[1], x1.coeffs[0]) > (x2.coeffs[1], x2.coeffs[0]) else x2
        return (x1.coeffs[1], x1.coeffs[0]) > (x2.coeffs[1], x2.coeffs[0])
    return None

for a in range(10):
    for b in range(10):
        x = FQ2([a, b]) ** 1
        print(f"({a}, {b}) - {x} - {FQ2([-x.coeffs[0], -x.coeffs[1]])}")
        y = modular_squareroot(x)
        print(f"x1 > x2 - {y}")

Prints

(0, 0) - (0, 0) - (0, 0)
x1 > x2 - None
(0, 1) - (0, 1) - (0, -1)
x2: (-2973677408986561043442465346520108879172042883009249989176415018091420807192182638567116318576472649347015917690530, -1028732146235106349975324479215795277384839936929757896155643118032610843298655225875571310552543014690878354869257)
x1 > x2 - True
(0, 2) - (0, 2) - (0, -2)
x2: (-1, -1)
x1 > x2 - True
(0, 3) - (0, 3) - (0, -3)
x2: (-3276987779355557806373750845889189077968193680945900882488430245112227261734582917308126182608378421639625449960584, -3276987779355557806373750845889189077968193680945900882488430245112227261734582917308126182608378421639625449960584)
x1 > x2 - True
(0, 4) - (0, 4) - (0, -4)
x2: (-2057464292470212699950648958431590554769679873859515792311286236065221686597310451751142621105086029381756709738514, -1944945262751454693467140867304313601787202946079492093020771900058809963893527412691545008023929634656137562821273)
x1 > x2 - True
(0, 5) - (0, 5) - (0, -5)
x2: (-3812379096704106451088749511229787994157686022564966047825536927476905165389684039502953357607714022383577266845174, -3812379096704106451088749511229787994157686022564966047825536927476905165389684039502953357607714022383577266845174)
x1 > x2 - True

...
...

As this is implementation specific and not spec related, I'm closing this.