summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--lib/cbrt.go56
-rw-r--r--lib/cbrt_test.go14
2 files changed, 70 insertions, 0 deletions
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
+ }
+}