diff --git a/include/ack/fe.hpp b/include/ack/fe.hpp index 41c5531..ec3ebb4 100644 --- a/include/ack/fe.hpp +++ b/include/ack/fe.hpp @@ -193,7 +193,7 @@ namespace ack { return underlying().mul( x ); } - /** + /** * Calculates square of this element. * @return R = this^2. */ @@ -202,6 +202,15 @@ namespace ack { return underlying().sqr(); } + /** + * Calculates square root of this element. + * @return R = sqrt(this). + */ + [[nodiscard]] Derived sqrt() const + { + return underlying().sqrt(); + } + /** * Calculates division of this element with another. * @param x - Element to divide. diff --git a/include/ack/fp.hpp b/include/ack/fp.hpp index 347ba15..240a597 100644 --- a/include/ack/fp.hpp +++ b/include/ack/fp.hpp @@ -92,6 +92,69 @@ namespace ack { { return ( a * b.modinv( p ) ) % p; } + + template + [[nodiscard]] inline bigint fp_sqrt(const bigint& n, const bigint& p) + { + // Check that n is indeed a square: (n | p) must be = 1 + // The Legendre symbol (n | p) denotes the value of n^(p-1)/2 (mod p). + if ( n.jacobi( p ) != 1 ) { + return 0; + } + + if (( p % 4 ) == 3 ) { + // sqrt(a) = a^((p + 1) / 4) (mod p) + return n.modexp( ( p + 1 ) / 4, p ); + } + + /* Fallback to general Tonelli-Shanks algorithm */ + + // Find q and s such that p - 1 = q2^s with q odd + auto q = p - 1; + std::size_t s = 0; + while ( q.is_even() ) { + s++; + q >>= 1; + } + + // Find non-square z such that (z | p) = -1 + bigint z = 2; + while ( z.jacobi( p ) != -1 ) { // TODO: Possible infinite loop + ++z; + } + + auto c = z.modexp( q, p ); + auto r = n.modexp(( q - 1 ) / 2, p ); + auto t = ( r.sqr() % p ) * n % p; + r = n * r % p; + + while ( t != 1 ) + { + std::size_t m = 0; + z = t; + do + { + m++; + z = z.sqr() % p; + if ( m == s ) { + return 0; + } + } + while ( z != 1 ); + + auto b = c; + for ( std::size_t i = 0; i < ( s - m - 1 ); i++ ) { + b = b.sqr() % p; + } + + c = b.sqr() % p; + s = m; + r = r * b % p; + t = t * c % p; + } + + return r; + } } /** Checks if finite field element a is valid. @@ -321,18 +384,40 @@ namespace ack { * Out of range values are not checked and will produce wrong results. * * @tparam BufferA - Type of the finite field element. - * @tparam BufferC - Type of the prime modulus. + * @tparam BufferB - Type of the prime modulus. * * @param a - Finite field element. * @param p - Prime modulus of the finite field. * @return a * a mod p. */ - template - [[nodiscard]] inline bigint fp_sqr(const bigint& a, const bigint& p) + template + [[nodiscard]] inline bigint fp_sqr(const bigint& a, const bigint& p) { return fp_mul( a, a, p ); } + /** + * Calculates the modular square root of an integer 'n' modulo odd prime 'p'. + * The function returns the first root, where r^2 == n (mod p), and not r^2 == -n (mod p). + * The caller is responsible for verifying if the returned root 'r' satisfies r^2 == n (mod p). + * If not, the caller should calculate r = p - r to get the correct root. + * + * @note This function has an optimization for the case where 'p' == 3 (mod 4). + * In this case, the square root is efficiently computed. + * If 'p' is not congruent to 3 modulo 4, the function uses a general + * Tonelli-Shanks algorithm to find the modular square root. + * + * @param n - The integer for which the modular square root is calculated. + * @param p - The odd prime modulus. + * @return The first modular square root of 'n' modulo 'p', or 0 if no square root exists. + * + */ + template + [[nodiscard]] inline bigint fp_sqrt(const bigint& n, const bigint& p) + { + return detail::fp_sqrt( n, p ); + } + /** * Calculates division of prime finite field elements a and b. * The result is a / b mod p. Internally it's calculated as (a * b^-1) mod p. @@ -379,7 +464,7 @@ namespace ack { } /** - * Constructs a one finite field element. + * Constructs a finite field element of value 1. * @note This instance doesn't have a valid modulus, * so it can't be used for any operations. * The returned instance is invalid and can be used only for comparison. @@ -563,7 +648,7 @@ namespace ack { } /** - * Calculate modular addition of this finite field element and another one. + * Calculates modular addition of this finite field element and another one. * @param x - Finite field element to add. * @return (this + x) % modulus. */ @@ -573,7 +658,7 @@ namespace ack { } /** - * Calculate modular addition of this finite field element and big integer. + * Calculates modular addition of this finite field element and big integer. * @tparam BufferT - Buffer type of big integer. * @param x - Big integer to add. * @return (this + x) % modulus. @@ -585,7 +670,7 @@ namespace ack { } /** - * Calculate modular addition of this finite field element and integer. + * Calculates modular addition of this finite field element and integer. * @tparam IntU - Small integer type. * @param x - Integer to add. * @return (this + x) % modulus. @@ -597,7 +682,7 @@ namespace ack { } /** - * Calculate modular subtraction of this finite field element and another one. + * Calculates modular subtraction of this finite field element and another one. * @param x - Finite field element to subtract. * @return (this - x) % modulus. */ @@ -607,7 +692,7 @@ namespace ack { } /** - * Calculate modular subtraction of this finite field element and big integer. + * Calculates modular subtraction of this finite field element and big integer. * @tparam BufferT - Buffer type of big integer. * @param x - Big integer to subtract. * @return (this - x) % modulus. @@ -619,7 +704,7 @@ namespace ack { } /** - * Calculate modular subtraction of this finite field element and integer. + * Calculates modular subtraction of this finite field element and integer. * @tparam IntU - Small integer type. * @param x - Integer to subtract. * @return (this - x) % modulus. @@ -631,7 +716,7 @@ namespace ack { } /** - * Calculate modular multiplication of this finite field element and another one. + * Calculates modular multiplication of this finite field element and another one. * @param x - Finite field element to multiply. * @return (this * x) % modulus. */ @@ -641,7 +726,7 @@ namespace ack { } /** - * Calculate modular multiplication of this finite field element and big integer. + * Calculates modular multiplication of this finite field element and big integer. * @tparam BufferT - Buffer type of big integer. * @param x - Big integer to multiply. * @return (this * x) % modulus. @@ -659,7 +744,7 @@ namespace ack { } /** - * Calculate modular multiplication of this finite field element and integer. + * Calculates modular multiplication of this finite field element and integer. * @tparam IntU - Small integer type. * @param x - Integer to multiply. * @return (this * x) % modulus. @@ -677,7 +762,7 @@ namespace ack { } /** - * Calculate modular square of this finite field element. + * Calculates modular square of this finite field element. * @return (this^2) % modulus. */ [[nodiscard]] fp_element sqr() const @@ -689,7 +774,23 @@ namespace ack { } /** - * Calculate modular division of this finite field element and another one. + * Calculates the modular square root of this finite field element. + * The function returns the first root, where r^2 == n (mod p), and not r^2 == -n (mod p). + * The caller is responsible for verifying if the returned root 'r' satisfies r^2 == n (mod p). + * If not, the caller should calculate r = p - r to get the correct root. + * + * @return The first modular square root of this mod modulus, or 0 if no square root exists. + */ + [[nodiscard]] fp_element sqrt() const + { + if ( v_.is_one() ) { + return *this; + } + return fp_element( fp_sqrt( v_, *pm_ ), *pm_ ); + } + + /** + * Calculates modular division of this finite field element and another one. * @param x - Finite field element to divide. * @return (this / x) % modulus. */ @@ -699,7 +800,7 @@ namespace ack { } /** - * Calculate modular division of this finite field element and big integer. + * Calculates modular division of this finite field element and big integer. * @tparam BufferT - Buffer type of big integer. * @param x - Big integer to divide. * @return (this / x) % modulus. diff --git a/tests/include/ack/tests/fp_test.hpp b/tests/include/ack/tests/fp_test.hpp index b2bc39d..2dfa9b3 100644 --- a/tests/include/ack/tests/fp_test.hpp +++ b/tests/include/ack/tests/fp_test.hpp @@ -383,6 +383,175 @@ namespace ack::tests { REQUIRE_EQUAL( fp_mul( bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11DCB29F14BACCFF196CE3F0AD2" ), -bn_t( "6AA11DCB29F14BACCFF196CE3F0AD2" ), bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), -bn_t( "1A945B1E1A80601352171738E860E7922A6EBAAA6") ) } + // fp_sqr + { + REQUIRE_EQUAL( fp_sqr( bn_t( 0 ), bn_t( 10 ) ), 0 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 1 ), bn_t( 10 ) ), 1 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 2 ), bn_t( 10 ) ), 4 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 3 ), bn_t( 10 ) ), 9 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 4 ), bn_t( 10 ) ), 6 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 5 ), bn_t( 10 ) ), 5 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 6 ), bn_t( 10 ) ), 6 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 7 ), bn_t( 10 ) ), 9 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 8 ), bn_t( 10 ) ), 4 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 9 ), bn_t( 10 ) ), 1 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 10 ), bn_t( 10 ) ), 0 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 11 ), bn_t( 10 ) ), 1 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 12 ), bn_t( 10 ) ), 4 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 13 ), bn_t( 10 ) ), 9 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 14 ), bn_t( 10 ) ), 6 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 15 ), bn_t( 10 ) ), 5 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 16 ), bn_t( 10 ) ), 6 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 17 ), bn_t( 10 ) ), 9 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 18 ), bn_t( 10 ) ), 4 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 19 ), bn_t( 10 ) ), 1 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 20 ), bn_t( 10 ) ), 0 ) + REQUIRE_EQUAL( fp_sqr( bn_t( 21 ), bn_t( 10 ) ), 1 ) + + REQUIRE_EQUAL( fp_sqr( bn_t( -1 ), bn_t( 10 ) ), 1 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -2 ), bn_t( 10 ) ), 4 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -3 ), bn_t( 10 ) ), 9 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -4 ), bn_t( 10 ) ), 6 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -5 ), bn_t( 10 ) ), 5 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -6 ), bn_t( 10 ) ), 6 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -7 ), bn_t( 10 ) ), 9 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -8 ), bn_t( 10 ) ), 4 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -9 ), bn_t( 10 ) ), 1 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -10 ), bn_t( 10 ) ), 0 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -11 ), bn_t( 10 ) ), 1 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -12 ), bn_t( 10 ) ), 4 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -13 ), bn_t( 10 ) ), 9 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -14 ), bn_t( 10 ) ), 6 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -15 ), bn_t( 10 ) ), 5 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -16 ), bn_t( 10 ) ), 6 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -17 ), bn_t( 10 ) ), 9 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -18 ), bn_t( 10 ) ), 4 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -19 ), bn_t( 10 ) ), 1 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -20 ), bn_t( 10 ) ), 0 ) + REQUIRE_EQUAL( fp_sqr( bn_t( -21 ), bn_t( 10 ) ), 1 ) + + // Big numbers + REQUIRE_EQUAL( fp_sqr( bn_t( "6A9E8EA23CB63C228F31" ), bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), bn_t( "2C67B052BA023CE9E8FDB4FD5866C975511BC761") ) + REQUIRE_EQUAL( fp_sqr( bn_t( "6AA11DCB29F14BACCFF1" ), bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), bn_t( "2C69D212BBBE3CE8824379FD4E8F6C414ABFA0E1") ) + + REQUIRE_EQUAL( fp_sqr( bn_t( "6A9E8EA23CB63C228F31AD9268B4097F1" ), bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), bn_t( "39FEFE5D80C16D765627970811925B885C5EE895B") ) + REQUIRE_EQUAL( fp_sqr( bn_t( "6AA11DCB29F14BACCFF196CE3F0AD2" ) , bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), bn_t( "61A073DF655912F0BD29C66FCD1FA82267512A555") ) + + REQUIRE_EQUAL( fp_sqr( bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11DCB29F14BACCFF196CE3F0AD2" ), bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), bn_t( "63221D108A0978A6D626C73774B45F2CBFF75E686") ) + REQUIRE_EQUAL( fp_sqr( bn_t( "6AA11DCB29F14BACCFF196CE3F0AD2" ), bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), bn_t( "61A073DF655912F0BD29C66FCD1FA82267512A555") ) + + REQUIRE_EQUAL( fp_sqr( -bn_t( "6A9E8EA23CB63C228F31AD9268B4097F1" ), bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), bn_t( "39FEFE5D80C16D765627970811925B885C5EE895B") ) + REQUIRE_EQUAL( fp_sqr( -bn_t( "6AA11DCB29F14BACCFF196CE3F0AD2" ) , bn_t( "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11" ) ), bn_t( "61A073DF655912F0BD29C66FCD1FA82267512A555") ) + } + + // fp_sqrt + { + const auto test_sqrt = [](const bn_t& a, const bn_t& p) { + REQUIRE_EQUAL( fp_sqrt( fp_sqr( a, p ), p ), a ) + }; + + const auto test_sqrt_neg = [](const bn_t& a, const bn_t& p) { + REQUIRE_EQUAL( p - fp_sqrt( fp_sqr( a, p ), p ), a ) + }; + + // Test cases where p % 4 != 3 + // i.e.: general Tonelli–Shanks algorithm + { + bn_t p = 17; + test_sqrt( 1, p ); + test_sqrt( 2, p ); + test_sqrt( 4, p ); + test_sqrt( 6, p ); + test_sqrt( 7, p ); + test_sqrt( 8, p ); + test_sqrt( 12, p ); + test_sqrt( 14, p ); + + // special cases (negated roots, invalid roots) + test_sqrt( 0, p ); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 3 ), p ), p ), 14 ) + test_sqrt_neg( 3, p); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 5 ), p ), p ), 12 ) + test_sqrt_neg( 5, p); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 9 ), p ), p ), 8 ) + test_sqrt_neg( 9, p); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 10 ), p ), p ), 7 ) + test_sqrt_neg( 10, p); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 11 ), p ), p ), 6 ) + test_sqrt_neg( 11, p); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 13 ), p ), p ), 4 ) + test_sqrt_neg( 13, p); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 15 ), p ), p ), 2 ) + test_sqrt_neg( 15, p); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 16 ), p ), p ), 1 ) + test_sqrt_neg( 16, p); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 17 ), p ), p ), 0 ) + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 18 ), p ), p ), 1 ) + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 19 ), p ), p ), 2 ) + + // Additional tests + test_sqrt( 7, 13 ); + test_sqrt( 37, 101 ); + test_sqrt( 651, 1009 ); + test_sqrt( 174, 1009 ); + REQUIRE_EQUAL( fp_sqrt( bn_t( 1032 ) , bn_t( 1009 )), 0 ) + + test_sqrt( 378633312, 1000000009 ); + test_sqrt( 378633312, 1000000009 ); + test_sqrt( "15f73ba6b33d469cb81975ea86c422a33be005070c", "446c3b15f9926687d2c40534fdb564000000000241" ); + } + + // Test cases where p % 4 == 3 + { + bn_t p = 11; + test_sqrt( 1, 11 ); + test_sqrt( 3, 11 ); + test_sqrt( 4, 11 ); + test_sqrt( 5, 11 ); + test_sqrt( 9, 11 ); + + // special cases (negated roots, invalid roots) + test_sqrt( 0, 11 ); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 2 ), p ), p ), 9 ) + test_sqrt_neg( 2, p ); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 6 ), p ), p ), 5 ) + test_sqrt_neg( 6, p ); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 7 ), p ), p ), 4 ) + test_sqrt_neg( 7, p ); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 8 ), p ), p ), 3 ) + test_sqrt_neg( 8, p ); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 10 ), p ), p ), 1 ) + test_sqrt_neg( 10, p ); + + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 11 ), p ), p ), 0 ) + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 12 ), p ), p ), 1 ) + REQUIRE_EQUAL( fp_sqrt( fp_sqr( bn_t( 13 ), p ), p ), 9 ) + + // Additional tests + test_sqrt( 791399408049ULL, 1000000000039ULL ); + + // secp256r1 G.y and p + test_sqrt( "4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5", "ffffffff00000001000000000000000000000000ffffffffffffffffffffffff" ); + + // secp256k1 G.y and p + test_sqrt( "483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f" ); + } + } + // fp_div tests { REQUIRE_EQUAL( fp_div( bn_t( 0 ), bn_t( 1 ), bn_t( 10 ) ), 0 ) @@ -2270,7 +2439,6 @@ namespace ack::tests { fe *= 65536; REQUIRE_EQUAL( fe, 8 ) - REQUIRE_EQUAL( fpe_t( -1, modulus ).mul( fpe_t( -1, modulus ) ), 1 ) REQUIRE_EQUAL( fpe_t( -1, modulus ).mul( bn_t( -1 ) ) , 1 ) REQUIRE_EQUAL( fpe_t( -1, modulus ).mul( -1 ) , 1 ) @@ -2767,6 +2935,181 @@ namespace ack::tests { REQUIRE_EQUAL( fe, -bn_t( "1A945B1E1A80601352171738E860E7922A6EBAAA6") ) } + // fp_sqr + { + bn_t modulus = 10; + + REQUIRE_EQUAL( fpe_t( bn_t(0), modulus ).sqr(), 0 ); + REQUIRE_EQUAL( fpe_t( bn_t(1), modulus ).sqr(), 1 ); + REQUIRE_EQUAL( fpe_t( bn_t(2), modulus ).sqr(), 4 ); + REQUIRE_EQUAL( fpe_t( bn_t(3), modulus ).sqr(), 9 ); + REQUIRE_EQUAL( fpe_t( bn_t(4), modulus ).sqr(), 6 ); + REQUIRE_EQUAL( fpe_t( bn_t(5), modulus ).sqr(), 5 ); + REQUIRE_EQUAL( fpe_t( bn_t(6), modulus ).sqr(), 6 ); + REQUIRE_EQUAL( fpe_t( bn_t(7), modulus ).sqr(), 9 ); + REQUIRE_EQUAL( fpe_t( bn_t(8), modulus ).sqr(), 4 ); + REQUIRE_EQUAL( fpe_t( bn_t(9), modulus ).sqr(), 1 ); + REQUIRE_EQUAL( fpe_t( bn_t(10), modulus ).sqr(), 0 ); + REQUIRE_EQUAL( fpe_t( bn_t(11), modulus ).sqr(), 1 ); + REQUIRE_EQUAL( fpe_t( bn_t(12), modulus ).sqr(), 4 ); + REQUIRE_EQUAL( fpe_t( bn_t(13), modulus ).sqr(), 9 ); + REQUIRE_EQUAL( fpe_t( bn_t(14), modulus ).sqr(), 6 ); + REQUIRE_EQUAL( fpe_t( bn_t(15), modulus ).sqr(), 5 ); + REQUIRE_EQUAL( fpe_t( bn_t(16), modulus ).sqr(), 6 ); + REQUIRE_EQUAL( fpe_t( bn_t(17), modulus ).sqr(), 9 ); + REQUIRE_EQUAL( fpe_t( bn_t(18), modulus ).sqr(), 4 ); + REQUIRE_EQUAL( fpe_t( bn_t(19), modulus ).sqr(), 1 ); + REQUIRE_EQUAL( fpe_t( bn_t(20), modulus ).sqr(), 0 ); + REQUIRE_EQUAL( fpe_t( bn_t(21), modulus ).sqr(), 1 ); + + REQUIRE_EQUAL( fpe_t( bn_t(-1), modulus ).sqr(), 1 ); + REQUIRE_EQUAL( fpe_t( bn_t(-2), modulus ).sqr(), 4 ); + REQUIRE_EQUAL( fpe_t( bn_t(-3), modulus ).sqr(), 9 ); + REQUIRE_EQUAL( fpe_t( bn_t(-4), modulus ).sqr(), 6 ); + REQUIRE_EQUAL( fpe_t( bn_t(-5), modulus ).sqr(), 5 ); + REQUIRE_EQUAL( fpe_t( bn_t(-6), modulus ).sqr(), 6 ); + REQUIRE_EQUAL( fpe_t( bn_t(-7), modulus ).sqr(), 9 ); + REQUIRE_EQUAL( fpe_t( bn_t(-8), modulus ).sqr(), 4 ); + REQUIRE_EQUAL( fpe_t( bn_t(-9), modulus ).sqr(), 1 ); + REQUIRE_EQUAL( fpe_t( bn_t(-10), modulus ).sqr(), 0 ); + REQUIRE_EQUAL( fpe_t( bn_t(-11), modulus ).sqr(), 1 ); + REQUIRE_EQUAL( fpe_t( bn_t(-12), modulus ).sqr(), 4 ); + REQUIRE_EQUAL( fpe_t( bn_t(-13), modulus ).sqr(), 9 ); + REQUIRE_EQUAL( fpe_t( bn_t(-14), modulus ).sqr(), 6 ); + REQUIRE_EQUAL( fpe_t( bn_t(-15), modulus ).sqr(), 5 ); + REQUIRE_EQUAL( fpe_t( bn_t(-16), modulus ).sqr(), 6 ); + REQUIRE_EQUAL( fpe_t( bn_t(-17), modulus ).sqr(), 9 ); + REQUIRE_EQUAL( fpe_t( bn_t(-18), modulus ).sqr(), 4 ); + REQUIRE_EQUAL( fpe_t( bn_t(-19), modulus ).sqr(), 1 ); + REQUIRE_EQUAL( fpe_t( bn_t(-20), modulus ).sqr(), 0 ); + REQUIRE_EQUAL( fpe_t( bn_t(-21), modulus ).sqr(), 1 ); + + // Big numbers + modulus = "6A9E8EA23CB63C228F31AD9268B4097F156D6AA11"; + + REQUIRE_EQUAL(fpe_t( bn_t("6A9E8EA23CB63C228F31"), modulus ).sqr(), bn_t("2C67B052BA023CE9E8FDB4FD5866C975511BC761")) + REQUIRE_EQUAL(fpe_t( bn_t("6AA11DCB29F14BACCFF1"), modulus ).sqr(), bn_t("2C69D212BBBE3CE8824379FD4E8F6C414ABFA0E1")) + + REQUIRE_EQUAL(fpe_t( bn_t("6A9E8EA23CB63C228F31AD9268B4097F1"), modulus ).sqr(), bn_t("39FEFE5D80C16D765627970811925B885C5EE895B")) + REQUIRE_EQUAL(fpe_t( bn_t("6AA11DCB29F14BACCFF196CE3F0AD2"), modulus ).sqr(), bn_t("61A073DF655912F0BD29C66FCD1FA82267512A555")) + + REQUIRE_EQUAL(fpe_t( bn_t("6A9E8EA23CB63C228F31AD9268B4097F156D6AA11DCB29F14BACCFF196CE3F0AD2" ), modulus).sqr(), bn_t("63221D108A0978A6D626C73774B45F2CBFF75E686")) + REQUIRE_EQUAL(fpe_t( bn_t("6AA11DCB29F14BACCFF196CE3F0AD2"), modulus ).sqr(), bn_t("61A073DF655912F0BD29C66FCD1FA82267512A555")) + + REQUIRE_EQUAL(fpe_t( -bn_t("6A9E8EA23CB63C228F31AD9268B4097F1"), modulus ).sqr(), bn_t("39FEFE5D80C16D765627970811925B885C5EE895B")) + REQUIRE_EQUAL(fpe_t( -bn_t("6AA11DCB29F14BACCFF196CE3F0AD2"), modulus ).sqr(), bn_t("61A073DF655912F0BD29C66FCD1FA82267512A555")) + } + + // fp_sqrt + { + const auto test_sqrt = [](const bn_t& a, const bn_t& p) { + REQUIRE_EQUAL( fpe_t( a, p ).sqr().sqrt(), a ) + }; + + const auto test_sqrt_neg = [](const bn_t& a, const bn_t& p) { + REQUIRE_EQUAL( -fpe_t( a, p ).sqr().sqrt(), a ) + }; + + // Test cases where p % 4 != 3 + // i.e.: general Tonelli–Shanks algorithm + { + bn_t p = 17; + test_sqrt( 1, p ); + test_sqrt( 2, p ); + test_sqrt( 4, p ); + test_sqrt( 6, p ); + test_sqrt( 7, p ); + test_sqrt( 8, p ); + test_sqrt( 12, p ); + test_sqrt( 14, p ); + + // special cases (negated roots, invalid roots) + test_sqrt( 0, p ); + + REQUIRE_EQUAL(fpe_t(bn_t(3), p).sqr().sqrt(), 14); + test_sqrt_neg(3, p); + + REQUIRE_EQUAL(fpe_t(bn_t(5), p).sqr().sqrt(), 12); + test_sqrt_neg(5, p); + + REQUIRE_EQUAL(fpe_t(bn_t(9), p).sqr().sqrt(), 8); + test_sqrt_neg(9, p); + + REQUIRE_EQUAL(fpe_t(bn_t(10), p).sqr().sqrt(), 7); + test_sqrt_neg(10, p); + + REQUIRE_EQUAL(fpe_t(bn_t(11), p).sqr().sqrt(), 6); + test_sqrt_neg(11, p); + + REQUIRE_EQUAL(fpe_t(bn_t(13), p).sqr().sqrt(), 4); + test_sqrt_neg(13, p); + + REQUIRE_EQUAL(fpe_t(bn_t(15), p).sqr().sqrt(), 2); + test_sqrt_neg(15, p); + + REQUIRE_EQUAL(fpe_t(bn_t(16), p).sqr().sqrt(), 1); + test_sqrt_neg(16, p); + + REQUIRE_EQUAL(fpe_t(bn_t(17), p).sqr().sqrt(), 0); + REQUIRE_EQUAL(fpe_t(bn_t(18), p).sqr().sqrt(), 1); + REQUIRE_EQUAL(fpe_t(bn_t(19), p).sqr().sqrt(), 2); + + // Additional tests + test_sqrt( 7, 13 ); + test_sqrt( 37, 101 ); + test_sqrt( 651, 1009 ); + test_sqrt( 174, 1009 ); + + p = bn_t( 1009 ); + REQUIRE_EQUAL( fpe_t( bn_t( 1032 ) , p).sqrt(), 0 ) + + test_sqrt( 378633312, 1000000009 ); + test_sqrt( 378633312, 1000000009 ); + test_sqrt( "15f73ba6b33d469cb81975ea86c422a33be005070c", "446c3b15f9926687d2c40534fdb564000000000241" ); + } + + // Test cases where p % 4 == 3 + { + bn_t p = 11; + test_sqrt( 1, 11 ); + test_sqrt( 3, 11 ); + test_sqrt( 4, 11 ); + test_sqrt( 5, 11 ); + test_sqrt( 9, 11 ); + + // special cases (negated roots, invalid roots) + test_sqrt( 0, 11 ); + + REQUIRE_EQUAL(fpe_t(bn_t(2), p).sqr().sqrt(), 9); + test_sqrt_neg(2, p); + + REQUIRE_EQUAL(fpe_t(bn_t(6), p).sqr().sqrt(), 5); + test_sqrt_neg(6, p); + + REQUIRE_EQUAL(fpe_t(bn_t(7), p).sqr().sqrt(), 4); + test_sqrt_neg(7, p); + + REQUIRE_EQUAL(fpe_t(bn_t(8), p).sqr().sqrt(), 3); + test_sqrt_neg(8, p); + + REQUIRE_EQUAL(fpe_t(bn_t(10), p).sqr().sqrt(), 1); + test_sqrt_neg(10, p); + + REQUIRE_EQUAL(fpe_t(bn_t(11), p).sqr().sqrt(), 0); + REQUIRE_EQUAL(fpe_t(bn_t(12), p).sqr().sqrt(), 1); + REQUIRE_EQUAL(fpe_t(bn_t(13), p).sqr().sqrt(), 9); + + // Additional tests + test_sqrt( 791399408049ULL, 1000000000039ULL ); + + // secp256r1 G.y and p + test_sqrt( "4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5", "ffffffff00000001000000000000000000000000ffffffffffffffffffffffff" ); + + // secp256k1 G.y and p + test_sqrt( "483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f" ); + } + } + // Division tests { bn_t modulus = 10;