JavaScript: Math, 80–109

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15const rad = (deg) => (deg * Math.PI) / 180;
const deg = (r) => (r * 180) / Math.PI;
const round12 = (x) => Math.round(x * 1e12) / 1e12;

const angles = [0, 30, 45, 60, 90, 180, 360];

console.log('30 degrees in radians:', rad(30), 'and Math.PI / 6 is', Math.PI / 6);
console.log('sin 30               :', Math.sin(rad(30)), '<- not 0.5');
console.log('cos 60               :', Math.cos(rad(60)), '<- not 0.5 either');
console.log('tan 45               :', Math.tan(rad(45)), '<- not 1');
console.log('rounded to 12 places :', [Math.sin(rad(30)), Math.cos(rad(60)), Math.tan(rad(45))].map(round12));
console.log('sin 90, cos 90       :', Math.sin(rad(90)), Math.cos(rad(90)));
console.log('sin PI               :', Math.sin(Math.PI));
console.log('tan PI / 2           :', Math.tan(Math.PI / 2), 'finite:', Number.isFinite(Math.tan(Math.PI / 2)));
console.log('round trip           :', angles.map((d) => deg(rad(d))));

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19// JavaScript has no modulo operator. % is a truncated remainder: it is
// whatever is left after integer division that rounds TOWARDS ZERO.
const remainder = (a, n) => a % n;

// A real modulo is never negative for positive n. Two % calls get there.
const mod = (a, n) => ((a % n) + n) % n;

// Twelve-hour clock. Hour 0 is really hour 12.
const hoursBefore = (hour, back) => mod(hour - back, 12) || 12;

console.log('7 % 3, -7 % 3   :', 7 % 3, -7 % 3);
console.log('7 % -3, -7 % -3 :', 7 % -3, -7 % -3);
console.log('-6 % 3          :', -6 % 3, 'is it minus zero?', Object.is(-6 % 3, -0));
console.log('5.5 % 2         :', 5.5 % 2, '<- % is not integers-only');
console.log('7 % 0           :', 7 % 0);
console.log('Infinity % 3    :', Infinity % 3, 'and 3 % Infinity:', 3 % Infinity);
console.log('mod(-7, 3)      :', mod(-7, 3), 'against -7 % 3 =', -7 % 3);
console.log('mod(-1, 12)     :', mod(-1, 12), '<- one hour before midnight');
console.log('three before one:', hoursBefore(1, 3));

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15const mod = (a, n) => ((a % n) + n) % n;

const idivTrunc = (a, b) => Math.trunc(a / b);
const idivFloor = (a, b) => Math.floor(a / b);
const idivBits = (a, b) => (a / b) | 0;
const idivShift = (a, b) => (a / b) >> 0;

console.log('7 / 3  :', idivTrunc(7, 3), idivFloor(7, 3), idivBits(7, 3), idivShift(7, 3));
console.log('-7 / 3 :', idivTrunc(-7, 3), idivFloor(-7, 3), idivBits(-7, 3), idivShift(-7, 3));
console.log('trunc pairs with %      :', idivTrunc(-7, 3) * 3 + (-7 % 3));
console.log('floor pairs with mod    :', idivFloor(-7, 3) * 3 + mod(-7, 3));
console.log('trunc(-0.5), floor(-0.5):', Math.trunc(-0.5), Math.floor(-0.5));
console.log('3000000000 | 0          :', 3000000000 | 0, 'against', Math.trunc(3000000000));
console.log('1e10 | 0                :', 1e10 | 0);
console.log('NaN | 0, Infinity | 0   :', NaN | 0, Infinity | 0);

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
33const gcd = (a, b) => {
  a = Math.abs(a);
  b = Math.abs(b);
  while (b) [a, b] = [b, a % b];
  return a;
};

// Divide first: a * b overflows the safe range long before the lcm does.
const lcm = (a, b) => (a === 0 || b === 0 ? 0 : Math.abs(a / gcd(a, b)) * Math.abs(b));
const lcmNaive = (a, b) => Math.abs(a * b) / gcd(a, b);

const reduceFraction = (n, d) => {
  const g = gcd(n, d);
  return [n / g, d / g];
};

// BigInt versions, for when the answer is bigger than a double can hold exactly.
const gcdBig = (a, b) => {
  if (a < 0n) a = -a;
  if (b < 0n) b = -b;
  while (b) [a, b] = [b, a % b];
  return a;
};
const lcmBig = (a, b) => (a / gcdBig(a, b)) * b;

console.log('gcd(48, 18)      :', gcd(48, 18));
console.log('gcd with 0 or -  :', gcd(0, 5), gcd(5, 0), gcd(0, 0), gcd(-48, 18));
console.log('lcm(4, 6)        :', lcm(4, 6), 'lcm(21, 6):', lcm(21, 6));
console.log('gcd * lcm = a * b:', gcd(21, 6) * lcm(21, 6), 21 * 6);
console.log('84 / 126 reduces :', reduceFraction(84, 126));
console.log('lists            :', [12, 18, 24].reduce(gcd), [4, 6, 10].reduce(lcm));
console.log('big pair, double :', lcm(123456789, 987654321), lcmNaive(123456789, 987654321));
console.log('big pair, BigInt :', String(lcmBig(123456789n, 987654321n)));

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20const radixResult = (r) => {
  try {
    return (255).toString(r);
  } catch (error) {
    return error.constructor.name;
  }
};

const toBinary8 = (n) => n.toString(2).padStart(8, '0');
const roundTrip = (n, radix) => parseInt(n.toString(radix), radix);

console.log('255 in four bases:', (255).toString(16), (255).toString(2), (255).toString(8), (255).toString(36));
console.log('-255 in hex      :', (-255).toString(16), '<- sign and magnitude');
console.log('parseInt with 16 :', parseInt('ff', 16), parseInt('0xff', 16), parseInt('0xff'));
console.log('parseInt bare ff :', parseInt('ff'));
console.log('partly readable  :', parseInt('102', 2), parseInt('12abc', 16), parseInt('  42abc'));
console.log('1e3 two ways     :', parseInt('1e3'), Number('1e3'));
console.log('bad radices      :', radixResult(1), radixResult(37));
console.log('binary byte      :', toBinary8(5), toBinary8(200));
console.log('0.1 in binary    :', (0.1).toString(2), 'length', (0.1).toString(2).length);

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// A shift count is taken modulo 32, which is why 1 << 32 is 1 again.
const isPowerOfTwo = (n) => n > 0 && (n & (n - 1)) === 0;

const popcount = (n) => {
  let bits = n >>> 0;               // unsigned, so negatives count all 32
  let count = 0;
  while (bits) {
    count += bits & 1;
    bits >>>= 1;
  }
  return count;
};

const xorSwap = (a, b) => {
  a ^= b;
  b ^= a;
  a ^= b;
  return [a, b];
};

// Permission flags, the one place bitwise code reads well.
const READ = 1;
const WRITE = 2;
const ADMIN = 4;
const grant = READ | ADMIN;
const can = (flags, flag) => (flags & flag) !== 0;

