muller_recurrence.go

// Muller's Recurrence in Go.
//
// https://latkin.org/blog/2014/11/22/mullers-recurrence-roundoff-gone-wrong/
package main

import (
	"fmt"
	"math/big"
	
	"github.com/cockroachdb/apd/v2"
)

func main() {
	fmt.Println("== float64 ==")
	runFloat(21)
	fmt.Println("")
	fmt.Println("== *big.Rat ==")
	runRat(101)
	fmt.Println("")
	fmt.Println("== *apd.Decimal ==")
	runDec(101)
}

func runFloat(ct int) {
	values := make([]float64, ct)
	values[0] = 4
	values[1] = 4.25
	for i := 2; i < len(values); i++ {
		v := fFloat(values[i-1], values[i-2])
		values[i] = v

		fmt.Printf("% 5d: %f\n", i, v)
	}
}

func fFloat(y, z float64) float64 {
	return 108 - (815-1500/z)/y
}

func runRat(ct int) {
	values := make([]big.Rat, ct)
	x0 := big.NewRat(4, 1)
	x1 := big.NewRat(425, 100)
	values[0] = *x0
	values[1] = *x1

	for i := 2; i < len(values); i++ {
		v := fRat(&values[i], &values[i-1], &values[i-2])

		fmt.Printf("% 5d: %s\n", i, v.FloatString(200))
	}
}

var (
	r1 = big.NewRat(108, 1)
	r2 = big.NewRat(815, 1)
	r3 = big.NewRat(1500, 1)
)

func fRat(x, y, z *big.Rat) *big.Rat {
	// return 108 - (815-1500/z)/y
	// return r1 - (r2-r3/z)/y
	x.Quo(r3, z)
	x.Sub(r2, x)
	x.Quo(x, y)
	x.Sub(r1, x)
	return x
}

func runDec(ct int) {
	values := make([]apd.Decimal, ct)
	x0 := apd.New(4, 0)
	x1 := apd.New(425, -2)
	values[0] = *x0
	values[1] = *x1
	
	fmt.Printf("dec x0: %f\n", x0)
	fmt.Printf("dec x1: %f\n", x1)

	for i := 2; i < len(values); i++ {
		v := fDec(&values[i], &values[i-1], &values[i-2])

		fmt.Printf("% 5d: %f\n", i, v)
	}
}

var (
	d1 = apd.New(108, 0)
	d2 = apd.New(815, 0)
	d3 = apd.New(1500, 0)
	dc = apd.BaseContext.WithPrecision(200)
)

func fDec(x, y, z *apd.Decimal) *apd.Decimal {
	// return 108 - (815-1500/z)/y
	// return d1 - (d2-d3/z)/y
	dc.Quo(x, d3, z)
	dc.Sub(x, d2, x)
	dc.Quo(x, x, y)
	dc.Sub(x, d1, x)
	return x
}