1module charpol;
2
3% Author: R. Liska.
4
5% Redistribution and use in source and binary forms, with or without
6% modification, are permitted provided that the following conditions are met:
7%
8%    * Redistributions of source code must retain the relevant copyright
9%      notice, this list of conditions and the following disclaimer.
10%    * Redistributions in binary form must reproduce the above copyright
11%      notice, this list of conditions and the following disclaimer in the
12%      documentation and/or other materials provided with the distribution.
13%
14% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
15% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
16% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
17% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
18% CONTRIBUTORS
19% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
20% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
21% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
22% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
23% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
24% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
25% POSSIBILITY OF SUCH DAMAGE.
26%
27
28
29% Version REDUCE 3.6     05/1991.
30
31fluid '(!*exp !*gcd !*prfourmat)$
32
33switch prfourmat$
34!*prfourmat:=t$
35
36procedure coefc1 uu$
37begin
38  scalar lco,l,u,v,a$
39  u:=car uu$
40  v:=cadr uu$
41  a:=caddr uu$
42  lco:=aeval list('coeff,u,v)$
43  lco:=cdr lco$
44  l:=length lco - 1$
45  for i:=0:l do
46    <<setel(list(a,i),car lco)$
47      lco:=cdr lco >>$
48  return (l . 1)
49end$
50
51deflist('((coefc1 coefc1)),'simpfn)$
52
53global '(cursym!* coords!* icoords!* unvars!*)$
54
55icoords!*:='(i j k l m n i1 j1 k1 l1 m1 n1)$
56
57flag('(tcon unit charmat ampmat denotemat),'matflg)$
58
59put('unit,'rtypefn,'getrtypecar)$
60put('charmat,'rtypefn,'getrtypecar)$
61put('ampmat,'rtypefn,'getrtypecar)$
62put('denotemat,'rtypefn,'getrtypecar)$
63
64procedure unit u$
65generateident length matsm u$
66
67procedure charmat u$
68matsm list('difference,list('times,'lam,list('unit,u)),u)$
69
70procedure charpol u$
71begin
72  scalar x,complexx;
73  complexx:=!*complex;
74  algebraic on complex;
75  x:=simp list('det,list('charmat,carx(u,'charpol)))$
76  if null complexx then algebraic off complex;
77  return x
78end;
79
80put('charpol,'simpfn,'charpol)$
81
82algebraic$
83korder lam$
84
85procedure re(x)$
86sub(i=0,x)$
87
88procedure im(x)$
89(x-re(x))/i$
90
91procedure con(x)$
92sub(i=-i,x)$
93
94procedure complexpol x$
95begin
96  scalar y$
97  y:=troot1 x$
98  return if im y=0 then y
99           else y*con y
100end$
101
102procedure troot1 x$
103begin
104  scalar y$
105  y:=x$
106  while not(sub(lam=0,y)=y) and sub(lam=1,y)=0 do y:=y/(lam-1)$
107  x:=sub(lam=1,y)$
108  if not(numberp y or (numberp num y and numberp den y)) then
109      write " If  ",re x,"  =  0  and  ",im x,
110            "  =  0  , a root of the polynomial is equal to 1."$
111  return y
112end$
113
114procedure hurw(x)$
115% X is a polynomial in LAM, all its roots are |LAMI|<1 <=> for all roots
116% of the polynomial HURW(X) holds RE(LAMI)<0.
117(lam-1)**deg(num x,lam)*sub(lam=(lam+1)/(lam-1),x)$
118
119symbolic$
120
121procedure unfunc u$
122<<unvars!*:=u$
123  for each a in u do put(a,'simpfn,'simpiden) >>$
124
125put('unfunc,'stat,'rlis)$
126
127global '(denotation!* denotid!*)$
128denotation!*:=nil$
129denotid!*:='a$
130
131procedure denotid u$
132<<denotid!*:=car u$nil>>$
133
134put('denotid,'stat,'rlis)$
135
136procedure cleardenot$
137denotation!*:=nil$
138
139put('cleardenot,'stat,'endstat)$
140flag('(cleardenot),'eval)$
141
142algebraic$
143array cofpol!*(20)$
144
145procedure denotepol u$
146begin
147  scalar nco,dco$
148  dco:=den u$
149  u:=num u$
150  nco:=coefc1 (u,lam,cofpol!*)$
151  for j:=0:nco do cofpol!*(j):=cofpol!*(j)/dco$
152  denotear nco$
153  u:=for j:=0:nco sum lam**j*cofpol!*(j)$
154  return u
155end$
156
157symbolic$
158
159put('denotear,'simpfn,'denotear)$
160
161procedure denotear u$
162begin
163  scalar nco,x$
164  nco:=car u$
165  for i:=0:nco do
166    <<x:=list('cofpol!*,i)$
167      setel(x,mk!*sq denote(getel x,0,i)) >>$
168  return (nil .1)
169end$
170
171procedure denotemat u$
172begin
173  scalar i,j,x$
174  i:=0$
175  x:=for each a in  matsm u collect
176       <<i:=i+1$
177         j:=0$
178         for each b in a collect
179           <<j:=j+1$
180             denote(mk!*sq b,i,j) >> >>$
181  return x
182end$
183
184procedure denote(u,i,j)$
185% U is prefix form, I,J are integers
186begin
187  scalar reu,imu,ireu,iimu,eij,fgcd$
188  if atom u then return simp u$
189  fgcd:=!*gcd$
190  !*gcd:=t$
191  reu:=simp!* list('re,u)$
192  imu:=simp!* list('im,u)$
193  !*gcd:=fgcd$
194  eij:=append(explode i,explode j)$
195  ireu:=intern compress append(append(explode denotid!* ,'(r)),eij)$
196  iimu:=intern compress append(append(explode denotid!* ,'(i)),eij)$
197  if car reu then insdenot(ireu,reu)$
198  if car imu then insdenot(iimu,imu)$
199  return simp list('plus,
200                   if car reu then ireu
201                     else 0,
202                   list('times,
203                        'i,
204                        if car imu then iimu
205                          else 0))
206end$
207
208procedure insdenot(iden,u)$
209denotation!*:=(u . iden) . denotation!*$
210
211procedure prdenot$
212for each a in reverse denotation!* do
213  assgnpri(list('!*sq,car a,t),list cdr a,'only)$
214
215put('prdenot,'stat,'endstat)$
216flag('(prdenot),'eval)$
217
218procedure ampmat u$
219begin
220  scalar x,i,h1,h0,un,rh1,rh0,ru,ph1,ph0,!*exp,!*gcd,complexx$
221  complexx:=!*complex;
222  !*exp:=t$
223  fouriersubs()$
224  u:=car matsm u$
225  x:=for each a in coords!* collect if a='t then 0
226                                      else list('times,
227                                                tcar get(a,'index),
228                                                get(a,'wave),
229                                                get(a,'step))$
230  x:=list('expp,'plus . x)$
231  x:=simp x$
232  u:=for each a in u collect resimp quotsq(a,x)$
233  gonsubs()$
234  algebraic on complex;
235  u:=for each a in u collect resimp a$
236  remfourier()$
237a:if null u then go to d$
238  ru:=caar u$
239  un:=unvars!*$
240  i:=1$
241b:if un then go to c$
242  rh1:=reverse rh1$
243  rh0:=reverse rh0$
244  h1:=rh1 . h1$
245  h0:=rh0 . h0$
246  rh0:=rh1:=nil$
247  u:=cdr u$
248  go to a$
249c:rh1:=coefck(ru,list('u1!*,i)) . rh1$
250  rh0:=negsq coefck(ru,list('u0!*,i)) . rh0$
251  un:=cdr un$
252  i:=i+1$
253  go to b$
254d:h1:=reverse h1$
255  h0:=reverse h0$
256  if !*prfourmat then
257      <<ph1:=for each a in h1 collect
258               for each b in a collect prepsq!* b$
259        setmatpri('h1,nil . ph1)$
260        ph1:=nil$
261        ph0:=for each a in h0 collect
262               for each b in a collect prepsq!* b$
263        setmatpri('h0,nil . ph0)$
264        ph0:=nil >>$
265  !*gcd:=t;
266  x:=if length h1=1 then list list quotsq(caar h0,caar h1)
267       else lnrsolve(h1,h0)$
268  if null complexx then algebraic off complex;
269  return x
270end$
271
272procedure coefck(x,y)$
273% X is standard form, Y is prefix form, returns coefficient of Y
274% appearing in X, i.e. X contains COEFCK(X,Y)*Y
275begin
276  scalar ky,xs$
277  ky:=!*a2k y$
278  xs:=car subf(x,list(ky . 0))$
279  xs:=addf(x,negf xs)$
280  if null xs then return(nil . 1)$
281  xs:=quotf1(xs,!*k2f ky)$
282  return if null xs then <<msgpri
283    ("COEFCK failed on ",list('sq!*,x . 1,nil)," with ",y,'hold)$
284                           xread nil$
285                           !*f2q xs>>
286           else !*f2q xs
287end$
288
289procedure simpfour u$
290begin
291  scalar nrunv,x,ex,arg,mv,cor,incr,lcor$
292  nrunv:=get(car u,'nrunvar)$
293a:u:=cdr u$
294  if null u then go to r$
295  arg:=simp car u$
296  mv:=mvar car arg$
297  if not atom mv or not numberp cdr arg then return msgpri
298      ("Bad index ",car u,nil,nil,'hold)$
299  cor:=tcar get(mv,'coord)$
300  if not(cor member coords!*) then return msgpri
301      ("Term ",car u," contains non-coordinate ",mv,'hold)$
302  if cor member lcor then return msgpri
303      ("Term ",car u," means second appearance of coordinate ",cor,
304       'hold)$
305  if not(cor='t) and cdr arg>get(cor,'maxden)
306    then put(cor,'maxden,cdr arg)$
307  lcor:=cor . lcor$
308  incr:=addsq(arg,negsq !*k2q mv)$
309  if not flagp(cor,'uniform) then return lprie
310      ("Non-uniform grids not yet supported")$
311  if cor='t then go to ti$
312  ex:=list('times,car u,get(cor,'step),get(cor,'wave)) . ex$
313  go to a$
314ti:if null car incr then x:=list('u0!*,nrunv)
315    else if incr= 1 . 1 then x:=list('u1!*,nrunv)
316    else return lprie "Scheme is not twostep in time"$
317  go to a$
318r:for each a in setdiff(coords!*,lcor) do
319    if a='t then x:=list('u0!*,nrunv)
320      else ex:=list('times,tcar get(a,'index),get(a,'step),get(a,'wave))
321                  . ex$
322  return simp list('times,x,list('expp,'plus . ex))
323end$
324
325procedure fouriersubs$
326begin
327  scalar x,i$
328  for each a in '(expp u1!* u0!*) do put(a,'simpfn,'simpiden)$
329  x:=unvars!*$
330  i:=1$
331a:if null x then go to b$
332  put(car x,'nrunvar,i)$
333  i:=i+1$
334  x:=cdr x$
335  go to a$
336b:flag(unvars!*,'full)$
337  for each a in unvars!* do put(a,'simpfn,'simpfour)$
338  for each a in coords!* do
339    if not(a='t) then
340    <<x:=intern compress append(explode 'h,explode a)$
341      put(a,'step,x)$
342      if not flagp(a,'uniform) then put(x,'simpfn,'simpiden)$
343      x:=intern compress append(explode 'k,explode a)$
344      put(a,'wave,x)$
345      x:=intern compress append(explode 'a,explode a)$
346      put(a,'angle,x)$
347      put(a,'maxden,1) >>$
348  algebraic for all z,y,v let
349    expp((z+y)/v)=expp(z/v)*expp(y/v),
350    expp(z+y)=expp z*expp y
351end$
352
353procedure gonsubs$
354begin
355  scalar xx$
356  algebraic for all z,y,v clear expp((z+y)/v),expp(z+y)$
357  for each a in coords!* do
358    if not(a='t) then
359      <<xx:=list('quotient,
360                 list('times,
361                      get(a,'maxden),
362                      get(a,'angle)),
363                 get(a,'step))$
364        setk(get(a,'wave),aeval xx)$
365        xx:=list('times,
366                 get(a,'wave),
367                 get(a,'step))$
368        mathprint list('setq,
369                       get(a,'angle),
370                       if get(a,'maxden)=1 then xx
371                         else list('quotient,
372                                   xx,
373                                   get(a,'maxden))) >>$
374  algebraic for all x let expp x=cos x+i*sin x$
375  algebraic for all x,n such that numberp n and n>1 let
376    sin(n*x)=sin x*cos((n-1)*x)+cos x*sin((n-1)*x),
377    cos(n*x)=cos x*cos((n-1)*x)-sin x*sin((n-1)*x)$
378  for each a in unvars!* do
379    <<put(a,'simpfn,'simpiden)$
380      remprop(a,'nrunvar) >>
381end$
382
383procedure remfourier$
384<<algebraic for all x clear expp x$
385  algebraic for all x,n such that numberp n and n>1 clear
386    sin(n*x),cos(n*x)$
387  for each a in coords!* do
388    if not(a='t) then
389    <<remprop(a,'step)$
390      remprop(remprop(a,'wave),'avalue)$
391      remprop(a,'maxden) >> >>$
392
393operator numberp$
394
395endmodule;
396
397end;
398