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