The built-in math of the scripting language: elementary functions, rounding, the gamma function, sums and statistics, interpolation, polynomial fitting and number-to-text conversion.
How the precision is kept
Values are 64-bit numbers in Float64Array tensors. Trigonometric, exponential and logarithm functions come from field/detmath.js, which the code comment says is built from + − × ÷ and square root so that a script gives the same bits in every engine; their accuracy is that module's. Several functions here are built on top of those with formulas that lose digits in part of their range, and the gamma function is a Lanczos approximation.
Open item. Give log1p, expm1, asinh, acosh and atanh forms that keep 15 digits for small arguments, let gamma(171) return the factorial of 170, and evaluate the gamma function for large non-whole arguments through logarithms so it does not overflow early.
Calculations under this policy
18 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.
log1p and expm1 keep three series terms below 1e-4, where the next term is up to 2.5e-13 (log1p) or 4e-14 (expm1) of the result; above the switch they form 1 + x or exp(x) − 1, which rounds away part of a small x. asinh, acosh and atanh take a logarithm of a value near 1 for small arguments and lose digits in proportion. Degree functions multiply by Math.PI ÷ 180 first, so their argument carries the rounding of π. round() rounds halves away from zero.
cs-lib-core.js
lines 267 to 278
PartCraft 0.8.1
const MATH1 = {
// the transcendental ones are field/detmath.js's (+ − × ÷ and sqrt only), so a script gives the same bits in Node, Chrome, Firefox and Safari
abs: Math.abs, sqrt: Math.sqrt, exp: dexp, log: dlog, log2: dlog2, log10: dlog10,
log1p: (x) => (Math.abs(x) < 1e-4 ? x - x * x / 2 + x * x * x / 3 : dlog(1 + x)), expm1: (x) => (Math.abs(x) < 1e-4 ? x + x * x / 2 + x * x * x / 6 : dexp(x) - 1),
sin: dsin, cos: dcos, tan: dtan, asin: dasin, acos: dacos, atan: datan, sinh: dsinh, cosh: dcosh, tanh: dtanh,
asinh: (x) => (x < 0 ? -1 : 1) * dlog(Math.abs(x) + Math.sqrt(x * x + 1)), acosh: (x) => (x < 1 ? NaN : dlog(x + Math.sqrt(x * x - 1))), atanh: (x) => (x > 1 || x < -1 ? NaN : 0.5 * dlog((1 + x) / (1 - x))),
floor: Math.floor, ceil: Math.ceil, fix: Math.trunc,
round: (x) => (x < 0 ? -Math.round(-x) : Math.round(x)), sign: Math.sign, cbrt: dcbrt,
sind: (x) => dsin(x * Math.PI / 180), cosd: (x) => dcos(x * Math.PI / 180), tand: (x) => dtan(x * Math.PI / 180),
deg2rad: (x) => x * Math.PI / 180, rad2deg: (x) => x * 180 / Math.PI, gamma, factorial: (x) => (x < 0 || !Number.isInteger(x) ? NaN : gamma(x + 1)),
real: (x) => x, conj: (x) => x, imag: () => 0,
};
Rounding to a number of decimals
64-bit
round(x, d): rounds to d decimal places.
round(x × 10^d) ÷ 10^d
The product x × 10^d is rounded in binary before the rounding to a whole number, so a value that is a half in decimal (such as 1.005 with d = 2) can round either way. The result is the 64-bit number nearest the decimal value.
cs-lib-core.js
line 281
PartCraft 0.8.1
if (name === 'round' && a[1] !== undefined) { const d = I(a[1], t('script.core.what.round_digits')), m = 10 ** d; return K.map(T(a[0], 'round'), (x) => MATH1.round(x * m) / m); }
Gamma function and factorial
64-bit
Γ(x), and factorial(n) = Γ(n + 1).
whole x < 171: (x − 1)! by product; x < 0.5: π ÷ (sin(πx) Γ(1 − x)); else √(2π) t^(x − ½) e^(−t) Σ c_i ÷ (x − 1 + i), t = x + 6.5
For whole arguments the running product is exact up to the factorial of 22 (its odd part still fits in 53 bits) and is rounded at each multiplication after that. The whole-number branch stops at x < 171, so gamma(171) returns Infinity although its value, the factorial of 170, is about 7.26e306 and representable. Other arguments use a Lanczos sum with nine coefficients written to 17 digits; the code states no error bound for it. The power t^(x − ½) overflows before the result does for non-whole x above about 142.
cs-lib-core.js
lines 284 to 291
PartCraft 0.8.1
function gamma(x) {
if (Number.isInteger(x)) { if (x <= 0) return Infinity; if (x < 171) { let r = 1; for (let i = 2; i < x; i++) r *= i; return r; } return Infinity; }
if (x < 0.5) return Math.PI / (dsin(Math.PI * x) * gamma(1 - x));
const g = 7, c = [0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
x -= 1; let s = c[0]; for (let i = 1; i < g + 2; i++) s += c[i] / (x + i);
const t = x + g + 0.5;
return Math.sqrt(2 * Math.PI) * dpow(t, x + 0.5) * dexp(-t) * s;
}
Binomial coefficient
64-bit
nchoosek(n, k): the number of ways to choose k of n.
r ← r × (n − m + i) ÷ i for i = 1 … m, m = min(k, n − k); result = round(r)
Each step multiplies and divides in 64-bit and the final value is rounded to a whole number. The result is exact while the running product stays below 2⁵³; beyond that it is the nearest representable value at best.
cs-lib-core.js
line 305
PartCraft 0.8.1
def('nchoosek', (a) => { const n = I(a[0], 'nchoosek: n'), k = I(a[1], 'nchoosek: k'); if (k < 0 || k > n) return 0; let r = 1; for (let i = 1; i <= Math.min(k, n - k); i++) r = r * (n - Math.min(k, n - k) + i) / i; return Math.round(r); }, 'nchoosek(n, k)');
Spacing of numbers
64-bit
eps(x): the gap between x and the next larger 64-bit number.
2^(floor(log₂|x|) − 52); eps(0) = 5e-324
The exponent is taken from Math.log2, whose rounding just below a power of two can give the next exponent up. For very small x the power underflows to 0 rather than the smallest positive number.
cs-lib-core.js
line 75
PartCraft 0.8.1
def('eps', (a) => { if (!a.length) return Number.EPSILON; return K.map(T(a[0], 'eps'), (x) => { x = Math.abs(x); if (!Number.isFinite(x)) return NaN; if (x === 0) return 5e-324; const e = Math.floor(Math.log2(x)); return 2 ** (e - 52); }); }, 'eps | eps(x)'); // i18n-ignore: call signature (code)
Evenly spaced values
64-bit
linspace(a, b, n): n values from a to b.
value i = a + (b − a) × i ÷ (n − 1)
Each value is computed from the ends, so spacing errors do not accumulate. The last value is a + (b − a), which can differ from b in its last bit; only the one-point case returns b itself.
cs-lib-core.js
line 92
PartCraft 0.8.1
def('linspace', (a) => { const lo = N(a[0], t('script.core.what.linspace_start')), hi = N(a[1], t('script.core.what.linspace_end')), n = a[2] === undefined ? 100 : Math.max(0, Math.floor(N(a[2], t('script.core.what.linspace_count')))); const d = new Float64Array(n); for (let i = 0; i < n; i++) d[i] = n === 1 ? hi : lo + (hi - lo) * i / (n - 1); return n === 1 ? d[0] : tensor([1, n], d); }, 'linspace(a, b, n)');
Logarithmically spaced values
64-bit
logspace(a, b, n): n values from 10^a to 10^b.
value i = 10 ^ (a + (b − a) × i ÷ (n − 1))
Uses JavaScript's built-in power operator, not the deterministic power from detmath.js that the other functions use, so the last bit can differ between engines.
cs-lib-core.js
line 93
PartCraft 0.8.1
def('logspace', (a) => { const lo = N(a[0]), hi = N(a[1]), n = a[2] === undefined ? 50 : I(a[2], 'logspace: n'); const d = new Float64Array(n); for (let i = 0; i < n; i++) d[i] = 10 ** (n === 1 ? hi : lo + (hi - lo) * i / (n - 1)); return tensor([1, n], d); }, 'logspace(a, b, n)');
Running sum and product
64-bit
cumsum and cumprod along a dimension.
acc ← acc + x_k (or acc × x_k), stored at each step
A plain running total in 64-bit, added in index order with no compensation, so rounding accumulates along a long vector.
cs-lib-core.js
lines 354 to 361
PartCraft 0.8.1
for (const [name, prod] of [['cumsum', false], ['cumprod', true]]) def(name, (a) => {
const t = T(a[0], name), dim = a[1] === undefined ? K.defaultDim(t.shape) : dimArg(a[1], name);
const s = t.shape.slice(); while (s.length <= dim) s.push(1);
let inner = 1; for (let k = 0; k < dim; k++) inner *= s[k];
const len = s[dim], outer = t.data.length / (inner * len || 1), o = new Float64Array(t.data.length);
for (let q = 0; q < outer; q++) for (let i = 0; i < inner; i++) { let acc = prod ? 1 : 0; for (let k = 0; k < len; k++) { const p = q * inner * len + k * inner + i; acc = prod ? acc * t.data[p] : acc + t.data[p]; o[p] = acc; } }
return out(t.shape, o);
}, `${name}(A, dim)`); // i18n-ignore: call signature (code)
Variance and standard deviation
64-bit
var and std of a set of values.
m = Σx ÷ n; var = Σ(x − m)² ÷ (n − 1), or ÷ n with weight 1; std = √var
Two passes: the mean first, then squared deviations from it, which avoids the cancellation of the Σx² − n m² form. Sums are added in order. One value or none gives 0.
cs-lib-core.js
line 365
PartCraft 0.8.1
function variance(x, w) { const n = x.length; if (n <= 1) return 0; let m = 0; for (const v of x) m += v; m /= n; let s = 0; for (const v of x) s += (v - m) ** 2; const ww = w === undefined || (isTensor(w) && !w.data.length) ? 0 : N(w); return s / (ww === 1 ? n : n - 1); }
Median
64-bit
The middle value of a set, or the mean of the two middle values.
Sorted in a Float64Array. For an even count the sum of the two middle values is formed before halving, which can overflow for values near the largest representable number.
cs-lib-core.js
line 362
PartCraft 0.8.1
def('median', (a) => statAlong(a, 'median', (x) => { const s = Float64Array.from(x).sort(); const n = s.length; return n === 0 ? NaN : n % 2 ? s[(n - 1) / 2] : (s[n / 2 - 1] + s[n / 2]) / 2; }), 'median(A, dim)');
Dot product
64-bit
dot(a, b): the sum of products of two vectors.
Σ a_i b_i
Products added in index order with no compensation. For nearly perpendicular long vectors the terms cancel and the result carries the rounding of the largest terms.
cs-lib-core.js
line 375
PartCraft 0.8.1
def('dot', (a) => { const x = T(a[0], 'dot'), y = T(a[1], 'dot'); if (x.data.length !== y.data.length) throw err(t('script.core.dot_len', { a: x.data.length, b: y.data.length })); let s = 0; for (let i = 0; i < x.data.length; i++) s += x.data[i] * y.data[i]; return s; }, 'dot(a, b)');
Trapezoid integration
64-bit
trapz(y) or trapz(x, y): the area under sampled data.
Σ (x_i − x_{i−1}) (y_i + y_{i−1}) ÷ 2
The arithmetic is exact to rounding; the answer's accuracy is that of the trapezoid rule on the given samples.
cs-lib-core.js
line 390
PartCraft 0.8.1
def('trapz', (a) => { if (a.length === 1) { const y = T(a[0], 'trapz').data; let s = 0; for (let i = 1; i < y.length; i++) s += (y[i] + y[i - 1]) / 2; return s; } const x = T(a[0], 'trapz').data, y = T(a[1], 'trapz').data; if (x.length !== y.length) throw err(t('script.core.trapz_len')); let s = 0; for (let i = 1; i < y.length; i++) s += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2; return s; }, 'trapz(y) | trapz(x, y)');
One-dimensional interpolation
64-bit
interp1(x, y, xi): linear, nearest or previous-value interpolation in a table.
The interval is found by bisection on the sorted x values, which must be strictly increasing. Outside the table the result is NaN unless extrapolation or a fill value is asked for; extrapolation continues the end interval's line.
cs-lib-core.js
lines 502 to 517
PartCraft 0.8.1
def('interp1', (a) => {
const x = T(a[0], 'interp1').data, y = T(a[1], 'interp1').data, xi = T(a[2], 'interp1');
const method = isStr(a[3]) ? a[3] : 'linear', extrap = a[a.length - 1] === 'extrap' ? true : (typeof a[a.length - 1] === 'number' && a.length >= 5 ? a[a.length - 1] : null);
if (x.length !== y.length || x.length < 2) throw err(t('script.core.interp1_len'));
for (let i = 1; i < x.length; i++) if (!(x[i] > x[i - 1])) throw err(t('script.core.interp1_increasing'));
const o = K.map(xi, (q) => {
if (q < x[0] || q > x[x.length - 1]) { if (extrap === null) return NaN; if (typeof extrap === 'number') return extrap; }
let lo = 0, hi = x.length - 1;
while (hi - lo > 1) { const m = (lo + hi) >> 1; if (x[m] <= q) lo = m; else hi = m; }
if (q < x[0]) { lo = 0; hi = 1; } else if (q > x[x.length - 1]) { lo = x.length - 2; hi = x.length - 1; }
if (method === 'nearest') return q - x[lo] <= x[hi] - q ? y[lo] : y[hi];
if (method === 'previous') return q >= x[hi] ? y[hi] : y[lo];
return y[lo] + (y[hi] - y[lo]) * (q - x[lo]) / (x[hi] - x[lo]);
});
return o;
}, 'interp1(x, y, xi, method, \'extrap\')');
Polynomial evaluation
64-bit
polyval(p, x): the value of a polynomial.
s ← s × x + c for each coefficient (Horner)
Horner's rule, one multiply and one add per coefficient. Near a root of the polynomial the result is a small difference of large terms.
cs-lib-core.js
line 518
PartCraft 0.8.1
def('polyval', (a) => { const p = T(a[0], 'polyval').data; return K.map(dense(a[1]), (x) => { let s = 0; for (const c of p) s = s * x + c; return s; }); }, 'polyval(p, x)');
Polynomial fit
64-bit
polyfit(x, y, n): the least-squares polynomial of degree n through data.
V_ij = x_i^(n − j); solve V p ≈ y
Builds the Vandermonde matrix from raw powers of x with the built-in power operator and passes it to the solver in cs-tensor.js. The x values are not centered or scaled, so for a high degree or x far from zero the matrix is poorly conditioned and coefficients lose digits.
cs-lib-core.js
line 519
PartCraft 0.8.1
def('polyfit', (a) => { const x = T(a[0], 'polyfit').data, y = T(a[1], 'polyfit').data, n = I(a[2], t('script.core.what.polyfit_degree')); const m = x.length; const V = new Float64Array(m * (n + 1)); for (let i = 0; i < m; i++) for (let j = 0; j <= n; j++) V[i + j * m] = x[i] ** (n - j); const p = K.solve(tensor([m, n + 1], V), tensor([m, 1], Float64Array.from(y))); return K.transpose(typeof p === 'number' ? tensor([1, 1], Float64Array.of(p)) : p); }, 'polyfit(x, y, n)');
Matrix to text
64-bit
mat2str(A): text that can be typed back in to recreate a matrix.
each value written with toPrecision(15) by default
The default of 15 significant digits is not always enough to read back the same 64-bit number, which can need 17. A precision argument changes the digit count. Whole numbers are written in full.
cs-lib-core.js
lines 645 to 647
PartCraft 0.8.1
def('mat2str', (a) => {
const v = dense(a[0]); const p = a[1] === undefined ? 15 : I(a[1], t('script.core.what.mat2str_precision'));
const f = (x) => (Number.isInteger(x) ? String(x) : !Number.isFinite(x) ? (x !== x ? 'NaN' : x > 0 ? 'Inf' : '-Inf') : String(+x.toPrecision(p)));
A display conversion: non-whole numbers are shortened to four digits beyond the whole part unless a precision or format is given. Text made this way does not carry the full value.
cs-lib-core.js
lines 626 to 636
PartCraft 0.8.1
def('num2str', (a) => {
const v = dense(a[0]);
if (isStr(v)) return v;
if (isStr(a[1])) return sprintf(a[1], [v]);
const fmt1 = (x) => {
if (a[1] !== undefined) return fmtG(x, I(a[1], t('script.core.what.num2str_precision')), false, false);
if (!Number.isFinite(x)) return x !== x ? 'NaN' : x > 0 ? 'Inf' : '-Inf';
if (Number.isInteger(x)) return String(x);
const d = Math.max(Math.ceil(Math.log10(Math.abs(x))), 1) + 4;
return fmtG(x, d, false, false);
};
Text to number
64-bit
str2double(s): reads a number from text.
Number(text) for text matching a decimal pattern, else NaN
The text is checked against a pattern and then converted by Number(), which returns the nearest 64-bit number to the decimal written.
cs-lib-core.js
line 656
PartCraft 0.8.1
def('str2double', (a) => { const v = a[0]; const one = (s) => { if (!isStr(s)) return NaN; const t = s.trim().replace(/^\+/, ''); if (/^-?inf$/i.test(t)) return t[0] === '-' ? -Infinity : Infinity; if (/^nan$/i.test(t)) return NaN; return /^[-]?(\d+\.?\d*|\.\d+)([eE][-+]?\d+)?$/.test(t) ? Number(t) : NaN; }; if (v && v.t === 'cell') return out(sizeOf(v), Float64Array.from(v.items, one)); return one(v); }, 'str2double(s)');
In detail
Method
Elementwise maps over tensors. log1p and expm1: three-term series below |x| = 1e-4, else log(1 + x) and exp(x) − 1. Inverse hyperbolic functions from logarithms. Gamma: exact product for whole numbers below 171, reflection below 0.5, Lanczos (g = 7, 9 coefficients) otherwise. Variance by two passes. Polynomial fit by a Vandermonde matrix.
Number types
64-bit numbers and Float64Array; logical results in Uint8Array.
Tolerances
log1p and expm1 switch formulas at |x| = 1e-4. mat2str writes 15 significant digits by default; num2str writes 4 decimals beyond the whole part.
Limits
Near |x| = 1e-4, log1p and expm1 carry about 12 to 13 digits (the dropped series term is x⁴ ÷ 4, about 2.5e-13 of the result; just above the switch, forming 1 + x rounds x at about 1e-12 relative). asinh and atanh have no small-argument form, so their relative error grows as the argument shrinks. For non-whole arguments above about 142 the Lanczos factor t^(x − ½) overflows to Infinity although the gamma value itself is still representable. gamma(171), and so factorial(170), returns Infinity although the factorial of 170 (about 7.26e306) is representable.
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.