// Copyright 2014 The Go Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. // This file is adapted from https://github.com/robpike/ivy/blob/master/value/loop.go. package apd import ( "fmt" "math" ) type loop struct { c *Context name string // The name of the function we are evaluating. i uint64 // Loop count. precision int32 maxIterations uint64 // When to give up. arg *Decimal // original argument to function; only used for diagnostic. prevZ Decimal // Result from the previous iteration. delta Decimal // |Change| from previous iteration. } const digitsToBitsRatio = math.Ln10 / math.Ln2 // newLoop returns a new loop checker. Arguments: // - name: name of the function being calculated (for error messages) // - arg: argument to the function (for error messages) // - precision: desired precision; the loop ends when consecutive estimates // differ less than the desired precision. Note that typically // the inner computations in an iteration need higher precision, // so this is normally lower than the precision in the context. // - maxItersPerDigit: after this many iterations per digit of precision, the // loop ends in error. func (c *Context) newLoop(name string, arg *Decimal, precision uint32, maxItersPerDigit int) *loop { return &loop{ c: c, name: name, arg: new(Decimal).Set(arg), precision: int32(precision), maxIterations: 10 + uint64(maxItersPerDigit*int(precision)), } } // done reports whether the loop is done. If it does not converge // after the maximum number of iterations, it returns an error. func (l *loop) done(z *Decimal) (bool, error) { if _, err := l.c.Sub(&l.delta, &l.prevZ, z); err != nil { return false, err } sign := l.delta.Sign() if sign == 0 { return true, nil } if sign < 0 { // Convergence can oscillate when the calculation is nearly // done and we're running out of bits. This stops that. // See next comment. l.delta.Neg(&l.delta) } // We stop if the delta is smaller than a change of 1 in the // (l.precision)-th digit of z. Examples: // // p = 4 // z = 12345.678 = 12345678 * 10^-3 // eps = 10.000 = 10^(-4+8-3) // // p = 3 // z = 0.001234 = 1234 * 10^-6 // eps = 0.00001 = 10^(-3+4-6) var eps Decimal eps.Coeff.Set(bigOne) eps.Exponent = -l.precision + int32(z.NumDigits()) + z.Exponent if l.delta.Cmp(&eps) <= 0 { return true, nil } l.i++ if l.i == l.maxIterations { return false, fmt.Errorf( "%s %s: did not converge after %d iterations; prev,last result %s,%s delta %s precision: %d", l.name, l.arg.String(), l.maxIterations, z.String(), l.prevZ.String(), l.delta.String(), l.precision, ) } l.prevZ.Set(z) return false, nil }