kill(tc,onepower)$
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
 to get a few large ones, and not all the intermediate ones */

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$

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)))$

/* ToCheby2 is similar, but this one converts an analytic function
  lambda([x]...expression in x) to Chebyshev series by evaluation.
  The expression f may be a polynomial or function in other variables, too.
   */

ToCheby2(f,n):=  /* this one takes a function, e,g, lambda([x],..) */
 block([in:1.0d0/n,p:?pi/n],
  apply('chebseries,makelist(AJS1x(f,j,n,in,p,0),j,0,n-1)))$

AJS1(%%f,j,n,in,p,h):=  /* compute a[j]s, by eval at gauss-lobatto points */
  2*in*(2-(if (j=n) then 0 else 1))*
  sum( (h:p*(k+.5d0),
       cos(h*j)*%%f(cos(h))),k,0,n-1)$


FromCheby(a,x):=  
/* If x is a symbol, this will convert Chebyshev series to an expression
   explicitly written out 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. (assuming typical situation: 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)$

/* We end up generating series that have very small terms in them,
like 1e-16, which are painful to look at so we can remove them with
this program: */

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)$

/* now we are ready for a test  or two*/

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

globrelerr:1.0d-15$
tcs(a):=TrimCS(a,globrelerr)$
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)))$

/* Multiplying chebyshev series. works. We could do much more of this,
  and also provide a kind of polymorphic arithmetic or generic arithmetic
  surface to this.*/
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)) */
 tcs( 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)))$


/* simple transformation examples */

(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)))$


/* 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):=  /* (slow version) of discrete cosine transform on a list u. */
 block([n:length(u), L:[],q],
  q:?pi/n,
for s from 1 thru n do
L: cons(sum(u[r]*cos(q*(r-0.5d0)*(s-1)),r, 1, n),L),reverse(sqrt(float(1/n))*L)),

DCTs(u):=  /* symbolic(slow version) of discrete cosine transform on a list u. */
 block([n:length(u), L:[]],
for i: 0 thru n-1 do
L: cons(sum(u[r]*cos(%pi*(r-1/2)*(n-1-i)/n),r, 1, n)/sqrt(n),L),L)$

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

iDCTs(u):=  /* symbolic (slow version) of Inverse discrete cosine transform on a list u. */
 block([n:length(u), L:[]],
for i: 0 thru n-1 do
L: cons((u[1]+2*sum(u[r]*cos(%pi*(r-1)*(n-1/2-i)/n),r, 2, n))/sqrt(n),L),L)$


/* [slow] fourier transform, just to check up on FFT values :)*/
DFT (a):=block([A,n:length(a),B,h,d],
local(A,B),d:1.0d0/n,
 A[i]:=a[i+1],
 B[i]:=0,
 for j:0 thru n-1 do (h: exp(2*j*%i*?pi/n),
           B[j]:d*sum(A[k]*h^k,k,0,n-1)),
 makelist(B[j],j,0,n-1))$	

(DFTe (a):=block([A,n:length(a),B,h,d],
 /*exact fourier, leave pi in, uses cos+i*sin instead of exp, for fun.*/
d:1/n,
local(A,B),
 A[i]:=a[i+1],
 Y[i]:=0,
 for k:0 thru n-1 do (h: cos(2*k*%pi/n)+%i*sin(2*k*%pi/n),
           Y[k]:d*sum(A[j]*h^j,j,0,n-1)),
makelist(Y[i],i,0,n-1)),	

DCTe(x):=block([X,Y,n:length(x),d], /*exact cosine, any length*/
 local(X,Y),
 d:sqrt(1/n),
 X[i]:=x[i+1],
 Y[i]:=0,
 for k:0 thru n-1 do Y[k]: d*sum(X[j]*cos(pi*(j+1/2)*k/n), j,0,n-1),
 makelist(Y[i],i,0,n-1)))$

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

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

 fastDCT(x):=
  block([n:length(x),a:fft(drv(x))],
   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))$

ToCheby2s(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)))$

/* how fast is this?  let g(x):= cos(30*x).
  computing 128 terms in the series takes
   .32 sec in compiled lisp, n^2 version
   .61 sec in Maxima code,   n^2 version
   .04 sec in Maxima code calling the FFT program.

We can now write a program that guesses at the number of terms necessary to encode a given function f.  If this is judged to be not enough, we can double the
number of terms.  If this is deemed to be too many, we can trim off the extra.



*/