From 734416dde97d6a7ed91bd5bbb2f1dcb68edcb80b Mon Sep 17 00:00:00 2001
From: Lumi Pakkanen <lumi.pakkanen@gmail.com>
Date: Mon, 25 Dec 2023 10:05:47 +0200
Subject: [PATCH] Implement radical approximation

---
 src/__tests__/approximation.spec.ts | 15 +++++++
 src/approximation.ts                | 61 +++++++++++++++++++++++++++++
 2 files changed, 76 insertions(+)

diff --git a/src/__tests__/approximation.spec.ts b/src/__tests__/approximation.spec.ts
index 2e75780..23e41f7 100644
--- a/src/__tests__/approximation.spec.ts
+++ b/src/__tests__/approximation.spec.ts
@@ -3,6 +3,7 @@ import {
   approximateOddLimit,
   approximatePrimeLimit,
   approximatePrimeLimitWithErrors,
+  approximateRadical,
   getConvergents,
 } from '../approximation';
 import {valueToCents} from '../conversion';
@@ -114,3 +115,17 @@ describe('Convergent calculator', () => {
     expect(semiconvergents[4].equals('13/4')).toBeTruthy();
   });
 });
+
+describe('Radical approximator', () => {
+  it("finds Ramanujan's approximation to pi", () => {
+    const {index, radicant} = approximateRadical(Math.PI);
+    expect(index).toBe(4);
+    expect(radicant.toFraction()).toBe('2143/22');
+  });
+
+  it('works with a random value without crashing', () => {
+    const value = Math.random() * 1000 - 100;
+    const {index, radicant} = approximateRadical(value);
+    expect(radicant.valueOf() ** (1 / index) / value).toBeCloseTo(1);
+  });
+});
diff --git a/src/approximation.ts b/src/approximation.ts
index fea696c..169640e 100644
--- a/src/approximation.ts
+++ b/src/approximation.ts
@@ -269,3 +269,64 @@ export function approximatePrimeLimit(
     maxLength
   ).map(result => result[0]);
 }
+
+/**
+ * Calculate an array of continued fraction elements representing the value.
+ * https://en.wikipedia.org/wiki/Continued_fraction
+ * @param value Value to turn into a continued fraction.
+ * @returns An array of continued fraction elements.
+ */
+export function continuedFraction(value: number) {
+  const result: number[] = [];
+  let coeff = Math.floor(value);
+  result.push(coeff);
+  value -= coeff;
+  while (value && coeff < 1e12 && result.length < 32) {
+    value = 1 / value;
+    coeff = Math.floor(value);
+    result.push(coeff);
+    value -= coeff;
+  }
+  return result;
+}
+
+/**
+ * Approximate a value with a radical expression.
+ * @param value Value to approximate.
+ * @param maxIndex Maximum index of the radical. 2 means square root, 3 means cube root, etc.
+ * @param maxHeight Maximum Benedetti height of the radicant in the approximation.
+ * @returns Object with index of the radical and the radicant. Result is "index'th root or radicant".
+ */
+export function approximateRadical(
+  value: number,
+  maxIndex = 5,
+  maxHeight = 50000
+) {
+  let index = 1;
+  let radicant = new Fraction(1);
+  let bestError = Math.abs(value - 1);
+  for (let i = 1; i <= maxIndex; ++i) {
+    const cf = continuedFraction(value ** i);
+    let candidate: Fraction | undefined;
+    for (let j = 0; j < cf.length - 1; ++j) {
+      let convergent = new Fraction(cf[j]);
+      for (let k = j - 1; k >= 0; --k) {
+        convergent = convergent.inverse().add(cf[k]);
+      }
+      if (convergent.n * convergent.d > maxHeight) {
+        break;
+      }
+      candidate = convergent;
+    }
+    if (candidate !== undefined) {
+      const error = Math.abs(candidate.valueOf() ** (1 / i) - value);
+      if (error < bestError) {
+        index = i;
+        radicant = candidate;
+        bestError = error;
+      }
+    }
+  }
+
+  return {index, radicant};
+}