JavaScript: Math, 110–114

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// Terms: whole part, subtract, invert, repeat.
function terms(x, count) {
  const out = [];
  let v = x;
  for (let i = 0; i < count; i++) {
    const a = Math.floor(v);
    out.push(a);
    const frac = v - a;
    if (frac === 0) break;
    v = 1 / frac;
  }
  return out;
}

// Convergents: p / q, each built from the two before it.
function convergents(list) {
  let pPrev = 1;
  let p = list[0];
  let qPrev = 0;
  let q = 1;
  const out = [[p, q]];
  for (let i = 1; i < list.length; i++) {
    const a = list[i];
    const pNext = a * p + pPrev;
    const qNext = a * q + qPrev;
    pPrev = p;
    p = pNext;
    qPrev = q;
    q = qNext;
    out.push([p, q]);
  }
  return out;
}

// The best fraction with denominator at most maxQ, found by looking at all of them.
function bestFraction(x, maxQ) {
  let best = [1, 1];
  let bestError = Infinity;
  for (let q = 1; q <= maxQ; q++) {
    const p = Math.round(x * q);
    const error = Math.abs(p / q - x);
    if (error < bestError) {
      bestError = error;
      best = [p, q];
    }
  }
  return best;
}

// Backwards: grow convergents until one is within tol of x.
function toFraction(x, tol) {
  const list = terms(x, 40);
  const all = convergents(list);
  for (const [p, q] of all) {
    if (Math.abs(p / q - x) <= tol) return [p, q];
  }
  return all[all.length - 1];
}

const piTerms = terms(Math.PI, 12);
const piConvergents = convergents(piTerms).slice(0, 5);
const errors = piConvergents.map(([p, q]) => p / q - Math.PI);
const golden = convergents(terms((1 + Math.sqrt(5)) / 2, 10)).map(([p, q]) => p + '/' + q);
const root2 = terms(Math.SQRT2, 8);
const exactValue = terms(2.75, 4);

console.log('pi terms      :', piTerms);
console.log('convergents   :', piConvergents.map(([p, q]) => p + '/' + q));
console.log('errors        :', errors);
console.log('best q <= 100 :', bestFraction(Math.PI, 100));
console.log('best q <= 1000:', bestFraction(Math.PI, 1000));
console.log('0.1 back      :', toFraction(0.1, 1e-12), 'and 0.07 gives', toFraction(0.07, 1e-12));
console.log('the double 0.1:', 3602879701896397, '/ 2 ** 55');
console.log('golden        :', golden);
console.log('root 2 terms  :', root2);
console.log('2.75 terms    :', exactValue, '<- the last one is not real');

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
75const shape = (m) => [m.length, m[0].length];

function multiply(a, b) {
  const [rowsA, colsA] = shape(a);
  const [rowsB, colsB] = shape(b);
  if (colsA !== rowsB) {
    throw new RangeError('cannot multiply ' + rowsA + 'x' + colsA + ' by ' + rowsB + 'x' + colsB);
  }
  const out = [];
  for (let i = 0; i < rowsA; i++) {
    const row = [];
    for (let j = 0; j < colsB; j++) {
      let total = 0;
      for (let k = 0; k < colsA; k++) total += a[i][k] * b[k][j];
      row.push(total);
    }
    out.push(row);
  }
  return out;
}

const identity = (n) =>
  Array.from({ length: n }, (_, i) => Array.from({ length: n }, (_, j) => (i === j ? 1 : 0)));

// Drop one row and one column.
const minor = (m, row, col) =>
  m.filter((_, i) => i !== row).map((r) => r.filter((_, j) => j !== col));

function det(m) {
  const n = m.length;
  if (n === 1) return m[0][0];
  if (n === 2) return m[0][0] * m[1][1] - m[0][1] * m[1][0];
  let total = 0;
  for (let j = 0; j < n; j++) total += (j % 2 ? -1 : 1) * m[0][j] * det(minor(m, 0, j));
  return total;
}

function inverse2(m) {
  const d = det(m);
  if (d === 0) return null;
  return [[m[1][1] / d, -m[0][1] / d], [-m[1][0] / d, m[0][0] / d]];
}