console.log('5 & 3, 5 | 3, 5 ^ 3 :', 5 & 3, 5 | 3, 5 ^ 3);
console.log('~5, 5 << 1, 5 >> 1  :', ~5, 5 << 1, 5 >> 1);
console.log('-1 >> 0 vs -1 >>> 0 :', -1 >> 0, -1 >>> 0);
console.log('1 << 30, 31, 32     :', 1 << 30, 1 << 31, 1 << 32);
console.log('2 ** 31 | 0         :', 2 ** 31 | 0, 'and 2 ** 32 | 0:', 2 ** 32 | 0);
console.log('~~3.7, ~~-3.7       :', ~~3.7, ~~-3.7);
console.log('powers of two       :', [1, 2, 3, 4, 1024, 0, -4].map(isPowerOfTwo));
console.log('popcount            :', popcount(255), popcount(1024), popcount(-1));
console.log('grant               :', grant, can(grant, READ), can(grant, WRITE), can(grant, ADMIN));

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17const closeAbsolute = (a, b) => Math.abs(a - b) <= Number.EPSILON;
const closeRelative = (a, b) => Math.abs(a - b) <= Number.EPSILON * Math.max(Math.abs(a), Math.abs(b));
const closeFixed = (a, b) => Math.abs(a - b) < 1e-9;

// The same arithmetic, a billion times larger.
const scaledA = (0.1 + 0.2) * 1e9;
const scaledB = 1e9 * 0.3;

console.log('EPSILON            :', Number.EPSILON, 'is 2 ** -52:', Number.EPSILON === 2 ** -52);
console.log('0.1 + 0.2          :', 0.1 + 0.2, 'off by', Math.abs(0.1 + 0.2 - 0.3));
console.log('4500 + 0.1 + 0.2   :', 4500 + 0.1 + 0.2, 'exactly 4500.3?', 4500 + 0.1 + 0.2 === 4500.3);
console.log('1 + EPSILON / 2    :', 1 + Number.EPSILON / 2, 'still 1?', 1 + Number.EPSILON / 2 === 1);
console.log('scaled pair        :', scaledA, scaledB, 'gap', Math.abs(scaledA - scaledB));
console.log('absolute test      :', closeAbsolute(0.1 + 0.2, 0.3), closeAbsolute(scaledA, scaledB));
console.log('relative test      :', closeRelative(0.1 + 0.2, 0.3), closeRelative(scaledA, scaledB));
console.log('fixed 1e-9 test    :', closeFixed(1e-10, 2e-10), '<- calls these equal');
console.log('round then compare :', Number((0.1 + 0.2).toFixed(10)));

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20const add = (a, b) => a + b;
const halves = [0.5, 1.5, 2.5, 3.5, -0.5, -1.5, -2.5];

// Round half to even, the accounting default.
const roundHalfEven = (x) => {
  const down = Math.floor(x);
  const frac = x - down;
  if (frac > 0.5) return down + 1;
  if (frac < 0.5) return down;
  return down % 2 === 0 ? down : down + 1;
};

console.log('typeof toFixed     :', typeof (1.5).toFixed(2), (1.5).toFixed(2));
console.log('1.005 to 2 places  :', (1.005).toFixed(2), 'because it is', (1.005).toPrecision(20));
console.log('2.675 and 8.575    :', (2.675).toFixed(2), (8.575).toFixed(2));
console.log('real halves        :', halves.map((x) => x.toFixed(0)));
console.log('Math.round         :', halves.map(Math.round));
console.log('half to even       :', halves.map(roundHalfEven));
console.log('sums               :', [0.5, 1.5, 2.5, 3.5].map(Math.round).reduce(add, 0), [0.5, 1.5, 2.5, 3.5].map(roundHalfEven).reduce(add, 0), 0.5 + 1.5 + 2.5 + 3.5);
console.log('minus zero, 1e21   :', (-0).toFixed(2), (1e21).toFixed(2));

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
26const mixError = () => {
  try {
    return 1n + 1;
  } catch (error) {
    return error.constructor.name;
  }
};

const stringifyError = () => {
  try {
    return JSON.stringify({ id: 1n });
  } catch (error) {
    return error.constructor.name;
  }
};

console.log('MAX_SAFE_INTEGER   :', Number.MAX_SAFE_INTEGER, '= 2 ** 53 - 1:', Number.MAX_SAFE_INTEGER === 2 ** 53 - 1);
console.log('2 ** 53 + 1        :', 2 ** 53 + 1, 'same as 2 ** 53?', 2 ** 53 === 2 ** 53 + 1);
console.log('the literal        :', 9007199254740993);
console.log('isInteger vs safe  :', Number.isInteger(2 ** 53), Number.isSafeInteger(2 ** 53));
console.log('parseInt of an id  :', parseInt('9007199254740993'));
console.log('through JSON       :', JSON.parse('{"id":9007199254740993}').id);
console.log('as a JSON string   :', JSON.parse('{"id":"9007199254740993"}').id);
console.log('BigInt             :', String(9007199254740993n + 2n), String(2n ** 64n));
console.log('1n + 1             :', mixError(), 'JSON.stringify(1n):', stringifyError());
console.log('past MAX_VALUE     :', Number.MAX_VALUE * 2);

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
28const factorial = (n) => {
  let result = 1;
  for (let i = 2; i <= n; i++) result *= i;
  return result;
};

const factorialRec = (n) => (n <= 1 ? 1 : n * factorialRec(n - 1));

const factorialBig = (n) => {
  let result = 1n;
  for (let i = 2n; i <= BigInt(n); i++) result *= i;
  return result;
};

// The empty loop returns 1 for any n below 2, including nonsense. Say so.
const factorialGuarded = (n) => (Number.isInteger(n) && n >= 0 ? factorial(n) : NaN);

const exactAt = (n) => BigInt(factorial(n)) === factorialBig(n);

console.log('18!, 19!        :', factorial(18), factorial(19));
console.log('safe?           :', Number.isSafeInteger(factorial(18)), Number.isSafeInteger(factorial(19)));
console.log('exact 19 to 23  :', [19, 20, 21, 22, 23].map(exactAt));
console.log('23! as a double :', String(BigInt(factorial(23))));
console.log('23! really      :', String(factorialBig(23)));
console.log('170!, 171!      :', factorial(170), factorial(171));
console.log('25! in BigInt   :', String(factorialBig(25)));
console.log('unguarded edges :', factorial(-3), factorial(2.5));
console.log('20! / 18!       :', factorial(20) / factorial(18));

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
42const add = (a, b) => a + b;

const permutations = (n, r) => {
  let result = 1;
  for (let i = 0; i < r; i++) result *= n - i;
  return result;
};

// Multiply then divide in step: the running value is a binomial coefficient at
// every stage, so it is always an integer and never grows past the answer.
const combinations = (n, k) => {
  const m = Math.min(k, n - k);
  let result = 1;
  for (let i = 1; i <= m; i++) result = (result * (n - m + i)) / i;
  return result;
};

const factorial = (n) => {
  let result = 1;
  for (let i = 2; i <= n; i++) result *= i;
  return result;
};
const combinationsByFactorial = (n, k) => factorial(n) / (factorial(k) * factorial(n - k));

