The elementwise operators, matrix product, backslash, ranges and sums that a script's own formulas run on.
How the precision is kept
Numeric arrays are Float64Array and logical results are Uint8Array; every operator is one 64-bit operation per element. Squares and square roots take a direct path, other powers and atan2 go through field/detmath.js so every engine gives the same bits. Sums are added in index order without compensation.
Open item. Give the square backslash a scaled singularity test like the one luFactor in cs-lib-linalg.js uses, so a nearly singular matrix is reported instead of solved.
Calculations under this policy
9 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.
Adds, subtracts, multiplies, divides, compares or takes the remainder of two numbers.
mod(x, y) = x − ⌊x ÷ y⌋ × y; rem(x, y) = x − trunc(x ÷ y) × y
Each case is one 64-bit expression. Comparisons are exact (=== and <), with no tolerance. mod and rem subtract a product from x, so when x ÷ y is very large the remainder has few correct digits; a zero divisor returns x.
cs-tensor.js
lines 34 to 46
PartCraft 0.8.1
function scalarOp(op, x, y) {
switch (op) {
case '+': return x + y; case '-': return x - y; case '.*': case '*': return x * y;
case './': case '/': return x / y; case '.\\': case '\\': return y / x;
case '.^': case '^': return pow(x, y);
case '==': return x === y; case '~=': return x !== y; case '<': return x < y; case '<=': return x <= y;
case '>': return x > y; case '>=': return x >= y; case '&': return x !== 0 && y !== 0; case '|': return x !== 0 || y !== 0;
case 'min': return x < y || Number.isNaN(y) ? x : y; case 'max': return x > y || Number.isNaN(y) ? x : y;
case 'mod': return y === 0 ? x : x - Math.floor(x / y) * y; case 'rem': return y === 0 ? x : x - Math.trunc(x / y) * y;
case 'atan2': return datan2(x, y); case 'hypot': return hypot2(x, y);
default: throw new CsError(t('script.tensor.unknown_op', { op }));
}
}
Hypotenuse
64-bit
√(x² + y²) for two numbers.
m × √((x ÷ m)² + (y ÷ m)²), m = max(|x|, |y|)
Both values are divided by the larger one before squaring, so the squares cannot overflow or underflow. Built from + − × ÷ and sqrt only; the code's comment gives the reason that Math.hypot may differ between engines.
cs-tensor.js
lines 47 to 56
PartCraft 0.8.1
/** √(x²+y²) without overflow, from + − × ÷ and sqrt only (Math.hypot is allowed to differ between engines). */
export function hypot2(x, y) {
const a = Math.abs(x), b = Math.abs(y);
if (a === Infinity || b === Infinity) return Infinity;
if (a !== a || b !== b) return NaN;
const m = a > b ? a : b;
if (m === 0) return 0;
const p = a / m, q = b / m;
return m * Math.sqrt(p * p + q * q);
}
Power
64-bit
x raised to the power y.
An exponent of exactly 2 is computed as x × x and exactly 0.5 as Math.sqrt(x), both correctly rounded. Every other exponent calls dpow in field/detmath.js. A negative base with a non-integer exponent raises an error instead of returning NaN.
cs-tensor.js
lines 57 to 64
PartCraft 0.8.1
/** x^y as MATLAB for real results; a negative base with a non-integer exponent is complex → NaN with a clear error. */
function pow(x, y) {
if (y === 2) return x * x;
if (y === 0.5) return Math.sqrt(x);
const r = dpow(x, y);
if (Number.isNaN(r) && x < 0 && !Number.isInteger(y) && Number.isFinite(y)) throw new CsError(t('script.tensor.complex_power', { x, y }));
return r;
}
Elementwise loops
64-bit
Applies an operator to every element of one or two arrays.
One 64-bit operation per element, written into a Float64Array (or a Uint8Array of 0 and 1 for comparisons). max and min return the other operand when one of the two is NaN.
cs-tensor.js
lines 66 to 115
PartCraft 0.8.1
/** Loop kernels: same-length arrays, or one side scalar. */
function kernel(op, a, b, x, y, n, out) {
// a/b typed arrays or null (then x/y scalars)
if (a && b) {
switch (op) {
case '+': for (let i = 0; i < n; i++) out[i] = a[i] + b[i]; return;
case '-': for (let i = 0; i < n; i++) out[i] = a[i] - b[i]; return;
case '.*': for (let i = 0; i < n; i++) out[i] = a[i] * b[i]; return;
case './': for (let i = 0; i < n; i++) out[i] = a[i] / b[i]; return;
case '.\\': for (let i = 0; i < n; i++) out[i] = b[i] / a[i]; return;
case '==': for (let i = 0; i < n; i++) out[i] = a[i] === b[i] ? 1 : 0; return;
case '~=': for (let i = 0; i < n; i++) out[i] = a[i] !== b[i] ? 1 : 0; return;
case '<': for (let i = 0; i < n; i++) out[i] = a[i] < b[i] ? 1 : 0; return;
case '<=': for (let i = 0; i < n; i++) out[i] = a[i] <= b[i] ? 1 : 0; return;
case '>': for (let i = 0; i < n; i++) out[i] = a[i] > b[i] ? 1 : 0; return;
case '>=': for (let i = 0; i < n; i++) out[i] = a[i] >= b[i] ? 1 : 0; return;
case 'max': for (let i = 0; i < n; i++) { const p = a[i], q = b[i]; out[i] = p > q || q !== q ? p : q; } return;
case 'min': for (let i = 0; i < n; i++) { const p = a[i], q = b[i]; out[i] = p < q || q !== q ? p : q; } return;
default: for (let i = 0; i < n; i++) out[i] = +scalarOp(op, a[i], b[i]); return;
}
}
if (a) { // array op scalar
switch (op) {
case '+': for (let i = 0; i < n; i++) out[i] = a[i] + y; return;
case '-': for (let i = 0; i < n; i++) out[i] = a[i] - y; return;
case '.*': for (let i = 0; i < n; i++) out[i] = a[i] * y; return;
case './': for (let i = 0; i < n; i++) out[i] = a[i] / y; return;
case '.^': if (y === 2) { for (let i = 0; i < n; i++) out[i] = a[i] * a[i]; return; } if (y === 1) { for (let i = 0; i < n; i++) out[i] = a[i]; return; }
for (let i = 0; i < n; i++) out[i] = pow(a[i], y); return;
case '==': for (let i = 0; i < n; i++) out[i] = a[i] === y ? 1 : 0; return;
case '~=': for (let i = 0; i < n; i++) out[i] = a[i] !== y ? 1 : 0; return;
case '<': for (let i = 0; i < n; i++) out[i] = a[i] < y ? 1 : 0; return;
case '<=': for (let i = 0; i < n; i++) out[i] = a[i] <= y ? 1 : 0; return;
case '>': for (let i = 0; i < n; i++) out[i] = a[i] > y ? 1 : 0; return;
case '>=': for (let i = 0; i < n; i++) out[i] = a[i] >= y ? 1 : 0; return;
case 'max': if (y === y) { for (let i = 0; i < n; i++) { const p = a[i]; out[i] = p > y ? p : y; } } else for (let i = 0; i < n; i++) out[i] = a[i]; return;
case 'min': if (y === y) { for (let i = 0; i < n; i++) { const p = a[i]; out[i] = p < y ? p : y; } } else for (let i = 0; i < n; i++) out[i] = a[i]; return;
default: for (let i = 0; i < n; i++) out[i] = +scalarOp(op, a[i], y); return;
}
}
switch (op) { // scalar op array
case '+': for (let i = 0; i < n; i++) out[i] = x + b[i]; return;
case '-': for (let i = 0; i < n; i++) out[i] = x - b[i]; return;
case '.*': for (let i = 0; i < n; i++) out[i] = x * b[i]; return;
case './': for (let i = 0; i < n; i++) out[i] = x / b[i]; return;
case 'max': if (x === x) { for (let i = 0; i < n; i++) { const q = b[i]; out[i] = x > q || q !== q ? x : q; } } else for (let i = 0; i < n; i++) out[i] = b[i]; return;
case 'min': if (x === x) { for (let i = 0; i < n; i++) { const q = b[i]; out[i] = x < q || q !== q ? x : q; } } else for (let i = 0; i < n; i++) out[i] = b[i]; return;
default: for (let i = 0; i < n; i++) out[i] = +scalarOp(op, x, b[i]); return;
}
}
Matrix product
64-bit
The product of two dense matrices in a script expression.
C_ij = Σ_p A_ip × B_pj
A plain 64-bit sum of products in order of p; terms where the right-hand entry is exactly zero are skipped. No compensated summation.
cs-tensor.js
lines 212 to 235
PartCraft 0.8.1
export function matmul(a, b) {
if ((typeof a === 'number' || typeof a === 'boolean') || (typeof b === 'number' || typeof b === 'boolean')) return binop('.*', a, b);
const A = view(a, '*'), B = view(b, '*');
if (!A.d || !B.d) return binop('.*', a, b);
if (A.s.length > 2 || B.s.length > 2) throw new CsError(t('script.tensor.mtimes_2d'));
const [m, k] = A.s, [k2, n] = B.s;
if (k !== k2) throw new CsError(t('script.tensor.mtimes_inner', { a: shapeText(A.s), b: shapeText(B.s) }));
checkSize(m * n);
const out = new Float64Array(m * n), ad = A.d, bd = B.d;
if (m === 1) { // row × matrix: dot products down columns
for (let j = 0; j < n; j++) { let s = 0; const o = j * k; for (let p = 0; p < k; p++) s += ad[p] * bd[o + p]; out[j] = s; }
} else {
for (let j = 0; j < n; j++) {
const oc = j * m;
for (let p = 0; p < k; p++) {
const bv = bd[p + j * k];
if (bv === 0) continue;
const ac = p * m;
for (let i = 0; i < m; i++) out[oc + i] += ad[ac + i] * bv;
}
}
}
return m === 1 && n === 1 ? out[0] : tensor([m, n], out);
}
Backslash for dense matrices
64-bit
Solves A x = b, or the least-squares problem when A is not square.
Square: LU with partial pivoting; the only singularity test is a pivot that is exactly zero, so a nearly singular matrix is solved without a warning. More rows than columns: Householder QR with the reflector sign chosen against the diagonal; a rank error is raised only when a column norm is exactly zero. Fewer rows than columns: solved through A Aᵀ, which squares the condition number.
cs-tensor.js
lines 237 to 283
PartCraft 0.8.1
/** A\B for dense A: LU with partial pivoting (square), else least squares by Householder QR. */
export function solve(a, b) {
const A = view(a, '\\'), B = view(b, '\\');
if (!A.d && A.x !== undefined && A.s[0] === 1 && A.s[1] === 1) return binop('./', b, A.x);
const [m, n] = A.s;
if (A.s.length > 2 || B.s.length > 2) throw new CsError(t('script.tensor.mldivide_2d'));
const bd = B.d || Float64Array.of(B.x), nr = B.s[0], nc = B.d ? B.s[1] : 1;
if (nr !== m) throw new CsError(t('script.tensor.mldivide_rows', { a: shapeText(A.s), b: shapeText(B.s) }));
if (m === n) {
const L = Float64Array.from(A.d), piv = new Int32Array(n);
for (let k = 0; k < n; k++) {
let p = k, best = Math.abs(L[k + k * n]);
for (let i = k + 1; i < n; i++) { const v = Math.abs(L[i + k * n]); if (v > best) { best = v; p = i; } }
piv[k] = p;
if (best === 0) throw new CsError(t('script.tensor.singular'));
if (p !== k) for (let j = 0; j < n; j++) { const t = L[k + j * n]; L[k + j * n] = L[p + j * n]; L[p + j * n] = t; }
const d = L[k + k * n];
for (let i = k + 1; i < n; i++) L[i + k * n] /= d;
for (let j = k + 1; j < n; j++) { const f = L[k + j * n]; if (f !== 0) for (let i = k + 1; i < n; i++) L[i + j * n] -= L[i + k * n] * f; }
}
const X = Float64Array.from(bd);
for (let c = 0; c < nc; c++) {
const x = X.subarray(c * n, (c + 1) * n);
for (let k = 0; k < n; k++) { const p = piv[k]; if (p !== k) { const t = x[k]; x[k] = x[p]; x[p] = t; } }
for (let k = 0; k < n; k++) { const v = x[k]; if (v !== 0) for (let i = k + 1; i < n; i++) x[i] -= L[i + k * n] * v; }
for (let k = n - 1; k >= 0; k--) { x[k] /= L[k + k * n]; const v = x[k]; if (v !== 0) for (let i = 0; i < k; i++) x[i] -= L[i + k * n] * v; }
}
return nc === 1 && n === 1 ? X[0] : tensor([n, nc], X);
}
// least squares (m > n) or minimum-norm-ish (m < n, via normal equations of A') — QR for m ≥ n
if (m < n) { const At = transpose(a); const AAt = matmul(a, At); return matmul(At, solve(AAt, b)); }
const R = Float64Array.from(A.d), X = Float64Array.from(bd);
for (let k = 0; k < n; k++) {
let s = 0; for (let i = k; i < m; i++) s += R[i + k * m] ** 2;
let nrm = Math.sqrt(s); if (nrm === 0) throw new CsError(t('script.tensor.rank'));
if (R[k + k * m] > 0) nrm = -nrm;
const v0 = R[k + k * m] - nrm; R[k + k * m] = nrm;
// apply H = I - 2 v v'/(v'v) with v = [v0; R(k+1:m,k)]
let vtv = v0 * v0; for (let i = k + 1; i < m; i++) vtv += R[i + k * m] ** 2;
if (vtv === 0) continue;
for (let j = k + 1; j < n; j++) { let d = v0 * R[k + j * m]; for (let i = k + 1; i < m; i++) d += R[i + k * m] * R[i + j * m]; const f = 2 * d / vtv; R[k + j * m] -= f * v0; for (let i = k + 1; i < m; i++) R[i + j * m] -= f * R[i + k * m]; }
for (let c = 0; c < nc; c++) { const o = c * m; let d = v0 * X[o + k]; for (let i = k + 1; i < m; i++) d += R[i + k * m] * X[o + i]; const f = 2 * d / vtv; X[o + k] -= f * v0; for (let i = k + 1; i < m; i++) X[o + i] -= f * R[i + k * m]; }
}
const out = new Float64Array(n * nc);
for (let c = 0; c < nc; c++) for (let k = n - 1; k >= 0; k--) { let s = X[c * m + k]; for (let j = k + 1; j < n; j++) s -= R[k + j * m] * out[c * n + j]; out[c * n + k] = s / R[k + k * m]; }
return nc === 1 && n === 1 ? out[0] : tensor([n, nc], out);
}
Number of elements in a range
64-bit
How many values a:s:b has.
n = ⌊(b − a) ÷ s + max(3ε × max(|a|, |b|) ÷ |s|, 1e-10)⌋ + 1
A small allowance is added before rounding down so that a range such as 0:0.1:1 keeps its last value. The allowance is at least 1e-10 of a step, so an end point that falls short of a whole step by less than that is still counted.
cs-tensor.js
lines 285 to 292
PartCraft 0.8.1
/** a:s:b as MATLAB (1×n row; empty 1×0 when it runs the wrong way). */
export function rangeCount(a, s, b) {
if (!Number.isFinite(a) || !Number.isFinite(s) || !Number.isFinite(b)) { if (Number.isNaN(a) || Number.isNaN(s) || Number.isNaN(b)) throw new CsError(t('script.tensor.range_nan')); throw new CsError(t('script.tensor.range_inf')); }
if (s === 0 || (b - a) / s < 0) return 0;
const q = (b - a) / s;
const tol = 3 * Number.EPSILON * Math.max(Math.abs(a), Math.abs(b)) / Math.abs(s);
return Math.floor(q + Math.max(tol, 1e-10)) + 1;
}
Values of a range
64-bit
The values a, a + s, a + 2s … up to b.
x_k = a + k × s
Each value is computed from a and k directly, not by repeated addition, so rounding does not build up along the range. For non-integer ranges the last value is replaced by b when it is within 1e-10 × max(1, |b|) of it.
cs-tensor.js
lines 293 to 300
PartCraft 0.8.1
export function range(a, s, b) {
const n = rangeCount(a, s, b);
checkSize(n);
const out = new Float64Array(n);
for (let k = 0; k < n; k++) out[k] = a + k * s;
if (n > 0 && Number.isInteger(a) && Number.isInteger(s)) { /* exact */ } else if (n > 1) { const last = a + (n - 1) * s; if (Math.abs(last - b) < 1e-10 * Math.max(1, Math.abs(b))) out[n - 1] = b; }
return tensor([1, n], out);
}
Sum, product, mean, maximum and minimum
64-bit
Reduces an array along one dimension.
mean = (Σ xₖ) ÷ n
Sums and products run in index order in 64-bit without compensation, so a long sum of values with mixed signs or very different sizes can lose low digits. Maximum and minimum compare exactly and skip NaN.
cs-tensor.js
lines 664 to 716
PartCraft 0.8.1
/**
* Reduce along dim (0-based; null = first non-singleton). kind: sum prod max min any all mean.
* Returns { value, index } (index = 1-based positions for max/min).
*/
export function reduce(kind, v, dim = null) {
const T = typeof v === 'number' || typeof v === 'boolean' ? tensor([1, 1], Float64Array.of(+v)) : typeof v === 'string' ? charCodes(v) : v;
if (!isTensor(T)) throw new CsError(t('script.tensor.fn_value', { fn: kind, got: describe(v) }));
let shape = T.shape;
if (dim === null) {
if (shape.length === 2 && shape[0] === 0 && shape[1] === 0) { // sum([]) = 0, max([]) = []
if (kind === 'sum' || kind === 'any') return { value: kind === 'any' ? false : 0 };
if (kind === 'prod' || kind === 'all') return { value: kind === 'all' ? true : 1 };
if (kind === 'mean') return { value: NaN };
return { value: empty(), index: empty() };
}
dim = defaultDim(shape);
}
const full = shape.slice(); while (full.length <= dim) full.push(1);
const len = full[dim];
let inner = 1; for (let k = 0; k < dim; k++) inner *= full[k];
let outer = 1; for (let k = dim + 1; k < full.length; k++) outer *= full[k];
const outShape = full.slice(); outShape[dim] = (kind === 'max' || kind === 'min') && len === 0 ? 0 : 1;
const n = inner * outer;
const d = T.data;
const logical = kind === 'any' || kind === 'all';
const out = logical ? new Uint8Array(n) : new Float64Array(n);
const idx = kind === 'max' || kind === 'min' ? new Float64Array(n) : null;
for (let o = 0; o < outer; o++) {
for (let i = 0; i < inner; i++) {
const base = o * inner * len + i, oi = o * inner + i;
switch (kind) {
case 'sum': case 'mean': { let s = 0; for (let k = 0; k < len; k++) s += d[base + k * inner]; out[oi] = kind === 'mean' ? s / len : s; break; }
case 'prod': { let s = 1; for (let k = 0; k < len; k++) s *= d[base + k * inner]; out[oi] = s; break; }
case 'any': { let s = 0; for (let k = 0; k < len; k++) if (d[base + k * inner] !== 0 && d[base + k * inner] === d[base + k * inner]) { s = 1; break; } out[oi] = s; break; }
case 'all': { let s = 1; for (let k = 0; k < len; k++) if (d[base + k * inner] === 0) { s = 0; break; } out[oi] = s; break; }
case 'max': case 'min': {
let best = NaN, bi = 0;
for (let k = 0; k < len; k++) {
const x = d[base + k * inner];
if (x !== x) continue;
if (best !== best || (kind === 'max' ? x > best : x < best)) { best = x; bi = k; }
}
out[oi] = best; idx[oi] = bi + 1;
break;
}
default: throw new CsError(t('script.tensor.unknown_reduction', { fn: kind }));
}
}
}
if (len === 0 && (kind === 'max' || kind === 'min')) return { value: tensor(outShape, new Float64Array(0)), index: tensor(outShape, new Float64Array(0)) };
const wrap = (a) => (a.length === 1 && outShape.length <= 2 ? (logical ? a[0] !== 0 : a[0]) : tensor(outShape, a));
return { value: wrap(out), index: idx ? wrap(idx) : undefined };
}
In detail
Method
Loop kernels per operator; LU with partial pivoting for a square backslash, Householder QR for more rows than columns, normal equations of the transpose for fewer rows than columns.
Number types
64-bit numbers, Float64Array; Uint8Array for logical results; Int32Array for pivots and indexes.
Tolerances
Range count: the larger of 3 × 2^-52 × max(|a|, |b|) ÷ |s| and 1e-10 is added before rounding down. Range end: the last element is set to b when it is within 1e-10 × max(1, |b|). Backslash: a pivot is singular only when exactly zero.
Limits
The square backslash has no near-singular warning. The underdetermined case forms A Aᵀ, which squares the condition number.
Approach: 64-bit, direct
Every value is a 64-bit number from input to result. Formulas are closed form or a direct solve, and every comparison uses a stated tolerance.