Algorithm for Horner's Method with base points in Julia

jakekap7 0
Viewed 17

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.

Thank you.

numerical-methods julia
1 Answers

  • 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
userm7ua0
Related