Contents of the file c:\Maple\polynomials\polynomials.p

#=============================================================================
# Polynomial conversion functions
#=============================================================================
# Marc A. Murison
# Astronomical Applications Dept.
# U.S. Naval Observatory
# 3450 Massachusetts Ave., NW
# Washignton, DC 20392
# murisonATusno.navy.mil
# http://www.alpheratz.net/murison/
#=============================================================================


#---------------------------------------------------------------
# Convert a Taylor series in x to orthogonal polynomials in x.
#
# Usage: convert( p, ptype, x ), where
# p = a polynomial
# ptype = one of the polynomial types
# x = the independent variable (i.e., p = p(x))
#
# Polynomial types are:
# T Chebyshev type 1
# U Chebyshev type 2
# P Legendre
# H Hermite
# L Laguerre
#---------------------------------------------------------------
polynomials[`convert/T`] := proc( poly::{`+`,series}, x::name )
global T, __T;
T := 'T';
`polynomials/convert_terms`(poly,x,T,args[3..nargs]);
subs( __T=T, % );
end:

polynomials[`convert/U`] := proc( poly::{`+`,series}, x::name )
global U, __T;
U := 'U';
`polynomials/convert_terms`(poly,x,U,args[3..nargs]);
subs( __T=U, % );
end:

polynomials[`convert/P`] := proc( poly::{`+`,series}, x::name )
global P, __T;
P := 'P';
`polynomials/convert_terms`(poly,x,P,args[3..nargs]);
subs( __T=P, % );
end:

polynomials[`convert/H`] := proc( poly::{`+`,series}, x::name )
global H, __T;
H := 'H';
`polynomials/convert_terms`(poly,x,H,args[3..nargs]);
subs( __T=H, % );
end:

polynomials[`convert/L`] := proc( poly::{`+`,series}, x::name )
global L, __T;
L := 'L';
`polynomials/convert_terms`(poly,x,L,args[3..nargs]);
subs( __T=L, % );
end:

polynomials[`convert/G`] := proc( poly::{`+`,series}, x::name, a::algebraic )
global G, __T;
G := 'G';
`polynomials/convert_terms`(poly,x,G,a,args[4..nargs]);
subs( __T=G, % );
end:


#---------------------------------------------------------------
# Convert an orthogonal polynomial in x to a Taylor series in x.
#
# Usage: convert( p, taylor, x ), where
# p = a polynomial
# x = the independent variable (i.e., p = p(x))
#
# Polynomial types are:
# T Chebyshev type 1
# U Chebyshev type 2
# P Legendre
# H Hermite
# L Laguerre
#---------------------------------------------------------------
polynomials[`convert/taylor`] := proc( poly::`+`, x::name )
orthosubs( args );
collect( %, x );
end:



#---------------------------------------------------------------
# Substitute for T, U, P, H, L, or G, in the polynomial poly,
# the corresponding orthogonal polynomial terms. For example,
#
# > randpoly(x);
#
# 5 4 3 2
# x - 47 x - 91 x - 47 x - 61 x + 41
#
# > convert(%,P,x);
#
# 1222 239 1618 4031 376
# - ---- P[2](x) + --- - ---- P[3](x) - ---- P[1](x) - --- P[4](x)
# 21 15 45 35 35
#
# + 8/63 P[5](x)
#
# > orthosubs(%,x);
#
# 5 4 3 2
# x - 47 x - 91 x - 47 x - 61 x + 41
#
#---------------------------------------------------------------
polynomials[orthosubs] := proc( poly::`+`, x::name )

local N, k, p, pargs, ptype, xpart;

#detect the orthogonal polynomial type Q
xpart := select( has, poly, x );
if type(xpart,`+`) then #xpart = A*Pn(x) + B*Pm(x) + ...
op(1,select(has,poly,x)); #grab a term
op(select(has,[op(%)],x)); #isolate Q[n](x)
else #xpart = A*Pn(x)
select( has, xpart, x ); #isolate Q[n](x)
fi;
ptype := op([0,0],%); #select Q
debug_print(procname,cat("polynomial type = ",ptype),5);

