Polynomial Operations

/* A polynomial is represented by a list of its coefficients, high-order   */
/*   first.  e.g. x^3 + 3*x^2 - 2*x + 4 is represented by [1,3,-2,4].      */

/* add(Xs, Ys, Zs) is true if adding the polynomial Xs to the polynomial   */
/*   Ys gives the polynomial Zs.                                           */
/* e.g. add([1,2,3], [4,4,4,4], [4,5,6,7]).                                */
add(Xs, Ys, Zs):-
  reverse(Xs, ReversedXs),
  reverse(Ys, ReversedYs),
  add_1(ReversedXs, ReversedYs, ReversedZs),
  reverse(ReversedZs, Zs0),
  trim_left(Zs0, 0, Zs).

add_1([], Ys, Ys).
add_1([X|Xs], [], [X|Xs]):-!.
add_1([X|Xs], [Y|Ys], [Z|Zs]):-
  Z is X + Y,
  add_1(Xs, Ys, Zs).

/* subtract(Xs, Ys, Zs) is true if subtracting the polynomial Ys from the  */
/*   polynomial Xs gives the polynomial Zs.                                */
/* e.g. subtract([4,5,6,7], [1,2,3], [4,4,4,4]).                           */
subtract(Xs, Ys, Zs):-
  reverse(Xs, ReversedXs),
  reverse(Ys, ReversedYs),
  subtract_1(ReversedXs, ReversedYs, ReversedZs),
  reverse(ReversedZs, Zs0),
  trim_left(Zs0, 0, Zs).

subtract_1([], Ys, Ys).
subtract_1([X|Xs], [], [X|Xs]).
subtract_1([X|Xs], [Y|Ys], [Z|Zs]):-
  Z is X - Y,
  subtract_1(Xs, Ys, Zs).

/* multiply(Xs, Ys, Zs) is true if multiplying the polynomial Xs by the    */
/*   polynomial Ys gives the polynomial Zs.                                */
/* e.g. multiply([1,2,3], [4,5,6,7], [4,13,28,34,32,21]).                  */
multiply(Xs, Ys, Zs):-
  reverse(Ys, ReversedYs),
  multiply_1(Xs, ReversedYs, [], ReversedZs),
  reverse(ReversedZs, Zs).
  
multiply_1([], _, Zs, Zs).  
multiply_1([X|Xs], Ys, Zs0, Zs):-
  multiply_2(Ys, X, Zs1),
  add_1([0|Zs0], Zs1, Zs2),
  multiply_1(Xs, Ys, Zs2, Zs).
  
multiply_2([], _, []).  
multiply_2([Y|Ys], X, [Z|Zs]):-
  Z is X * Y,
  multiply_2(Ys, X, Zs).
  
/* divide(Ns, Ds, Qs, Rs) is true if Qs and Rs are the quotient and        */
/*   remainder when dividing the polynomial Ns by the polynomial Ds.       */
/* e.g. divide([4,5,6,7], [1,2,4], [4,-3] ,[-4,19]).                       */
divide(_, [0|_], [0], [0]):-!, 
  write(`Error - attempted division by zero`), nl, fail.
divide(Ns, [D|Ds], Qs, Rs):-
  divide_1(Ns, D, Ds, Qs, Rs).

divide_1([N|Ns], D, Ds, [Q|Qs], Rs):-
  Q is N / D,
  multiply_2(Ds, Q, Ds1),
  subt(Ds1, Ns, Ns1), !,
  divide_1(Ns1, D, Ds, Qs, Rs).
divide_1(Ns, _, _, [], Ns). 
 
/* e.g. subt([8,16], [5,6,7], [-3,-10,7]). */
subt([], Ys, Ys).
subt([X|Xs], [Y|Ys], [Z|Zs]):-
  Z is Y - X,
  subt(Xs, Ys, Zs).

/* 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).

/* trim_left(Xs, Y, Zs) is true if Zs is the list Xs with all leading      */
/*   occurrences of Y removed.                                             */
trim_left([Y|Xs], Y, Zs):-!, trim_left(Xs, Y, Zs).
trim_left(Xs, _, Xs).

LPA Index     Home Page