const factorialBig = (n) => {
  let result = 1n;
  for (let i = 2n; i <= n; i++) result *= i;
  return result;
};
const combinationsBig = (n, k) => factorialBig(n) / (factorialBig(k) * factorialBig(n - k));

const pascalRow = (n) => Array.from({ length: n + 1 }, (unused, k) => combinations(n, k));

console.log('P(5, 3), P(30, 5)  :', permutations(5, 3), permutations(30, 5));
console.log('C(5, 3), C(52, 5)  :', combinations(5, 3), combinations(52, 5));
console.log('C(49, 6)           :', combinations(49, 6));
console.log('C(30, 15) two ways :', combinationsByFactorial(30, 15), combinations(30, 15));
console.log('integer?           :', Number.isInteger(combinationsByFactorial(30, 15)), Number.isInteger(combinations(30, 15)));
console.log('C(60, 30) printed  :', combinations(60, 30), 'exact value', String(BigInt(combinations(60, 30))));
console.log('C(100, 50) double  :', String(BigInt(combinations(100, 50))));
console.log('C(100, 50) BigInt  :', String(combinationsBig(100n, 50n)));
console.log('Pascal row 6       :', pascalRow(6), 'sums to', pascalRow(6).reduce(add, 0));

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
56const digitSum = (n) => {
  let x = Math.abs(n);
  let total = 0;
  while (x > 0) {
    total += x % 10;
    x = Math.floor(x / 10);
  }
  return total;
};

const digitSumByString = (n) => String(Math.abs(n)).split('').reduce((s, d) => s + Number(d), 0);

const digitalRootByLoop = (n) => {
  let x = Math.abs(n);
  while (x > 9) x = digitSum(x);
  return x;
};

// Casting out nines: the root is n mod 9, with 9 in place of 0.
const digitalRoot = (n) => (n === 0 ? 0 : 1 + ((Math.abs(n) - 1) % 9));

const agreeUpTo = (limit) => {
  for (let n = 0; n <= limit; n++) if (digitalRootByLoop(n) !== digitalRoot(n)) return false;
  return true;
};

const persistence = (n) => {
  let x = Math.abs(n);
  let steps = 0;
  while (x > 9) {
    x = digitSum(x);
    steps++;
  }
  return steps;
};

const reverseDigits = (n) => {
  let x = Math.abs(n);
  let out = 0;
  while (x > 0) {
    out = out * 10 + (x % 10);
    x = Math.floor(x / 10);
  }
  return n < 0 ? -out : out;
};
const isPalindrome = (n) => n >= 0 && reverseDigits(n) === n;

console.log('digit sums       :', [0, 9, 18, 493193, 12345].map(digitSum));
console.log('by string        :', [0, 9, 18, 493193, 12345].map(digitSumByString));
console.log('1e21 prints as   :', String(1e21), '-> by string:', digitSumByString(1e21), 'by arithmetic:', digitSum(1e21));
console.log('roots by loop    :', [0, 9, 18, 493193, 12345, 999999999].map(digitalRootByLoop));
console.log('roots by formula :', [0, 9, 18, 493193, 12345, 999999999].map(digitalRoot));
console.log('27 and 28 mod 9  :', digitalRoot(27), 27 % 9, digitalRoot(28), 28 % 9);
console.log('persistence      :', persistence(9), persistence(199), persistence(999999999999));
console.log('reversed         :', reverseDigits(12345), reverseDigits(-120));
console.log('palindromes      :', [121, 12321, 123, 10].map(isPalindrome));

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
58const NUMERALS = [
  [1000, 'M'], [900, 'CM'], [500, 'D'], [400, 'CD'],
  [100, 'C'], [90, 'XC'], [50, 'L'], [40, 'XL'],
  [10, 'X'], [9, 'IX'], [5, 'V'], [4, 'IV'], [1, 'I'],
];
const PLAIN = [[1000, 'M'], [500, 'D'], [100, 'C'], [50, 'L'], [10, 'X'], [5, 'V'], [1, 'I']];

const greedy = (table) => (n) => {
  let out = '';
  let left = n;
  for (const [value, letters] of table) {
    while (left >= value) {
      out += letters;
      left -= value;
    }
  }
  return out;
};

const toRoman = greedy(NUMERALS);
const toRomanPlain = greedy(PLAIN);

const LETTER = { I: 1, V: 5, X: 10, L: 50, C: 100, D: 500, M: 1000 };

const fromRoman = (s) => {
  let total = 0;
  for (let i = 0; i < s.length; i++) {
    const here = LETTER[s[i]];
    const next = i + 1 < s.length ? LETTER[s[i + 1]] : 0;
    total += here < next ? -here : here;
  }
  return total;
};

const isCanonical = (s) => s.length > 0 && toRoman(fromRoman(s)) === s;

const roundTripHolds = (limit) => {
  for (let n = 1; n <= limit; n++) if (fromRoman(toRoman(n)) !== n) return false;
  return true;
};

const longestUnder4000 = (() => {
  let best = '';
  for (let n = 1; n < 4000; n++) {
    const r = toRoman(n);
    if (r.length > best.length) best = r;
  }
  return best;
})();

console.log('small numbers :', [1, 4, 9, 14, 40, 90, 400].map(toRoman));
console.log('1994, 3999    :', toRoman(1994), toRoman(3999));
console.log('without pairs :', toRomanPlain(9), toRomanPlain(1994));
console.log('read back     :', fromRoman('MCMXCIV'), fromRoman('MMMCMXCIX'));
console.log('nonsense      :', fromRoman('IIII'), fromRoman('IC'));
console.log('canonical?    :', isCanonical('MCMXCIV'), isCanonical('IIII'), isCanonical('IC'));
console.log('zero, negative:', JSON.stringify(toRoman(0)), JSON.stringify(toRoman(-5)));
console.log('longest       :', longestUnder4000, longestUnder4000.length, 'which is', fromRoman(longestUnder4000));

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
35const add = (a, b) => a + b;

const mean = (xs) => xs.reduce(add, 0) / xs.length;

const median = (xs) => {
  const sorted = [...xs].sort((a, b) => a - b);
  const mid = sorted.length >> 1;
  return sorted.length % 2 === 1 ? sorted[mid] : (sorted[mid - 1] + sorted[mid]) / 2;
};

const modes = (xs) => {
  const counts = new Map();
  for (const x of xs) counts.set(x, (counts.get(x) ?? 0) + 1);
  const top = Math.max(...counts.values());
  return [...counts.keys()].filter((k) => counts.get(k) === top).sort((a, b) => a - b);
};

const scores = [12, 7, 3, 7, 9, 12, 15, 7];
const withOutlier = [12, 7, 3, 7, 9, 12, 1500, 7];

// Proof that sort mutates.
const untouched = [10, 9, 100, 1];
const copySorted = [...untouched].sort((a, b) => a - b);
const inPlace = [10, 9, 100, 1];
inPlace.sort((a, b) => a - b);