if ptype='G' then
if nargs=2 then
pargs := 0,x;
else
pargs := args[3],x;
fi;
else
pargs := x;
fi;

#determine the polynomial order
N := 0;
for p in poly do
if has(p,ptype) then
N := max( N, op([0,1],op(select(has,[op(p)],x))) );
fi;
od;

#substitute for Q[k](x) the orthogonal polynomial term in x
subs( seq( ptype[k](x)=orthopoly[ptype](k,pargs), k=0..N ), poly );

end:


#---------------------------------------------------------------
# A plotting procedure for demonstration purposes.
#---------------------------------------------------------------
polynomials[plotpoly] := proc( intlist::list(posint), ptype::name )

local p, k, xmin, xmax, ymin, ymax, ytext, ttext, pargs;

if ptype='G' then
pargs := args[3], x;
else
pargs := x;
fi;

xmin := -1;
xmax := 1;
ymin := -1;
ymax := 1;
if ptype='L' then
xmin := 0;
xmax := 5;
ymin := -5;
ymax := 5;
elif ptype='H' then
xmin := 0;
xmax := 3;
ymin := -2;
ymax := 8;
elif ptype='U' then
ymin := -5;
ymax := 5;
elif ptype='G' and args[3] > 1/2 then
ymin := -5;
ymax := 5;
fi;

p := [];
for k from 1 to nops(intlist) do
if ptype='H' then
p := [ op(p),
plot( orthopoly[ptype](intlist[k],pargs)/intlist[k]^3,
x=xmin..xmax, color=blue ) ];
else
p := [ op(p),
plot( orthopoly[ptype](intlist[k],pargs),
x=xmin..xmax, color=blue ) ];
fi;
od;

ytext := cat(ptype,"[n](x)");
if ptype='P' then
ttext := "Legendre";
elif ptype='T' then
ttext := "Chebyshev type 1";
elif ptype='U' then
ttext := "Chebyshev type 2";
elif ptype='H' then
ttext := "Hermite";
ytext := cat(ptype,"[n](x)/n^3");
elif ptype='L' then
ttext := "Laguerre";
elif ptype='G' then
ttext := cat("Gegenbauer (a=",convert(args[3],string),")");
fi;

plots[display]( [op(p),plot(0,x=xmin..xmax,color=black)],
view=[xmin..xmax,ymin..ymax],
title=ttext, labels=["x",ytext] );

end:



#---------------------------------------------------------------
# Initialization function, where we define the workhorse
# routines used by the public polynomials package routines.
#---------------------------------------------------------------
polynomials[init] := proc()

global `polynomials/convert_terms`, `polynomials/get_horner_coeffs`;

#---------------------------------------------------------------
# The polynomial conversion workhorse routine.
# Not intended for public use.
#---------------------------------------------------------------
`polynomials/convert_terms` := proc( poly::{`+`,series}, x::name, ptype::name )

global __T, `polynomials/get_horner_coeffs`;
local n, #order of the polynomial
Q, #local storage of the polynomial
r, #structure to hold the recursion formulae for x*T[n]
r1, #structure to hold the recursion formulae for x (usually T[1])
k, c, p, j, tmp, m, S, a, pargs, astart, i, nterms;

if not type(poly,{polynom(anything,x),series}) then
ERROR(cat("not a polynomial or series in ",x));
fi;

if type(poly,series) then
S := convert(poly,polynom);
else
S := poly;
fi;

__T := '__T';
if ptype='G' then
if nargs=3 then
a := 0;
astart := 4;
else
a := args[4];
astart := 5;
fi;
if a=0 then
ERROR("Computation of G[n](x,0) is erroneous in the polynomials package.");
fi;
pargs := __T, a;
else
astart := 4;
pargs := __T;
fi;