// The shoelace formula, written as a sum of 2 by 2 determinants.
function area(points) {
  let total = 0;
  for (let i = 0; i < points.length; i++) {
    const p = points[i];
    const q = points[(i + 1) % points.length];
    total += det([[p[0], q[0]], [p[1], q[1]]]);
  }
  return total / 2;
}

const A = [[1, 2, 3], [4, 5, 6], [7, 8, 10]];
const B = [[2, 0, 1], [1, 3, 0], [0, 1, 4]];
const singular = [[1, 2, 3], [4, 5, 6], [7, 8, 9]];
const fourByFour = [[1, 0, 2, -1], [3, 0, 0, 5], [2, 1, 4, -3], [1, 0, 5, 0]];
const M = [[4, 7], [2, 6]];
const F = [[0.1, 0.2], [0.3, 0.4]];
const G = [[0.5, 0.6], [0.7, 0.8]];
const plot = [[0, 0], [40, 0], [40, 25], [15, 40], [0, 25]];

let mismatch = 'no error';
try { multiply(A, M); } catch (error) { mismatch = error.constructor.name + ': ' + error.message; }

console.log('A B          :', multiply(A, B));
console.log('B A          :', multiply(B, A), '<- a different matrix');
console.log('dets         :', det(A), det(B), 'product', det(A) * det(B), 'det of AB', det(multiply(A, B)));
console.log('4 by 4       :', det(fourByFour), 'transpose of A', det([0, 1, 2].map((j) => A.map((r) => r[j]))));
console.log('singular     :', det(singular), 'inverse of a singular 2 by 2', inverse2([[1, 2], [2, 4]]));
console.log('M times M^-1 :', multiply(M, inverse2(M)), '<- not quite the identity');
console.log('float dets   :', det(multiply(F, G)), det(F) * det(G));
console.log('shape error  :', mismatch);
console.log('plot         :', area(plot), 'square metres, reversed', area(plot.slice().reverse()));

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
47const residual = (a, b, c, x) => a * x * x + b * x + c;

// The formula as taught.
function naiveRoots(a, b, c) {
  const d = b * b - 4 * a * c;
  const root = Math.sqrt(d);
  return [(-b + root) / (2 * a), (-b - root) / (2 * a)];
}

// The stable formula. Add the square root to b with b's own sign, so nothing
// cancels, then get the second root from the product of the roots, c / a.
function stableRoots(a, b, c) {
  const d = b * b - 4 * a * c;
  if (d < 0) return [];
  const root = Math.sqrt(d);
  const q = -(b + Math.sign(b || 1) * root) / 2;
  return [q / a, q === 0 ? 0 : c / q];
}

// Divide through by the largest coefficient first, so b * b cannot overflow.
function scaledRoots(a, b, c) {
  const s = Math.max(Math.abs(a), Math.abs(b), Math.abs(c));
  return stableRoots(a / s, b / s, c / s);
}

const naiveBig = naiveRoots(1, 1e8, 1);
const stableBig = stableRoots(1, 1e8, 1);
const naiveResiduals = naiveBig.map((x) => residual(1, 1e8, 1, x));
const stableResiduals = stableBig.map((x) => residual(1, 1e8, 1, x));

// Vieta: the roots must add to -b / a and multiply to c / a.
const vietaNaive = naiveBig[0] * naiveBig[1] - 1;
const vietaStable = stableBig[0] * stableBig[1] - 1;

const overflowed = naiveRoots(1, 1e200, 1);
const stableOverflowed = stableRoots(1, 1e200, 1);
const rescued = scaledRoots(1, 1e200, 1);

