1
0
Fork 0
mirror of https://github.com/LadybirdBrowser/ladybird.git synced 2025-06-08 05:27:14 +09:00

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.
This commit is contained in:
Manuel Zahariev 2025-03-20 12:58:18 -07:00 committed by Ali Mohammad Pur
parent d2ea77c099
commit 4ed8e9e596
Notes: github-actions[bot] 2025-03-23 18:34:20 +00:00
2 changed files with 91 additions and 5 deletions

View file

@ -1,5 +1,6 @@
/*
* Copyright (c) 2022, Lucas Chollet <lucas.chollet@free.fr>
* Copyright (c) 2025, Manuel Zahariev <manuel@duck.com>
*
* SPDX-License-Identifier: BSD-2-Clause
*/
@ -8,6 +9,7 @@
#include <AK/ByteString.h>
#include <AK/Math.h>
#include <AK/StringBuilder.h>
#include <LibCrypto/BigInt/UnsignedBigInteger.h>
#include <LibCrypto/NumberTheory/ModularFunctions.h>
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<double>;
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;
};

View file

@ -1,5 +1,6 @@
/*
* Copyright (c) 2024, Tim Ledbetter <timledbetter@gmail.com>
* Copyright (c) 2025, Manuel Zahariev <manuel@duck.com>
*
* SPDX-License-Identifier: BSD-2-Clause
*/
@ -7,6 +8,18 @@
#include <LibCrypto/BigFraction/BigFraction.h>
#include <LibTest/TestCase.h>
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);
}