diff options
Diffstat (limited to 'lib')
| -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 +	} +} | 
