/* Chebyshev series: a small set of programs to get you started */
ToChebyt(e,x,n):= 
/* Convert analytic expression, depending on x, to Chebyshev series,
  via Taylor series  */

  block([ tay: taylor(e,x,0,n),q,r,ans:[]],
  power2cheby(makelist(ratcoef(tay,x,i),i,0,n)))$


/* Convert analytic function lambda([x]...expression in x) to Chebyshev series
  The expression f may be a polynomial or function in other variables, too.
  The result will be a list of the Chebyshev coefficients.
  Here we try to make them all floats.
  We are going to do this by a Discrete Cosine Transform, (DCT) and also
  by a Fast DCT.
   */


/* if u is a list of n values of a function f at selected points ... 
{ f(cos( pi*(k+1/2)/n)), k=0.. n-1 }
s1: DCT(u) is a list of the Chebyshev series coefficients. */

DCT(u):= 
 /* numeric (slow version) of discrete cosine transform on a list u. */
 block([n:length(u), L:[],q,in,nm1],
  q:?pi/n, in:sqrt(1.0d0/n), nm1:n-1,
  for i: 0 thru nm1 do
  L: cons(sum(u[r]*cos(q*(r-1/2)*(nm1-i)),r, 1, n)*in,L),
  L)$

( /*this fastDCT works if u is a list of numbers of even length.
  It uses whatever FFT is in maxima, for good or ill.*/

load(fft)$  /* needs this file loaded */

drv(m):=append(space2(m),space2(reverse(m))),
space2(k):= if k=[] then [] else cons(0,cons( first(k),space2(rest(k)))),

fastDCT(u):=
  block([n:length(u),a:fft(drv(u))],
   2*sqrt(float(n))* makelist(a[i],i,1,n)),

/* gauss-lobatto points*/
glp(n):= makelist(cos(?pi*(i+0.5d0)/n),i,0,n-1)$
glpf(%%f,n):= map(%%f, glp(n)))$

