#### Roots of a Polynomial

This uses Newton's method to find a root of a polynomial.
A polynomial is represented by a list of its coefficients,
high-order first.
For example, x^3 + 3*x^2 - 2*x + 4 is represented by [1,3,-2,4].

/* solve_polynomial(Polynomial, Delta, X0, X) is true if X is a root of */
/* the Polynomial, X0 is an initial guess for the root, and Delta is the */
/* absolute value of the difference between successive estimates of the */
/* root which causes iteration to terminate. */
/* e.g. solve_polynomial([1000, -30, -1071], 0.001, 1, X) gives X = 1.05 */
solve_polynomial(As, Delta, X0, X):-
differential(As, Bs),
solve_polynomial_1(X0, As, Bs, Delta, X).
solve_polynomial_1(X0, As, Bs, Delta, X):-
evaluate_polynomial(As, X0, F),
evaluate_polynomial(Bs, X0, Fprime),
X1 is X0 - F / Fprime,
solve_polynomial_2(X1, X0, As, Bs, Delta, X).
solve_polynomial_2(X1, X0, As, Bs, Delta, X):-
abs(X1 - X0) > Delta,
!,
solve_polynomial_1(X1, As, Bs, Delta, X).
solve_polynomial_2(X, _, _, _, _, X).
/* differential(Polynomial, Differential) is true if Differential is the */
/* result of differentiating the Polynomial. */
/* e.g. differential([95, -3, -3, -103], [285, -6, -3]). */
differential([A|As], Bs):-
length(As, Length),
differential_1(As, A, Length, Bs).
differential_1([], _, _, []).
differential_1([A|As], A0, I, [B|Bs]):-
B is I * A0,
I1 is I - 1,
differential_1(As, A, I1, Bs).
/* evaluate_polynomial(Polynomial, X, Y) is true if Y is the value of the */
/* Polynomial evaluated for the value X. */
/* e.g. evaluate_polynomial([2,3,4,5], 2, 41). */
evaluate_polynomial([A|As], X, Y):-
evaluate_polynomial_1(As, X, A, Y).
evaluate_polynomial_1([], _, Y, Y).
evaluate_polynomial_1([A|As], X, Y0, Y):-
Y1 is Y0 * X + A,
evaluate_polynomial_1(As, X, Y1, Y).
/* length(Xs, L) is true if L is the number of elements in the list Xs. */
%length(Xs, L):-length_1(Xs, 0, L).
%length_1([], L, L).
%length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).

LPA Index
Home Page