console.log('mean, median, modes :', mean(scores), median(scores), modes(scores));
console.log('odd length median   :', median([5, 1, 3]));
console.log('two modes           :', modes([1, 2, 2, 3, 3, 4]), 'all distinct:', modes([1, 2, 3]));
console.log('default sort        :', [10, 9, 100, 1].sort());
console.log('after copying       :', untouched, copySorted);
console.log('after sorting in place:', inPlace);
console.log('empty               :', mean([]), median([]));
console.log('one outlier         :', mean(withOutlier), median(withOutlier));
console.log('is the median in the data?', scores.includes(median(scores)));

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
48const add = (a, b) => a + b;
const mean = (xs) => xs.reduce(add, 0) / xs.length;

// Two passes: find the mean, then sum the squared deviations from it.
const populationVariance = (xs) => {
  const m = mean(xs);
  return xs.reduce((s, x) => s + (x - m) ** 2, 0) / xs.length;
};
const sampleVariance = (xs) => {
  const m = mean(xs);
  return xs.reduce((s, x) => s + (x - m) ** 2, 0) / (xs.length - 1);
};

// One pass, textbook, and dangerous.
const onePassVariance = (xs) => {
  const n = xs.length;
  const sum = xs.reduce(add, 0);
  const sumSquares = xs.reduce((s, x) => s + x * x, 0);
  return sumSquares / n - (sum / n) ** 2;
};

// One pass, done properly: update the mean as you go.
const welfordVariance = (xs) => {
  let n = 0;
  let m = 0;
  let m2 = 0;
  for (const x of xs) {
    n++;
    const delta = x - m;
    m += delta / n;
    m2 += delta * (x - m);
  }
  return m2 / n;
};

const sample = [2, 4, 4, 4, 5, 5, 7, 9];
const shifted = sample.map((x) => x + 1e9);
const scaled = sample.map((x) => x * 10);

console.log('mean                  :', mean(sample));
console.log('population var, sd    :', populationVariance(sample), Math.sqrt(populationVariance(sample)));
console.log('sample var, sd        :', sampleVariance(sample), Math.sqrt(sampleVariance(sample)));
console.log('one pass, small data  :', onePassVariance(sample), 'Welford:', welfordVariance(sample));
console.log('shifted by 1e9, two-pass:', populationVariance(shifted));
console.log('shifted, one pass     :', onePassVariance(shifted), '<- the spread did not change');
console.log('shifted, Welford      :', welfordVariance(shifted));
console.log('scaled by 10          :', populationVariance(scaled));
console.log('a single value        :', populationVariance([5]), sampleVariance([5]));

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
32const compound = (principal, rate, perYear, years) =>
  principal * (1 + rate / perYear) ** (perYear * years);

const simple = (principal, rate, years) => principal * (1 + rate * years);
const continuous = (principal, rate, years) => principal * Math.exp(rate * years);

// Effective annual rate: what 12% nominal, compounded monthly, really pays.
const effectiveAnnual = (rate, perYear) => (1 + rate / perYear) ** perYear - 1;

// The same five years in integer kobo, rounded to a whole kobo every month.
const monthByMonthKobo = (startKobo, rate, months) => {
  let balance = startKobo;
  for (let m = 0; m < months; m++) balance = Math.round(balance * (1 + rate / 12));
  return balance;
};
const closedFormKobo = (startKobo, rate, months) =>
  Math.round(startKobo * (1 + rate / 12) ** months);

const naira = 250000;
const kobo = naira * 100;
const rate = 0.12;

console.log('simple, 5 years   :', simple(naira, rate, 5));
console.log('annual            :', compound(naira, rate, 1, 5), '->', compound(naira, rate, 1, 5).toFixed(2));
console.log('monthly           :', compound(naira, rate, 12, 5), '->', compound(naira, rate, 12, 5).toFixed(2));
console.log('daily             :', compound(naira, rate, 365, 5));
console.log('continuous        :', continuous(naira, rate, 5));
console.log('effective annual  :', effectiveAnnual(rate, 12));
console.log('kobo, closed form :', closedFormKobo(kobo, rate, 60));
console.log('kobo, month by month:', monthByMonthKobo(kobo, rate, 60));
console.log('difference        :', monthByMonthKobo(kobo, rate, 60) - closedFormKobo(kobo, rate, 60), 'kobo');
console.log('doubling time     :', Math.log(2) / Math.log(1 + rate), 'rule of 72 says', 72 / 12);

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18const percentChange = (from, to) => ((to - from) / from) * 100;

// The fall that exactly undoes a given rise, as a percentage.
const reverseChange = (p) => 100 / (1 + p / 100) - 100;

// VAT in Nigeria is 7.5%. Held in kobo, it stays an integer.
const vatKobo = (k) => Math.round(k * 0.075);

console.log('4000 -> 5000, and back:', percentChange(4000, 5000), percentChange(5000, 4000));
console.log('up 20 then down 20    :', 4000 * 1.2 * 0.8, 'started at 4000');
console.log('to undo +25 you need  :', reverseChange(25), 'and to undo -20:', reverseChange(-20));
console.log('two 10% rises         :', 1.1 * 1.1, 'against 1.21');
console.log('20% off then 10% off  :', 1 - 0.8 * 0.9, 'against 0.3');
console.log('VAT on 12500 naira    :', 12500 * 0.075, 'total', 12500 * 1.075);
console.log('backing VAT out again :', 13437.5 / 1.075);
console.log('from zero             :', percentChange(0, 5), 'and zero to zero:', percentChange(0, 0));
console.log('from -100 to -50      :', percentChange(-100, -50), '<- that is a rise');
console.log('VAT in kobo           :', vatKobo(1250000), 'total', 1250000 + vatKobo(1250000));

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
30const heron = (a, b, c) => {
  const s = (a + b + c) / 2;
  return Math.sqrt(s * (s - a) * (s - b) * (s - c));
};

// Kahan's rearrangement. Same algebra, different grouping, and the sort matters.
const kahan = (a, b, c) => {
  const [x, y, z] = [a, b, c].slice().sort((p, q) => q - p);
  return Math.sqrt((x + (y + z)) * (z - (x - y)) * (z + (x - y)) * (x + (y - z))) / 4;
};

const isTriangle = (a, b, c) =>
  a > 0 && b > 0 && c > 0 && a + b > c && a + c > b && b + c > a;

const needle = [100000, 99999.99979, 0.00029];

// Computed offline in exact rational arithmetic, for the triangle whose sides
// are the doubles nearest those three decimals: 10.0000000770210383232996...
const exactNeedleArea = 10.000000077021038;

console.log('3-4-5      :', heron(3, 4, 5), kahan(3, 4, 5));
console.log('13-14-15   :', heron(13, 14, 15), kahan(13, 14, 15));
console.log('equilateral:', heron(2, 2, 2), 'against sqrt(3)', Math.sqrt(3));
console.log('needle, Heron:', heron(needle[0], needle[1], needle[2]));
console.log('needle, Kahan:', kahan(needle[0], needle[1], needle[2]));
console.log('the exact area:', exactNeedleArea);
console.log('Heron is out by:', Math.abs(heron(needle[0], needle[1], needle[2]) - exactNeedleArea));
console.log('degenerate 1-2-3:', heron(1, 2, 3), kahan(1, 2, 3));
console.log('impossible 1-2-10:', heron(1, 2, 10), kahan(1, 2, 10));
console.log('valid?     :', isTriangle(3, 4, 5), isTriangle(1, 2, 3), isTriangle(1, 2, 10));

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
37const shoelace = (points) => {
  let sum = 0;
  for (let i = 0; i < points.length; i++) {
    const [x1, y1] = points[i];
    const [x2, y2] = points[(i + 1) % points.length];
    sum += x1 * y2 - x2 * y1;
  }
  return sum / 2;
};