ToCheby2(f,n):=  /* n need not be a power of 2, not especially fast */
  block([dc:DCT(glpf(f,n)), in:sqrt(1.0d0/n)],
 apply('chebseries,makelist(dc[i]*2*in*(2- (if(j=n) then 0 else 1)),i,1,n)))$

ToCheby2f(f,n):=  /* n MUST be a power of 2, should be faster */
  block([dc:fastDCT(glpf(f,n)), in:sqrt(1.0d0/n)],
 apply('chebseries,makelist(dc[i]*2*in*(2- (if(j=n) then 0 else 1)),i,1,n)))$

/* if you don't like ToCheby2(lambda([x],x^2*sin(x)),8), you can do
   TC(x^2*sin(x),8, where TC is ...*/
TC(%%f,%%x,n):=ToCheby2s(buildq([f:%%f,x:%%x],lambda([x],f)),n)$



FromCheby(a,x):=  
/* If x is a symbol, this will convert Chebyshev series to an expression
   in the power basis, i.e. polynomial in x.
   If x is a number, this will evaluate the series at that point and
   provide a number. (if a is just numbers).
   This algorithm is based on recurrence in Clenshaw alg. */

  block([N:length(a)-1,keepfloat:true,
         bj2:0,bj1:0,bj0:0, z:reverse(a)],
 for j: 1 thru N do
  (bj0:2*x*bj1-bj2+first(z), 
  bj2:bj1, bj1:bj0,z:rest(z)),
  x*bj1+first(z)/2 -bj2)$
/* a test */

small_to_zero(expr,eps):=
 /* Replaces all numbers smaller in absolute value than EPS by zero . */
   fullmap(lambda([x], if numberp(x) and abs(x)<eps then 0 else x), 
             expr)$
stz(z):=small_to_zero(z,1.0d-9)$

	

/* Given a Chebyshev series, return another one in which we may have removed
 terms that don't matter much; that is, trailing coeffs were small relative to
large coeffs. */

TrimCS(a,relerr):=
 block([max: apply(max,args(a)),
        min: apply(min,args(a)),
        r: reverse(a),
        bigmaginv],
 bigmaginv: 1/max(abs(max),abs(min)),
 while (abs(first(r))*bigmaginv < relerr) do r:rest(r),
 return(reverse(r)))$


/* EXTRAS: Some simple transformations,  examples */


stz(FromCheby(ToCheby2(lambda([x],a*x^5+b*x^2+45),6),rat(x)));
FromCheby(ToChebyt(a*x^5+b*x^2+45,x,6),rat(x));


kill(tc,tc2)$
tc[0](x):=1$
tc[1](x):=rat(x)$
tc[n](x):=2*x*tc[n-1](x)-tc[n-2](x)$

/*here's another way to generate cheby polys, if you are in a hurry */
tc2[0](x):=1$
tc2[1](x):=x$
tc2[n](x):= if evenp(n) then 2*tc2[n/2](x)^2-1
  else 2*tc2[(n+1)/2](x)*tc2[(n-1)/2](x)-x$


/* Multiplying chebyshev series. works*/
ChebTimes(a,b) := 
/* assuming matching variable, e.g. x, and range e.g. [-1,1] */
block([la:length(a),lb:length(b),ans,bb,h],
local(ans),
array(ans, la+lb-1),
for i: 0 thru la+lb-1 do ans[i]:0,
 for i:0 thru la-1 do
 (bb:b,
  for j:0 thru lb-1 do
    ( term: first(a)*first(bb)/2,
       ans[i+j]: ans[i+j]+term,
       h:abs(i-j),
       ans[h]:ans[h]+term,
       bb:rest(bb)),
     a:rest(a)),
  ans[0]:2*ans[0],
 apply('chebseries, makelist(ans[i],i,0,la+lb-1)))$


/* a simple formula (Thacher, 1964) sufficient to convert from x^r to Cheby.
    onepowcs[r] is the chebseries for x^r */
  
(kill(onepowcs), onepowcs[r]:=apply('chebseries,makelist(
   if evenp(r-q)then 2^(1-r)*binomial(r,(r-q)/2) else 0,
                     q,0,r)))$

/*if powers(a,b,c,d) is an encoding of a+b*x+c*x^2+d*x^3,   with "x" understood, then converting
from power basis [exactly!] to chebshev basis is */

power2cheby(p):=
 block([pa,ans,n:length(p)],
  local(pa,ans),
  pa[i]:=inpart(p,i+1),
  ans[i]:=0,
  for r:0 thru n-1 do /* for each power r*/
   for q:0 thru r do (
     if evenp(r-q)then 
           ans[q]:ans[q]+pa[r]*2^(1-r)*binomial(r,(r-q)/2)),
  apply ('chebseries, makelist(ans[i],i,0,n-1)))$

cheby2power(c):=  /* convert to power basis list, not so efficient perhaps */
block([f:FromCheby(c,rat(zz))],
apply('powers,makelist(ratcoeff(f,zz,i),i,0,length(c)-1)))$



(matchdeclare(atrue,true),
defrule(r1,x^atrue,onepowcs[atrue]),
s:a+b*x+c*x^3,
TTT:apply1(s,r1));  /* note that x is treated as x^1 */

(chebseriesp(x):= is (not(atom(x))and inpart(x,0)=chebseries),
matchdeclare(cs,chebseriesp),
defrule(r2,cs,FromCheby(cs,xx)),
AllFromCheby(s,xx):=apply1(s,r2),
AllFromCheby(TTT,rat(z)))$


iDCT(u):=  
/* (slow version) of Inverse discrete cosine transform on a list u. */
  block([n:length(u), L:[],q,nh,ni],
   q:?pi/n,ni:sqrt(1.0d0/n),nh:n-0.5d0,
   for i: 0 thru n-1 do
    L: cons((u[1]+2*sum(u[r]*cos(q*(r-1)*(nh-i)),r, 2, n))*ni,L),
   L)$