#assign the recursion formulae for r=x*T[n](x) and r1=x for each polynomial type
r[T] := (P,n) -> (P[n+1] + P[n-1])/2;
r1[T] := (P) -> P[1];
r[U] := (P,n) -> (P[n+1] + P[n-1])/2;
r1[U] := (P) -> P[1]/2;
r[P] := (P,n) -> ((n+1)*P[n+1] + n*P[n-1])/(2*n+1);
r1[P] := (P) -> P[1];
r[H] := (P,n) -> P[n+1]/2 + n*P[n-1];
r1[H] := (P) -> P[1]/2;
r[L] := (P,n) -> (n+1)*(P[n] - P[n+1]) + n*(P[n] - P[n-1]);
r1[L] := (P) -> 1-P[1];
r[G] := proc(P,a,n)
((n+1)*P[n+1] + (n+2*a-1)*P[n-1])/(2*(n+a));
end;
r1[G] := proc(P,a)
if a=0 then
P[1]/2;
else
P[1]/(2*a);
fi;
end;

#get the horner-form coefficients
c := `polynomials/get_horner_coeffs`(S,x);
n := nops(c) - 1;
debug_print(procname,"horner coefficients",5,c);

#starting with the innermost horner nesting, replace x and T[n]*x with
#equivalents from the recursion formula for the polynomial type T
Q := 0;
for k from 1 to n do
if k=1 then
Q := c[n-k+1] + c[n-k+2]*r1[ptype](pargs);
else
Q := c[n-k+1] + Q*x;
Q := collect( Q, [seq(__T[j],j=0..k-1)] );
if type(Q,`+`) then
nterms := nops(Q);
else
nterms := 1;
fi;
tmp := Q;
for i from 1 to nterms do
if nterms=1 then
p := tmp;
else
p := op(i,tmp);
fi;
if has(p,__T) then
m := op(1,select(has,p,__T)); #order of T (i.e., T's subscript)
Q := subs( p=algsubs( __T[m]*x=r[ptype](pargs,m), p ), Q );
Q := subs(__T[0]=1,Q); #T[0]=1 for all orthogonal polynomials
else
Q := subs( p=subs(x=r1[ptype](pargs),p), Q );
fi;
od;
fi;
od;
Q := collect( Q, [seq(__T[j],j=0..n)], args[astart..nargs] );

#replace T[n] with T[n](x)
for p in Q do
m := op(1,select(has,p,__T));
Q := subs( p=subs(__T[m]=__T[m](x),p), Q );
od;

Q;

end:

#---------------------------------------------------------------
# Get the coefficients of the horner form of a polynomial in x.
#
# For example, 1 + 2*x - 3*x^2 becomes 1 + (2 - 3*x)*x, and the
# coefficients are [1,2,-3].
#
# Not intended for public use.
#---------------------------------------------------------------
`polynomials/get_horner_coeffs` := proc( poly::{`+`,series}, x::name )

local S, C, k, n, c, j, m, xpart;

if not type(poly,{polynom(anything,x),series}) then
ERROR(cat("not a polynomial or series in ",x));
fi;

if type(poly,series) then
S := convert(poly,polynom);
else
S := poly;
fi;

S := convert( S, horner, x );
c := [];
n := degree( collect(S,x), x );
for k from 1 to n do
if type(S,`+`) then #S = a + (b + ...)*x
c := [ op(c), remove(has,S,x) ];
else #S = (a + ...)*x
c := [ op(c), 0 ];
fi;
xpart := select(has,S,x);
if xpart=x then
S := 1;
else
if type(xpart,`^`) or nops(xpart) > 2 then #xpart=a*b*...*x^m
m := degree( xpart, x );
else #xpart=(a + ...)*x^m
m := degree( op(2,xpart), x );
fi;
if m > 1 then
for j from 2 to m do
c := [ op(c), 0 ];
xpart := xpart/x;
k := k + 1;
od;
fi;
S := xpart/x;
fi;
od;
[ op(c), S ];

end:

end:


save( polynomials, cat(POLY_PATH,"polynomials.m") );