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