new package iterate

doc.go

@@ -86,6 +86,7 @@
 //	2. About Accuracy                                    accuracy
 //	3. Interpolation                                     interp
 //	4. Curve Fitting                                     fit
+//	5. Iteration                                         iterate
 //	7. Julian Day                                        julian
 //	8. Date of Easter                                    easter
 //	9. Jewish and Moslem Calendars                       jm

iterate/iterate.go

@@ -0,0 +1,77 @@
+// Copyright 2013 Sonia Keys
+// License MIT: http://www.opensource.org/licenses/MIT
+
+// Iterate: Chapter 5, Iteration.
+//
+// This package is best considered illustrative.  While the functions are
+// usable, they are minimal in showing the points of the chapter text.  More
+// robust functions would handle more cases of overflow, loss of precision,
+// and divergence.
+package iterate
+
+import (
+	"errors"
+	"math"
+)
+
+// BetterFunc is a convience type definition.
+type BetterFunc func(float64) float64
+
+// DecimalPlaces iterates to a fixed number of decimal places.
+//
+// Inputs are an improvement function, a starting value, the number of
+// decimal places desired in the result, and an iteration limit.
+func DecimalPlaces(better BetterFunc, start float64, places, maxIterations int) (float64, error) {
+	d := math.Pow(10, float64(-places))
+	for i := 0; i < maxIterations; i++ {
+		n := better(start)
+		if math.Abs(n-start) < d {
+			return n, nil
+		}
+		start = n
+	}
+	return 0, errors.New("Maximum iterations reached")
+}
+
+// FullPrecison iterates to (nearly) the full precision of a float64.
+//
+// To allow for a little bit of floating point jitter, FullPrecision iterates
+// to 15 significant figures, which is the maximum number of full significant
+// figures representable in a float64, but still a couple of bits shy of the
+// full representable precision.
+func FullPrecision(better BetterFunc, start float64, maxIterations int) (float64, error) {
+	for i := 0; i < maxIterations; i++ {
+		n := better(start)
+		if math.Abs((n-start)/n) < 1e-15 {
+			return n, nil
+		}
+		start = n
+	}
+	return 0, errors.New("Maximum iterations reached")
+}
+
+// RootFunc is a convience type definition.
+type RootFunc func(float64) float64
+
+// BinaryRoot finds a root between given bounds by binary search.
+//
+// Inputs are a function on x and the bounds on x.  A root must exist between
+// the given bounds, otherwise the result is not meaningful.
+func BinaryRoot(f RootFunc, lower, upper float64) float64 {
+	yLower := f(lower)
+	var mid float64
+	for j := 0; j < 52; j++ {
+		mid = (lower + upper) / 2
+		yMid := f(mid)
+		if yMid == 0 {
+			break
+		}
+		if math.Signbit(yLower) == math.Signbit(yMid) {
+			lower = mid
+			yLower = yMid
+		} else {
+			upper = mid
+		}
+	}
+	return mid
+}

iterate/iterate_test.go

@@ -0,0 +1,98 @@
+// Copyright 2013 Sonia Keys
+// License MIT: http://www.opensource.org/licenses/MIT
+
+package iterate_test
+
+import (
+	"fmt"
+	"math"
+
+	"github.com/soniakeys/meeus/iterate"
+)
+
+func ExampleDecimalPlaces() {
+	// Example 5.a, p. 48.
+	betterSqrt := func(N float64) iterate.BetterFunc {
+		return func(n float64) float64 {
+			return (n + N/n) / 2
+		}
+	}
+	start := 12.
+	places := 8
+	maxIter := 20
+	n, err := iterate.DecimalPlaces(betterSqrt(159), start, places, maxIter)
+	if err != nil {
+		fmt.Println(err)
+		return
+	}
+	fmt.Printf("%.*f", places, n)
+	// Output:
+	// 12.60952021
+}
+
+func ExampleFullPrecision() {
+	// Example 5.b, p. 48.
+	betterRoot := func() iterate.BetterFunc {
+		return func(x float64) float64 {
+			return (8 - math.Pow(x, 5)) / 17
+		}
+	}
+	x, err := iterate.FullPrecision(betterRoot(), 0, 20)
+	if err != nil {
+		fmt.Println(err)
+		return
+	}
+	fmt.Printf("%.9f\n", x)
+	fmt.Println(x, "(full precision)")
+	// Output:
+	// 0.469249878
+	// 0.4692498784547387 (full precision)
+}
+
+func ExampleFullPrecision_diverging() {
+	// Example 5.c, p. 49.
+	betterRoot := func() iterate.BetterFunc {
+		return func(x float64) float64 {
+			return (8 - math.Pow(x, 5)) / 3
+		}
+	}
+	x, err := iterate.FullPrecision(betterRoot(), 0, 20)
+	if err != nil {
+		fmt.Println(err)
+		return
+	}
+	fmt.Printf("%.9f\n", x)
+	fmt.Println(x, "(full precision)")
+	// Output:
+	// Maximum iterations reached
+}
+
+func ExampleFullPrecision_converging() {
+	// Example 5.d, p.49.
+	betterRoot := func() iterate.BetterFunc {
+		return func(x float64) float64 {
+			return math.Pow(8-3*x, .2)
+		}
+	}
+	x, err := iterate.FullPrecision(betterRoot(), 0, 30)
+	if err != nil {
+		fmt.Println(err)
+		return
+	}
+	fmt.Printf("%.9f\n", x)
+	fmt.Println(x, "(full precision)")
+	// Output:
+	// 1.321785627
+	// 1.321785627117658 (full precision)
+}
+
+func ExampleBinaryRoot() {
+	// Example from p. 53.
+	f := func(x float64) float64 {
+		return math.Pow(x, 5) + 17*x - 8
+	}
+	x := iterate.BinaryRoot(f, 0, 1)
+	fmt.Println(x)
+	// Output:
+	// 0.46924987845473876
+}

readme.md

@@ -17,6 +17,7 @@ Cross reference from chapters to package names
 	<tr><td>2.  About Accuracy</td><td>accuracy</td><td>complete</td></tr>
     <tr><td>3.  Interpolation</td><td>interp</td><td>under renovation</td></tr>
     <tr><td>4.  Curve Fitting</td><td>fit</td><td>complete</td></tr>
+    <tr><td>5.  Iteration</td><td>iterate</td><td>complete</td></tr>
     <tr><td>7.  Julian Day</td><td>julian</td><td>complete</td></tr>
     <tr><td>8.  Date of Easter</td><td>easter</td><td>complete</td></tr>
     <tr><td>9.  Jewish and Moslem Calendars</td><td>jm</td><td>complete</td></tr>