From 4ed8e9e596d4508330639d01cc2f678ea11b7b64 Mon Sep 17 00:00:00 2001 From: Manuel Zahariev Date: Thu, 20 Mar 2025 12:58:18 -0700 Subject: [PATCH] LibCrypto: Improve precision of Crypto::BigFraction::to_double() Before: - FIXME: very naive implementation - was preventing passing some Temporal tests - https://github.com/tc39/test262 - https://github.com/LadybirdBrowser/libjs-test262 Bonus: Unrelated formatting change (Line 249) that unblocks the CI lint check. --- .../LibCrypto/BigFraction/BigFraction.cpp | 59 +++++++++++++++++-- Tests/LibCrypto/TestBigFraction.cpp | 37 ++++++++++++ 2 files changed, 91 insertions(+), 5 deletions(-) diff --git a/Libraries/LibCrypto/BigFraction/BigFraction.cpp b/Libraries/LibCrypto/BigFraction/BigFraction.cpp index d51c64bf758..971c7df56e7 100644 --- a/Libraries/LibCrypto/BigFraction/BigFraction.cpp +++ b/Libraries/LibCrypto/BigFraction/BigFraction.cpp @@ -1,5 +1,6 @@ /* * Copyright (c) 2022, Lucas Chollet + * Copyright (c) 2025, Manuel Zahariev * * SPDX-License-Identifier: BSD-2-Clause */ @@ -8,6 +9,7 @@ #include #include #include +#include #include namespace Crypto { @@ -134,10 +136,53 @@ BigFraction::BigFraction(double d) m_numerator.set_to(negative ? (m_numerator.negated_value()) : m_numerator); } +/* + * Complexity O(N^2), where N = number of words in the larger of denominator, numerator. + * - shifts: O(N); two copies + * - division: O(N^2): Knuth's D algorithm (UnsignedBigInteger::divided_by) + * - conversion to double: constant (64-bit quotient) + */ double BigFraction::to_double() const { - // FIXME: very naive implementation - return m_numerator.to_double() / m_denominator.to_double(); + bool const sign = m_numerator.is_negative(); + if (m_numerator.is_zero()) + return sign ? -0.0 : +0.0; + + UnsignedBigInteger numerator = m_numerator.unsigned_value(); // copy + UnsignedBigInteger const& denominator = m_denominator; + + size_t top_bit_numerator = numerator.one_based_index_of_highest_set_bit(); + size_t top_bit_denominator = denominator.one_based_index_of_highest_set_bit(); + size_t shift_left_numerator = 0; + + // 1. Shift numerator so that its most significant bit is exaclty 64 bits left tha than that of the denominator. + // NOTE: the precision of the result will be 63 bits (more than 53 bits necessary for the mantissa of a double). + if (top_bit_numerator < (top_bit_denominator + 64)) { + shift_left_numerator = top_bit_denominator + 64 - top_bit_numerator; + numerator = numerator.shift_left(shift_left_numerator); // copy + } + // NOTE: Do nothing if numerator already has more than 64 bits more than denominator. + + // 2. Divide [potentially shifted] numerator by the denominator. + auto division_result = numerator.divided_by(denominator); + if (!division_result.remainder.is_zero()) { + division_result.quotient = division_result.quotient.shift_left(1).plus(1); // Extend the quotient with a "fake 1". + // NOTE: Since the quotient has at least 63 bits, this will only affect the mantissa + // on rounding, and have the same effect on rounding as any fractional digits (from the remainder). + shift_left_numerator++; + } + + using Extractor = FloatExtractor; + Extractor double_extractor; + + // 3. Convert the quotient to_double using UnsignedBigInteger::to_double. + double_extractor.d = division_result.quotient.to_double(); + double_extractor.sign = sign; + + // 4. Shift the result back by the same number of bits as the numerator. + double_extractor.exponent -= shift_left_numerator; + + return double_extractor.d; } bool BigFraction::is_zero() const @@ -198,11 +243,15 @@ String BigFraction::to_string(unsigned rounding_threshold) const auto const number_of_digits = [](auto integer) { unsigned size = 1; - for (auto division_result = integer.divided_by(UnsignedBigInteger { 10 }); - division_result.remainder == UnsignedBigInteger { 0 } && division_result.quotient != UnsignedBigInteger { 0 }; - division_result = division_result.quotient.divided_by(UnsignedBigInteger { 10 })) { + UnsignedBigInteger const ten { 10 }; + + auto division_result = integer.divided_by(ten); + + while (division_result.remainder.is_zero() && !division_result.quotient.is_zero()) { + division_result = division_result.quotient.divided_by(ten); ++size; } + return size; }; diff --git a/Tests/LibCrypto/TestBigFraction.cpp b/Tests/LibCrypto/TestBigFraction.cpp index b0ce75fabe2..a3c19c3ec0e 100644 --- a/Tests/LibCrypto/TestBigFraction.cpp +++ b/Tests/LibCrypto/TestBigFraction.cpp @@ -1,5 +1,6 @@ /* * Copyright (c) 2024, Tim Ledbetter + * Copyright (c) 2025, Manuel Zahariev * * SPDX-License-Identifier: BSD-2-Clause */ @@ -7,6 +8,18 @@ #include #include +static Crypto::UnsignedBigInteger bigint_fibonacci(size_t n) +{ + Crypto::UnsignedBigInteger num1(0); + Crypto::UnsignedBigInteger num2(1); + for (size_t i = 0; i < n; ++i) { + Crypto::UnsignedBigInteger t = num1.plus(num2); + num2 = num1; + num1 = t; + } + return num1; +} + TEST_CASE(roundtrip_from_string) { Array valid_number_strings { @@ -26,3 +39,27 @@ TEST_CASE(roundtrip_from_string) EXPECT_EQ(result.to_string(precision), valid_number_string); } } + +TEST_CASE(big_fraction_to_double) +{ + // Golden ratio: + // - limit (inf) ratio of two consecutive fibonacci numbers + // - also ( 1 + sqrt( 5 ))/2 + Crypto::BigFraction phi(Crypto::SignedBigInteger { bigint_fibonacci(500) }, bigint_fibonacci(499)); + // Power 64 of golden ratio: + // - limit ratio of two 64-separated fibonacci numbers + // - also (23725150497407 + 10610209857723 * sqrt( 5 ))/2 + Crypto::BigFraction phi_64(Crypto::SignedBigInteger { bigint_fibonacci(564) }, bigint_fibonacci(500)); + + EXPECT_EQ(phi.to_double(), 1.618033988749895); // 1.6180339887498948482045868343656381177203091798057628621... (https://oeis.org/A001622) + EXPECT_EQ(phi_64.to_double(), 23725150497407); // 23725150497406.9999999999999578506361799772097881088769... (https://www.calculator.net/big-number-calculator.html) +} + +TEST_CASE(big_fraction_temporal_duration_precision_support) +{ + // https://github.com/tc39/test262/blob/main/test/built-ins/Temporal/Duration/prototype/total/precision-exact-mathematical-values-1.js + // Express 4000h and 1ns in hours, as a double + Crypto::BigFraction temporal_duration_precision_test = Crypto::BigFraction { Crypto::SignedBigInteger { "14400000000000001"_bigint }, "3600000000000"_bigint }; + + EXPECT_EQ(temporal_duration_precision_test.to_double(), 4000.0000000000005); +}