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

/* evaluate(Polynomial, X, Y) is true if Y is the value of the Polynomial    */
/*   evaluated for the value X.                                              */
/* e.g. evaluate([2,3,4,5], 2, 41).                                          */
evaluate([A|As], X, Y):-
  evaluate_1(As, X, A, Y).
  
evaluate_1([], _, Y, Y).
evaluate_1([A|As], X, Y0, Y):-
  Y1 is Y0 * X + A,
  evaluate_1(As, X, Y1, Y).

/* p_add(Xs, Ys, Zs) is true if adding the polynomial Xs to the polynomial   */
/*   Ys gives the polynomial Zs.                                             */
/* e.g. p_add([1,2,3], [4,5,6,7], [4,6,8,10]).                               */
p_add(Xs, Ys, Zs):-
  reverse(Xs, ReversedXs),
  reverse(Ys, ReversedYs),
  p_add_1(ReversedXs, ReversedYs, ReversedZs),
  reverse(ReversedZs, Zs0),
  trimleft(Zs0, Zs).

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

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

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

/* p_mult(Xs, Ys, Zs) is true if multiplying the polynomial Xs by the        */
/*   polynomial Ys gives the polynomial Zs.                                  */
/* e.g. p_mult([1,2,3], [4,5,6,7], [4,13,28,34,32,21]).                      */
p_mult(Xs, Ys, Zs):-
  reverse(Ys, ReversedYs),
  p_mult_1(Xs, ReversedYs, [], ReversedZs),
  reverse(ReversedZs, Zs).
  
p_mult_1([], _, Zs, Zs).  
p_mult_1([X|Xs], Ys, Zs0, Zs):-
  p_mult_2(Ys, X, Zs1),
  p_add_1([0|Zs0], Zs1, Zs2),
  p_mult_1(Xs, Ys, Zs2, Zs).
  
p_mult_2([], _, []):-!.  
p_mult_2([Y|Ys], X, [Z|Zs]):-
  Z is X * Y,
  p_mult_2(Ys, X, Zs).

/* p_div(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. p_div([10,27,45,32], [2,3,4], [5,6], [7,8]).                         */
p_div(_, [0|_], [0], [0]):-!, 
  write(`Error - attempted division by zero`), nl, fail.
p_div(Ns, [D|Ds], Qs, Rs):-
  p_div_1(Ns, D, Ds, Qs, Rs0),
  trimleft(Rs0, Rs).

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

/* pseudiv(Us, Vs, Qs, Rs) is true if the polynomials Us, Vs, Qs and Rs      */
/*   with integer coefficients are such that V^(M-N+1)Us=QsVs+Rs where V is  */
/*   the leading coefficient of Vs, M=deg(Us), N=deg(Vs) and M>=N>=0.        */
/*   [Knuth pseudo-division algorithm 4.6.1R]                                */
/* e.g. pseudiv([1,1,-1,2,3,-1,2], [2,2,-1,3], [8,0,-4,8], [28,4,8]).        */
/*      pseudiv([1,0,1,1], [3,1,1], [3,-1], [7,10]).                         */
pseudiv(Us, [V|Vs], Qs, Rs):- % deg(Vs)>=0
  V \= 0,                     % lc(Vs)!=0
  length(Us, Lu),             % deg(Us)+1
  length(Vs, Lv),             % deg(Vs)
  K is Lu - Lv - 1,           % deg(Us)-deg(Vs)
  K >= 0,                     % deg(Us)>=deg(Vs)
  power(V, K, Vk),            % V^K
  pseudiv_1(K, Vk, Us, V, Vs, Qs, Rs).

pseudiv_1(-1, _, Us, _, _, [], Rs):-!,
  trimleft(Us, Rs).
pseudiv_1(K, Vk, [U|Us], V, Vs, [Q|Qs], Rs):-
  Q is U * Vk,                % lc(Us)lc(Vs)^K
  p_mult_2(Us, V, VUs),       % lc(Vs)Us
  p_mult_2(Vs, U, UVs),       % lc(Us)Vs
  p_subt_1(VUs, UVs, Us1),    % lc(Vs)Us-lc(Us)Vs
  Vk1 is Vk // V,             % V^(K-1)
  K1 is K - 1,
  pseudiv_1(K1, Vk1, Us1, V, Vs, Qs, Rs).

/* power(X, N, Y) is true if Y is the result of raising X to the power N.    */
power(X, N, Y):-N >= 0, power_1(N, X, 1, Y).

power_1(0, _, Y, Y):-!.
power_1(1, X, Y0, Y):-!,
  Y is Y0 * X.
power_1(N, X, Y0, Y):-
  1 =:= N mod 2, !,
  N1 is N // 2,
  X1 is X * X,
  Y1 is Y0 * X,
  power_1(N1, X1, Y1, Y).
power_1(N, X, Y0, Y):-
  N1 is N // 2,
  X1 is X * X,
  power_1(N1, X1, Y0, Y).

/* trimleft(Xs, Zs) is true if Zs is the list Xs with all leading zeroes     */
/*   removed, except that [] is replaced by [0].                             */
trimleft([], [0]):-!.
trimleft(Xs, Zs):-
  trimleft_1(Xs, Zs).

trimleft_1([0|Xs], Zs):-
  Xs=[_|_], !,
  trimleft_1(Xs, Zs).
trimleft_1(Xs, Xs).

% Display all the polynomials in the list
% e.g. show_polys([[1,-1],[1,1,1],[1,1,1,1,1],[1,-1,0,1,-1,1,0,-1,1]]) displays
%   x-1
%   x^2+x+1
%   x^4+x^3+x^2+x+1
%   x^8-x^7+x^5-x^4+x^3-x+1
show_polys([]).
show_polys([P|Ps]):-
  show_poly(P),
  show_polys(Ps).

% Display one polynomial
show_poly([V]):-!,
  show_sign_1(V), AbsV is abs(V), write(AbsV), nl.
show_poly([V|Vs]):-
  length(Vs, L), 
  show_sign_1(V), 
  show_value(V), 
  show_exponent(L), 
  L1 is L - 1, 
  show_poly(Vs, L1),
  nl.
  
show_poly([], _):-!.
show_poly([0], _):-!.
show_poly([V], _):-!, 
  show_sign(V), AbsV is abs(V), write(AbsV).
show_poly([0|Vs], L):-!,
  L1 is L - 1, show_poly(Vs, L1).
show_poly([V|Vs], L):-
  show_sign(V), show_value(V), show_exponent(L), L1 is L-1, show_poly(Vs, L1).

show_sign_1(V):-V >= 0, !.  
show_sign_1(_):-write(`-`).

show_sign(0):-!.  
show_sign(V):-V > 0, write(`+`), !.  
show_sign(_):-write(`-`).
  
show_value(1):-!.
show_value(-1):-!.
show_value(V):-AbsV is abs(V), write(AbsV).
  
show_exponent(0):-!.
show_exponent(1):-!, write(`x`).
show_exponent(L):-write(`x^`), write(L).  
  
/* 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(Xs, L0, L) is true if L is equal to L0 plus the number of        */
/*   elements in the list Xs.                                                */
%length_1([], L, L).
%length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).

/* reverse(Xs, Ys) is true if Ys is the result of reversing the order of the */
/*   elements in the list Xs.                                                */
%reverse(Xs, Ys):-reverse_1(Xs, [], Ys).

%reverse_1([], As, As).
%reverse_1([X|Xs], As, Ys):-reverse_1(Xs, [X|As], Ys).

LPA Index     Home Page

Valid HTML 4.01 Strict