Sparse Polynomial Operations

A sparse polynomial is represented by a list of terms, each term being of the form t(PowerOfX, Coefficient), where PowerOfX is not necessarily an integer value, and the powers are in increasing order of sequence (at least for those operations using directly or indirectly add/3 or subtract/3). No terms are included for coefficients of value 0.

/* Calculates Goldbach numbers, that is, the number of representations of  */
/*   even numbers (greater than 4) as the sum of two primes.               */
/* (With this data, the values are correct for even numbers up to 100.)    */
goldbach(X):-  
  S = [t( 3,1),t( 5,1),t( 7,1),t(11,1),t(13,1),t(17,1),t(19,1),t(23,1),
       t(29,1),t(31,1),t(37,1),t(41,1),t(43,1),t(47,1),t(53,1),t(59,1),
       t(61,1),t(67,1),t(71,1),t(73,1),t(79,1),t(83,1),t(89,1),t(97,1)],
  multiply(S, S, X).

/* multiply(Xs, Ys, Zs) is true if Zs is the result of multiplying the     */
/*   polynomial Ys by the polynomial Xs.                                   */
multiply(Xs, Ys, Zs):-
  multiply_1(Xs, Ys, [], Zs).
  
multiply_1([], _, Zs, Zs).  
multiply_1([X|Xs], Ys, Zs0, Zs):-
  multiply_2(Ys, X, Zs1),
  add(Zs0, Zs1, Zs2),
  multiply_1(Xs, Ys, Zs2, Zs).
  
multiply_2([], _, []).
multiply_2([t(N,B)|Ys], t(M,A), [t(P,C)|Zs]):-
  P is N + M,
  C is B * A,
  multiply_2(Ys, t(M,A), Zs).

/* add(Xs, Ys, Zs) is true if Zs is the result of adding the polynomial Ys */
/*   to the polynomial Xs.                                                 */
add([], Ys, Ys). 
add([X|Xs], [], [X|Xs]). 
add([t(M,A)|Xs], [t(M,B)|Ys], [t(M,C)|Zs]):-
  C is A + B, C =\= 0, !, add(Xs, Ys, Zs).
add([t(M,_)|Xs], [t(M,_)|Ys], Zs):-
  !, add(Xs, Ys, Zs).
add([t(M,A)|Xs], [t(N,B)|Ys], [t(M,A)|Zs]):-
  M < N, !, add(Xs, [t(N,B)|Ys], Zs).
add([X|Xs], [Y|Ys], [Y|Zs]):-
  add([X|Xs], Ys, Zs).

/* subtract(Xs, Ys, Zs) is true if Zs is the result of subtracting the     */
/*   polynomial Ys from the polynomial Xs.                                 */
subtract([], [], []). 
subtract([X|Xs], [], [X|Xs]). 
subtract([t(M,A)|Xs], [t(M,B)|Ys], [t(M,C)|Zs]):-
  C is A - B, C =\= 0, !, subtract(Xs, Ys, Zs).
subtract([t(M,_)|Xs], [t(M,_)|Ys], Zs):-
  !, subtract(Xs, Ys, Zs).
subtract([t(M,A)|Xs], [t(N,B)|Ys], [t(M,A)|Zs]):-
  M < N, !, subtract(Xs, [t(N,B)|Ys], Zs).
subtract(Xs, [t(N,B)|Ys], [t(N,B1)|Zs]):-
  B1 is -B, subtract(Xs, Ys, Zs).
  
/* evaluate_polynomial(Polynomial, X, Y) is true if Y is the value of the  */
/*   Polynomial evaluated for the value X.                                 */
/* e.g. evaluate_polynomial([t(0,5),t(1,4),t(2,3),t(3,2)], 2, 41).         */
/*                               5 +  4*x + 3*x^2 +2*x^3                   */
evaluate_polynomial(As, X, Y):-
  evaluate_polynomial_1(As, X, 0, Y).
  
evaluate_polynomial_1([], _, Y, Y).
evaluate_polynomial_1([t(N,A)|As], X, Y0, Y):-
  Y1 is Y0 + A * aln(N * ln(X)),
  evaluate_polynomial_1(As, X, Y1, Y).

/* differentiate_polynomial(Polynomial, Differential) is true if           */
/*   Differential is the result of differentiating the Polynomial.         */
/* e.g. differentiate_polynomial([t(0,5),t(1,4),t(3,2)], [t(0,4),t(2,6)]). */
differentiate_polynomial([], []).
differentiate_polynomial([t(0,_)|As], Bs):-
  !,
  differentiate_polynomial(As, Bs).
differentiate_polynomial([t(M,A)|As], [t(N,B)|Bs]):-
  N is M - 1,
  B is M * A,
  differentiate_polynomial(As, Bs).
  
/* 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([t(0,-1071),t(1,-30),t(2,1000)], 0.001, 1, 1.05). */
solve_polynomial(As, Delta, X0, X):-
  differentiate_polynomial(As, Bs),
  solve_polynomial_1(As, Bs, Delta, X0, X).

solve_polynomial_1(As, Bs, Delta, X0, X):-
  evaluate_polynomial(As, X0, F),
  evaluate_polynomial(Bs, X0, Fprime),
  X1 is X0 - F / Fprime,
  solve_polynomial_2(As, Bs, Delta, X0, X1, X).

solve_polynomial_2(As, Bs, Delta, X0, X1, X):-
  abs(X0 - X1) > Delta,
  !,
  solve_polynomial_1(As, Bs, Delta, X1, X).
solve_polynomial_2(_, _, _, _, X, X).

LPA Index     Home Page