1%******************************************************************** 2module shortening$ 3%******************************************************************** 4% Routines for algebraically combining de's to reduce their length 5% Author: Thomas Wolf 6% Jan 1998 7% 8 9% BSDlicense: ***************************************************************** 10% * 11% Redistribution and use in source and binary forms, with or without * 12% modification, are permitted provided that the following conditions are met: * 13% * 14% * Redistributions of source code must retain the relevant copyright * 15% notice, this list of conditions and the following disclaimer. * 16% * Redistributions in binary form must reproduce the above copyright * 17% notice, this list of conditions and the following disclaimer in the * 18% documentation and/or other materials provided with the distribution. * 19% * 20% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * 21% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * 22% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * 23% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR CONTRIBUTORS BE * 24% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR * 25% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF * 26% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS * 27% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN * 28% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * 29% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * 30% POSSIBILITY OF SUCH DAMAGE. * 31%****************************************************************************** 32 33symbolic procedure alg_length_reduction(arglist)$ 34% Do one length-reducing combination of two equations 35% If cadddr arglist(=:p) is =pdes or =nil then shorten anyone with anyone 36% else shorten only elements (i.e. equations) of cadddr arglist 37 38begin scalar pdes,l,l1,p,cpu,gc$ 39 cpu:=time()$ gc:=gctime()$ 40 pdes:=car arglist$ 41 42 if expert_mode then <<l:=selectpdes(pdes,2); p:=nil>> 43 else << 44 l:=pdes; p:=cadddr arglist; % which may be nil 45 if l eq p then p:=nil; % then no selective reduction of only one PDE 46 % i.e. shorten any one with any one 47 >>; 48 !*rational_bak := cons(!*rational,!*rational_bak)$ 49 if !*rational then algebraic(off rational)$ 50 if struc_eqn then << 51 while l do << 52 if is_algebraic(car l) then l1:=cons(car l,l1); 53 l:=cdr l 54 >>$ 55 l:=reverse l1 56 >>$ 57 if l and cdr l and (l1:=err_catch_short(pdes,l,p)) then 58 <<for each a in cdr l1 do pdes:=drop_pde(a,pdes,nil)$ 59 for each a in car l1 do 60 if a then pdes:=eqinsert(a,pdes)$ 61 % checking whether a new equation contains a function which occurs 62 % only there, then not decoupling this equation anymore: 63 for each a in car l1 do 64 if a then dec_fct_check(a,pdes)$ 65 l:=nil; 66 for each a in car l1 do if a then l:=cons(a,l); 67 l:=list(pdes,cadr arglist,l)>> 68 else l:=nil$ 69 70 if print_ and !*time then << 71 write " time : ", time() - cpu, 72 " ms GC time : ", gctime() - gc," ms. " 73 >>$ 74 75 if !*rational neq car !*rational_bak then 76 if !*rational then algebraic(off rational) else algebraic(on rational)$ 77 !*rational_bak:= cdr !*rational_bak$ 78 79 return l$ 80end$ 81 82%------------------- 83 84symbolic procedure err_catch_short(allpdes,pdes,p2)$ 85% allpdes is the complete list of pdes as needed in mkeqSQ 86% pdes are the ones to be used for shortening 87% if p2 neq nil then only p2 shall be reduced else any one in pdes 88begin scalar g,h,bak,kernlist!*bak,kord!*bak,mi,newp,p1,bakup_bak,s; 89 90 if t then << % wrap up in an errorset shell 91 bak:=max_gc_counter$ max_gc_counter:=my_gc_counter+max_gc_short; 92 bakup_bak:=backup_; backup_:='max_gc_short$ 93 kernlist!*bak:=kernlist!*$ 94 kord!*bak:=kord!*$ 95 h:=errorset({'shorten_pdes,mkquote pdes,mkquote p2},nil,nil) 96 where !*protfg=t; 97 kord!*:=kord!*bak; 98 kernlist!*:=kernlist!*bak$ 99 erfg!*:=nil; 100 backup_:=bakup_bak; 101 max_gc_counter:=bak; 102 if (errorp h) or (car h=nil) then return nil; 103 h:=car h 104 >> else << 105 h:=shorten_pdes(pdes,p2)$ 106 if null h then return nil 107 >>; 108 109 mi:=car h; % mi ={list of 2 used eqns} 110 newp:=cdr h; % newp=({list of new values} . {list of eqns to be dropped}) 111 p1:=0; 112 for each pc in cdr newp do p1:=p1+get(pc,'terms); 113 h:=for each pc in car newp collect 114 if sqzerop pc then <<nequ_:=add1 nequ_;nil>> 115 else mkeqSQ(pc,nil,nil,fctsort union(get(car mi,'fcts), 116 get(cadr mi,'fcts)), 117 union(get(car mi,'vars), 118 get(cadr mi,'vars)),allflags_,t,list(0),nil,allpdes); 119 for each pc in h do if pc then p1:=p1-get(pc,'terms); 120 s:=(h . cdr newp); % = ({list of new eqn names} . {list of eqns to be dropped}) 121 122if nil then 123 if print_ then << 124 if tr_short then << 125 for each g in cdr newp do <<write g,": "$typeeq g>>$ 126 for each g in h do if null g then 127 <<write "This gives identity 0=0."$terpri()>> 128 else 129 <<write g,": "$typeeq g>>$ 130 >>$ 131 if p1>=0 then write "shortening by ",p1," term" else 132 if p1< 0 then write "increasing by ",-p1," term"$ 133 if (p1 neq 1) and (p1 neq -1) then write"s"$ 134 135 % h is the list of new equations: if there is exactly one eqn then 136 if h and car h and null cdr h then << 137 write" to now ",g:=get(car h,'terms)," term"$ 138 if g neq 1 then write"s"$ 139 write"."$ terpri()$ 140 >> else 141 if null h then << 142 write" to 0=0 ."$ terpri() 143 >> 144 145 >>; 146 147 if p1>0 then for each pc in cdr newp do drop_pde(pc,nil,nil) 148 else for each pc in h do if pc then << 149 % so that no cycle happens and the equation is 150 % back length reduced with the same equations 151 add_rl_with(car mi,pc); 152 add_rl_with(cadr mi,pc) 153 >>$ 154 155 return s 156end$ 157 158%------------------- 159 160symbolic procedure alg_length_reduce_1(arglist)$ 161% Do one length-reducing combination of a single equation wrt all others 162begin scalar p,q,h,p_to_eval,hcp$ %,cpu,gc$ %,toevalcp 163% cpu:=time()$ gc:=gctime()$ 164 if !*batch_mode then return nil; 165 if print_ or null old_history then << 166 terpri()$ 167 change_prompt_to ""$ 168 write"All equations:"$terpri()$ 169 listprint(car arglist)$terpri()$terpri()$ 170 write"Which equation should be length reduced: "$ 171 >>$ 172 p:=termread()$ 173 174% if print_ or null old_history then << 175% write"Specify a list of equations to be used to reduce ",p$terpri()$ 176% write"Input `;' to select all equations apart from ",p,": "$ 177% >>$ 178% l:=termlistread()$ 179% if l=nil then l:=car arglist; 180% l:=delete(p,l); 181% toevalcp:=for each h in l collect (h . flagp(h,'to_eval))$ 182% for each h in l do remflag1(h,'to_eval)$ 183 p_to_eval:=flagp(p,'to_eval)$ 184 flag1(p,'to_eval)$ 185 186 % the complete list of pdes is needed in 187 % alg_length_reduction() because mkeqSQ() needs it 188 189 190 repeat << 191 h:=alg_length_reduction({car arglist,cadr arglist,caddr arglist,p})$ 192 193 if h then << 194 q:=setdiff(car h,cons(p,car arglist)); 195 if q then <<arglist:=h; p:=car q>>; 196 hcp:=h; 197 >> 198 >> until null h or freeof(car h,p); 199% for each q in toevalcp do if cdr q then flag1(car q,'to_eval)$ 200 if not p_to_eval and h and not freeof(car h,p) then remflag1(p,'to_eval)$ 201 restore_interactive_prompt()$ 202 return hcp 203end$ 204 205%------------------- 206 207symbolic procedure is_algebraic(p)$ 208if get(p,'no_derivs) then t else nil$ 209 210%------------------- 211 212symbolic procedure shorten_pdes(des,p2)$ 213% uses global variable largest_fully_shortened which is the name of the 214% equation furthest in pdes (largest number of terms) that is fully 215% simplified wrt. all earlier equations 216% returns (mi . output_of_shorten) 217if not pairp des or not pairp cdr des then (nil . nil) else 218begin scalar p2rl,p,pc,pcc,newp,ineq_pre$ %,p2_to_eval$ 219 for each p1 in ineq_ do 220 if one_termpSF(numr p1) then ineq_pre:=cons(prepsq p1,ineq_pre)$ 221 if alg_poly then ineq_pre:=sort_according_to(ineq_pre,ftem_)$ 222 223 % find the pair of pdes not yet reduced with each other 224 % with the lowest product of their number of terms % printlength's 225 % mi is a list of 2 integers, composing a rational number which is 226 % supposed to be minimal for the 2 equations to be shortened. 227 228 if p2 then << 229 if print_ and tr_short then << 230 write"Trying to shorten the following equation: "$%terpri() 231 >>$ 232 % try only to shorten p2 wrt. to equations with at most twice as many terms. 233 p2rl:=2*get(p2,'terms)$ 234 for each p in des do 235 if (p neq p2) and get(p,'terms)<=p2rl then newp:=cons(p,newp); 236 if null newp then return nil; 237 des:=reversip newp$ 238 pcc:={p2}$ 239 des:=nconc(des,pcc); 240 pcc:=cons(1,pcc); 241 newp:=nil; 242 p2rl:=nil; 243 >> else << 244 if print_ and tr_short then << 245 write"Trying to shorten the following equations: "$%terpri() 246 >>$ 247 pcc:=des; 248 if largest_fully_shortened then << 249 while pcc and largest_fully_shortened neq car pcc do pcc:=cdr pcc; 250 if null pcc or null cdr pcc then return nil 251 >> 252 >>$ 253 if null cdr pcc then return nil; 254 255 repeat << 256 pcc:=cdr pcc$ 257 p2:=car pcc$ 258 if print_ and tr_short then write p2,"(",get(p2,'terms),") "$ 259% p2_to_eval:=flagp(p2,'to_eval)$ 260 p2rl:=get(p2,'rl_with)$ 261 262 pc:=des; 263 while null newp and not(pc eq pcc) do << 264 if (not member(car pc, p2rl)) % and 265 % We do not test (not member(p2,get(car pc,'rl_with))) because both 266 % ways of reduction are tested and in add_rl_with() are stored. 267 then <<if (backup_='max_gc_short) and (my_gc_counter > max_gc_counter) then 268 rederr "Stop in shorten_pdes()."$ 269 newp:=shorten(car pc,p2,ineq_pre); 270 if null newp then <<add_rl_with(car pc,p2);pc:=cdr pc>> 271 >> 272 else pc:=cdr pc 273 >>; 274 if null newp then largest_fully_shortened:=car pcc 275 >> until newp or null cdr pcc$ 276 277 return if newp then ({car pc,p2} . newp) 278 else nil 279end$ 280 281%------------------- 282 283symbolic procedure partition_3(de,AnB,A)$ 284% l is an equation in SQ-form, 285% returning {(#_of_terms_in_1 . ( 1 . inhomog )), 286% (#_of_terms_in_SFcoeff_1 . (flin_1 . SFcoeff_1)),...} 287% where flin_1,.. are in AnB (in A and in B, i.e. intersection of A and B) 288% and the leading term of inhomog must not be in A. 289% It is assumed that l in linear on all variables of A. 290% 291begin scalar l,l1,mv; 292 l:=reorder numr get(de,'sqval); 293 while l and not domainp l and member(mv:=mvar l,AnB) do << 294 l1:=cons(cons(no_of_tm_sf lc l,cons(mv,lc l)), l1)$ 295 l:=red l 296 >>; 297 while l and not domainp l and member(mvar l,A) do l:=red l; 298 return if l then cons(cons(no_of_tm_sf l,cons(1,l)), l1) % inhom case 299 else l1 % homog case 300end$ 301 302%------------------- 303 304symbolic procedure setdiff_according_to(a,b,ft)$ 305% Computes a\b where a and b are subsets of ft and are 306% sorted according to ft. 307% exception: partial differential equations, when elements of a,b 308% include derivatives, then ft is not used and execution is slower 309% which does not matter for fewer functions 310begin scalar amb,noft,bcp$ 311 while a and b do 312 if car a eq car b then <<a:=cdr a;b:=cdr b;ft:=cdr ft>> else << 313 314 if null noft and (pairp car a or pairp car b) then noft:=t; 315 if noft then << 316 bcp:=b$ 317 while bcp and (car a neq car bcp) do bcp:=cdr bcp; 318 if null bcp then amb:=cons(car a,amb); 319 a:=cdr a 320 >> else << 321 while car ft neq car a and 322 car ft neq car b do ft:=cdr ft; 323 324 if car ft eq car a then << 325 amb:=cons(car a,amb); 326 a:=cdr a 327 >> else b:=cdr b; 328 ft:=cdr ft 329 >> 330 >>$ 331 return 332 if a then nconc(reversip amb,a) 333 else reversip amb 334end$ 335 336%------------------- 337 338symbolic procedure add_fct_kern(de,ineq_pre)$ 339if alg_poly then 340 if lin_problem then (get(de,'fcts) . nil) % algebraic linear 341 else 342 if flin_ then begin scalar l,n; % algebraic partially linear 343 l:=get(de,'fct_kern_lin)$ 344 n:=get(de,'fct_kern_nli)$ 345 if null l and null n then << 346 n:=setdiff_according_to(l:=get(de,'fcts),flin_,ftem_); 347 l:=setdiff_according_to(l,n,ftem_); 348 n:=setdiff_according_to(n,ineq_pre,ftem_); 349 put(de,'fct_kern_lin,l); 350 put(de,'fct_kern_nli,n) 351 >>; 352 return (l . n) 353 end else begin scalar n; % algebraic fully non-linear 354 n:=get(de,'fct_kern_nli)$ 355 if null n then << 356 n:=setdiff_according_to(get(de,'fcts),ineq_pre,ftem_); 357 put(de,'fct_kern_nli,n); 358 >>; 359 return (nil . n) 360 end 361else % not alg_poly 362 if lin_problem then begin scalar l,r,s; % differential linear 363 l:=get(de,'fct_kern_lin); 364 if null l then << 365 r:=get(de,'fcts); 366 for each s in get(de,'kern) do 367 if not freeoflist(s,r) then l:=cons(s,l); 368 put(de,'fct_kern_lin,l); 369 >>; 370 return (l . nil) 371 end else 372 if flin_ then begin scalar l,n,r,s; % differential partially linear 373 l:=get(de,'fct_kern_lin)$ 374 n:=get(de,'fct_kern_nli)$ 375 if null l and null n then << 376 r:=get(de,'fcts); 377 for each s in get(de,'kern) do 378 if not freeoflist(s,r) then 379 if not freeoflist(s,flin_) then l:=cons(s,l) 380 else if not member(s,ineq_pre) then n:=cons(s,n); 381 put(de,'fct_kern_lin,l); 382 put(de,'fct_kern_nli,n) 383 >>; 384 return (l . n) 385 end else begin scalar n,r,s; % differential fully non-linear 386 n:=get(de,'fct_kern_nli)$ 387 if null n then << 388 r:=get(de,'fcts); 389 for each s in get(de,'kern) do 390 if not freeoflist(s,r) and not member(s,ineq_pre) then n:=cons(s,n); 391 put(de,'fct_kern_nli,n) 392 >>; 393 return (nil . n) 394 end$ 395 396%------------------- 397 398symbolic procedure shorten(de1,de2,ineq_pre)$ 399% shorten the two pdes with each other 400% returns a dotted pair, where car is a list of the values of new pdes 401% (currently only one) and cdr is a list of names of pdes to be dropped 402% (currently also only one) 403begin scalar a,b,r,l1,l2,nl1,nl2,l1ul2,l1ml2,l2ml1,l1il2,oldorder, 404 de1p,de2p,termsof1,termsof2,flip,m1,m2,n1,n2,ql,ql1,ql2,maxcancel, 405 take_first,tr_short_local,no_factors_x_1,no_factors_x_2,homo, 406 changed_korder,DoNotReplace1,DoNotReplace2,tryboth,replace1,bestq; 407 % ,max_success; 408 409 % take_first:=t; 410 % =t is not so radical, --> resulting eqn.s may be longer in total 411 % so the complete crack computation tends to be slower. 412 413 %tr_short_local:=t; 414 if tr_short_local then deprint list({'!*sq,get(de1,'sqval),t}, 415 {'!*sq,get(de2,'sqval),t} )$ 416 417 % 15.12.07: The homogeneity test is replaced by the more general concept of 418 % strictly linearly occuring functions flin_, so the following is commented out 419 420 if fhom_ and (1=car get(de1,'hom_deg) ) 421 and (1=car get(de2,'hom_deg) ) 422 % and ((cadr get(de1,'hom_deg)) neq 423 % (cadr get(de2,'hom_deg)) ) % changed when it didn't work 424 % for two factorizable equations with same flin_ factor and same degree 425 %non-flin_factors 426 then homo:=t; 427 428 % Which equation can be replaced? 429 430 termsof1:=get(de1,'terms)$ 431 termsof2:=get(de2,'terms)$ 432 433 if null flagp(de1,'to_eval) or 434 (homo and (cadr get(de1,'hom_deg)<cadr get(de2,'hom_deg))) or 435 (2*termsof1<termsof2) then DoNotReplace1:=t; 436 437 if null flagp(de2,'to_eval) or 438 (homo and (cadr get(de2,'hom_deg)<cadr get(de1,'hom_deg))) or 439 (2*termsof2<termsof1) then DoNotReplace2:=t; 440 441 if DoNotReplace1 and DoNotReplace2 then return nil$ 442 443 %---- determination of l1,l2 nl1,nl2 444 % In the following l1 and l2 are all functions and their derivatives in 445 % de1, resp. de2 which definitely turn up linearly in both equations 446 % (maybe there are even more) and which can not be multipliers such that it 447 % is possible to partition both equations. This is the case if the problem 448 % is linear or in the case of non null flin_ which shall stay linear. 449 % nl1, nl2 are the kernels in d1,d2 which do involve ftem_ but not flin_ 450 % and not ineq_ kernels. 451 452 nl1:=add_fct_kern(de1,ineq_pre)$ 453 nl2:=add_fct_kern(de2,ineq_pre)$ 454 l1:=sort_according_to(car nl1,ftem_)$ nl1:=cdr nl1$ 455 l2:=sort_according_to(car nl2,ftem_)$ nl2:=cdr nl2$ 456 % NOTE: nl1,nl2 are currently NOT sorted as they are currently not used below. 457 458 if l1 and l2 then << % fully linear or at least partially linear with flin_ 459 460 %---- partitioning both equations into subexpressions 461 % which have a chance to cancel each other 462 l1ml2:=setdiff_according_to(l1,l2,ftem_); % l1 - l2 463 l1il2:=setdiff_according_to(l1,l1ml2,ftem_); % intersection 464 if null l1il2 then return nil$ 465 l2ml1:=setdiff_according_to(l2,l1,ftem_); % l2 - l1 466 l1ul2:=union(l1,l2); % union 467 468 %---- writing both equations as polynomials in elements of l1ul2 469 oldorder:=setkorder append(l1il2,append(l1ml2,l2ml1)); changed_korder:=t$ 470 de1p:=partition_3(de1,l1il2,l1ml2)$ 471 de2p:=partition_3(de2,l1il2,l2ml1)$ 472 setkorder oldorder; 473 474 if (null de1p) or (null de2p) then return nil; 475 %---- drop inhomogeneity part if only one equation is inhomogeneous 476 %---- there is an inhom. part if cadar de?p = 1 477 if (cadar de1p = 1) and (cadar de2p neq 1) then de1p:=cdr de1p; 478 if (cadar de2p = 1) and (cadar de1p neq 1) then de2p:=cdr de2p; 479 480 % Now comes a stronger restriction: The maximum that can be canceled 481 % is sum of the minima of terms of each pair of the de1p,de2p sublists 482 % corresponding to the coefficients of different ftem functions/deriv. 483 484 a:=de1p; b:=de2p; n2:=nil; 485 while a and b do 486 if (cadar a) neq (cadar b) then << 487 % This case should not happen but it can happen if flin_ is not nil 488 % but at least one equation is non-linear and then a function can appear 489 % in an equation but will only turn up as a coefficient of another function 490 % and not as cadar a or cadar b and then one must step forward in one of 491 % the two lists a,b to get (cadar a) = (cadar b) 492 r:=b; while r and ((cadar r) neq (cadar a)) do r:=cdr r; 493 if r then b:=r 494 else << 495 r:=a; while r and ((cadar r) neq (cadar b)) do r:=cdr r; 496 a:=r 497 >> 498 >> else << 499 n1:=if (caar a)<(caar b) then caar a else caar b; 500 % n1 is min of terms of the coefficients of the same ftem function/der. 501 n2:=cons(2*n1,n2); 502 a:=cdr a; b:=cdr b; 503 >>$ 504 505 % maxcancel is the maximal number of cancellations in all the remaining 506 % runs of short() (which of course decreases from run to run). 507 508 maxcancel:=list(0); 509 n1:=0; 510 while n2 do << 511 n1:=n1+car n2; 512 n2:=cdr n2; 513 maxcancel:=cons(n1,maxcancel); 514 >>; 515 516 if tr_short_local then << 517 write"-----"$terpri()$ 518 write"flin_=",flin_$terpri()$ 519 write"flin_\ftem_=",setdiff_according_to(flin_,ftem_,ftem_)$ 520 write"ftem_\flin_=",setdiff_according_to(ftem_,flin_,ftem_)$ 521 write"get(de1,'fcts)=",get(de1,'fcts)$terpri()$ 522 write"get(de2,'fcts)=",get(de2,'fcts)$terpri()$ 523 write"l1=",l1$terpri()$ 524 write"nl1=",nl1$terpri()$ 525 write"l2=",l2$terpri()$ 526 write"nl2=",nl2$terpri()$ 527 write"l1ml2=",l1ml2$terpri()$ 528 write"l2ml1=",l2ml1$terpri()$ 529 write"l1il2=",l1il2$terpri()$ 530 write"length de1p=",length de1p$terpri()$ 531 write"length de2p=",length de2p$terpri()$ 532 write"de1p=",de1p$terpri()$ 533 write"de2p=",de2p$terpri()$ 534 write"---------"$terpri()$ 535 write"de1:",de1," with ",termsof1," terms"$terpri()$ 536 a:=de1p; 537 while a do << 538 write "caar a=", caar a;terpri()$ 539 write "cadar a=",cadar a;terpri()$ 540 write "cddar a=",cddar a;terpri()$ 541 write "cddar a="$algebraic write lisp {'!*sq,(cddar a . 1),t}; 542 a:=cdr a; 543 >>;terpri()$ 544 write"de2:",de2," with ",termsof2," terms"$terpri()$ 545 a:=de2p; 546 while a do << 547 write "caar a=", caar a;terpri()$ 548 write "cadar a=",cadar a;terpri()$ 549 write "cddar a=",cddar a;terpri()$ 550 write "cddar a="$algebraic write lisp {'!*sq,(cddar a . 1),t}; 551 a:=cdr a; 552 >>;terpri()$ 553 >> 554 >> else << % fully non-linear, i.e. at least one equation has no flin_ 555 de1p:=list cons(termsof1,cons(1,numr get(de1,'sqval)))$ 556 de2p:=list cons(termsof2,cons(1,numr get(de2,'sqval)))$ 557 maxcancel:=if termsof1<termsof2 then {2*termsof1,0} 558 else {2*termsof2,0} 559 >>$ 560 561 if ((car maxcancel)<termsof1) then if DoNotReplace1 then return nil 562 else DoNotReplace2:=t$ 563 if ((car maxcancel)<termsof2) then if DoNotReplace2 then return nil 564 else DoNotReplace1:=t$ 565 566 % flip is determined, so that after a flip (if necessary) de2 is replaced, 567 % i.e. if n1 is the number of terms of de1 then at least n1 terms must be 568 % canceled so that the sum a*de1+b*de2 has not more terms than de2. 569 % This would be achieved if at least n1/2 pair-cancellations would occur. 570 571 % Although nl1,2 are not used further below, the use of nconc instead of 572 % append let once to a wrong 'fcts and 'fct_kern_nli property of an equation 573 if DoNotReplace2 or (termsof2<termsof1) then << 574 flip:=t; 575 a:=de1p; de1p:=de2p; de2p:=a; 576 no_factors_x_2:=if l1 or null l2 then nl2 577 else append(nl2,setdiff_according_to(l2,ineq_pre,ftem_)); 578 no_factors_x_1:=if l2 or null l1 then nl1 579 else append(nl1,setdiff_according_to(l1,ineq_pre,ftem_)); 580 n1:=termsof2; 581 n2:=termsof1 582 >> else << 583 no_factors_x_1:=if l1 or null l2 then nl2 584 else append(nl2,setdiff_according_to(l2,ineq_pre,ftem_)); 585 no_factors_x_2:=if l2 or null l1 then nl1 586 else append(nl1,setdiff_according_to(l1,ineq_pre,ftem_)); 587 n1:=termsof1; 588 n2:=termsof2 589 >>; 590 591 if null DoNotReplace1 and null DoNotReplace2 then tryboth:=t; 592% max_success:=if tryboth then if termsof1<termsof2 then 2*termsof1 593% else 2*termsof2 594% else 2*termsof1$ 595 596 % To allow the generation of new *length increased* equations 597 % which have to be *added* and which should not be shortened again 598 % (not to get into a cycle) , i.e. 'to_eval=nil, one has to lower 599 % n1 here. n1 is the length of the equation to be kept. 600 % n1:=n1-4$ 601 602 if length_inc_alg neq 1.0 then << 603 m1:=reval algebraic(ceiling lisp {'quotient,n1,length_inc_alg}); 604 m2:=reval algebraic(ceiling lisp {'quotient,n2,length_inc_alg}) 605 >> else << 606 m1:=n1; 607 m2:=n2 608 >>; 609 610 repeat << % one shortening 611 b:=short(ql1,ql2,cddar de1p,cddar de2p,m1,m2,2*(caar de1p), 612 car maxcancel-cadr maxcancel,cadr maxcancel,%max_success, 613 take_first,no_factors_x_1,no_factors_x_2,tryboth)$ 614 % take_first:=car b; b:=cdr b; % to activate see end of short() 615 if b then << 616 ql1:=car b; 617 ql2:=cadr b; 618 a:=cddr b; 619 if a and take_first then << % the result 620 de1p:=car a; % numerator of a successful quotient 621 de2p:=cdr a; % denominator of a successful quotient 622 >> else << % the next pairing 623 de1p:=cdr de1p; 624 de2p:=cdr de2p; 625 >>; 626 maxcancel:=cdr maxcancel; 627 >> else a:=nil; 628 >> until (null b) or % failure 629 (a and take_first) or % early success 630 null de1p$ % end of all possible run --> needs check 631 632 if b and (null take_first) then << 633 % search of the best shortening 634 % car find..: highest number of saved terms in the quotient list 635 % de1p: numerator of the best quotient so far 636 % de2p: denominator of the best quotient so far 637 638 if null ql2 and (null tryboth or null ql1) then << 639 write"##### SOMETHING IS WRONG ###"$terpri()$ 640 write" short() should have recognized that there is no successful quotient."$ 641 terpri() 642 >>$ 643 644%write"ql2="$prettyprint ql2$ 645%write"ql1="$prettyprint ql1$ 646 647 ql2:=find_best_quotient(ql2)$ % works also if ql2=nil 648 ql1:=find_best_quotient(ql1)$ 649 % ql1 and ql2 are now of different structure, see find_best_quotient() 650 651%write"n1=",n1$terpri()$ 652%write"n2=",n2$terpri()$ 653%write"flip=",flip$terpri()$ 654%write"de1=",de1," de2=",de2$terpri()$ 655%write"car ql1=",car ql1$terpri()$ 656%write"car ql2=",car ql2$terpri()$ 657%write"replace1-1 = ",replace1$terpri()$ 658 659 if (caar ql1 - n2) > % at least n2 terms would have to be dropped 660 (caar ql2 - n1) then << % at least n1 terms would have to be dropped 661 replace1:=t$ 662 ql:=ql1 663 >> else ql:=ql2; 664 bestq:=car ql; 665 r:=n1+n2-car bestq; % the new number of terms 666 de1p:=cadr bestq; 667 de2p:=cddr bestq; 668 ql:=cdr ql % remaining multipliers 669 670 >>; 671%write"replace1-2 = ",replace1$terpri()$ 672 673 return 674 if null b or 675 (take_first and null a) then nil else % num and denom are de1p, de2p 676% if t then 677<< 678 %--- computing the shorter new equation 679 % (local variables l1,l2,nl1,nl2,m1,m2 are re-used) 680 % nl2 shall be the name and l2 the value of the equation to be replaced 681 % nl1 shall be the name and l1 the value of the other equation 682 % m1 shal be the sum of rational multipliers to l1 683 nl2:=if replace1 then if flip then de2 else de1 684 else if flip then de1 else de2; 685%write"nl2=",nl2$ 686 nl1:=if nl2=de2 then de1 else de2; 687 l1:=get(nl1,'sqval); l2:=get(nl2,'sqval); 688 689 m1:=if nl2=de2 then if flip then (de1p . de2p) else (de2p . de1p) % m1 is the multiplier to equation nl1 690 else if flip then (de2p . de1p) else (de1p . de2p); 691 l2:=subtrsq(l2,multsq(m1,l1)); 692 693%%%% TEST: 694%if nil then 695if null take_first and (r neq (m2:=no_of_tm_sf numr l2)) then << 696 write"##############################################################################"$ 697 write"Wrong expected length value:"$terpri()$ 698 write"expected value: r=",r$terpri()$ 699 write"real no of terms: ",m2$terpri()$ 700 701 algebraic<< 702 off nat$ 703 write"de1:=",lisp {'!*sq,get(de1,'sqval),t}$ 704 write"de2:=",lisp {'!*sq,get(de2,'sqval),t}$ 705 write"m1 =",lisp {'!*sq,m1,t}$ 706 on nat 707 >>; 708 write"ql="$terpri()$ 709 prettyprint ql 710>>; 711 712 if take_first then r:=no_of_tm_sf numr l2 713 else 714 while ql do << % r is the current length and ql is a list of further multipliers 715 if bestq neq car ql then << 716 m2:=if nl2=de2 then if flip then ((cadar ql) . (cddar ql)) else ((cddar ql) . (cadar ql)) 717 else if flip then ((cddar ql) . (cadar ql)) else ((cadar ql) . (cddar ql)); 718 a:=subtrsq(l2,multsq(m2,l1)); 719 b:=no_of_tm_sf numr a; 720 if b<r then <<l2:=a; r:=b; m1:=addsq(m1,m2)>> 721 >>; 722 ql:=cdr ql; 723 >>$ 724 if in_cycle({11,stepcounter_, 725 get(l1,'printlength),length get(l1,'fcts),numr m1, 726 get(l2,'printlength),length get(l2,'fcts),denr m1 }) then << 727 if print_ and tr_short then << 728 write"To avoid a loop, a possible shortening was not performed."$ 729 terpri() 730 >>; 731 nil 732 >> else << 733 if print_ then << 734 a:=get(nl2,'terms); b:=a-r; 735 if null tr_short then write"Shortening by ",b,if b=1 then " term" else " terms", 736 " to now ",r,if r=1 then " term." else " terms." 737 else << 738 terpri()$ 739 if r=0 then << 740 write"Equation ",nl2,"(",a,") is deleted as it is a consequence of ",nl1,"(", 741 get(nl1,'terms),") :"$ 742 algebraic(write " 0 = ",lisp {'!*sq,(numr subtrsq(mksq(nl2,1),multsq(m1,mksq(nl1,1))) . 1),t}) 743 744 >> else << 745 if null car recycle_eqns then m2:=mkid(eqname_,nequ_) 746 else m2:=caar recycle_eqns$ 747 % m2 is the name of the new equation. The same name will be generated 748 % in err_catch_short() which calls mkeqSQ() which calls new_pde(). 749 write"Equation ",nl2,"(",a,") is shortened by ",b,if b=1 then " term" else " terms", 750 " using ",nl1,"(",get(nl1,'terms),") ","and ","is ","replaced by:"$ 751 algebraic(write lisp m2,"(",r,") = ", 752 lisp {'!*sq,(numr subtrsq(mksq(nl2,1),multsq(m1,mksq(nl1,1))) . 1),t}) 753 >> 754 >> 755 >>; 756 ({l2} . {nl2}) 757 >> 758 >> 759% else << 760% if flip then <<a:=get(de2,'sqval); b:=get(de1,'sqval)>> 761% else <<a:=get(de1,'sqval); b:=get(de2,'sqval)>>$ 762% if print_ and tr_short then << 763% if null car recycle_eqns then n1:=mkid(eqname_,nequ_) 764% else n1:=caar recycle_eqns$ 765% algebraic write "The new equation ",n1," = ", 766% lisp {'difference,{'times,!*f2a de2p,if flip then de2 else de1}, 767% {'times,!*f2a de1p,if flip then de1 else de2} }, 768% " replaces " 769% >>$ 770% a:=subtrsq(if de2p=1 then a 771% else multsq((de2p . 1),a), 772% if de1p=1 then b 773% else multsq((de1p . 1),b) )$ 774% 775% if in_cycle(cons(11,cons(stepcounter_,if flip then { 776% get(de2,'printlength),length get(de2,'fcts),de2p, 777% get(de1,'printlength),length get(de1,'fcts),de1p } 778% else { 779% get(de1,'printlength),length get(de1,'fcts),de1p, 780% get(de2,'printlength),length get(de2,'fcts),de2p }))) 781% then nil 782% else (list a . list(if replace1 then if flip then de2 else de1 783% else if flip then de1 else de2)) 784% >> 785end$ % of shorten 786 787%------------------- 788 789symbolic procedure clean_num(qc,j)$ 790begin 791 scalar qc1,nall$ 792 return 793 if 2*(cdaar qc)<=j then t else << 794 qc1:=car qc; % the representative and list to proportional factors 795 nall:=cdar qc1; 796 while cdr qc1 do 797 if (cdadr qc1)+nall<=j then rplacd(qc1,cddr qc1) 798 else qc1:=cdr qc1; 799 if qc1=car qc then t else nil % whether empty or not after cleaning 800 >> 801end$ 802 803%-------------------- 804 805symbolic procedure clean_den(qc,j)$ 806begin 807 scalar qcc$ 808 qcc:=qc$ 809 while cdr qc do 810 if clean_num(cdr qc,j) then rplacd(qc,cddr qc) 811 else qc:=cdr qc$ 812 return null cdr qcc % Are there any numerators left? 813end$ 814 815%------------------- 816 817symbolic procedure find_best_quotient(ql)$ 818% - ql is a list of denominator classes: (dcl1 dcl2 dcl3 ...) 819% - each denominator class dcli is a dotted pair (di . nclist) where 820% - di is the denominator and 821% - nclist is a list of numerator classes. 822% Each numerator class is a list with 823% - first element: (ncl . n) where ncl is the numerator 824% up to a rational numerical factor and n is the number of 825% occurences of ncl (up to a rational numerical factor) 826% - further elements: (nfi . ni) where nfi is the numerical 827% proportionality factor and ni the number of occurences 828% of this factor 829begin scalar %r, 830 s,cp,n,m,qlist,bestnu,bestq; 831% r:=0; 832 bestq:=(0 . (nil . nil)); 833 while ql do << 834%write"denumerator=",caar ql$terpri()$ 835 s:=cdar ql; % the numerator list of the first denominator class 836 while s do << % for each numerator class 837 cp:=car s; s:=cdr s; % cp is the first numerator class 838 n:=cdar cp; % n is the sum of occurences of all numerators in this class 839%write"numerator=",caar cp$terpri()$ 840%write"numerator frequency=",n$terpri()$ 841%for each m in cdr cp do write"num fac=",car m," frequency=",cdr m$terpri()$ 842 843 % finding the best numerical factor bestnu in the numerator class 844 bestnu:=cdr cp; 845%write"bestnu1=",bestnu$terpri()$ 846 while cdr bestnu do 847 if cdar bestnu <= cdadr bestnu then bestnu:=cdr bestnu 848 else rplacd(bestnu,cddr bestnu); 849%write"bestnu2=",bestnu$terpri()$ 850 % multiplying 851 % the numerator of the numerical factor to the numerator of the quotient and 852 % the denomiator of the numerical factor to the denomiator of the quotient 853 854 m:=n + cdar bestnu; 855%write"m=",m$terpri()$ 856% qlist:=cons((m . % saving 857% ((if caaar bestnu=1 then caar cp 858% else multf(caar cp,caaadr cp)) . % numerator 859% (if cdaar bestnu=1 then caar ql 860% else multf(caar ql,cdaadr cp)) )), % denominator 861% qlist); 862 863 qlist:=cons((m . % saving 864 ((if caaar bestnu=1 then caar cp 865 else multf(caar cp,caaar bestnu)) . % numerator 866 (if cdaar bestnu=1 then caar ql 867 else multf(caar ql,cdaar bestnu)) )), % denominator 868 qlist); 869 870% if m>r then <<r:=m;bestq:=cdar qlist>> 871%write"qlist=",qlist$terpri()$ 872%write"bestq1=",bestq$terpri()$ 873 if m>car bestq then bestq:=car qlist 874%;write"bestq2=",bestq$terpri()$ 875 >>; 876 ql:=cdr ql 877 >>; 878%write"r=",r," bestq=",bestq$terpri()$ 879 return (bestq . qlist) 880end$ 881 882%-------------------- 883 884symbolic procedure short(ql1,ql2,d1,d2,n1,n2,n1_now,max_save_now,max_save_later, 885 %max_success, 886 take_first,no_factors_x_1,no_factors_x_2,tryboth)$ 887 % - ql1,ql2 are lists of quotients so far determined where de1 will be 888 % multiplied with the denominator of `the' quotient and de2 is multiplied 889 % with the numerator of that quotient. if tryboth=nil then ql1=nil. 890 % Numerators of ql2 do not depend on no_factors_x_2 so that these quotients 891 % are used to replace de2 and denominators of ql1 do not depend on 892 % no_factors_x_1 so that these quotients are used to replace de1. 893 % - d1,d2 are two subexpressions of de1, de2. 894 % - The full equations de1,de2 have n1,n2 terms. n1 many must be saved in total 895 % if de2 is to be replaced and n2 many must be saved if de1 is to be replaced. 896 % - n1_now is 2 times the number of terms of the input expression d1. 897 % - max_save_now is the maximum number of cancellations in this 898 % run of short() based on 2*min(terms of d1, terms of d2) 899 % - max_save_later is the maximal number of terms that could be saved in all 900 % the remaining short() runs under this shorten() call. 901 % - if no_factors_x_i <> nil if non-linear equations --> must check that di 902 % is not multiplied with a factor from no_factors_x_i 903 % - if tryboth=t then try to shorten any one of both equations, else try only 904 % to shorten de2, i.e. the mother of d2. 905begin 906 scalar LastQSuccess,d2cop,j1,j2,e1,e2,q,dcl,nu,ldcl,lnu,h,repl1,repl2,allowedq 907 ,max_coeff_len_reached 908%,bynow1,bynow2 909%,qisnumber,qcount 910; 911 % tr_short_local:=t; 912 % - n1_now is the maximum number of terms cancelling each other 913 % in this run of short, based on 2*(number of remaining terms of d1 914 % still to check). Initially n1_now=max_save_now or n1_now>max_save_now 915 % but as the computation goes on, n1_now decreases, so eventually it will 916 % be less than max_save_now. 917 918%n1:=reval n1; 919%n2:=reval n2; 920 j1:=max_save_later + max_save_now$ 921 j2:=n1-j1$ 922 j1:=n2-j1$ 923 % - The following j2-value is the minimal number of cancellations that 924 % should already have been found to have a chance to replace de2. 925 % - The following j1-value is the minimal number of cancellations that 926 % should already have been found to have a chance to replace de1. 927 % - LastQSuccess = t iff the last quotient gives a length reduction. 928 % This is the case if >n1 terms are killed for a deletion of de2 or >n2 929 % terms are killed for a deletion of de1. 930%bynow1:=0; bynow2:=0; qcount:=0; 931 repeat << % for each term e1 of d1 932 n1_now:=n1_now-2; 933 e1:=first_term_SF d1; d1:=subtrf(d1,e1); 934 %---- divide each term e1 of d1 through all terms e2 of d2 935 d2cop:=d2; 936 %---- continue as long as possible and un-successful 937 while d2cop and not(take_first and LastQSuccess) do << % correct version, 3% slower than: 938% while d2cop and null LastQSuccess do << % The first quotient which proves to be 939 % beneficial is nearly always the best 940 e2:=first_term_SF d2cop; d2cop:=subtrf(d2cop,e2); 941 q:=cancel(e1 ./ e2); % otherwise already successful 942 943%qcount:=add1 qcount$ 944%if q neq resimp q then << 945%write"###### q <> q !!!. "$terpri()$ 946%write"qold=",q$terpri()$ 947%write"qnew=",resimp quotsq((e1 . 1),(e2 . 1))$terpri()$ 948%write"qnew=",simp reval {'!*sq,quotsq((e1 . 1),(e2 . 1)),nil}$terpri()$ 949%write"re1=",reval {'!*sq,q,nil}$terpri()$ 950%write"re2=",reval {'!*sq,(e1 . 1),(e2 . 1),nil}$terpri()$ 951%>>$ 952 953%if (denr q = 1) and (domainp numr q) then << 954% qisnumber:=t; 955% bynow1:=add1 bynow1; 956% if 5184 = car q then bynow2:=add1 bynow2; 957%>> else qisnumber:=nil; 958%%if qisnumber then write"t " else write "n "$ 959 960 %---- Determination of 961 % ldcl: the numerical factor of the denominator 962 % dcl: the rest of the denominator 963 dcl:=cdr q; % dcl is the denominator of the current quotient 964 repl1:=tryboth; 965% if numberp dcl or (pairp dcl and car dcl= '!:GI!:) then <<ldcl:=dcl;dcl:=1>> 966 if domainp dcl then <<ldcl:=dcl;dcl:=1>> 967 else << 968 if repl1 and no_factors_x_1 and not freeoflist(dcl,no_factors_x_1) then repl1:=nil; 969 h:=dcl; 970% while <<ldcl:=lc h; not(numberp ldcl or car ldcl = '!:GI!:)>> do h:=ldcl; 971 while not domainp(ldcl:=lc h) do h:=ldcl; 972 rplacd(car h,1)$ % This seems to work reliably, i.e. we do not need here: 973 % if ldcl neq 1 then dcl:=numr cancel(dcl ./ ldcl) 974 >>; 975 976 %---- Determination of 977 % lnu: the numerical factor of the numerator 978 % nu: the rest of the numerator 979 nu:=car q; % nu is the numerator of the current quotient 980% if numberp nu or (pairp nu and car nu= '!:GI!:) then <<lnu:=nu;nu:=1;repl2:=t;allowedq:=t>> 981 if domainp nu then <<lnu:=nu;nu:=1;repl2:=t;allowedq:=t>> 982 else << 983 if no_factors_x_2 and not freeoflist(nu,no_factors_x_2) then << 984 allowedq:=repl1; 985 repl2:=nil 986 >> else <<repl2:=t;allowedq:=t>>$ 987 if allowedq then << 988 h:=nu; 989% while <<lnu:=lc h; not(numberp lnu or car lnu = '!:GI!:)>> do h:=lnu; 990 while not domainp(lnu:=lc h) do h:=lnu; 991 % rplacd(car h,1)$ % is not possible as it might change e1 (as it has occured) 992 if lnu neq 1 then nu:=numr cancel(nu ./ lnu) 993 >> 994 >>; 995 996 %---- Do numerical coefficients get too long? 997 if allowedq and 998 ((numberp lnu and (( lnu>max_coeff_len) or ( lnu<-max_coeff_len))) or 999 (pairp lnu and (car lnu='!:gi!:) and 1000 ((cadr lnu>max_coeff_len) or (cadr lnu<-max_coeff_len) or 1001 (cddr lnu>max_coeff_len) or (cddr lnu<-max_coeff_len) ) ) or 1002 (numberp ldcl and ((ldcl>max_coeff_len) or (ldcl<-max_coeff_len))) or 1003 (pairp ldcl and (car ldcl='!:gi!:) and 1004 ((cadr ldcl>max_coeff_len) or (cadr ldcl<-max_coeff_len) or 1005 (cddr ldcl>max_coeff_len) or (cddr ldcl<-max_coeff_len) ) ) ) then << 1006 if tr_short and null max_coeff_len_reached then << 1007 write"### Num. factors grew too large in shortening."$ 1008 terpri()$ 1009% write"N" 1010 >>; 1011 max_coeff_len_reached:=t 1012 >> else 1013 1014 %---- Record the quotient, check whether it is successful 1015 if null allowedq then LastQSuccess:=nil 1016 else 1017 if repl1 and (null repl2 or (n1>n2)) then << % to replace de1 1018%if qisnumber then write"###### wrongly adding to ql1!"$ 1019 h:=add_quotient(ql1,nu,lnu,dcl,ldcl,j1)$ 1020 if car h > n2 then LastQSuccess:=t 1021 else LastQSuccess:=nil$ 1022 ql1:=cdr h 1023 >> else << % to replace de1 1024 h:=add_quotient(ql2,nu,lnu,dcl,ldcl,j2)$ 1025 if car h > n1 then LastQSuccess:=t 1026 else LastQSuccess:=nil$ 1027 ql2:=cdr h 1028 >> 1029 1030 >>$ % all terms of d2 1031 1032 % Computing updated values of j1,j2 as n1_now has 1033 % been lowered and may now be < max_save_now 1034 if n1_now<max_save_now then << 1035 j1:=max_save_later + n1_now; 1036 j2:=n1-j1$ 1037 j1:=n2-j1$ 1038 1039 % Deleting quotients that have no chance. 1040 if j1>0 then ql1:=clean_ql(ql1,j1)$ 1041 if j2>0 then ql2:=clean_ql(ql2,j2)$ 1042 >>$ 1043 1044 >> % all terms of d1 1045 until (null d1 ) or % everything divided 1046 (take_first and LastQSuccess ) or % successful: saving > cost 1047 (((null tryboth ) or 1048 ((j1 > 0) and (null ql1)) ) and 1049 ( j2 > 0) and (null ql2) )$ % all quotients are too rare 1050 1051% write"bynow1=",bynow1," bynow2=",bynow2," qcount=",qcount$terpri()$ 1052 1053 return 1054 if ((null tryboth ) or 1055 ((j1 > 0) and (null ql1)) ) and 1056 ( j2 > 0) and (null ql2 ) then nil else 1057 if null d1 then (ql1 . (ql2 . nil)) else 1058 (ql1 . (ql2 . q )); 1059 % take_first and LastQSuccess 1060 1061end$ % of short 1062 1063%-------------------- 1064 1065symbolic procedure clean_ql(ql,j)$ 1066begin scalar qc; 1067 while ql and clean_den(car ql,j) do ql:=cdr ql; 1068 if ql then << 1069 qc:=ql; 1070 while cdr qc do 1071 if clean_den(cadr qc,j) then rplacd(qc,cddr qc) 1072 else qc:=cdr qc 1073 >>; 1074 return ql 1075end$ 1076 1077%-------------------- 1078 1079symbolic procedure add_quotient(ql,nu,lnu,dcl,ldcl,j)$ 1080% - ql is a list of denominator classes: (dcl1 dcl2 dcl3 ...) 1081% - each denominator class dcli is a dotted pair (di . nclist) where 1082% - di is the denominator and 1083% - nclist is a list of numerator classes. 1084% Each numerator class is a list with 1085% - first element: (ncl . n) where ncl is the numerator 1086% up to a rational numerical factor and n is the number of 1087% occurences of ncl (up to a rational numerical factor) 1088% - further elements: (nfi . ni) where nfi is the numerical 1089% proportionality factor and ni the number of occurences 1090% of this factor 1091% - lnu is the numerical factor of the numerator and 1092% nu is the rest of the numerator 1093% - ldcl is the numerical factor of the denominator 1094% dcl is the rest of the denominator 1095% - The j-value is the minimal number of cancellations that should 1096% already have been found (including cancellations from this quotient) 1097% to have a chance for shortening. 1098% 1099% It returns (number_of_saved_terms_by_this_quotient_as_known_so_far . 1100% quotient_list) 1101 1102begin scalar nall,m,qc,qq,preqc$ 1103 1104 if null !*complex then << % This case handles rational although 1105 % they are not supposed to come up 1106 if pairp ldcl then << 1107 write"####### ldcl=",ldcl," is not an integer!"$ terpri() 1108 >>$ 1109 if pairp lnu and car lnu neq '!:rn!: then << 1110 write"####### lnu=",lnu," is not rational!"$ terpri() 1111 >>$ 1112 if fixp lnu then qq:=(lnu . ldcl) else 1113 % all following under the assumption: car lnu=!:RN!: 1114 if ldcl=1 then qq:=cdr lnu 1115 else <<m:=cancel (cadr lnu ./ ldcl); 1116 qq:=(car m . (cddr lnu * cdr m))>> 1117 >> else % ON COMPLEX, lnu and ldcl could both be Gaussian 1118 % integers but without real common integer factor 1119 % --> no need to do cancel( ./ ) 1120 if fixp ldcl then qq:=(lnu . ldcl) 1121 else qq:=quotsq((lnu . 1),(ldcl . 1)); 1122 1123%if domainp ldcl and (ldcl neq 1) then << 1124%write"***** ldcl=",ldcl," lnu=",lnu$ 1125%algebraic(on rational); 1126%m:=quotsq((lnu . 1),(ldcl . 1)); 1127%lnu:=numr m; 1128%ldcl:=denr m; m:=nil; 1129%write" | ldcl=",ldcl," lnu=",lnu$terpri()$ 1130%>>; 1131 1132 %---- search for the denominator class 1133 qc:=ql; 1134 while qc and (dcl neq caar qc) do qc:=cdr qc; 1135 1136 1137 if null qc then % denominator class not found 1138 if j <= 0 then << % add denominator class 1139 nall:=2; % to give a sum of 2 1140 ql:=cons((dcl . list(list((nu . 1),(qq . 1)))), ql) 1141 >> else nall:=0 % too late to add this denominator 1142 else << % denominator class has been found 1143 1144 %---- now search of the numerator class 1145 qc:=cdar qc; % qc is the list of numerator classes nclist 1146 while qc and (nu neq caaar qc) do <<preqc:=qc; qc:=cdr qc>>; 1147 1148 if null qc then % numerator class not found 1149 if j <= 0 then << % add numerator class %@@@@ 1150 nall:=2; 1151 rplacd(preqc,list(list((nu . 1),(qq . 1))) ) 1152 >> else nall:=0 % too late to add this numerator 1153 else <<% numerator class found 1154 nall:=add1 cdaar qc; % increasing the total number of occur. 1155 rplacd(caar qc,nall); 1156 1157 %---- now search for the numerical factor 1158 qc:=cdar qc; 1159 while qc and (qq neq caar qc) do <<preqc:=qc;qc:=cdr qc>>; 1160 if null qc then << % numerical factor not found 1161 rplacd(preqc,list((qq . 1))); 1162%;if (nu=1) and (dcl=1) then 1163%if nall neq bynow1 or ((lnu=5184) and (bynow2 neq 1)) then << 1164%write"&&&&&& bynow1=",bynow1," nall=",nall," bynow2=",bynow2," m=",1," q=",q$terpri() 1165%>>$ 1166 nall:=add1 nall 1167 >> else << 1168 m:=add1 cdar qc$ 1169 rplacd(car qc,m); 1170%;if (nu=1) and (dcl=1) then 1171%if nall neq bynow1 or ((lnu=5184) and (bynow2 neq m)) then << 1172%write"&&&&&& bynow1=",bynow1," nall=",nall," bynow2=",bynow2," m=",m," q=",q$terpri() 1173%>>$ 1174 nall:=nall+m 1175 >> 1176 >> % numerator class found 1177 >>$ % denominator class found 1178 return (nall . ql) 1179end$ % of add_quotient 1180 1181%-------------------- 1182 1183symbolic procedure drop_lin_dep(arglist)$ 1184% drops linear dependent equations 1185begin scalar pdes,tr_drop,p,cp,incre,newpdes,m,h,s,r,a,v, 1186 vli,indp,indli,conli,mli,success$ 1187 % the pdes are assumed to be sorted by the number of terms, 1188 % shortest come first 1189 % vli is the list of all `independent variables' v in this lin. algebra 1190 % computation, i.e. a list of all different products of powers of 1191 % derivatives of ftem functions and constants 1192 % format: ((product1, v1, sum1),(product2, v2, sum2),...) 1193 % where sumi is the sum of all terms of all equations multiplied 1194 % with the multiplier of that equation 1195 % indli is a list marking whether equations are necessarily lin 1196 % indep. because they involve a `variable' v not encountered yet 1197 % mli is the list of multipliers of the equations 1198 pdes:=car arglist$ 1199 % tr_drop:=t$ 1200 if pdes and cdr pdes then << 1201 while pdes do << 1202 p:=car pdes; pdes:=cdr pdes; newpdes:=cons(p,newpdes); 1203 m:=gensym()$ 1204 a:=sort_partition(p,nil,get(p,'fcts),nil); 1205 if tr_drop then <<write "new eqn:"$prettyprint a; 1206 write"multiplier for this equation: m=",m$terpri() 1207 >>$ 1208 indp:=nil; 1209 for each h in a do << 1210 s:=car h; 1211 % Does s occur in vli? 1212 % Adding the terms multiplied with the multiplier to the corresp. sum 1213 cp:=vli; 1214 while cp and (s neq caar cp) do cp:=cdr cp; 1215 if tr_drop then << 1216 write"searched for: s=",s$terpri(); 1217 if cp then write"found: car cp=",car cp$terpri()$ 1218 >>$ 1219 if cp then <<r:=cadar cp; 1220 incre:=reval {'quotient, 1221 {'times,m,r, 1222 if cadr h>1 then cons('plus,caddr h) 1223 else caaddr h}, 1224 s}; 1225 rplaca(cddar cp,cons(incre,caddar cp)) 1226 >> 1227 else <<r:=if (pairp s) or (numberp s) then gensym() else s; 1228 indp:=s; 1229 incre:=reval {'quotient, 1230 {'times,m,r, 1231 if cadr h>1 then cons('plus,caddr h) 1232 else caaddr h}, 1233 s}; 1234 vli:=cons({s,r,{incre}},vli) 1235 >>; 1236 if tr_drop then << 1237 write"corresponding symbol: r=",r$terpri()$ 1238 1239 write"upd: incre=",incre$terpri()$ 1240 write"vli="$prettyprint vli 1241 >>$ 1242 >>; 1243 mli:=cons(m,mli); 1244 indli:=cons(indp,indli) 1245 >>$ 1246 1247 % Formulating a list of conditions 1248 while vli do << 1249 v:=caddar vli; vli:=cdr vli; 1250 conli:=cons(if cdr v then cons('plus,v) 1251 else car v, 1252 conli) 1253 >>; 1254 1255 % Now the investigation of linear independence 1256 pdes:=nil; % the new list of lin. indep. PDEs 1257 while cdr newpdes do << 1258 if tr_drop then << 1259 terpri()$ 1260 if car indli then write"lin. indep. without search of ",car newpdes, 1261 " due to the occurence of ",car indli 1262 else write"lin. indep. investigation for ",car newpdes$ 1263 >>; 1264 if car indli then pdes:=cons(car newpdes,pdes) 1265 else << 1266 s:=cdr solveeval {cons('list,subst(1,car mli,conli)),cons('list,cdr mli)}; 1267 if s then <<drop_pde(car newpdes,nil,nil)$ % lin. dep. 1268 success:=t$ 1269 if print_ then << 1270 terpri()$ 1271 write"Eqn. ",car newpdes, 1272 " has been dropped due to linear dependence."$ 1273 >> 1274 >> 1275 else <<pdes:=cons(car newpdes,pdes); % lin. indep. 1276 if tr_drop then << 1277 terpri()$ 1278 write"Eqn. ",car newpdes," is lin. indep."$ 1279 >> 1280 >>; 1281 >>; 1282 newpdes:=cdr newpdes; 1283 indli:=cdr indli; 1284 conli:=subst(0,car mli,conli); 1285 mli:=cdr mli 1286 >>; 1287 pdes:=cons(car newpdes,pdes) 1288 >>; 1289 return if success then list(pdes,cadr arglist) 1290 else nil 1291end$ 1292 1293%-------------------- 1294 1295symbolic procedure find_1_term_eqn(arglist)$ 1296% checks whether a linear combination of the equations can produce 1297% an equation with only one term 1298if not lin_problem then nil else 1299begin scalar pdes,tr_drop,p,cp,incre,m,h,s,r,a,v, 1300 vli,indp,indli,conli,mli,mpli,success, 1301 sli,slilen,maxlen,newconli,newpdes,newp,fl,vl$ 1302%tr_drop:=t$ 1303 if tr_drop then terpri()$ 1304 pdes:=car arglist$ 1305 newpdes:=pdes$ 1306%--------------------------------- 1307% if struc_eqn then << 1308% cp:=pdes; 1309% while cp do << 1310% if is_algebraic(car cp) then r:=cons(car cp,r) 1311% else s:=cons(car cp,s); 1312% cp:=cdr cp 1313% >>; 1314% r:=nil; 1315% s:=nil; 1316% >>$ 1317% Drop all PDEs which have at least two derivs which no other have 1318%--------------------------------- 1319 if pdes and cdr pdes then << 1320 while pdes do << 1321 p:=car pdes; pdes:=cdr pdes; 1322 m:=gensym()$ 1323 if tr_drop then <<terpri()$write"multiplier m=",m$terpri()>>$ 1324 a:=sort_partition(p,nil,get(p,'fcts),nil); 1325 for each h in a do << 1326 s:=car h; 1327 % Does s occur in vli? 1328 % Adding the terms multiplied with the multiplier to the corresp. sum 1329 cp:=vli; 1330 while cp and (s neq caar cp) do cp:=cdr cp; 1331 if tr_drop then << 1332 write"searched for: s=",s$terpri(); 1333 if cp then <<write"found: car cp=",car cp$terpri()$>> 1334 >>$ 1335 if cp then <<r:=cadar cp; 1336 incre:=reval {'quotient, 1337 {'times,m,%r, 1338 if cadr h>1 then cons('plus,caddr h) 1339 else caaddr h}, 1340 s}; 1341 rplaca(cddar cp,cons(incre,caddar cp)) 1342 >> 1343 else <<r:=if (pairp s) or (numberp s) then gensym() else s; 1344 indp:=s; 1345 incre:=reval {'quotient, 1346 {'times,m,%r, 1347 if cadr h>1 then cons('plus,caddr h) 1348 else caaddr h}, 1349 s}; 1350 vli:=cons({s,r,{incre}},vli) 1351 >>; 1352 if tr_drop then << 1353 write"corresponding symbol: r=",r$terpri()$ 1354 1355 write"upd: incre=",incre$terpri()$ 1356 write"vli="$prettyprint vli 1357 >>$ 1358 >>; 1359 mli:=cons(m,mli); 1360 mpli:=cons((m . p),mpli); 1361 indli:=cons(indp,indli) 1362 >>$ 1363 1364 % Formulating a list of conditions 1365 while vli do << 1366 sli:=cons(caar vli,sli); 1367 v:=caddar vli; vli:=cdr vli; 1368 conli:=cons(if cdr v then cons('plus,v) 1369 else car v, 1370 conli) 1371 >>; 1372 1373 % Now the investigation of the existence of equations with only one 1374 % term 1375 slilen:=length sli; 1376 mli :=cons('list, mli); 1377 conli:=cons('list,conli); 1378 if tr_drop then << 1379 write"sli=",sli$terpri()$ 1380 algebraic(write"mli=",mli)$ 1381 algebraic(write"conli=",conli)$ 1382 write"mpli=",mpli$terpri()$ 1383 >>; 1384 for h:=1:slilen do << % for each possible single term 1385 newp:=car sli;sli:=cdr sli; 1386 pdes:=newpdes; 1387 while pdes and 1388 ((get(car pdes,'terms)>1) or 1389 (not zerop reval {'difference,get(car pdes,'val),newp})) do 1390 pdes:=cdr pdes; 1391 if null pdes then << 1392 cp:=conli; 1393 for s:=1:h do cp:=cdr cp; 1394 rplaca(cp,reval {'plus,1,car cp}); 1395 if tr_drop then << 1396 write"h=",h$terpri()$ 1397 algebraic(write"new conli=",conli)$ 1398 >>; 1399 1400 s:=cdr solveeval {conli,mli}; 1401 if (null s) and tr_drop then <<write"no success"$terpri()>>$ 1402 if s then << % found 1-term equation 1403 if null success then 1404 for each p in newpdes do << 1405 fl:=union(get(p,'fcts),fl); 1406 vl:=union(get(p,'vars),vl) 1407 >>$ 1408 success:=t$ 1409 maxlen:=0$ 1410 s:=cdar s; % first solution (lin. system), dropping 'list 1411 1412 % now find the equation to be replaced by the 1-term-equation 1413 % find the non-vanishing m in s, such that the corresponding p in 1414 % mpli has a maximum number of terms 1415 while s do << 1416 if caddar s neq 0 then << 1417 r:=cadar s; 1418 cp:=mpli; 1419 while caar cp neq r do cp:=cdr cp; 1420 if get(cdar cp,'terms)>maxlen then << 1421 p:=cdar cp; % p will be the equation to be replaced 1422 m:=r; 1423 maxlen:=get(p,'terms); 1424 >> 1425 >>$ 1426 s:=cdr s 1427 >>$ 1428 1429 % Now replacing the old equation p by the new 1-term equation in conli: 1430 r:=0; 1431 newconli:=nil$ 1432 while conli do << 1433 v:=subst(0,m,car conli)$ 1434 conli:=cdr conli$ 1435 if r=h then << 1436 v:=reval {'plus,{'times,m,newp},v}$ 1437 >>$ 1438 newconli:=cons(v,newconli); 1439 r:=add1(r) 1440 >>$ 1441 conli:=reverse newconli$ 1442 1443 % the new equation: 1444 newp:=mkeqSQ(nil,nil,newp,fl,vl,allflags_,t,list(0),nil,newpdes); 1445 % last argument is nil as no new inequalities can result 1446 % if new equation has only one term 1447 % --> has been changed as this may change, e.g. ineq_or 1448 newpdes:=cons(newp,newpdes); 1449 if print_ then << 1450 terpri()$ 1451 write"The new equation ",newp$ typeeq newp$ 1452 write" replaces ",p$ typeeq p$ 1453 >>; 1454 drop_pde(p,nil,nil)$ 1455 newpdes:=delete(p,newpdes); 1456 1457 % update of mpli: 1458 mpli:=subst(newp,p,mpli)$ 1459 1460 if tr_drop then << 1461 write"mpli=",mpli$terpri()$ 1462 >>; 1463 1464 >>; % end of successful find 1465 cp:=conli; 1466 for s:=1:h do cp:=cdr cp; 1467 rplaca(cp,reval {'plus,-1,car cp}); 1468 >> % if the 1-term PDE is not already known 1469 >>$ % for each possible single term 1470 1471 >>; 1472 return if success then list(newpdes,cadr arglist) 1473 else nil 1474end$ 1475 1476endmodule$ 1477 1478end$ 1479 1480----------------------------------------------------------------------------------- 1481Possible Improvements: 1482 1483 - auch subtrahieren, wenn 0 Gewinn (Zyklus!) 1484 - kann Zyklus mit decoupling geben 1485 - Erweiterung auf mehrere Gleichungen 1486 - counter example: 1487 0 = +a+b+c 1488 0 = -b +d+e 1489 0 = -c-d +f 1490 0 = -a -e-f 1491 combining any 2 gives a longer one, summing 3 an equally long one and 1492 the sum of all 4 is even zero. 1493 - In order not to run into a cycle with decouple one could use 1494 dec_hist_list but that costs memory. 1495 --------- 1496 1497 DONE 1498 - when a ftem_ function becomes non-zero then drop all equations from the rl_with lists 1499 which involve this function 1500 1501 DONE: 1502 - Another problem with allowing some non-zero multipliers and others not: 1503 0 = a + b + c + d + g 1504 0 = e*a + e*b + e*c + f*d 1505 Here the SHORTER one could and should be replaced by 0 = f*d - e*d - e*g using 1506 the longer first one. 1507 1508 - Parallelize the shortening 1509 - It could be beneficial to produce and ADD new equations if they are short, eg 1510 0 = a*g + a*h + b*c \____ b**2*c - a**2*c = c*(a-b)*(a+b) = 0 1511 0 = b*g + b*h + a*c / 1512 Disadvantage: Checking ALL quotients takes much longer. 1513 - If a new equation is slightly longer than the one to be replaced (de2) then one 1514 could ADD this new equation and mark it as redundant (and check whether it has 1515 any useful properties, e.g. factorizability, involving only few ftem_, or being 1516 an ODE when de1,de2 were PDEs. 1517 - perhaps return in shorten() and shorten_pdes() only a single new and to delete PDE, 1518 not lists with each only one element 1519 - One probably could speed up the shortening of homogeneous polynomials 1520 if the degree of both polynomials are equal. 1521 This would need a modification of partition_3 and inputting all 'fcts as 1522 2nd parameter. 1523 - test adding some extra length allowance to the new generated equation 1524 - try take_first=t 1525 - stop when take_first=nil and max_reduction = 2*min(n1,n2) is achieved 1526 1527 DONE: 1528 - When one equation has no flin_ and the other has then treat both as having 1529 no flin_. 1530 1531 - How can the knowledge of factors be used for shortening, especially if they 1532 have more than one term? 1533 E.g. if both equations have common factors, or one is factorizable and the other 1534 not, or both are factorizable. 1535 - Write procedures that plot the size of equations that were shortened over the step number, 1536 and `number of terms in all equations' over 'length interval of equations' 1537 - produce plots for different degrees of overdetermination. 1538 1539 - About rounding: 1540 "also im lisp mode macht quotient einfach division ohne rest 1541 und Ruecksicht auf Umstehende.Da ist ceiling auch simpel 1542 gestrickt.Das macht das LISP system. 1543 1544 Im Algebraic mode ist das anders. Da machen quotient 1545 und ceiling schon mehr. 1546 1547 Die Wandlung machst Du am einfachsten mit 1548 1549 simpquot '(quotient 4 5); 1550 1551 bzw round!* falls Du die nackte float haben willst. 1552 1553 round!* '(!:rd!: . 0.8); 1554 " 1555 - The memory to store which equation has already been length reduced wrt to which 1556 equation grows quadratically with the number of equations and becomes the more 1557 serious the more equations one has. Perhaps one could delete all rl_with entries 1558 of all equations coming in the list of equations before largest_fully_shortened. 1559 This has the consequence that one has to look up the rl_with information at the 1560 larger of both equations. But then, why not just storing this information only 1561 for the larger of both equations? 1562 1563 DONE: 1564 - When shortening a single equation wrt all others, then consider also longer equations 1565 which are less than 2 times as long as the equation to be shortened. 1566 1567 - Because the order of pairing equations for shortening does depend on the number of terms 1568 and the place of the equation in the list pdes of equations, it should be fixed some 1569 time that when multi-term factors are dropped and the number of terms of the equation 1570 is dropped than they should get a new place in pdes, moving more to the start. 1571 1572 - load trysub22 1573then do 1574 1575 s 20 l 11 10 30 30 30 11 11 11 11 11 11 11 11 l 11 10 td 30) 1576 1577and get: 1578 1579next: tr_decouple is now on. 1580next: 30 1581 1582Step 41352: 1583*** Garbage collection starting 1584*** GC 679: 13-Jan-2008 02:09:54 (~ 11358250 ms cpu time, gc : 3 %) 1585*** time 1340 ms, 56164584 occupied, 318491626 recovered, 318491666 free 1586 1587 first pde e_1659: 6 terms, 28 factors, with derivatives 1588 of functions of all variables: 1589 1590 3 2 2 7 1591 (p9 ,p9 ,p9,p24 ,p3 ,p3) 1592 1593is not differentiated, second pde e_837: 9 terms, 44 factors, with derivatives 1594 of functions of all variables: 1595 1596 3 2 2 8 5 2 1597 (p9 ,p9 ,p9,p24 ,p3 ,p3 ,p3 ) 1598 1599is not differentiated, to eliminate 1600p9 1601 1602 1603e_1659 (resp its derivative) is multiplied with 1604***** Non-numeric argument in arithmetic 1605 1606 - A serious problem is that numerical factors may grow too large. 1607 Maybe this is unavoidable. One should have a constant for the 1608 largest numerical factor so that one can change this easily. 1609 1610 - Is there a bug when equations are shortened by 0 terms? 1611 Shall one ignore it, or use it for replacing the previous equation or add it? 1612 1613 - Is it useful to add many short equations in order to reduce larger ones better? 1614 1615 - maybe one should inspect new equations and stop shortening them it they have 1616 too long numerical coefficients? That would be better than trying to shorten 1617 this equation when there is not a chance anyway? 1618 1619 - Or should one have a counter of how often a quotient is disregarded because of too 1620 large numerical coefficients and when this counter gets too large, then stop 1621 shortening this equation? 1622 1623 - Larger equations are reduced for a while but then not anymore because their 1624 Num factors grow too large. But then they do not shrink and can not become 1625 very short to shorten other and/or to provide a breakthrough with very short 1626 equations. 1627 1628 - When shortening a single equation, ak how often maximally. 1629 1630----------------------------------------------------------------------------------- 1631alg_length_reduction % interface to crack 1632 err_catch_short % wrap up 1633 mkeqSQ 1634 shorten_pdes % find a pair of equations 1635 shorten % shorten two pdes given by names with each other 1636 % returns a dotted pair, where car is a list of the values 1637 % of new pdes (only one) and cdr is a list of names of pdes 1638 % to be dropped (also only one) 1639 sort_partition % in crint.red 1640 partition_1 1641 partition_2 1642 short 1643 clean_den(qc,j)$ 1644 clean_num(qc,j)$ 1645 1646 1647 1648alg_length_reduce_1 % Do one length-reducing combination of a single equation wrt all others 1649 1650 alg_length_reduction 1651 1652is_algebraic % very short test 1653 1654 1655drop_lin_dep(arglist)$ % module 12 1656% drops linear dependent equations 1657 Calls: drop_pde neq reval solveeval sort_partition !~prettyprint 1658 1659find_1_term_eqn(arglist)$ % module 13 1660% checks whether a linear combination of the equations can produce 1661% an equation with only one term 1662 Calls: aeval aeval!* assgnpri drop_pde mkeq neq reval 1663 solveeval sort_partition typeeq union !~prettyprint 1664 1665 1666tr alg_length_reduction 1667tr err_catch_short 1668tr mkeqSQ 1669tr shorten_pdes 1670tr shorten 1671tr sort_partition 1672tr short 1673tr clean_den 1674tr clean_num 1675 1676