Newton's algorithm for a polynomial of arbitary degree

Question

Improvement over the last one I posted (now deleted, had no answers and was only cubics). Uses Horner's algorithm this time...

package main

import (
    "fmt"
    "math"
)

const ε float64 = 0.01

type polynomial struct {
    coeffs []float64
    x      float64
}

// Both of these functions use Horner's algorithm

func evalPoly(p polynomial, n int) float64 {
    var result float64 = 0
    for i := 0; i < n; i++ {
        result = result*p.x + p.coeffs[i]
    }
    return result
}

func evalPolyDeriv(p polynomial, n int) float64 {
    var result float64 = 0
    for i := 0; i < n-1; i++ {
        result = result*p.x + float64(n-i-1)*p.coeffs[i]
    }
    return result
}

func f(p polynomial) float64 {
    return evalPoly(p, len(p.coeffs)) / evalPolyDeriv(p, len(p.coeffs))
}

func solveNewton(p polynomial) float64 {
    h := f(p)
    for math.Abs(h) >= ε {
        p.x -= h
        h = f(p)
    }
    return p.x - h
}

func main() {
    guess := float64(2)
    coeffs := []float64{1, -2, 1, -4} // Can be of arbitary length
    fmt.Println(solveNewton(polynomial{coeffs: coeffs, x: guess})) // 2.3146...
}

Answer 1

In general, LGTM. Few notes:

• The name f is meaningless. newtonDelta perhaps?

• It feels right to compute a derivative's coefficients once, and reuse evalPoly for both the function and its derivative. After all, the polynomial's derivative is also polynomial. There is an usual space-time tradeoff, but in the case of polynomials DRY rule casts a deciding vote.

• Newton's algorithm does not necessarily converge. You should be prepared to handle the divergent case.

• Hardcoding ε is dubious. I recommend to solveNewton have it as a parameter.

• Using Horner schedule is a definite improvement.

Answer 2

type polynomial struct {
    coeffs []float64
    x      float64
}

The name is misleading: it would be appropriate for

type polynomial struct {
    coeffs []float64
}

Also, this very much needs either a comment or some Hungarian notation to indicate the endianness of the coefficient array.

---

func evalPoly(p polynomial, n int) float64 {
    var result float64 = 0
    for i := 0; i < n; i++ {
        result = result*p.x + p.coeffs[i]
    }
    return result
}

Why is n a parameter? Shouldn't it be

func evalPoly(p polynomial, x float64) float64 {
    var result float64 = 0
    for i := 0; i < len(p.coeffs); i++ {
        result = result*x + p.coeffs[i]
    }
    return result
}

(maybe pulling out len(p.coeffs) to a local variable)?

---

func solveNewton(p polynomial) float64 {
    ...
        p.x -= h
    ...
    return p.x - h
}

This reinforces my first point about the polynomial type: the result is returned as a return value, but a very close approximation of it is also returned as a side-effect, modifying the argument! There is a case to be made that in computer algebra systems the majority (if not all) of the objects should be immutable, and the functions pure.

---

I will not repeat the suggestions vnp has made, but I agree with them.

Answer 3

vnp and Peter Taylors suggestions are very good. One could go a step further and make nice methods on the Polynomial type: Evaluate, Derivative and SolveNewton for instance:

package main

import (
    "fmt"
    "math"
)

// Polynomial represents a polynomial
type Polynomial struct {
    // Coeffs are stored so that Coeffs[len(Coeffs)-1] is the constant
    Coeffs []float64
}

// Evaluate p using Horner's algorithm
func (p Polynomial) Evaluate(x float64) (result float64) {
    for _, coeff := range p.Coeffs {
        result = result*x + coeff
    }
    return result
}

// Derivative returns the derivative of the polynomial
func (p Polynomial) Derivative() (d Polynomial) {
    l := len(p.Coeffs)
    if l <= 1 {
        return d
    }
    d.Coeffs = make([]float64, l-1)
    for i, coeff := range p.Coeffs[:l-1] {
        d.Coeffs[i] = float64(l-i-1) * coeff
    }
    return d
}

// SolveNewton finds a zero using Newton's algorithm, starting at firstGuess until precision is reached
func (p Polynomial) SolveNewton(firstGuess, precision float64) float64 {
    d := p.Derivative()

    newtonDelta := func(x float64) float64 {
        return p.Evaluate(x) / d.Evaluate(x)
    }
    guess := firstGuess
    h := newtonDelta(guess)
    for math.Abs(h) >= precision {
        guess -= h
        h = newtonDelta(guess)
    }
    return guess - h
}

func main() {
    poly := Polynomial{Coeffs: []float64{1, -2, 1, -4}}
    fmt.Println(poly.SolveNewton(float64(2), float64(0.01))) // 2.3146...
}