Works out sin, cos, exp, log, powers and inverse trigonometric functions from + − × ÷ and √ only, so every browser gives the same bits.
How the precision is kept
Each function reduces its argument and sums a truncated series in 64-bit; the series length and the argument reduction set the accuracy. The file's header states about 1e-15 for sine and cosine for |x| < 1e6 and a few units in the last place for the others. Only basic IEEE 754 operations are used, which every engine rounds the same way.
Open item. Split 2π so that n × TWO_PI_HI is exact (a high part with trailing zero bits), or state the accuracy of sine and cosine as about |x| × 1.1e-16.
Calculations under this policy
21 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
n! is built by repeated multiplication and each coefficient is one division. Factorials from 23! up are no longer exact in 64-bit; the terms that use them are below 2e-18 for |r| ≤ π/2.
field/detmath.js
lines 14 to 17
PartCraft 0.8.1
let f = 1; // 1/n! built by exact divisions, not typed in
const inv = [1];
for (let n = 1; n <= 26; n++) { f *= n; inv.push(1 / f); }
for (let n = 0; n <= 12; n++) { SIN_C.push((n % 2 ? -1 : 1) * inv[2 * n + 1]); COS_C.push((n % 2 ? -1 : 1) * inv[2 * n]); }
Angle reduction
64-bit
Brings an angle into [−π, π].
n = floor(x ÷ 2π + ½); r = (x − n × 2π_hi) − n × 2π_lo
TWO_PI_HI is the 64-bit number nearest 2π and uses nearly all of its 53 bits, so n × TWO_PI_HI is rounded for all but the smallest n, by up to half a unit in the last place of x. That error stays in r: the reduced angle is good to about |x| × 1.1e-16, not to a fixed 1e-15. TWO_PI_LO then corrects for the difference between TWO_PI_HI and 2π.
field/detmath.js
lines 19 to 22
PartCraft 0.8.1
function reduce(x) { // → r in [−π, π]
const n = Math.floor(x * INV_TWO_PI + 0.5);
return (x - n * TWO_PI_HI) - n * TWO_PI_LO;
}
Sine series
64-bit
The sine of a reduced angle.
sin r = r(c₀ + c₁r² + … + c₁₂r²⁴)
Horner's rule in r², 13 coefficients (through r²⁵). The first term left out is below 2e-23 for |r| ≤ π/2.
field/detmath.js
line 23
PartCraft 0.8.1
function sinPoly(r) { const r2 = r * r; let s = SIN_C[12]; for (let n = 11; n >= 0; n--) s = s * r2 + SIN_C[n]; return s * r; }
Cosine series
64-bit
The cosine of a reduced angle.
cos r = c₀ + c₁r² + … + c₁₂r²⁴
Horner's rule in r², 13 coefficients (through r²⁴). The first term left out is below 4e-22 for |r| ≤ π/2.
field/detmath.js
line 24
PartCraft 0.8.1
function cosPoly(r) { const r2 = r * r; let s = COS_C[12]; for (let n = 11; n >= 0; n--) s = s * r2 + COS_C[n]; return s; }
Sine
64-bit
sin x for any finite x.
sin r = sin(π − r) for r > π/2; sin r = sin(−π − r) for r < −π/2
After reduction the angle is folded into [−π/2, π/2]. π here is TWO_PI_HI ÷ 2, the 64-bit π, so the fold adds up to about 1.2e-16 to the folded angle.
field/detmath.js
lines 25 to 29
PartCraft 0.8.1
export function dsin(x) {
let r = reduce(x);
if (r > HALF_PI) r = PI - r; else if (r < -HALF_PI) r = -PI - r;
return sinPoly(r);
}
Cosine
64-bit
cos x for any finite x.
cos r = −cos(π − |r|) for |r| > π/2
The reduced angle's absolute value is taken first, then folded with the 64-bit π as in the sine.
field/detmath.js
lines 30 to 33
PartCraft 0.8.1
export function dcos(x) {
const r = Math.abs(reduce(x));
return r > HALF_PI ? -cosPoly(PI - r) : cosPoly(r);
}
Exponential
64-bit
The exponential eˣ of a number.
k = floor(x ÷ ln 2 + ½); r = (x − k × ln2_hi) − k × ln2_lo; eˣ = 2ᵏ × Σ rⁿ ÷ n! (n = 0 … 18)
LN2_HI is the split constant whose low 21 bits are zero, so k × LN2_HI is exact for every k that occurs. The 18 terms are added to 1 from the largest down. 2ᵏ is built from exponent bits, in two steps when k is outside −1022 … 1023. Above 709.78 the result is Infinity, below −745.13 it is 0.
field/detmath.js
lines 43 to 53
PartCraft 0.8.1
export function dexp(x) {
if (x !== x) return NaN;
if (x > 709.782712893384) return Infinity;
if (x < -745.1332191019411) return 0;
const k = Math.floor(x * INV_LN2 + 0.5), r = (x - k * LN2_HI) - k * LN2_LO;
let p = 1, term = 1;
for (let n = 1; n <= 18; n++) { term = term * r / n; p += term; }
if (k > 1023) return p * pow2(1023) * pow2(k - 1023);
if (k < -1022) return p * pow2(-1022) * pow2(k + 1022);
return p * pow2(k);
}
Natural logarithm
64-bit
The natural logarithm ln x of a number.
x = m × 2ᵉ, m in [√2/2, √2]; ln m = 2(s + s³/3 + … + s²⁹/29), s = (m − 1) ÷ (m + 1); ln x = e × ln2_hi + (ln m + e × ln2_lo)
The exponent and mantissa are read from the bits of the double, so the split is exact; subnormal inputs are scaled by 2⁵⁴ first. With |s| ≤ 0.172 the first term left out is below 2e-25. The small part e × ln2_lo is added to the series before the large part e × ln2_hi.
field/detmath.js
lines 54 to 69
PartCraft 0.8.1
export function dlog(x) {
if (x !== x || x < 0) return NaN;
if (x === 0) return -Infinity;
if (x === Infinity) return Infinity;
F64[0] = x;
let e = (U32[HI] >>> 20) & 0x7ff;
if (e === 0) { F64[0] = x * 18014398509481984; e = ((U32[HI] >>> 20) & 0x7ff) - 54; } // subnormal: scale by 2^54
e -= 1023;
U32[HI] = (U32[HI] & 0x000fffff) | 0x3ff00000; // mantissa m in [1, 2)
let m = F64[0];
if (m > 1.4142135623730951) { m *= 0.5; e += 1; } // m in [√2/2, √2]
const s = (m - 1) / (m + 1), s2 = s * s; // ln m = 2 atanh s = 2 (s + s³/3 + s⁵/5 + …), |s| ≤ 0.172
let acc = 1 / 29;
for (let n = 27; n >= 1; n -= 2) acc = acc * s2 + 1 / n;
return e * LN2_HI + (2 * s * acc + e * LN2_LO);
}
Base-2 logarithm
64-bit
The base-2 logarithm log₂ x of a number.
An exact power of two (all mantissa bits zero, a normal exponent) returns its exponent as a whole number. Otherwise ln x × 1.4426950408889634.
When the quotient is within 1e-9 of a whole number n between −22 and 22, and 10ⁿ built by repeated squaring equals x exactly, n itself is returned. Otherwise the quotient is returned as computed.
field/detmath.js
lines 75 to 80
PartCraft 0.8.1
export function dlog10(x) {
const r = dlog(x) / LN10;
const n = Math.round(r);
if (n >= -22 && n <= 22 && Math.abs(r - n) < 1e-9 && (n >= 0 ? ipow(10, n) === x : 1 / ipow(10, -n) === x)) return n;
return r;
}
Whole-number power
64-bit
xⁿ for a whole n ≥ 0.
Repeated squaring: each squaring and each multiplication into the result is rounded. Exact when every intermediate product can be held exactly.
field/detmath.js
line 81
PartCraft 0.8.1
function ipow(x, n) { let r = 1; while (n > 0) { if (n & 1) r *= x; x *= x; n = Math.floor(n / 2); } return r; }
Power
64-bit
A number x raised to a power y.
xʸ = exp(y × ln x) for a fractional y
A whole power up to ±1024 uses ipow, with one division for a negative power. A fractional power is dexp(y × dlog(x)): the rounding of the product y × ln x becomes a relative error of that size in the result, so large results lose digits. A negative base with a fractional power gives NaN.
field/detmath.js
lines 82 to 89
PartCraft 0.8.1
export function dpow(x, y) {
if (y === 0) return 1;
if (x !== x || y !== y) return Math.pow(x, y);
if (Number.isInteger(y) && Math.abs(y) <= 1024) { const r = ipow(x, Math.abs(y)); return y < 0 ? 1 / r : r; }
if (x === 1 || x === 0 || !Number.isFinite(x) || !Number.isFinite(y)) return Math.pow(x, y);
if (x < 0) return NaN; // a negative base with a fractional power
return dexp(y * dlog(x));
}
Arctangent
64-bit
The arctangent atan x of a number.
atan z = 2·atan(z ÷ (1 + √(1 + z²))), twice; then 4(z − z³/3 + … + z⁴¹/41)
Arguments above 1 use π/2 − atan(1/x). After two halvings |z| ≤ 0.199, where the first term left out is below 1e-31. π/2 is the 64-bit value.
field/detmath.js
lines 90 to 103
PartCraft 0.8.1
export function datan(x) {
if (x !== x) return NaN;
const a = Math.abs(x);
if (a === Infinity) return x > 0 ? HALF_PI : -HALF_PI;
let z = a > 1 ? 1 / a : a;
z = z / (1 + Math.sqrt(1 + z * z)); // atan z = 2 atan(z / (1 + √(1+z²)))
z = z / (1 + Math.sqrt(1 + z * z)); // |z| ≤ tan(π/16) ≈ 0.199 after two halvings
const z2 = z * z; // atan z = z − z³/3 + z⁵/5 − …
let acc = 0;
for (let n = 41; n >= 1; n -= 2) acc = acc * (-z2) + 1 / n;
let r = 4 * z * acc;
if (a > 1) r = HALF_PI - r;
return x < 0 ? -r : r;
}
Two-argument arctangent
64-bit
The angle of the point (x, y).
The quotient y ÷ x is formed first and can overflow to Infinity for a tiny x, which datan maps to ±π/2. For x < 0 the 64-bit π is added or subtracted. Both arguments zero defers to Math.atan2.
field/detmath.js
lines 104 to 110
PartCraft 0.8.1
export function datan2(y, x) {
if (x !== x || y !== y) return NaN;
if (x === 0 && y === 0) return Math.atan2(y, x);
if (x > 0) return datan(y / x);
if (x < 0) return y >= 0 && !(y === 0 && 1 / y < 0) ? datan(y / x) + PI : datan(y / x) - PI;
return y > 0 ? HALF_PI : -HALF_PI;
}
Arcsine
64-bit
The arcsine asin x of a number in [−1, 1].
asin x = atan2(x, √((1 − x)(1 + x)))
(1 − x)(1 + x) is used instead of 1 − x², which keeps digits near |x| = 1. Outside [−1, 1] the result is NaN; the argument is not clamped.
Two series and one division; the division is not guarded where the cosine is zero.
field/detmath.js
line 113
PartCraft 0.8.1
export const dtan = (x) => dsin(x) / dcos(x);
Hyperbolic sine
64-bit
The hyperbolic sine sinh x of a number.
sinh x = (eˣ − e⁻ˣ) ÷ 2
Below |x| = 0.1 a nested series through x¹¹ is used, which avoids the subtraction of nearly equal numbers; above it the formula with dexp and one reciprocal.
A sum of two positive terms from one dexp and one reciprocal; nothing cancels.
field/detmath.js
line 115
PartCraft 0.8.1
export function dcosh(x) { const e = dexp(x); return (e + 1 / e) / 2; }
Hyperbolic tangent
64-bit
The hyperbolic tangent tanh x of a number.
tanh x = (e²ˣ − 1) ÷ (e²ˣ + 1)
±1 beyond |x| = 20; dsinh ÷ dcosh below |x| = 0.1; the formula in between.
field/detmath.js
line 116
PartCraft 0.8.1
export function dtanh(x) { if (Math.abs(x) > 20) return x > 0 ? 1 : -1; if (Math.abs(x) < 0.1) return dsinh(x) / dcosh(x); const e = dexp(2 * x); return (e - 1) / (e + 1); }
Cube root
64-bit
The cube root ∛x of a number.
r₀ = exp(ln|x| ÷ 3); r ← r − (r³ − |x|) ÷ (3r²), twice
A first value from dexp and dlog, then exactly two Newton steps; the steps are not repeated to a tolerance. The sign of x is put back at the end.
field/detmath.js
line 117
PartCraft 0.8.1
export function dcbrt(x) { if (x === 0 || !Number.isFinite(x)) return x; const a = Math.abs(x); let r = dexp(dlog(a) / 3); r = r - (r * r * r - a) / (3 * r * r); r = r - (r * r * r - a) / (3 * r * r); return x < 0 ? -r : r; }
In detail
Method
Sine and cosine: reduction by 2π in two parts (the header calls it Cody–Waite), then 13 Taylor terms on [−π/2, π/2]. exp: reduction by ln 2 in two parts, 18 Taylor terms, scaling by 2ᵏ. log: mantissa moved to [√2/2, √2], series of 2·atanh((m − 1) ÷ (m + 1)). atan: two half-angle reductions, 21 series terms.
Number types
64-bit numbers; the bits of a double are read through a Float64Array and Uint32Array pair.
Tolerances
dlog10 returns a whole number n when the result is within 1e-9 of it, |n| ≤ 22, and 10ⁿ equals the input exactly. dsinh and dtanh switch to a series below |x| = 0.1; dtanh returns ±1 above |x| = 20.
Limits
In sine and cosine the product n × TWO_PI_HI is rounded, so the reduced angle is known to about |x| × 1.1e-16; the error grows with the size of the angle. A fractional power is exp(y × log x), so the rounding of y × log x is carried into the result. dasin and dacos return NaN outside [−1, 1]; nothing is clamped.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.