1module transfrm; 2 3% Cobasis transformations 4 5% Author: David Hartley 6 7% Redistribution and use in source and binary forms, with or without 8% modification, are permitted provided that the following conditions are met: 9% 10% * Redistributions of source code must retain the relevant copyright 11% notice, this list of conditions and the following disclaimer. 12% * Redistributions in binary form must reproduce the above copyright 13% notice, this list of conditions and the following disclaimer in the 14% documentation and/or other materials provided with the distribution. 15% 16% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 17% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 18% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 19% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 20% CONTRIBUTORS 21% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 22% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 23% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 24% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 25% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 26% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 27% POSSIBILITY OF SUCH DAMAGE. 28% 29 30 31COMMENT. Data structure: 32 33 xform ::= list of 1-form kernel . 1-form pf 34 35endcomment; 36 37 38fluid '(xvars!* kord!* subfg!*); 39global '(!*sqvar!*); 40 41fluid '(cfrmcob!* cfrmcrd!* cfrmdrv!* cfrmrsx!* !*edssloppy); 42 43 44% Type coersions 45 46 47symbolic procedure !*a2xform u; 48 % u:list of equation -> !*a2xform:xform 49 % Turn off subfg!* to stop let rules being applied. 50 % should remove x=x entries 51 for each j in getrlist u collect 52 if eqexpr j then !*a2k lhs j . xpartitop rhs j where subfg!* = nil 53 else rerror(eds,000,"Incorrectly formed transform"); 54 55 56symbolic procedure !*xform2map x; 57 % x:xform -> !*xform2map:map 58 foreach p in x collect 59 car p . mk!*sq !*pf2sq cdr p; 60 61 62symbolic procedure !*map2xform u; 63 % u:map -> !*map2xform:xform 64 % Turn off subfg!* to stop let rules being applied. 65 (for each x in u collect 66 car u . xpartitop cadr u) where subfg!* = nil; 67 68 69symbolic procedure !*xform2drv x; 70 % x:xform -> !*xform2drv:drv 71 % Turn off subfg!* to stop let rules being applied. 72 (foreach p in x collect 73 {'replaceby,car p,!*pf2a cdr p}) where subfg!* = nil; 74 75 76symbolic procedure !*xform2sys x; 77 % x:xform -> !*xform2sys:sys 78 foreach p in x collect addpf(!*k2pf car p,negpf cdr p); 79 80 81% Transform 82 83 84put('transform,'rtypefn,'getrtypecar); 85put('transform,'cfrmfn,'transformcfrm); 86put('transform,'edsfn,'transformeds); 87 88 89symbolic procedure transformcfrm(m,x); 90 % m:cfrm, x:list of equation -> transformcfrm:cfrm 91 % Transform m using map x. 92 checkcfrm xformcfrm(m,!*a2xform x); 93 94 95symbolic procedure transformeds(s,x); 96 % s:eds, x:list of equation -> transform:eds 97 % pulls back s using map x 98 xformeds(s,!*a2xform x); 99 100 101symbolic procedure xformcfrm(m,x); 102 % m:cfrm, x:xform -> xformcfrm:cfrm 103 % Apply transform x to m, where x may be either old=f(new,old) or 104 % new=f(old). 105 xformcfrm1(m,car u,cadr u,caddr u) 106 where u = getxform(x,cfrm_cob m); 107 108 109symbolic procedure xformcfrm1(m,x,y,new); 110 % m:cfrm, x,y:xform, new:cob -> xformcfrm1:cfrm 111 % Apply transform x to m, where x is old=f(new,old), y is x inverse, 112 % and new gives the new cobasis elements. 113 begin scalar p,z; 114 m := copycfrm m; 115 z := pair(foreach p in x collect car p,new); 116 cfrm_cob m := % replace old by new in-place 117 foreach k in cfrm_cob m collect % sublis here destroys kernels 118 if p := atsoc(k,z) then cdr p else k; 119 cfrm_crd m := % retain all old coordinates (may appear in eds) 120 reverse union(foreach k in new join if exact k then {cadr k}, 121 reverse cfrm_crd m); 122 cfrm_drv m := % add new differentials and structure equations 123 append(xformdrv(cfrm_drv m,x), 124 append(!*xform2drv foreach p in x join if exact car p then {p}, 125 structeqns(y,x))); 126 if !*edssloppy then m := updatersx m; % invxform may have added new 127 % rsx 128 m := purgecfrm m; 129 return m; 130 end; 131 132 133symbolic procedure xformcfrm0(m,x,new); 134 % m:cfrm, x:xform, new:cob -> xformcfrm0:cfrm 135 % Cut down version of xformcfrm1 which doesn't update structure 136 % equations. Useful when following operations are purely algebraic. 137 begin scalar p,z; 138 m := copycfrm m; 139 z := pair(foreach p in x collect car p,new); 140 cfrm_cob m := % replace old by new in-place 141 foreach k in cfrm_cob m collect % sublis here destroys kernels 142 if p := atsoc(k,z) then cdr p else k; 143 cfrm_crd m := % retain all old coordinates (may appear in eds) 144 reverse union(foreach k in new join if exact k then {cadr k}, 145 reverse cfrm_crd m); 146 if !*edssloppy then m := updatersx m; % invxform may have added new 147 % rsx 148 m := purgecfrm m; 149 return m; 150 end; 151 152 153symbolic procedure xformdrv(d,x); 154 % d:drv, x:xform -> xformdrv:drv 155 % Apply xform to drv. Must suppress active rules, for example if d a 156 % => d b is active and x = {d b => d a}, then after applying x, it 157 % will be undone immediately. 158 pullbackdrv(d,!*xform2map x) where subfg!* = nil; 159 160 161symbolic procedure updatersx m; 162 % m:cfrm -> updatersx:cfrm 163 % Reinstall restrictions in s from global variable, typically 164 % after solvepfsys when !*edssloppy is t. 165 begin 166 m := copycfrm m; 167 cfrm_rsx m := foreach f in purge cfrmrsx!* collect !*pf2a f; 168 return m; 169 end; 170 171 172symbolic procedure xformeds(s,x); 173 % s:eds, x:xform -> xformeds:eds 174 % Apply transform x to m, where x may be either old=f(new,old) or 175 % new=f(old). 176 % possibly changes kernel order 177 xformeds1(s,car u,cadr u,caddr u) 178 where u = getxform(x,edscob s); 179 180 181 182symbolic procedure xformeds1(s,x,y,new); 183 % s:eds, x,y:xform, new:cob -> xformeds1:eds 184 % Apply transform x to m, where x is old=f(new,old), y is x inverse, 185 % and new gives the new cobasis elements. Changes background 186 % coframing. 187 begin scalar k; 188 s := copyeds s; 189 % Transform coframing 190 eds_cfrm s := xformcfrm1(eds_cfrm s,x,y,new); 191 % Make sure old get eliminated (and add new to kord!* for safety) 192 k := updkordl append(foreach p in x collect car p,new); 193 x := !*xform2sys x; 194 % Transform rest of eds 195 eds_sys s := foreach f in eds_sys s collect xreduce(xreorder f,x); 196 eds_ind s := foreach f in eds_ind s collect xreduce(xreorder f,x); 197 remkrns s; 198 s := purgeeds!* s; 199 rempropeds(s,'jet0); 200 foreach f in {'solved,'reduced} do rempropeds(s,f); 201 setkorder k; 202 s := normaleds s; % Refine this a bit? 203 setcfrm eds_cfrm!* s; 204 return s; 205 end; 206 207 208symbolic procedure xformeds0(s,x,new); 209 % s:eds, x:xform, new:cob -> xformeds0:eds 210 % Cut down version of xformeds1 which doesn't care about structure 211 % equations (some are lost). Useful when following operations are 212 % purely algebraic. Changes background coframing. 213 begin scalar k; 214 s := copyeds s; 215 % Transform coframing (ignore structure equations) 216 eds_cfrm s := xformcfrm0(eds_cfrm s,x,new); 217 % Make sure old get eliminated (and add new to kord!* for safety) 218 k := updkordl append(foreach p in x collect car p,new); 219 x := !*xform2sys x; 220 % Transform rest of eds 221 eds_sys s := foreach f in eds_sys s collect xreduce(xreorder f,x); 222 eds_ind s := foreach f in eds_ind s collect xreduce(xreorder f,x); 223 remkrns s; 224 s := purgeeds!* s; 225 rempropeds(s,'jet0); 226 foreach f in {'solved,'reduced} do rempropeds(s,f); 227 setkorder k; 228 s := normaleds s; % Refine this a bit? 229 setcfrm eds_cfrm!* s; 230 return s; 231 end; 232 233 234symbolic procedure getxform(x,cob); 235 % x:xform, cob:cob -> getxform:{xform,xform,cob} 236 % Analyse transform x, which may be either old=f(new,old) or 237 % new=f(old). The sense is established by cob, which contains the old 238 % cobasis. Return value is {x,y,new} where x is in the sense old = 239 % f(new,old), and y is the inverse of x (ie new = f(old)). The 240 % inverse y is calculated only if x is old = f(new,old) and there 241 % are anholonomic forms in new. 242 begin scalar old,new,y; 243 foreach p in x do 244 << new := union(xpows cdr p,new); 245 old := car p . old >>; 246 if not xnp(old,cob) then % x is new=f(old), must invert 247 << y := x; x := invxform x; 248 new := old; old := foreach p in x collect car p >>; 249 new := sort(setdiff(new,cob),'termordp); 250 edsdebug("New cobasis elements...",new,'cob); 251 edsdebug("... replacing old cobasis elements",old,'cob); 252 if length new neq length old or not subsetp(old,cob) then 253 rerror(eds,000,"Bad cobasis transformation"); 254 if not allexact new and null y then 255 y := invxform x; % for structure equations 256 return {x,y,new}; 257 end; 258 259 260% Structure equations 261 262 263symbolic procedure xformdrveval u; 264 % u:{rlist,rlist} or {rlist} -> xformdrveval:rlist 265 begin scalar x,y,k; 266 y := !*a2xform car u; 267 x := if cdr u then !*a2xform cadr u else invxform y; 268 k := updkordl foreach p in x collect car p; 269 y := structeqns(y,x); 270 setkorder k; 271 return makelist y; 272 end; 273 274symbolic procedure xformdrveval u; 275 % u:{rlist,rlist} or {rlist} -> xformdrveval:rlist 276 begin scalar x,y,xvars!*; 277 y := !*a2xform car u; 278 x := if cdr u then !*a2xform cadr u else invxform y; 279 y := structeqns(y,x); 280 return makelist y; 281 end; 282 283 284symbolic procedure structeqns(y,x); 285 % y,x:xform -> structeqns:list of rule 286 % y is the inverse of x, and d lhs x are known. 287 % Returns rules for d lhs y. 288 begin scalar ok; 289 ok := updkordl foreach p in x collect car p; 290 x := !*xform2sys x; 291 y := foreach p in y join 292 if not exact car p then 293 {{'replaceby, 294 {'d,car p}, 295 !*pf2a xreduce(exdfpf cdr p,x)}}; 296 setkorder ok; 297 return y; 298 end; 299 300symbolic procedure structeqns(y,x); 301 % y,x:xform -> structeqns:list of rule 302 % y is the inverse of x, and d lhs x are known. 303 % Returns rules for d lhs y. 304 begin scalar ok; 305 ok := updkordl sort(foreach p in x collect car p,function ordop); 306 x := !*xform2sys x; 307 y := foreach p in y join 308 if not exact car p then 309 {{'replaceby, 310 {'d,car p}, 311 !*pf2a xreduce(exdfpf cdr p,x)}}; 312 setkorder ok; 313 return y; 314 end; 315 316% Inverting tranformations 317 318 319put('invert,'rtypefn,'quotelist); 320put('invert,'listfn,'inverteval); 321 322symbolic procedure inverteval(u,v); 323 % u:{prefix list of eqn} -> inverteval:{prefix list of eqn} 324 % u is unevaluated. 325 makelist foreach p in invxform !*a2xform(u := reval car u) collect 326 {'equal,car p,!*pf2a1(cdr p,v)}; 327 328 329symbolic procedure invxform x; 330 % x:xform -> invxform:xform 331 % Returns inverse transformation. Selects kernels to eliminate based 332 % on prevailing order 333 begin scalar old,y,k, subfg!*; 334 subfg!* := nil; 335 foreach p in x do old := union(xpows cdr p,old); 336 old := sort(old,'termordp); 337 % ensure old eliminated, and add new to kord!* for safety 338 k := updkordl append(old,foreach p in x collect car p); 339 edsdebug("Inverting transform",nil,nil); 340 y := solvepfsys1(!*xform2sys x,old); % invert transformation 341 if cadr y then 342 rerror(eds,000,"Cobasis transform could not be inverted"); 343 setkorder k; 344 return foreach f in car y collect 345 lpow f . negpf xreorder red f; 346 end; 347 348 349symbolic procedure tmpind s; 350 % s:eds -> tmpind:{eds,xform} 351 % Returns s with eds_ind s all kernels, transforming to a new 352 % cobasis if necessary. Second return value is nil if no change 353 % made, or the list of transformation relations. Structure 354 % equations are not transformed, so s should be closed first if 355 % necessary. NB. Background coframing changed. 356 begin scalar new,x; 357 if singleterms eds_ind s then return {s,nil}; 358 new := foreach f in eds_ind s collect mkform!*(intern gensym(),1); 359 x := invxform pair(new,eds_ind s); 360 updkordl foreach p in x collect car p; 361 return {xformeds0(s,x,new),x}; 362 end; 363 364endmodule; 365 366end; 367