From 73b0bec6c8c2aaf0f8b4e22a60a925881b202a15 Mon Sep 17 00:00:00 2001
From: siddharth ravikumar <s@ricketyspace.net>
Date: Sun, 20 Nov 2022 16:07:47 -0500
Subject: lib: add `BigIntCubeRoot`

---
 lib/cbrt.go      | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 lib/cbrt_test.go | 14 ++++++++++++++
 2 files changed, 70 insertions(+)

diff --git a/lib/cbrt.go b/lib/cbrt.go
index 6a856a3..122f860 100644
--- a/lib/cbrt.go
+++ b/lib/cbrt.go
@@ -63,3 +63,59 @@ func BigCubeRoot(a *big.Float) *big.Float {
 	}
 	return nil
 }
+
+// Returns cube root of a.
+//
+// Uses Newton's method.
+// https://en.wikipedia.org/wiki/Newton's_method
+func BigIntCubeRoot(a *big.Int) *big.Int {
+	// If x^3 = a, then our f(x) is:
+	//     f(x) = x^3 - a
+	fx := func(x *big.Int) *big.Int {
+		// x^3
+		e := big.NewInt(0)
+		e = e.Mul(x, x)
+		e = e.Mul(e, x)
+
+		// x^3 - a
+		z := big.NewInt(0).Sub(e, a)
+
+		return z
+	}
+
+	// f'(x) is:
+	//    f'(x) = 3 * x^2
+	fxPrime := func(x *big.Int) *big.Int {
+		// x^2
+		x2 := big.NewInt(0).Mul(x, x)
+
+		// 3 * x^2
+		z := big.NewInt(0).Mul(big.NewInt(3), x2)
+
+		return z
+	}
+
+	x0 := a // Initial guess.
+	zero := big.NewInt(0)
+	one := big.NewInt(1)
+	for i := 0; i < 1000000; i++ {
+		// f(x0) / f'(x0)
+		d := fx(x0)
+		d = d.Div(d, fxPrime(x0))
+
+		// x0 - ( f(x0) / f'(x0) )
+		x1 := big.NewInt(0).Set(x0)
+		x1 = x1.Sub(x1, d)
+
+		// x0 - x1
+		df := big.NewInt(0).Set(x0)
+		df = df.Sub(df, x1)
+		df = df.Abs(df)
+		if df.Cmp(zero) == 0 {
+			return x1.Sub(x1, one)
+		}
+
+		x0 = x1
+	}
+	return nil
+}
diff --git a/lib/cbrt_test.go b/lib/cbrt_test.go
index a5722c8..6b3c9f3 100644
--- a/lib/cbrt_test.go
+++ b/lib/cbrt_test.go
@@ -21,3 +21,17 @@ func TestBigCubeRoot(t *testing.T) {
 		return
 	}
 }
+
+func TestBigIntCubeRoot(t *testing.T) {
+	a := big.NewInt(19683)
+	acr := BigIntCubeRoot(a)
+	if acr == nil {
+		t.Errorf("Could not find cube root of %v\n", a)
+		return
+	}
+	expected := big.NewInt(27)
+	if big.NewInt(0).Sub(acr, expected).Cmp(big.NewInt(0)) != 0 {
+		t.Errorf("Could not find cube root of %v (%v)\n", a, acr)
+		return
+	}
+}
-- 
cgit v1.2.3