3
\$\begingroup\$

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...
}
\$\endgroup\$
5
\$\begingroup\$

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.

\$\endgroup\$
2
\$\begingroup\$
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.

\$\endgroup\$
  • \$\begingroup\$ 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) \$\endgroup\$ – oliverpool Nov 20 '17 at 8:02
0
\$\begingroup\$

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...
}
\$\endgroup\$
  • \$\begingroup\$ I actually ended up doing this, but in C# (ewww) because the friend I was explaining this for wanted it... \$\endgroup\$ – alchzh Nov 20 '17 at 14:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.