console.log('x2 - 5x + 6      :', naiveRoots(1, -5, 6), stableRoots(1, -5, 6));
console.log('x2 - 2x + 1      :', naiveRoots(1, -2, 1), 'discriminant', 4 - 4);
console.log('x2 + x + 1       :', naiveRoots(1, 1, 1), 'against a list of', stableRoots(1, 1, 1).length, 'real roots');
console.log('b = 1e8, naive   :', naiveBig, 'residuals', naiveResiduals);
console.log('b = 1e8, stable  :', stableBig, 'residuals', stableResiduals);
console.log('b = 1e9, naive   :', naiveRoots(1, 1e9, 1), '<- the small root is gone');
console.log('a = 0, naive     :', naiveRoots(0, 2, -4), 'stable', stableRoots(0, 2, -4));
console.log('b = 1e200        :', overflowed, 'discriminant', 1e200 * 1e200);
console.log('stable, scaled   :', stableOverflowed, rescued);

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
67function trapezium(f, a, b, n) {
  const h = (b - a) / n;
  let total = (f(a) + f(b)) / 2;
  for (let i = 1; i < n; i++) total += f(a + i * h);
  return total * h;
}

function midpoint(f, a, b, n) {
  const h = (b - a) / n;
  let total = 0;
  for (let i = 0; i < n; i++) total += f(a + (i + 0.5) * h);
  return total * h;
}

function simpson(f, a, b, n) {
  if (n % 2 !== 0) throw new RangeError('Simpson needs an even number of strips');
  const h = (b - a) / n;
  let total = f(a) + f(b);
  for (let i = 1; i < n; i++) total += f(a + i * h) * (i % 2 ? 4 : 2);
  return (total * h) / 3;
}

const square = (x) => x * x;
const cube = (x) => x ** 3;
const line = (x) => 2 * x + 1;
const reciprocal = (x) => 1 / x;
const spike = (x) => 1 / Math.sqrt(x);          // infinite at x = 0

const trapErrors = [1, 2, 6, 12].map((n) => trapezium(square, 0, 3, n) - 9);
const T = trapezium(square, 0, 3, 6);
const Mid = midpoint(square, 0, 3, 6);

const t10 = trapezium(reciprocal, 1, 2, 10);
const t20 = trapezium(reciprocal, 1, 2, 20);
const s20 = simpson(reciprocal, 1, 2, 20);
const richardson = (4 * t20 - t10) / 3;
const halvingRatio = (t10 - Math.log(2)) / (t20 - Math.log(2));

// How many trapezium strips does it take to match Simpson's twenty?
let crossover = 0;
for (let n = 2; n <= 600 && crossover === 0; n++) {
  if (Math.abs(trapezium(reciprocal, 1, 2, n) - Math.log(2)) <= Math.abs(s20 - Math.log(2))) crossover = n;
}

// Keep refining and watch the tenfold gain in accuracy decay.
const ratios = [];
let previous = 0;
for (const n of [100, 1000, 10000, 100000]) {
  const error = Math.abs(trapezium(reciprocal, 1, 2, n) - Math.log(2));
  if (previous !== 0) ratios.push(previous / error);
  previous = error;
}

let oddStrips = 'no error';
try { simpson(square, 0, 3, 3); } catch (error) { oddStrips = error.constructor.name; }

console.log('x2 trapezium  :', [1, 2, 6, 12].map((n) => trapezium(square, 0, 3, n)), 'exact 9');
console.log('errors        :', trapErrors, '<- each a quarter of the last');
console.log('midpoint      :', Mid, 'trapezium', T, '(T + 2M) / 3 =', (T + 2 * Mid) / 3);
console.log('Simpson x2, x3:', simpson(square, 0, 3, 2), simpson(cube, 0, 2, 2), 'exact 9 and 4');
console.log('a straight line:', trapezium(line, 0, 4, 1), trapezium(line, 0, 4, 7), trapezium(line, 0, 4, 8));
console.log('1/x, ln 2     :', t10, t20, s20, Math.log(2));
console.log('richardson    :', richardson, 'against Simpson', s20, 'gap', richardson - s20);
console.log('halving ratio :', halvingRatio, 'refinement ratios', ratios);
console.log('crossover     :', crossover, 'trapezium strips to match Simpson on 20');
console.log('spike         :', trapezium(spike, 0, 1, 100), 'midpoint', midpoint(spike, 0, 1, 100), 'exact 2');
console.log('odd strips    :', oddStrips);

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
86const digits = (s) => [...s].map(Number);