const area = (points) => Math.abs(shoelace(points));
const isCounterClockwise = (points) => shoelace(points) > 0;
const translate = (points, dx, dy) => points.map(([x, y]) => [x + dx, y + dy]);
const perimeter = (points) =>
  points.reduce((s, p, i) => {
    const q = points[(i + 1) % points.length];
    return s + Math.hypot(q[0] - p[0], q[1] - p[1]);
  }, 0);

const square = [[0, 0], [4, 0], [4, 4], [0, 4]];
const triangle = [[0, 0], [4, 0], [0, 3]];
const lShape = [[0, 0], [6, 0], [6, 2], [2, 2], [2, 5], [0, 5]];
const pentagon = [[3, 4], [5, 11], [12, 8], [9, 5], [5, 6]];
const hexagon = Array.from({ length: 6 }, (unused, i) => [
  Math.cos((i * Math.PI) / 3),
  Math.sin((i * Math.PI) / 3),
]);

console.log('square, both windings:', shoelace(square), shoelace([...square].reverse()));
console.log('triangle             :', shoelace(triangle), 'base times height over 2 =', (4 * 3) / 2);
console.log('L-shape              :', shoelace(lShape), 'perimeter', perimeter(lShape));
console.log('pentagon             :', shoelace(pentagon), '<- negative, so clockwise');
console.log('unit hexagon         :', shoelace(hexagon), 'against 3 root 3 over 2 =', (3 * Math.sqrt(3)) / 2);
console.log('repeating the first point:', shoelace([...square, [0, 0]]));
console.log('three collinear points   :', shoelace([[0, 0], [1, 1], [2, 2]]));
console.log('pentagon moved by 1000   :', shoelace(translate(pentagon, 1000, 1000)));
console.log('pentagon moved by 1e9    :', shoelace(translate(pentagon, 1e9, 1e9)), '<- was -30');

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
27const MILE_IN_KM = 1.609344;         // exact, by international agreement
const POUND_IN_KG = 0.45359237;      // also exact

const kmToMiles = (km) => km / MILE_IN_KM;
const milesToKm = (mi) => mi * MILE_IN_KM;
const kgToPounds = (kg) => kg / POUND_IN_KG;
const cToF = (c) => (c * 9) / 5 + 32;
const fToC = (f) => ((f - 32) * 5) / 9;
const cToK = (c) => c + 273.15;
const kToC = (k) => k - 273.15;

const failingCelsius = (lo, hi) => {
  const out = [];
  for (let c = lo; c <= hi; c++) if (fToC(cToF(c)) !== c) out.push(c);
  return out;
};

console.log('one mile          :', milesToKm(1), 'km, exactly');
console.log('100 km            :', kmToMiles(100), 'miles, back to', milesToKm(kmToMiles(100)));
console.log('7 km              :', kmToMiles(7), 'miles, back to', milesToKm(kmToMiles(7)));
console.log('a marathon        :', kmToMiles(42.195), 'miles');
console.log('body temperature  :', cToF(37), 'F, and back', fToC(98.6));
console.log('the crossing point:', cToF(-40));
console.log('0.1 C out and back:', fToC(cToF(0.1)));
console.log('integer C failures:', failingCelsius(-50, 50));
console.log('kelvin round trip :', kToC(cToK(25)), kToC(cToK(25)) === 25);
console.log('70 kg             :', kgToPounds(70), 'pounds');

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
52const add = (a, b) => a + b;

const factorise = (n) => {
  const factors = [];
  let x = n;
  for (let d = 2; d * d <= x; d++) {
    while (x % d === 0) {
      factors.push(d);
      x /= d;
    }
  }
  if (x > 1) factors.push(x);
  return factors;
};

const grouped = (n) => {
  const pairs = [];
  for (const p of factorise(n)) {
    const last = pairs[pairs.length - 1];
    if (last && last[0] === p) last[1]++;
    else pairs.push([p, 1]);
  }
  return pairs;
};

// Both of these read the answer off the exponents, listing nothing.
const divisorCount = (n) => grouped(n).reduce((c, pair) => c * (pair[1] + 1), 1);
const divisorSum = (n) =>
  grouped(n).reduce((s, pair) => s * ((pair[0] ** (pair[1] + 1) - 1) / (pair[0] - 1)), 1);

const divisors = (n) => {
  const out = [];
  for (let d = 1; d * d <= n; d++) {
    if (n % d === 0) {
      out.push(d);
      if (d !== n / d) out.push(n / d);
    }
  }
  return out.sort((a, b) => a - b);
};

const rebuilds = (n) => factorise(n).reduce((p, q) => p * q, 1) === n;

console.log('360        :', factorise(360));
console.log('97, 1, 2   :', factorise(97), factorise(1), factorise(2));
console.log('1024, 13195:', factorise(1024).length, 'twos and', factorise(13195));
console.log('grouped    :', grouped(360), grouped(5040));
console.log('9973 * 9973:', factorise(9973 * 9973), 'which is', 9973 * 9973);
console.log('divisors of 360 :', divisors(360));
console.log('count, from pairs:', divisorCount(360), 'from the list:', divisors(360).length);
console.log('sum, from pairs  :', divisorSum(360), 'from the list:', divisors(360).reduce(add, 0));
console.log('zero and negative:', factorise(0), factorise(-12), '<- the loop never starts');

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// Sum of proper divisors. Starts at 1 and never adds n itself.
const aliquotSum = (n) => {
  if (n < 2) return 0;
  let total = 1;
  for (let d = 2; d * d <= n; d++) {
    if (n % d === 0) {
      total += d;
      const pair = n / d;
      if (pair !== d) total += pair;
    }
  }
  return total;
};

const classify = (n) => {
  const s = aliquotSum(n);
  return s === n ? 'perfect' : s > n ? 'abundant' : 'deficient';
};

const isAmicablePair = (a, b) => a !== b && aliquotSum(a) === b && aliquotSum(b) === a;

const perfectUpTo = (limit) => {
  const out = [];
  for (let n = 2; n <= limit; n++) if (aliquotSum(n) === n) out.push(n);
  return out;
};

const amicableUnder = (limit) => {
  const out = [];
  for (let a = 2; a < limit; a++) {
    const b = aliquotSum(a);
    if (b > a && b < limit && aliquotSum(b) === a) out.push([a, b]);
  }
  return out;
};

const euclidPerfect = (p) => 2 ** (p - 1) * (2 ** p - 1);

// Computed once: each of these walks a lot of divisors.
const perfects = perfectUpTo(10000);
const amicable = amicableUnder(1300);

