/* Chebyshev related functions.  author: Richard Fateman. 2016
See also cheby.lisp for faster versions of some of this */
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 */
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)))$


/* another version, which has been written elswhere in lisp for maxima */

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

ToCheby2(f,n):=  /* this one takes a function, e,g, lambda([x],..) */
  apply('chebseries,makelist(AJS1(f,j,n),j,0,n-1))$

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

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

AJS1x(%%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), print([h,cos(h),%%f(cos(h))]),
       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
   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)$

FromChebyfloat(a,x):=  
/* x is a floating-point number.
   This algorithm is based on recurrence in Clenshaw alg. */

  block([N:length(a)-1,keepfloat:true,fz:0.0,
         bj2:0.0,bj1:0.0,bj0:0.0, z:reverse(a)],
modedeclare([bj0,bj1,bj2,x,fz],float,[j,N], fixnum),
 for j: 1 thru N do
  (fz:first(z),
 bj0:2.0*x*bj1-bj2+fz, 
  bj2:bj1, bj1:bj0,z:rest(z)),
  x*bj1+first(z)*0.5 -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)$


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

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


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

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

/* copy from wolfram documentation */
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.


*/


   
/* Thoughts about accurate evaluation of a polynomial.
 if you don't want to evaluate on -1,1, but f(x) on [a,b],
then translate the center (a+b)/2  to zero , a to -1, b to 1.*/

mapto_basep(f%,a,b,x):= subst(1/2*( (b-a)*x+b+a),x,f%)$ /* f is polynomial expression */

mapfrom_base(g%,a,b,x):=g%((2*y-b-a)/(b-a))$


/*to create a function pf
to  evaluate poly expression p, in x, accurately between a and b: */

accuratize(p,x,a,b):=buildq(
            [p,x,q%:x-float((a+b)/2),
             body:float(horner(ratsubst((a+b)/2+%y,x,p)))],
            block([%y:q%],body))$

/*usage */

/* define(pf(x),accuratize(x^4-1,x,3,4));
results in
pf(x):=block([%y:x-3.5],
    0.0625*(%y*(%y*(%y*(16.0*%y+224.0)+1176.0)+2744.0)+2385.0))


*/

/*For Chebyshev approx, we can put everything in one function also. */

accuratizeC(p,x,a,b):=
 (p:ratexpand(p),
buildq([p,x,
       q%:  float((2*x-b-a)/(b-a)),
       body:float(ToChebyt(mapto_basep(p,a,b,x),x,hipow(p,x)))],
         FromCheby(body, q%)))$

accuratizeCchop(p,x,a,b,deg):=
 (p:ratexpand(p),
buildq([p,x,
       q%:  float((2*x-b-a)/(b-a)),
       body:float(ToChebyt(mapto_basep(p,a,b,x),x,deg))],
         FromCheby(body, q%)))$

/*usage 

define(pf(x),accuratizeC(x^4-1,x,3,4));  
results in 

(x):=x^4-1 for example.  note. load("cheby.lisp"); load("cheby.mac");

pf(x):=FromCheby(chebseries(316.546875,87.0625,9.21875,0.4375,0.0078125),
   2.0*x-7.0)$
 
more stuff in fateman's c:/lisp/generic/cheby.lisp
*/