// Luhn: double every second digit counting from the right.
function luhnTotal(s) {
  const d = digits(s);
  let total = 0;
  for (let i = 0; i < d.length; i++) {
    const fromRight = d.length - 1 - i;
    let v = d[i];
    if (fromRight % 2 === 1) {
      v *= 2;
      if (v > 9) v -= 9;
    }
    total += v;
  }
  return total;
}
const luhnValid = (s) => luhnTotal(s) % 10 === 0;
const luhnDigit = (body) => (10 - (luhnTotal(body + '0') % 10)) % 10;

// NUBAN: weights 3, 7, 3 repeating over the bank code and the serial number.
const WEIGHTS = [3, 7, 3, 3, 7, 3, 3, 7, 3, 3, 7, 3];
function nubanDigit(bank, serial) {
  const d = digits(bank + serial);
  const total = d.reduce((sum, v, i) => sum + v * WEIGHTS[i], 0);
  return (10 - (total % 10)) % 10;
}

// ISBN-10: weights 10 down to 1, modulus 11, and a check value of 10 prints as X.
function isbnDigit(body) {
  const total = digits(body).reduce((sum, v, i) => sum + v * (10 - i), 0);
  const check = (11 - (total % 11)) % 11;
  return check === 10 ? 'X' : String(check);
}
const isbnValid = (s) =>
  [...s].map((c) => (c === 'X' ? 10 : Number(c))).reduce((sum, v, i) => sum + v * (10 - i), 0) % 11 === 0;

// The control: add the digits, no weights at all.
const sumDigit = (body) => (10 - (digits(body).reduce((a, b) => a + b, 0) % 10)) % 10;
const sumValid = (s) => sumDigit(s.slice(0, -1)) === Number(s[s.length - 1]);

const card = '4111111111111111';

// Every single-digit change to the card: 16 positions, 9 wrong digits each.
let missedSingle = 0;
for (let i = 0; i < card.length; i++) {
  for (let d = 0; d <= 9; d++) {
    if (String(d) === card[i]) continue;
    if (luhnValid(card.slice(0, i) + d + card.slice(i + 1))) missedSingle++;
  }
}

// Every ordered pair of distinct digits, swapped in the same two positions.
function missedSwaps(check, valid) {
  const missed = [];
  for (let a = 0; a <= 9; a++) {
    for (let b = 0; b <= 9; b++) {
      if (a === b) continue;
      const body = '0306' + a + b + '15';
      const full = body + check(body);
      const swapped = full.slice(0, 4) + b + a + full.slice(6);
      if (valid(swapped)) missed.push(a + '' + b);
    }
  }
  return missed;
}

const bank = '011';
const serial = '123456789';
const account = serial + nubanDigit(bank, serial);
// Which adjacent swaps inside the serial does NUBAN see? Only those across
// positions whose weights differ.
const nubanSwaps = [0, 1, 2, 3, 4, 5, 6, 7].map((i) => {
  const swapped = serial.slice(0, i) + serial[i + 1] + serial[i] + serial.slice(i + 2);
  return nubanDigit(bank, swapped) !== nubanDigit(bank, serial);
});

console.log('luhn 79927398713 :', luhnValid('79927398713'), 'total', luhnTotal('79927398713'));
console.log('luhn test card   :', luhnValid(card), 'check digit for 411111111111111 is', luhnDigit('411111111111111'));
console.log('single-digit miss:', missedSingle, 'of', card.length * 9);
console.log('luhn swaps missed:', missedSwaps(luhnDigit, luhnValid));
console.log('isbn swaps missed:', missedSwaps(isbnDigit, isbnValid));
console.log('sum swaps missed :', missedSwaps(sumDigit, sumValid).length, 'of 90');
console.log('nuban account    :', account, 'from bank', bank, 'and serial', serial);
console.log('nuban swaps seen :', nubanSwaps, 'weights', WEIGHTS.slice(3));
console.log('isbn             :', isbnDigit('030640615'), isbnDigit('043942089'), isbnValid('097522980X'));