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

check_determ
domains
  reals = real*
predicates
  solve_polynomial(reals,real,real,real)
  solve_polynomial_1(real,reals,reals,real,real)
  solve_polynomial_2(real,real,reals,reals,real,real)
  differential(reals,reals)
  differential_1(reals,real,integer,reals)
  evaluate_polynomial(reals,real,real)
  evaluate_polynomial_1(reals,real,real,real)
  length(reals,integer)
  length_1(reals,integer,integer)
clauses

/* 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 = 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 = I * A0,
  I1 = 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 = 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 = L0 + 1, length_1(Xs, L1, L).

Turbo Prolog Index     Home Page