%       sfun.s
%
%       Purpose:        Special Mathematical Functions
%                       1. Orthogonal Polynomials
%       Copyright:      Ilja Heitlager, Leiden University, 1997
%       Updates:
%	Literature: H. Bateman, Higher Transcendental Functions, Vol II,
%	            McGraw-Hill, 1953, A. Erdelyi (ed.) e.a.,
%	            Bateman manuscript project.

PolySimpStrategy := {[]}.

% Gegenbauer Polynomials.
%
%   L
%  C (X) = cln(X, L, N)
%   n
%
% Example: cln(x, 3, 7)
%          x * (x ^ 2 * (x ^ 2 * (4608 * x ^ 2 - 5376) + 1680) - 120)

cln(1, 0, 0) := 1.
cln(1, 0, N) := 2 / N.
cln(- X, L, N) := -1 ^ N * cln(X, L, N).
(cln(X, L, N) :=
 (iff is_integer(N) and N > 0 then twocln(X, L, N) : 1 else
  error("cln: integer argument expected")
 ))
.
	
twocln(X, L, 1) := [2 * L * X, 1].
(twocln(X, L, N) :=
 ([ (( 2 * (N + L - 1) / N * X * cln : 1
     - (N + 2 * L - 2) / N * cln : 2
     ) with
     PolySimpStrategy
    )
  , cln : 1 ] where
  cln = twocln(X, L, N - 1)
 ))
.

cln(X, L, N) latex "C^", L, "_", N, "(", X, ")".

% Legendre Polynomials.
%
% First Kind:
%
%  P (X) = pn(X,N)
%   n
%
% Second Kind:
%
%  Q (X) = qn(X,N)
%   n

pn(X, N) := cln(X, 1 / 2, N).

pn(X, N) latex "P_", N, "(", X, ")".

qn(X, 0) := 1 / 2 * log((X + 1) / (X - 1)).
(qn(X, N) :=
 (iff is_integer(N) and N > 0 then twoqn(X, N) : 1 else
  error("qn: integer argument expected")
 ))
.

qn(X, N) latex "Q_", N, "(", X, ")".

(twoqn(X, 1) :=
 [qn * X - 1, qn] where qn = 1 / 2 * log((X + 1) / (X - 1)))
.
(twoqn(X, N) :=
 ([ (( (2 * N - 1) / N * X * qn : 1
     - (N - 1) / N * qn : 2
     ) with
     PolySimpStrategy
    )
  , qn : 1 ] where
  qn = twoqn(X, N - 1)
 ))
.

% Chebychev Polynomials.
%
% First Kind:
%
%  T (X) => tn(X,N)
%   n
%
% Second Kind:
%
%  U (X) => un(X,N)
%   n

(tn(X, N) :=
 (iff is_integer(N) and N > 0 then twotn(X, N) : 1 else
  error("tn: integer argument expected")
 ))
.

tn(X, N) latex "T_", N, "(", X, ")".

un(X, N) := cln(X, 1, N).

un(X, N) latex "U_", N, "(", X, ")".

twotn(X, 1) := [X, 1].
(twotn(X, N) :=
 ([2 * X * tn : 1 - tn : 2 with PolySimpStrategy, tn : 1] where
  tn = twotn(X, N - 1)
 ))
.