diff options
Diffstat (limited to 'lib/cbrt.go')
-rw-r--r-- | lib/cbrt.go | 56 |
1 files changed, 56 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 +} |