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