1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
/*
* 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 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 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);
}
} // namespace JSC
|
