Achieving Precision Equal to or Greater Than R in JavaScript: Fundamentals of Double-Double Arithmetic and Its Usage in rStats.js

rStats.js is a statistics library for JavaScript and the browser. Its function names, arguments, default values, and return values ​​are modeled after R, and every function has been tested by comparing its output against that of R. It covers areas such as probability distributions, random number generation, various statistical tests, analysis of variance (ANOVA), linear models, mixed-effects models, and survival analysis. It has no external dependencies and is released under the MIT License.

• GitHub: https://github.com/numero-uno-jpn/rStats.js
• Demo: https://numero-uno-jpn.github.io/rStats.js/demo/ (You can compare the results of all functions against the values ​​produced by R.)
• npm: npm install @numero-uno-jpn/rstats

The background behind its creation and an overall introduction can be found in my previous article.

Why use "double-double" arithmetic?

In JavaScript, numbers are represented as 64-bit doubles, which provide a precision of approximately 16 decimal digits. When performing calculations using standard doubles, the discrepancy from the true value can sometimes exceed 1e-14—for instance, in the tails of probability distributions. rStats.js addresses this by employing "double-double" arithmetic (combining two doubles to achieve approximately 32 digits of precision) for these specific cases, ensuring accuracy equal to or better than that of R.

In this article, I will explain the fundamentals of double-double arithmetic using JavaScript code, describe where and how it is applied within rStats.js, and discuss how I optimized the performance of these otherwise computationally intensive operations.

When I asked an AI how to achieve precision equal to or better than that of R, it suggested "double-double" arithmetic. While I was aware of the technique, I had never used it before. Although the concept made sense immediately, I was concerned that it might be too computationally expensive. In fact, this became an issue (discussed later).

Calculations with double (double-precision floats) result in more errors than one might expect

Let’s start with a well-known example:

0.1 + 0.2; // 0.30000000000000004

When adding values ​​that differ significantly in magnitude, the smaller value is lost due to rounding.

let s = 0;
for (const x of [1e16, 1, -1e16]) s += x;
s;  // 0 (the correct result is 1)

As soon as 1 is added to 1e16, the 1 is lost due to rounding.

The same thing happens in R. While R’s sum function may utilize long double—which offers higher precision than double—if the environment supports it, the actual level of precision varies by environment. In the version of R (4.5.2) I used on Windows to write this article, .Machine$longdouble.digits was 53, meaning it had the same precision as a standard double. Consequently, the following sum is incorrect in this R environment:

sum(c(1e20, 1e3, 1e-10, -1e20, -1e3))  # -1000 (the correct result is 1e-10)

In contrast, the sum function in rStats.js returns 1e-10 for this calculation.

While developing rStats.js, I was surprised to discover quite a few instances where R’s own results were imprecise. I have compiled these findings in a document titled Differences from R. Furthermore, the fact that the precision of intermediate calculations in R can vary depending on the environment was, in itself, a surprise. WebR is an implementation of R compiled into WebAssembly, allowing it to run in browsers or within Node.js. Since WebAssembly instructions only support 32-bit and 64-bit floating-point numbers, WebR performs all calculations using 64-bit doubles. In the version of R for Windows that I examined, the precision of long double was also identical to that of double. Learning about these types of numerical errors convinced me of the utility of rStats.js.

Fundamentals of Double-Double Arithmetic

Double-double arithmetic represents a single number as the sum of two double values: hi + lo. The hi component holds the bulk of the value, while the lo component captures the remainder that cannot be represented by hi alone. This approach yields a precision of approximately 106 bits (roughly 32 decimal digits).

The foundation for this lies in the "error-free transformations" described by Dekker (1971) and Knuth (in The Art of Computer Programming, section 4.2.2). These methods allow for the precise decomposition of an arithmetic operation (such as addition or multiplication) into two parts: the rounded result and the error term discarded during the rounding process. ### Addition: twoSum

// Returns [s, e] such that a + b = s + e holds exactly
function twoSum(a, b) {
const s = a + b;
const bb = s - a;
const e = (a - (s - bb)) + (b - bb);
return [s, e];
}

twoSum(0.1, 0.2);  // [0.30000000000000004, -2.7755575615628914e-17]

Here, s is the result of a standard addition, and e represents the portion lost due to rounding. In the first example, the sum can be corrected by adding back the lost portion.

