Search code examples
julianumerical-methods

Algorithm for Horner's Method with base points in Julia


I'm attempting to write an algorithm in Julia for Horner's method with base points (b1, b2, etc.) to evaluate the polynomial: y = p(x) = c1 + c2(x-b1) + c3(x-b1)(x-b2) + c4(x-b1)(x-b2)(x-b3), where the coefficients are represented by "c" and the base points are represented by "b". The polynomial that I need to evaluate is: p(x) = 10 - (21/2)(x+2) + (45/6)(x+2)(x) - (7/3)(x+2)(x)(x-1).

Here is my code so far:

function hornerpolybase(x,c,b)
    
    # assign y to the first coefficient
    # make y the same type as input x
    y = c[1] * one(x) 
    
    n = length(c)
    
    for i = 2:n
        y = y + c[i] * (x - b[i-1])   
    end
    
    return y
end

My issue is that after the second iteration of the for loop, y = c1 + c2(x-b1) + c3(x-b2), whereas the equation should be y = c1 + c2(x-b1) + c3 * (x-b1)(x-b2). The problem repeats in the third iteration. I'm thinking that possibly an "if, then" conditional statement within the for loop should solve this issue, but I'm completely stuck. I would appreciate any help.


Solution

    • The code isn't fully type stable, despite the way y is initialized. E.g. try @code_warntype hornerpolybase(1, [1,1], [1.0]). You could fix it by including the type of b:

      # assign y to the first coefficient
      T = promote_type(typeof(x), eltype(c), eltype(b))
      y::T = c[1] # implicit conversion to T
      

      (evalpoly in Julia actually has the same problem...)

    • You're accessing c in reverse (if c[1] is supposed to correspond with c_1). You need to start with the last coefficient in Horner's method.

    • In each iteration it's y that should be multiplied, not the coefficient:

      y::T = c[end]
      
      n = length(c)
      for i = n-1:-1:1
          y = c[i] + y * (x - b[i])   
      end