Source/JavaScriptCore/runtime/FractionToDouble.cpp

/*
 * Copyright (C) 2025 Igalia, S.L. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY APPLE INC. AND ITS CONTRIBUTORS ``AS IS''
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL APPLE INC. OR ITS CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
 * THE POSSIBILITY OF SUCH DAMAGE.
 */

#include "config.h"
#include "FractionToDouble.h"

#include "MathCommon.h"

// The calculations here are based on algorithms from two sources. The second
// one builds on the first.
//
// Shewchuk (1997). Adaptive precision floating-point arithmetic and fast robust
//   geometric predicates. Discrete & Computational Geometry 18(3), pp. 305–363.
//   https://doi.org/10.1007/PL00009321
//
// Hida, Li, Bailey (2008). Library for double-double and quad-double
//   arithmetic. Manuscript. https://www.davidhbailey.com/dhbpapers/qd.pdf
//   and the accompanying QD library https://github.com/BL-highprecision/QD,
//   which is BSD-licensed.

namespace JSC {

// Double-double precision floating point number, represented as the unevaluated
// sum of two doubles. In other words, dd[0] is the double approximation term
// and dd[1] is the error term.
//
// There are many such representations, but only one is 'normalized' meaning the
// dd[0] term is the most accurate possible double-precision approximation of
// the double-double value.
using DD = std::array<double, 2>;

// Conversion of Int128 to double-double precision floating point. The
// calculations follow from the definition of hi and lo: hi is the closest
// double-precision approximation to the exact value (which itself will be a
// safe integer) and lo is the error term.
static DD NODELETE int128ToDD(const Int128& value)
{
    double hi = static_cast<double>(value);
    double lo = static_cast<double>(value - static_cast<Int128>(hi));
    return { hi, lo };
}

// Computes double-double precision a + b of two doubles a and b. This is the
// Two-Sum algorithm in theorem 7 of the Shewchuk paper.
static DD NODELETE ddSum(double a, double b)
{
    // First compute the double-precision approximation of the sum by regular
    // double addition.
    double sum = a + b;

    // Compute the error term.
    double bVirtual = sum - a;
    double aVirtual = sum - bVirtual;
    double bRoundoff = b - bVirtual;
    double aRoundoff = a - aVirtual;
    double error = aRoundoff + bRoundoff;

    return { sum, error };
}

// Computes double-double precision a * b of two doubles a and b. The
// optimization using std::fma() is suggested in section 2 of the Hida-Li-Bailey
// paper.
static DD ddProduct(double a, double b)
{
    // First compute the double-precision approximation of the product by
    // regular double multiplication.
    double product = a * b;

    // On armv8, this emits the fnmsub instruction.
    // On x86_64, this emits the vfmsub213sd instruction if compiled with SSE
    // instructions. If not, it calls libm's fma(), which is comparably fast to
    // using the Two-Product algorithm in theorem 18 of the Shewchuk paper.
    double error = std::fma(a, b, -product);

    return { product, error };
}

// Computes double-double precision numerator / denominator, where divisor is a
// double, and rounds the result to double precision. This is described in
// section 3.5 of the Hida-Li-Bailey paper.
static double fractionToDoubleSlow(const Int128& numerator, double denominator)
{
    DD ddNumerator = int128ToDD(numerator);

    // Compute a first approximation of the quotient by regular double division.
    double quotient0 = ddNumerator[0] / denominator;

    // Compute remainder, ddNumerator - quotient0 * denominator.
    DD product = ddProduct(quotient0, denominator);
    DD remainder = ddSum(ddNumerator[0], -product[0]);

    // Compute the next approximation term.
    double error = remainder[1] + ddNumerator[1] - product[1];
    double quotient1 = (remainder[0] + error) / denominator;

    // The result is DD { quotient0, quotient1 }. If we wanted double-double
    // precision here, we would have to use the Fast-Two-Sum algorithm from
    // theorem 6 of the Shewchuk paper to renormalize the two terms, but since
    // we only need double precision we can discard the error term.
    return quotient0 + quotient1;
}

double fractionToDouble(const Int128& numerator, double denominator)
{
    ASSERT(denominator > 0);
    ASSERT(isSafeInteger(denominator));

    if (!numerator)
        return 0;

    // When the denominator is 1, we are just calculating the double
    // approximation of the numerator.
    if (denominator == 1)
        return static_cast<double>(numerator);

    // When the numerator can be represented exactly as a double the algorithm
    // collapses to a simple double division.
    if (isSafeInteger(static_cast<double>(numerator))) [[likely]]
        return static_cast<double>(numerator) / denominator;

    // Otherwise use double-double precision to compute the result.
    return fractionToDoubleSlow(numerator, denominator);
}

// Extends section 3.5 of the Hida-Li-Bailey paper to an Int128 denominator.
static double fractionToDoubleSlow(const Int128& numerator, const Int128& denominator)
{
    // n0 + n1 = N exactly, d0 + d1 = D exactly.
    DD n = int128ToDD(numerator);
    DD d = int128ToDD(denominator);

    // We want Q = N / D, where Q = q0 + q1.
    // => (q0 + q1) * D = N
    // => residual = q1 * D = N - q0 * D
    // => q1 = (N - q0 * D) / D

    // Step 1: q0 = first approximation of Q.
    // n0 and d0 are the best doubles for N and D; n0/d0 is the best
    // double approximation of N/D using regular double division.
    // (Can't divide N/D directly — both may exceed 2^53 as Int128.)
    double q0 = n[0] / d[0];

    // Step 2: residual = N - q0 * D, split as r0 + error.
    //
    // Expand: N - q0 * D
    // => (n0 + n1) - q0 * (d0 + d1)
    // => (n0 - q0 * d0) + n1 - q0 * d1
    //
    DD p = ddProduct(q0, d[0]); // p0 + p1 = q0 * d0
    DD r = ddSum(n[0], -p[0]); // r0 + r1 = n0 - p0
    // Expand: n0 - q0 * d0
    // => n0 - p0 - p1
    // => r0 + r1 - p1
    //
    // Substituting back:
    //   residual = r0 + r1 - p1 + n1 - q0 * d1
    //            = r0 + error
    double error = r[1] + n[1] - p[1] - q0 * d[1];

    // Step 3: q1 = residual / D ≈ (r[0] + error) / d0.
    // Using d0 instead of D: error ≈ q1 * (d1 / d0) ≈ Q*2^-104, negligible.
    double q1 = (r[0] + error) / d[0];

    // Discard DD renormalization (Shewchuk Fast-Two-Sum, theorem 6):
    // single-precision result is sufficient here.
    return q0 + q1;
}

double fractionToDouble(const Int128& numerator, const Int128& denominator)
{
    ASSERT(denominator > 0);

    if (!numerator)
        return 0;
    if (denominator == 1)
        return static_cast<double>(numerator);
    if (isSafeInteger(static_cast<double>(numerator)) && isSafeInteger(static_cast<double>(denominator))) [[likely]]
        return static_cast<double>(numerator) / static_cast<double>(denominator);
    return fractionToDoubleSlow(numerator, denominator);
}

} // namespace JSC