Algorithm to calculate the Arithmetic Geometric Mean

Maths/ArithmeticGeometricMean.js

@@ -0,0 +1,28 @@
+/**
+ * @function agm
+ * @description This finds the Arithmetic-Geometric Mean between any 2 numbers.
+ * @param {Number} a - 1st number, also used to store Arithmetic Mean.
+ * @param {Number} g - 2nd number, also used to store Geometric Mean.
+ * @return {Number} - AGM of both numbers.
+ * @see [AGM](https://en.wikipedia.org/wiki/Arithmetic%E2%80%93geometric_mean)
+ */
+
+export const agm = (a, g) => {
+  if (a === Infinity && g === 0) return NaN
+  if (Object.is(a, -0) && !Object.is(g, -0)) return 0
+  if (a === g) return a // avoid rounding errors, and increase efficiency
+  let x // temp var
+  do {
+    [a, g, x] = [(a + g) / 2, Math.sqrt(a * g), a]
+  } while (a !== x && !isNaN(a))
+  /*
+  `x !== a` ensures the return value has full precision,
+  and prevents infinite loops caused by rounding differences between `div` and `sqrt` (no need for "epsilon").
+  If we were to compare `a` with `g`, some input combinations (not all) can cause an infinite loop,
+  because the rounding mode never changes at runtime.
+  Precision is not the same as accuracy, but they're related.
+  This function isn't always 100% accurate (round-errors), but at least is more than 95% accurate.
+  `!isNaN(x)` prevents infinite loops caused by invalid inputs like: negatives, NaNs and Infinities.
+  */
+  return a
+}

Maths/test/ArithmeticGeometricMean.test.js

@@ -0,0 +1,56 @@
+import { agm } from '../ArithmeticGeometricMean.js'
+
+describe('Tests for AGM', () => {
+  it('should be a function', () => {
+    expect(typeof agm).toEqual('function')
+  })
+
+  it('number of parameters should be 2', () => {
+    expect(agm.length).toEqual(2)
+  })
+
+  const m = 0x100 // scale for rand
+
+  it('should return NaN if any or all params has a negative argument', () => {
+    // I multiplied by minus one, because the sign inversion is more clearly visible
+    expect(agm(-1 * Math.random() * m, Math.random() * m)).toBe(NaN)
+    expect(agm(Math.random() * m, -1 * Math.random() * m)).toBe(NaN)
+    expect(agm(-1 * Math.random() * m, -1 * Math.random() * m)).toBe(NaN)
+  })
+
+  it('should return Infinity if any arg is Infinity and the other is not 0', () => {
+    expect(agm(Math.random() * m + 1, Infinity)).toEqual(Infinity)
+    expect(agm(Infinity, Math.random() * m + 1)).toEqual(Infinity)
+    expect(agm(Infinity, Infinity)).toEqual(Infinity)
+  })
+
+  it('should return NaN if some arg is Infinity and the other is 0', () => {
+    expect(agm(0, Infinity)).toBe(NaN)
+    expect(agm(Infinity, 0)).toBe(NaN)
+  })
+
+  it('should return +0 if any or all args are +0 or -0, and return -0 if all are -0', () => {
+    expect(agm(Math.random() * m, 0)).toBe(0)
+    expect(agm(0, Math.random() * m)).toBe(0)
+    expect(agm(Math.random() * m, -0)).toBe(0)
+    expect(agm(-0, Math.random() * m)).toBe(0)
+    expect(agm(0, -0)).toBe(0)
+    expect(agm(-0, 0)).toBe(0)
+    expect(agm(-0, -0)).toBe(-0)
+  })
+
+  it('should return NaN if any or all args are NaN', () => {
+    expect(agm(Math.random() * m, NaN)).toBe(NaN)
+    expect(agm(NaN, Math.random() * m)).toBe(NaN)
+    expect(agm(NaN, NaN)).toBe(NaN)
+  })
+
+  it('should return an accurate approximation of the AGM between 2 valid input args', () => {
+    // all the constants are provided by WolframAlpha
+    expect(agm(1, 2)).toBeCloseTo(1.4567910310469068)
+    expect(agm(2, 256)).toBeCloseTo(64.45940719438667)
+    expect(agm(55555, 34)).toBeCloseTo(9933.4047239552)
+    // test "unsafe" numbers
+    expect(agm(2 ** 48, 3 ** 27)).toBeCloseTo(88506556379265.7)
+  })
+})