diff options
-rw-r--r-- | lib/cbrt.go | 56 | ||||
-rw-r--r-- | lib/cbrt_test.go | 14 |
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 + } +} |