# Newton's algorithm for a polynomial of arbitary degree

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...
}


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.

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.

• And even further with the evalPoly: func (p polynomial) eval(x float64) float64 (and as @vnp suggested, maybe a func (p polynomial) derivative() polynomial method as well) – oliverpool Nov 20 '17 at 8:02

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...
}

• I actually ended up doing this, but in C# (ewww) because the friend I was explaining this for wanted it... – alchzh Nov 20 '17 at 14:44