console.log('aliquot sums   :', [1, 6, 8, 12, 28, 97].map(aliquotSum));
console.log('classified     :', [1, 6, 8, 12, 28, 97].map(classify));
console.log('perfect <= 1e4 :', perfects);
console.log('220 and 284    :', aliquotSum(220), aliquotSum(284));
console.log('amicable < 1300:', amicable);
console.log('Euclid-Euler   :', [2, 3, 5, 7, 13].map(euclidPerfect));
console.log('the fifth one  :', euclidPerfect(13), 'is perfect:', aliquotSum(33550336) === 33550336);

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
52const collatz = (n) => {
  const seq = [n];
  let x = n;
  while (x !== 1) {
    x = x % 2 === 0 ? x / 2 : 3 * x + 1;
    seq.push(x);
  }
  return seq;
};

const steps = (n) => collatz(n).length - 1;
const peak = (n) => collatz(n).reduce((m, x) => (x > m ? x : m), 0);

// The honest version: bounded, and it reports whether it finished.
const reachesOneWithin = (n, cap) => {
  let x = n;
  let used = 0;
  while (x !== 1 && used < cap) {
    x = x % 2 === 0 ? x / 2 : 3 * x + 1;
    used++;
  }
  return x === 1;
};

const allTerminate = (limit, cap) => {
  for (let n = 1; n <= limit; n++) if (!reachesOneWithin(n, cap)) return false;
  return true;
};

const longestUnder = (limit) => {
  let bestN = 1;
  let bestSteps = -1;
  for (let n = 1; n < limit; n++) {
    const s = steps(n);
    if (s > bestSteps) {
      bestSteps = s;
      bestN = n;
    }
  }
  return [bestN, bestSteps];
};

const longest = longestUnder(1000);

console.log('from 6        :', collatz(6));
console.log('from 1        :', collatz(1), '<- the loop never runs');
console.log('27            :', steps(27), 'steps, peaking at', peak(27));
console.log('97            :', steps(97), 'steps, peaking at', peak(97));
console.log('7             :', steps(7), 'steps, peaking at', peak(7));
console.log('powers of two :', steps(1), steps(2), steps(1024));
console.log('longest < 1000:', longest);
console.log('capped at 500 :', allTerminate(200, 500), 'capped at 50:', allTerminate(200, 50));

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
38const add = (a, b) => a + b;

const fib = (n) => {
  let a = 0;
  let b = 1;
  for (let i = 0; i < n; i++) {
    const next = a + b;
    a = b;
    b = next;
  }
  return a;
};

const fibBig = (n) => {
  let a = 0n;
  let b = 1n;
  for (let i = 0; i < n; i++) {
    const next = a + b;
    a = b;
    b = next;
  }
  return a;
};

const PHI = (1 + Math.sqrt(5)) / 2;
const binet = (n) => Math.round(PHI ** n / Math.sqrt(5));

const firstTen = Array.from({ length: 10 }, (unused, i) => fib(i));

console.log('the first ten  :', firstTen);
console.log('F(10), F(30)   :', fib(10), fib(30), 'F(50):', fib(50));
console.log('F(78)          :', fib(78), 'safe:', Number.isSafeInteger(fib(78)));
console.log('F(79) by loop  :', fib(79), 'safe:', Number.isSafeInteger(fib(79)));
console.log('F(79) truly    :', String(fibBig(79)));
console.log('Binet at 70, 71:', binet(70), binet(71), 'against', fib(70), fib(71));
console.log('F(300)         :', String(fibBig(300)));
console.log('ratio to PHI   :', fib(30) / fib(29), fib(50) / fib(49), PHI);
console.log('the identity   :', fib(10) ** 2 + fib(11) ** 2, 'should be F(21) =', fib(21));

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
69const add = (a, b) => a + b;

const sieve = (limit) => {
  const composite = new Uint8Array(limit + 1);
  const primes = [];
  for (let i = 2; i <= limit; i++) {
    if (!composite[i]) {
      primes.push(i);
      for (let j = i * i; j <= limit; j += i) composite[j] = 1;
    }
  }
  return primes;
};

const isPrime = (n) => {
  if (!Number.isInteger(n) || n < 2) return false;
  if (n % 2 === 0) return n === 2;
  for (let d = 3; d * d <= n; d += 2) if (n % d === 0) return false;
  return true;
};

const primes = sieve(1000);
const primeSet = new Set(primes);

const trialAgrees = (limit) => {
  for (let n = 0; n <= limit; n++) if (isPrime(n) !== primeSet.has(n)) return false;
  return true;
};

const twins = primes.filter((p, i) => primes[i + 1] === p + 2).map((p) => [p, p + 2]);

const largestGap = () => {
  let gap = 0;
  let after = 0;
  for (let i = 1; i < primes.length; i++) {
    const g = primes[i] - primes[i - 1];
    if (g > gap) {
      gap = g;
      after = primes[i - 1];
    }
  }
  return [gap, after];
};

// Returns the first even number that is NOT a sum of two primes, or 0 if none.
const goldbachFailsAt = (limit) => {
  for (let n = 4; n <= limit; n += 2) {
    let found = false;
    for (const p of primes) {
      if (p > n / 2) break;
      if (primeSet.has(n - p)) {
        found = true;
        break;
      }
    }
    if (!found) return n;
  }
  return 0;
};

console.log('how many      :', primes.length, 'the last one:', primes[primes.length - 1]);
console.log('the first ten :', primes.slice(0, 10));
console.log('their sum     :', primes.reduce(add, 0));
console.log('trial agrees  :', trialAgrees(1000));
console.log('twin pairs    :', twins.length, twins[0], 'up to', twins[twins.length - 1]);
console.log('largest gap   :', largestGap());
console.log('Goldbach to 200:', goldbachFailsAt(200), '<- zero means it held');
console.log('n / ln n says :', 1000 / Math.log(1000), 'against the real', primes.length);
console.log('tiny limits   :', sieve(10), sieve(2), sieve(1));

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
80const mod = (a, n) => ((a % n) + n) % n;

// Square and multiply. Nothing ever grows past m squared.
const powMod = (base, exponent, m) => {
  let result = 1;
  let b = mod(base, m);
  let e = exponent;
  while (e > 0) {
    if (e % 2 === 1) result = (result * b) % m;
    b = (b * b) % m;
    e = Math.floor(e / 2);
  }
  return result;
};

// Extended Euclid: returns [g, x, y] with a * x + b * y === g.
const egcd = (a, b) => {
  let [oldR, r] = [a, b];
  let [oldS, s] = [1, 0];
  let [oldT, t] = [0, 1];
  while (r !== 0) {
    const q = Math.floor(oldR / r);
    [oldR, r] = [r, oldR - q * r];
    [oldS, s] = [s, oldS - q * s];
    [oldT, t] = [t, oldT - q * t];
  }
  return [oldR, oldS, oldT];
};

const modInverse = (a, m) => {
  const [g, x] = egcd(mod(a, m), m);
  return g === 1 ? mod(x, m) : NaN;
};