let hi = 0, lo = 0;
for (const x of [1e16, 1, -1e16]) {
const [s, e] = twoSum(hi, x);
hi = s;
lo += e;
}
hi + lo;  // 1

Multiplication: twoProd

To extract the error from a multiplication, the simplest approach in many other languages ​​is to use FMA (an instruction that performs a fused multiply-add, with only a single rounding step). However, JavaScript lacks an FMA operation. Instead, we use Veltkamp’s method to split each number into two parts—a high part and a low part—and then calculate their products without rounding. Since the product of two 26-bit values ​​fits within 52 bits, it can be represented exactly using a double-precision float (which has a 53-bit significand).

const SPLITTER = 134217729;  // 2^27 + 1

// Splits 'a' into 'hi' and 'lo' components, each with a 26-bit significand (a = hi + lo)
function split(a) {
const t = SPLITTER * a;
const hi = t - (t - a);
return [hi, a - hi];
}

// Returns [p, e] such that a * b = p + e holds exactly
function twoProd(a, b) {
const p = a * b;
const [ah, al] = split(a);
const [bh, bl] = split(b);
const e = ((ah * bh - p) + ah * bl + al * bh) + al * bl;
return [p, e];
}

twoProd(0.1, 0.1); // [0.010000000000000002, -8.326672684688674e-19]

javascript

Addition, multiplication, and division for "double-double" numbers are constructed by combining these two operations. As an example, addition is implemented as follows—this follows the same procedure as ddAdd in rStats.js.

// A version of twoSum with fewer steps; valid when |a| >= |b|
function quickTwoSum(a, b) {
const s = a + b;
return [s, b - (s - a)];
}

// Addition of two double-double numbers: (ah + al) + (bh + bl)
function ddAdd([ah, al], [bh, bl]) {
let [s, e] = twoSum(ah, bh);
const [t, f] = twoSum(al, bl);
e += t;
[s, e] = quickTwoSum(s, e);
e += f;
return quickTwoSum(s, e);
}

// Sum of 1e16 + 1 and -1e16 + 0.5
ddAdd([1e16, 1], [-1e16, 0.5]);  // [1.5, 0]

This approach uses twoSum to add the high-order components (ah and bh) and the low-order components (al and bl) separately. It then incorporates the resulting error terms step-by-step, finally formatting the result back into a high-order (hi) and low-order (lo) pair.

The code provided here is a minimal implementation intended to illustrate the underlying mechanism. The implementation in rStats.js (src/special/dd.ts) includes additional handling for edge cases: if an overflow occurs, the error term is set to 0; if a value is so large that SPLITTER * a would overflow, the value is scaled down (by a power of 2) to bring it closer to 1 before the split is performed.

Usage in rStats.js

rStats.js does not use double-double arithmetic for every calculation. It is employed only in scenarios where standard double-precision floating-point numbers are insufficient to achieve the required precision of $1 \times 10^{-14}$ . For example, the following cases are relevant:

• Exponentials of small probabilities: When calculating probabilities in the form $\exp(-E)$ , the absolute error in $E$ translates directly into a relative error in the probability. As the probability approaches $10^{-300}$ , $E$ reaches several hundred, and the error using double precision reaches $10^{-13}$ . At $E = 700$ , the mere act of rounding $E$ introduces an error of $6 \times 10^{-14}$ .

• Continued fractions for the incomplete beta function: When parameters are large and $x$ is near the mean, the result becomes highly sensitive to the rounding of $x$ . For $a = 10^{5}$ and $x = 0.999$ , using only double precision results in a deviation of $6 \times 10^{-13}$ .

• Summations for statistical tests: Calculations for means, variances, and correlations are performed using double-double precision.

Conversely, in cases where the upper bound of the error indicates that standard double precision is sufficient, the double-double operations are omitted.

The impact of this can be verified by disabling double-double arithmetic internally. When disabled, only about half of the test results maintain a precision of $10^{-14}$ . In some cases, the degradation is significant: the beta function (beta) shows a deviation of up to $3 \times 10^{-5}$ , and the pt function (Student's t-distribution) deviates by up to $6 \times 10^{-7}$ .

The absence of FMA (Fused Multiply-Add) in JavaScript also influences the choice of R—the reference implementation used for comparison. The arm64 version of R utilizes FMA, which causes the final digit of the second value generated by set.seed(42); rnorm(3) to differ from that of the x86_64 version. rStats.js...