summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--lib/cbrt.go68
-rw-r--r--lib/cbrt_test.go23
2 files changed, 91 insertions, 0 deletions
diff --git a/lib/cbrt.go b/lib/cbrt.go
new file mode 100644
index 0000000..69ddc17
--- /dev/null
+++ b/lib/cbrt.go
@@ -0,0 +1,68 @@
+// Copyright © 2021 siddharth ravikumar <s@ricketyspace.net>
+// SPDX-License-Identifier: ISC
+
+package lib
+
+import (
+ "math/big"
+)
+
+// Cube root tolerance.
+var bigCubeRootTolerance = big.NewFloat(0.00001)
+
+// Returns cube root of a.
+//
+// Uses Newton's method.
+// https://en.wikipedia.org/wiki/Newton's_method
+func BigCubeRoot(a *big.Float) *big.Float {
+ // If x^3 = a, then our f(x) is:
+ // f(x) = x^3 - a
+ fx := func(x *big.Float) *big.Float {
+ // x^3
+ e := big.NewFloat(0)
+ e = e.Mul(x, x)
+ e = e.Mul(e, x)
+
+ // x^2 - a
+ z := big.NewFloat(0).Sub(e, a)
+
+ return z
+ }
+
+ // f'(x) is:
+ // f'(x) = 3 * x^2
+ fxPrime := func(x *big.Float) *big.Float {
+ // x^2
+ x2 := big.NewFloat(0).Mul(x, x)
+
+ // 3 * x^2
+ z := big.NewFloat(0).Mul(big.NewFloat(3), x2)
+
+ return z
+ }
+
+ x0 := a // Initial guess.
+ max := 1000 // Max iterations.
+ i := 0 // Current iteration.
+ for i < max {
+ // f(x0) / f'(x0)
+ d := fx(x0)
+ d = d.Quo(d, fxPrime(x0))
+
+ // x0 - ( f(x0) / f'(x0) )
+ x1 := big.NewFloat(0).Set(x0)
+ x1 = x1.Sub(x1, d)
+
+ // x0 - x1
+ df := big.NewFloat(0).Set(x0)
+ df = df.Sub(df, x1)
+ df = df.Abs(df)
+ if df.Cmp(bigCubeRootTolerance) == -1 {
+ return x1
+ }
+
+ i += 1
+ x0 = x1
+ }
+ return nil
+}
diff --git a/lib/cbrt_test.go b/lib/cbrt_test.go
new file mode 100644
index 0000000..a5722c8
--- /dev/null
+++ b/lib/cbrt_test.go
@@ -0,0 +1,23 @@
+// Copyright © 2021 siddharth ravikumar <s@ricketyspace.net>
+// SPDX-License-Identifier: ISC
+
+package lib
+
+import (
+ "math/big"
+ "testing"
+)
+
+func TestBigCubeRoot(t *testing.T) {
+ a := big.NewFloat(612)
+ acr := BigCubeRoot(a)
+ if acr == nil {
+ t.Errorf("Could not find cube root of %v\n", a)
+ return
+ }
+ expected := big.NewFloat(8.490184748)
+ if big.NewFloat(0).Sub(acr, expected).Cmp(bigCubeRootTolerance) != -1 {
+ t.Errorf("Could not find cube root of %v (%v)\n", a, acr)
+ return
+ }
+}