// ISBN-10: weights 10 down to 1, and the whole thing must be 0 mod 11.
const isbn10Weighted = (chars) => {
  let total = 0;
  for (let i = 0; i < chars.length; i++) {
    const ch = chars[i];
    total += (10 - i) * (ch === 'X' ? 10 : Number(ch));
  }
  return total;
};
const digitsOf = (s) => s.split('').filter((ch) => ch !== '-');
const isbn10CheckDigit = (first9) => {
  const c = mod(-isbn10Weighted(digitsOf(first9)), 11);
  return c === 10 ? 'X' : String(c);
};
const isbn10Valid = (s) => {
  const d = digitsOf(s);
  return d.length === 10 && isbn10Weighted(d) % 11 === 0;
};

// Luhn: mod 10, doubling every second digit from the right.
const luhn = (s) => {
  let sum = 0;
  let double = false;
  for (let i = s.length - 1; i >= 0; i--) {
    let d = Number(s[i]);
    if (double) {
      d *= 2;
      if (d > 9) d -= 9;
    }
    sum += d;
    double = !double;
  }
  return sum % 10 === 0;
};

console.log('2 ** 10 mod 1000  :', powMod(2, 10, 1000));
console.log('7 ** 100 mod 13   :', powMod(7, 100, 13), 'naively:', 7 ** 100 % 13, 'BigInt:', String(7n ** 100n % 13n));
console.log('7 ** 1000 naively :', 7 ** 1000, '->', 7 ** 1000 % 13);
console.log('Fermat mod 13     :', [2, 3, 5, 6, 7].map((a) => powMod(a, 12, 13)));
console.log('the same on 15    :', powMod(2, 14, 15), '<- 15 is composite');
console.log('egcd(240, 46)     :', egcd(240, 46), 'check:', 240 * -9 + 46 * 47);
console.log('inverses mod 11   :', Array.from({ length: 10 }, (unused, i) => modInverse(i + 1, 11)));
console.log('6 mod 9           :', modInverse(6, 9), '<- they share a factor of 3');
console.log('check digits      :', isbn10CheckDigit('030640615'), isbn10CheckDigit('043942089'));
console.log('ISBN transposed   :', isbn10Valid('0306406152'), isbn10Valid('0306406512'));
console.log('Luhn transposed   :', luhn('4091234567890121'), luhn('4901234567890121'));

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
60const toKobo = (amount) => Math.round(amount * 100);
const naira = (kobo) => '₦' + (kobo / 100).toFixed(2);

// Hand out `total` kobo in the given ratio. Floor every share, then give the
// unallocated kobo to the shares with the largest discarded fraction, largest
// first. This is the largest-remainder rule, and it never invents or loses a kobo.
function allocate(total, weights) {
  const weight = weights.reduce((a, b) => a + b, 0);
  const exact = weights.map((w) => (total * w) / weight);
  const shares = exact.map(Math.floor);
  const left = total - shares.reduce((a, b) => a + b, 0);
  const byFraction = exact
    .map((v, i) => ({ i, frac: v - Math.floor(v) }))
    .sort((a, b) => b.frac - a.frac || a.i - b.i);
  for (let k = 0; k < left; k++) shares[byFraction[k].i]++;
  return shares;
}

// The wrong way: round every share independently and hope.
const roundEach = (total, parts) => new Array(parts).fill(Math.round(total / parts));

// Two baskets of naira prices, priced the same way. Only one of them adds up.
const basketA = [1250.75, 899.99, 45.5, 12345.67];
const basketB = [12000.05, 3450.15, 89.9];
const floatA = basketA.reduce((a, b) => a + b, 0);
const floatB = basketB.reduce((a, b) => a + b, 0);
const koboA = basketA.map(toKobo).reduce((a, b) => a + b, 0);
const koboB = basketB.map(toKobo).reduce((a, b) => a + b, 0);

const naiveThirds = roundEach(10000, 3);              // ₦100.00 three ways
const thirds = allocate(10000, [1, 1, 1]);
const centres = allocate(999999, [3, 5, 7]);          // ₦9,999.99 as 3 : 5 : 7
const exactCentres = [3, 5, 7].map((w) => (999999 * w) / 15);
const drift = centres.map((share, i) => Math.abs(share - exactCentres[i]));

// 7.5% VAT taken back out of a VAT-inclusive price. Round the net, and let the
// VAT be whatever is left, so the two halves always reconstruct the price.
const RATE = 0.075;
const netOf = (gross) => Math.round(gross / (1 + RATE));
const vatOf = (gross) => gross - netOf(gross);

// What happens if you round both halves instead?
let bothRounded = 0;
let firstLoss = 0;
for (let gross = 1; gross <= 10000; gross++) {
  const net = netOf(gross);
  if (net + Math.round(net * RATE) !== gross) {
    bothRounded++;
    if (firstLoss === 0) firstLoss = gross;
  }
}

console.log('basket A     :', floatA, 'kobo total', koboA, 'agree:', floatA === koboA / 100);
console.log('basket B     :', floatB, 'kobo total', koboB, 'agree:', floatB === koboB / 100);
console.log('naive thirds :', naiveThirds, 'sum', naiveThirds.reduce((a, b) => a + b, 0));
console.log('allocated    :', thirds, 'sum', thirds.reduce((a, b) => a + b, 0));
console.log('3 : 5 : 7    :', centres, 'sum', centres.reduce((a, b) => a + b, 0));
console.log('exact shares :', exactCentres, 'drift', drift);
console.log('VAT out of   :', naira(5375), '-> net', naira(netOf(5375)), 'VAT', naira(vatOf(5375)));
console.log('both rounded :', bothRounded, 'of the first 10000 prices fail, from', firstLoss, 'kobo');

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
66const naive = (xs) => xs.reduce((a, b) => a + b, 0);

// Kahan: `c` carries what the last addition lost, and is removed from the next term.
function kahan(xs) {
  let sum = 0;
  let c = 0;
  for (const x of xs) {
    const y = x - c;
    const t = sum + y;
    c = t - sum - y;
    sum = t;
  }
  return sum;
}

// Neumaier: works out which operand was larger, accumulates the loss instead of
// folding it into the next term, and adds it back at the very end.
function neumaier(xs) {
  let sum = 0;
  let c = 0;
  for (const x of xs) {
    const t = sum + x;
    c += Math.abs(sum) >= Math.abs(x) ? sum - t + x : x - t + sum;
    sum = t;
  }
  return sum + c;
}

// Pairwise: halve, sum each half, add. No compensation, no extra state.
function pairwise(xs, lo, hi) {
  const from = lo === undefined ? 0 : lo;
  const to = hi === undefined ? xs.length : hi;
  if (to - from <= 8) {
    let s = 0;
    for (let i = from; i < to; i++) s += xs[i];
    return s;
  }
  const mid = (from + to) >> 1;
  return pairwise(xs, from, mid) + pairwise(xs, mid, to);
}

const N = 100000;
const tenths = new Array(N).fill(0.1);
const harmonic = Array.from({ length: N }, (_, i) => 1 / (i + 1));

// A reference for the harmonic sum, in exact integers: 30 decimal places of
// each term, added as BigInt, so the only error is the final conversion.
const SCALE = 10n ** 30n;
let scaled = 0n;
for (let k = 1n; k <= BigInt(N); k++) scaled += SCALE / k;
const reference = Number(scaled) / 1e30;

const ascending = harmonic.slice().reverse();          // smallest term first
const hostile = [1, 1e100, 1, -1e100];                 // true sum: 2
const cancelling = [1e16, 3, -1e16];                   // true sum: 3

console.log('100000 tenths, naive :', naive(tenths));
console.log('kahan / neumaier     :', kahan(tenths), neumaier(tenths));
console.log('pairwise             :', pairwise(tenths));
console.log('harmonic reference   :', reference);
console.log('naive, descending    :', naive(harmonic), naive(harmonic) - reference);
console.log('naive, ascending     :', naive(ascending), naive(ascending) - reference);
console.log('kahan / pairwise     :', kahan(harmonic), pairwise(harmonic));
console.log('hostile list         :', naive(hostile), kahan(hostile), neumaier(hostile), '<- want 2');
console.log('cancelling list      :', naive(cancelling), kahan(cancelling), neumaier(cancelling), '<- want 3');
console.log('a single addition    :', kahan([0.1, 0.2]), neumaier([0.1, 0.2]));

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
71function factorial(n) {
  let product = 1n;
  for (let k = 2n; k <= BigInt(n); k++) product *= k;
  return product;
}

function factorialDouble(n) {
  let product = 1;
  for (let k = 2; k <= n; k++) product *= k;
  return product;
}

function fib(n) {
  let a = 0n;
  let b = 1n;
  for (let i = 0; i < n; i++) {
    const next = a + b;
    a = b;
    b = next;
  }
  return a;
}

function fibDouble(n) {
  let a = 0;
  let b = 1;
  for (let i = 0; i < n; i++) {
    const next = a + b;
    a = b;
    b = next;
  }
  return a;
}

// Multiply then divide, alternately: the running value is a binomial coefficient
// at every step, so it is always whole and never larger than the answer.
function choose(n, k) {
  let result = 1n;
  for (let i = 1n; i <= BigInt(k); i++) result = (result * (BigInt(n) - BigInt(k) + i)) / i;
  return result;
}

const f100 = factorial(100);
const text = String(f100);
let trailingZeros = 0;
while (text[text.length - 1 - trailingZeros] === '0') trailingZeros++;
const legendre = Math.floor(100 / 5) + Math.floor(100 / 25);
const digitSum = [...text].reduce((a, d) => a + Number(d), 0);

// The first index where the double version stops telling the truth.
let firstWrong = 0;
for (let i = 1; i <= 100 && firstWrong === 0; i++) {
  if (BigInt(fibDouble(i)) !== fib(i)) firstWrong = i;
}

// The four refusals, caught so the panel survives them.
const refusals = [];
try { refusals.push(1n + 1); } catch (e) { refusals.push(e.constructor.name); }
try { refusals.push(JSON.stringify({ total: 1n })); } catch (e) { refusals.push(e.constructor.name); }
try { refusals.push(Math.sqrt(4n)); } catch (e) { refusals.push(e.constructor.name); }
try { refusals.push(BigInt(1.5)); } catch (e) { refusals.push(e.constructor.name); }

console.log('100!            :', text);
console.log('digits          :', text.length, 'trailing zeros', trailingZeros, 'digit sum', digitSum);
console.log('as a double     :', Number(f100), 'against the double loop', factorialDouble(100));
console.log('first wrong fib :', firstWrong, String(fib(firstWrong)), fibDouble(firstWrong));
console.log('fib 200         :', String(fib(200)));
console.log('C(100, 50)      :', String(choose(100, 50)), 'as a double', Number(choose(100, 50)));
console.log('refusals        :', refusals);
console.log('7n / 2n         :', String(7n / 2n), 'and -7n / 2n is', String(-7n / 2n));
console.log('lost before conversion:', String(BigInt(Number.MAX_SAFE_INTEGER + 2)));

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
82const step = (x, n) => (x + n / x) / 2;

// Iterate until the value stops changing, or repeats the one before it.
function heron(n) {
  if (n < 0) return { root: NaN, steps: 0 };
  if (n === 0) return { root: 0, steps: 0 };
  let x = n;
  let prev = -1;
  let steps = 0;
  while (steps < 100) {
    const next = step(x, n);
    steps++;
    if (next === x || next === prev) return { root: Math.min(x, next), steps };
    prev = x;
    x = next;
  }
  return { root: x, steps };
}

// Same iteration, started from a power of two near the root.
function heronGuessed(n) {
  if (n <= 0) return heron(n);
  let x = 2 ** Math.round(Math.log2(n) / 2);
  let prev = -1;
  let steps = 0;
  while (steps < 100) {
    const next = step(x, n);
    steps++;
    if (next === x || next === prev) return { root: Math.min(x, next), steps };
    prev = x;
    x = next;
  }
  return { root: x, steps };
}

// Two stopping rules, both capped so a panel cannot hang.
function stepsUntilAbsolute(n, cap) {
  let x = n;
  let steps = 0;
  while (Math.abs(x * x - n) > 1e-10 && steps < cap) { x = step(x, n); steps++; }
  return steps;
}
function stepsUntilRelative(n, cap) {
  let x = n;
  let steps = 0;
  while (Math.abs(x * x - n) > n * Number.EPSILON && steps < cap) { x = step(x, n); steps++; }
  return steps;
}

// No guard at all: five steps, whatever the input.
function unguarded(n, steps) {
  let x = n;
  for (let i = 0; i < steps; i++) x = step(x, n);
  return x;
}

const trace = [1];
for (let i = 0; i < 5; i++) trace.push(step(trace[trace.length - 1], 2));
const errors = trace.map((x) => x - Math.sqrt(2));
// Quadratic convergence says the next error is about the square of this one,
// divided by twice the root.
const predicted = errors[3] ** 2 / (2 * Math.sqrt(2));

let identical = 0;
let worstGap = 0;
for (let n = 1; n <= 10000; n++) {
  const mine = heron(n).root;
  const theirs = Math.sqrt(n);
  if (mine === theirs) identical++;
  const gap = Math.abs(mine - theirs) / theirs / Number.EPSILON;
  if (gap > worstGap) worstGap = gap;
}

console.log('trace from 1     :', trace);
console.log('errors           :', errors);
console.log('predicted 5th    :', predicted, 'actual', errors[4]);
console.log('heron(2)         :', heron(2).root, 'in', heron(2).steps, 'steps, against', Math.sqrt(2));
console.log('steps, 1e10      :', heron(1e10).steps, 'with a guess', heronGuessed(1e10).steps);
console.log('bit identical    :', identical, 'of 10000, worst gap', worstGap, 'x EPSILON');
console.log('absolute rule    :', stepsUntilAbsolute(2e20, 500), 'steps (capped at 500)');
console.log('relative rule    :', stepsUntilRelative(2e20, 500), 'steps');
console.log('unguarded 0      :', unguarded(0, 5), 'and -4 gives', unguarded(-4, 5));