1%******************************************************************** 2module simplifications$ 3%******************************************************************** 4% Routines for simplifications, contradiction testing 5% and substitution of functions 6% Author: Andreas Brand 1991 1993 1994 7% Thomas Wolf since 1996 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 signchange(g)$ % #### to be dropped? still used in crint 34% ensure, that the first term is positive 35if pairp g then 36 if (car g='minus) then cadr g 37 else if (car g='plus) and (pairp cadr g) and (caadr g='minus) 38 then reval list('minus,g) 39 else g 40else g$ 41 42%------------------------------- 43 44symbolic procedure raise_contradiction(g,text)$ 45<<contradiction_:=t$ 46 if print_ then 47 <<terpri()$if text then write text 48 else write "contradiction : "$ 49 deprint list g>> >>$ 50 51%---------------------------- 52 53symbolic procedure fcteval(p)$ 54% looks for a function which can be eliminated 55% if one is found, it is stored in p as (coeff_of_f . f) 56% 'to_eval neq nil iff not checked yet for substitution 57% or if subst. possible 58% i.e. 'to_eval=nil if checked and not possible 59% 'fcteval_lin includes subst. with coefficients that do not 60% include ftem functions/constants 61% 'fcteval_nca includes subst. with non-vanishing coefficients 62% and therefore no case distinctions (linearity) 63% 'fcteval_nli includes subst. with possibly vanishing coefficients 64% and therefore case distinctions (non-linearity) 65% 'fcteval_n2l =t if a function not in flin_ is to be replaced and 66% the equation contains allvarfcts in flin_, otherwise nil 67 68if flagp(p,'to_eval) then 69if pairp get(p,'fac) then remflag1(p,'to_eval) % factorizable equations should 70 % not be used for substitutions 71 else 72begin scalar ft,a,fl,li,nc,nl,f,cpf,fv,fc,f_is_in_flin,flin_sub_found$ 73 74 if (not get(p,'fcteval_lin)) and 75 (not get(p,'fcteval_nca)) and 76 (not get(p,'fcteval_nli)) then << % otherwise already determined 77 ft:=get(p,'allvarfcts)$ 78 % if flin_ and (fl:=intersection(ft,flin_)) then <<ft:=fl; fl:=nil>>$ 79 % This strict condition never to substitute a non-flin_ function by 80 % an expression including flin_ functions is eased below. 81 if null ft then << 82 ft:=setdiff(get(p,'rational),get(p,'not_to_eval))$ 83 % i.e. for example, no functions that replace a derivative 84 85 % drop all functions f from ft for which there is another 86 % function which is a function of all variables of f + at 87 % least one extra variable 88 if vl_ then << 89 for each f in ft do << 90 cpf:=get(p,'fcts)$ 91 fv:=fctargs f$ 92 while cpf and 93 (not_included(fv,fc:=fctargs car cpf) or 94 (length fv >= length fc) ) do 95 cpf:=cdr cpf; 96 if null cpf then fl:=cons(f,cpf) 97 >>$ 98 ft:=fl 99 >> 100 >>$ 101 li:=nc:=nl:=nil$ 102 for each f in ft do 103 if null member(f,get(p,'not_to_eval)) then % i.e. not a function just 104 % introduced through this equation 105 if alg_linear_fct(p,f) then % i.e. only linear algebr. fcts 106 if confirm_subst or (null flin_) or 107 <<f_is_in_flin:=member(f,flin_)$ 108 if flin_sub_found and null f_is_in_flin then nil 109 else << 110 % either f is in flin_ or no fnc in flin_ is found yet 111 if f_is_in_flin and null flin_sub_found then << 112 % f is the first fct in flin_ that is found 113 li:=nc:=nl:=nil$ 114 % That means that a case generating substitution of a 115 % flin_ function has a higher priority than a non-case 116 % generating substitution of a non-flin_ function in 117 % terms of flin_ functions 118 flin_sub_found:=t 119 >>$ 120 t 121 >> 122 >> then 123 <<a:=coeffn({'!*sq,get(p,'sqval),t},f,1)$ 124 a:=if pairp a and (car a='!*sq) then cadr a 125 else simp a; 126 if null (fl:=smemberl(delete(f,get(p,'fcts)),a)) then 127 li:=cons(cons(a,f),li) else 128 if can_not_become_zeroSQ(a,fl) then nc:=cons(cons(a,f),nc) 129 else nl:=cons(cons(a,f),nl) 130 >>$ 131 132 if flin_ and null flin_sub_found and vl_ and % the 'and vl_' test is 133 % made to avoid the more time expensive intersection test if possible 134 (fl:=intersection(ft,flin_)) then << 135 136 % We can not get to here for purely algebraic problems 137 % because a case generating substitution of a flin_function 138 % would have a higher priority than a substitution of a 139 % non-flin_ function. 140 141 % At this point it is not clear whether this substitution will be selected 142 % in get_subst, so it is not investigated now whether f is not actually 143 % occuring only linearly wrt. flin_. Also, such an investigation would 144 % potentially be repeated many times for the same function if done here. 145 % --> For informing get_subst() that this substitution is potentially 146 % replacing a non-flin_ function by an flin_ function : 147 put(p,'fcteval_n2l,t) 148 149 >>$ 150 if li then put(p,'fcteval_lin,reverse li); % else 151 if nc then put(p,'fcteval_nca,reverse nc); % else 152 if nl then put(p,'fcteval_nli,reverse nl); 153 if not (li or nc or nl) then remflag1(p,'to_eval) 154 >>$ 155 156 return (get(p,'fcteval_lin) or 157 get(p,'fcteval_nca) or 158 get(p,'fcteval_nli) ) 159end$ 160 161%---------------------------- 162 163symbolic procedure best_mdu(p)$ 164if get(p,'fcteval_lin) then 1 else 165if get(p,'fcteval_nca) then 2 else 166if get(p,'fct_nli_lin) then 3 else 167if get(p,'fct_nli_nca) then 4 else 168if get(p,'fct_nli_nli) then 5 else 169if get(p,'fct_nli_nus) then 6 else 7$ 170 171%---------------------------- 172 173symbolic procedure get_subst(pdes,l,length_limit,less_vars,no_df,no_cases)$ 174% --> new version to minimize no of terms substituted, old version to minimize 175% no of cases see at end 176% get the most simple pde from l which leads to a function substitution 177% if less_vars neq nil: the expr. to subst. has only fcts. of fewer vars 178% if no_df neq nil: the expr. to subst. has no derivatives 179 180% Requirements for substitutions are potentially contradicting: 181% - not to generate cases 182% - not to multiply with extra factors, i.e. not to have no-number coefficients 183% - to replace by as few terms as possible 184% - not to replace a function that is known to be non-vanishing 185% - not to replace a non-flin_ function/constant by an expression 186% involving flin_ functions 187% - to use equations involving the highest number of independent variables 188begin scalar p,q,h,l1,l2,m,ntms,mdu,ineq_cp,rtn,lcop,fcteval_cop,necount, 189 found_non_ineq_sub,found_non_n2l_sub,is_non_ineq_sub,is_non_n2l_sub$ 190 % mdu=(1:lin, 2:nca, 3:nli_lin, 4:nli_nca, 5:nli_nli, 6:nli_nus) 191 lcop:=l; 192 % next: if less_vars then pick only equations with one allvarfcts 193 % or no independent variables 194 if less_vars then << 195 while l do << 196 if get(car l,'nvars)=0 then l1:=cons(car l,l1) 197 else << 198 l2:=get(car l,'allvarfcts)$ 199 if l2 and null cdr l2 then l1:=cons(car l,l1)$ 200 >>; 201 l:=cdr l 202 >>$ 203 l:=reverse l1 204 >>$ 205 206 % next: drop all equations longer than length_limit 207 if length_limit then << 208 l1:=nil; 209 while l do 210 if (h:=get(car l,'length))>length_limit then 211 if h>2*length_limit then l:=nil 212 else l:=cdr l 213 214 else << 215 l1:=cons(car l,l1)$ 216 l:=cdr l 217 >>$ 218 l:=reverse l1 219 >>$ 220 % l is now the list of equations <= length_limit 221 222 % next: substitution only if no_df=nil or 223 % no derivative of any function occurs 224 if no_df then << 225 l1:=nil; 226 for each s in l do << 227 l2:=get(s,'derivs)$ 228 while l2 do 229 if pairp(cdaar l2) then 230 <<l2:=nil; 231 l1:=cons(s,l1)>> 232 else l2:=cdr l2>>$ 233 l:=setdiff(l,l1) 234 >>$ 235 236 % next: initialize 'fcteval_lin,... 237 % if no_cases then restrict to substitutions that don't lead to cases. 238 239 l1:=nil; 240 necount:=0; 241 for each s in l do 242 if fcteval(s) then 243 if get(s,'fcteval_lin) or get(s,'fcteval_nca) then l1:=cons(s,l1) else 244 if null no_cases then 245 if (h:=get(s,'fcteval_nli)) then << 246 if (null get(s,'fct_nli_lin)) and 247 (null get(s,'fct_nli_nca)) and 248 (null get(s,'fct_nli_nli)) and 249 (null get(s,'fct_nli_nus)) then << 250 ineq_cp:=ineq_; ineq_:=nil; 251 % partition of get(s,'fcteval_nli) into the above 4 cases 252 for each l2 in h do << 253 q:=mkeqSQ(car l2,nil,nil,get(s,'fcts),get(s,'vars),allflags_, 254 t,list(0),nil,nil); 255 % the pdes-argument in mkeqSQ() is nil to avoid lasting effect on pdes 256 necount:=add1 necount$ 257 fcteval(q)$ % less_vars 258 if get(q,'fcteval_lin) then 259 put(s,'fct_nli_lin,cons(l2,get(s,'fct_nli_lin))) else 260 if get(q,'fcteval_nca) then 261 put(s,'fct_nli_nca,cons(l2,get(s,'fct_nli_nca))) else 262 if get(q,'fcteval_nli) then 263 put(s,'fct_nli_nli,cons(l2,get(s,'fct_nli_nli))) else 264 put(s,'fct_nli_nus,cons(l2,get(s,'fct_nli_nus)))$ 265 drop_pde(q,nil,nil)$ 266 if necount>100 then << 267 clean_prop_list(pdes)$ 268 necount:=0 269 >> 270 >>$ 271 ineq_:=ineq_cp 272 >>$ 273 l1:=cons(s,l1)$ % Before the number of splitting was minimized. Now the 274 % number of indep. variables is max.ed and length min.ed. 275 >>$ 276 l:=l1$ 277 278 % next: find an equation with the following priority list: 279 % - as many as possible variables 280 % - few as possible terms for substitution 281 % - generating as few as possible cases 282 % - not substituting a non-flin_ function by flin_ functions 283 m:=-1$ mdu:=8$ 284 for each s in l do << 285 l1:=get(s,'nvars); 286 if get(s,'starde) then l1:=sub1 l1; 287 if l1>m then m:=l1$ 288 >>$ 289 290 while m>=0 do << 291 l1:=l$ 292 ntms:=10000000$ 293 found_non_ineq_sub:=nil$ % whether a possibly vanishing function can be substituted 294 found_non_n2l_sub :=nil$ % whether it can be avoided to substitute a non-flin_ 295 % function by an expression involving flin_ functions 296 while l1 do % check each suitable equation car l1 297 if ((get(car l1,'nvars) - 298 if get(car l1,'starde) then 1 299 else 0) = m ) and 300 (q:=fcteval(car l1)) and 301 ((is_non_ineq_sub:=null member(simp cdar q,ineq_)) or 302 null found_non_ineq_sub 303 % Try to find a substitution which does not replace a 304 % non-vanishing function. 305 % if one is found then is_non_ineq_sub <> nil (good) 306 % if none is found then is_non_ineq_sub = nil ( bad) 307 % Only the first of the possible substitutions of this equation is 308 % checked because non-vanishing functions have a low priority and 309 % should coming only later in q. 310 ) and 311 ((is_non_n2l_sub:=null get(car l1,'fcteval_n2l)) or 312 null found_non_n2l_sub 313 % Try to find a substitution which does not replace a 314 % non-flin_ function by an expression invloving allvarfcts 315 % flin_ functions 316 % if one is found then is_non_n2l_sub <> nil (good) 317 % if none is found then is_non_n2l_sub = nil ( bad) 318 ) and 319 ((get(car l1,'terms) < ntms ) or 320 ((get(car l1,'terms) = ntms) and 321 (best_mdu(car l1) < mdu ) ) ) then 322 <<p:=car l1$ l1:=cdr l1$ 323 ntms:=get(p,'terms)$ 324 mdu:=best_mdu(p)$ 325 if null found_non_ineq_sub then found_non_ineq_sub:=is_non_ineq_sub$ 326 if null found_non_n2l_sub then found_non_n2l_sub:=is_non_n2l_sub 327 >> else l1:=cdr l1$ 328 m:=if p then -1 329 else sub1 m 330 >>$ 331 332 if p then return << 333 334 fcteval_cop:=if mdu=1 then get(p,'fcteval_lin) else 335 if mdu=2 then get(p,'fcteval_nca) else 336 if mdu=3 then get(p,'fct_nli_lin) else 337 if mdu=4 then get(p,'fct_nli_nca) else 338 if mdu=5 then get(p,'fct_nli_nli) else 339 if mdu=6 then get(p,'fct_nli_nus); 340 rtn:={mdu,p,pick_fcteval(pdes,mdu,fcteval_cop)}; 341 % prevent the substitution of a function<>0 342 if rtn and fhom_ and 343 setdiff(for each h in get(p,'fcts) collect simp h,ineq_) and 344 cdr pdes and (get(p,'terms)>1) then << 345 % i.e. not all ftem_ have to be non-zero 346 % and it is not the last pde 347 348 % n0f:=setdiff(ineq_,setdiff(ineq_,for each h in ftem_ collect simp h)); 349 % if freeof(n0f,cdaddr rtn) then rtn 350 if not member(simp cdaddr rtn,ineq_) then rtn 351 else 352 if null cdr fcteval_cop then % rtn was the only substitution of this eqn. 353 if cdr lcop then % there are other eqn.s to choose from 354 <<h:=get_subst(pdes,delete(p,lcop),length_limit,less_vars,no_df,no_cases); 355 if null h then rtn else h 356 >> else rtn % nil % no substitution --> changed to rtn 357 else << 358 fcteval_cop:=delete(caddr rtn,fcteval_cop); 359 {mdu,p,pick_fcteval(pdes,mdu,fcteval_cop)} 360 >> 361 >> else rtn 362 >> 363end$ 364 365symbolic procedure total_deg_of_first_tm(sf)$ 366if domainp sf then 0 else 367ldeg sf + total_deg_of_first_tm lc sf$ 368 369symbolic procedure nco_SF(sf)$ 370% returns the numerical coefficient of the first term 371<< while pairp sf and not domainp sf and not domainp lc sf do sf:=lc sf; sf>>$ 372 373symbolic procedure pick_fcteval(pdes,mdu,fctli)$ 374if fctli then 375if (not expert_mode) or (length fctli = 1) then 376% automatic pick of all the possible substitutions 377if null cdr fctli then car fctli else 378if mdu<3 then begin % substitute the function coming first in ftem_ 379 scalar best, besth1, besth2, besth3, h1, h2, h3$ 380 381 besth1:=100000000$ 382 while fctli do << 383 if besth1 > (h1:=no_of_tm_sf numr caar fctli) then << 384 best:=car fctli; 385 besth1:=h1; 386 if h1=1 then << 387 besth2:=total_deg_of_first_tm numr caar fctli$ 388 besth3:=nco_SF numr caar fctli 389 >> 390 >> else 391 if besth1=h1 then % else car fctli is not better than best 392 if h1>1 then % besth2 and besth3 are not determined 393 if which_first(cdr best,cdar fctli,ftem_) neq cdr best then best:=car fctli else 394 else % besth1=h1=1, then besth2, besth3 are determined 395 if besth2>(h2:=total_deg_of_first_tm numr caar fctli) then << 396 best:=car fctli; 397 besth2:=h2; 398 besth3:=nco_SF numr caar fctli 399 >> else 400 if besth2=h2 then % else car fctli is not better than best 401 if (fixp(h3:=nco_SF numr caar fctli) and not fixp besth3) or 402 (fixp h3 and fixp besth3 and (abs h3 < abs besth3) ) then << 403 best:=car fctli; 404 besth3:=h3 405 >>; 406 fctli:=cdr fctli 407 >>; 408 return best 409end else 410begin scalar co,minfinco,minnofinco,finco,nofinco,fctlilen, 411 n,maxnopdes,nopdes,f,bestn$ 412 % 1. find a substitution where the coefficient involves as few as possible functions 413 fctlilen:=length fctli$ 414 minnofinco:=10000$ 415 for n:=1:fctlilen do << 416 co:=nth(fctli,n)$ 417 finco:=smemberl(ftem_,car co); 418 nofinco:=length finco; 419 if nofinco<minnofinco then <<minfinco:=list(cons(n,finco)); 420 minnofinco:=nofinco>> else 421 if nofinco=minnofinco then minfinco:=cons(cons(n,finco),minfinco); 422 >>$ 423 if (length minfinco=1) or (minnofinco>1) 424 % if there is only one substitution where the coefficient has a 425 % minimal number of ftem_ functions or 426 % if the minimal number of functions in any coefficient is >1 427 then return nth(fctli,caar minfinco) % return any ony one of the minimal ones 428 else return << % find the one with the ftem_ function that occurs in the 429 % fewest equations, to complicate as few as possible equations 430 maxnopdes:=1000000; 431 for each su in minfinco do << 432 f:=cadr su; 433 nopdes:=0; 434 for each p in pdes do if not freeof(get(p,'fcts),f) then nopdes:=add1 nopdes; 435 if nopdes<maxnopdes then <<maxnopdes:=nopdes;bestn:=car su>>; 436 >>$ 437 nth(fctli,bestn) 438 >> 439end else 440begin scalar fl,a,h; 441 fl:=for each a in fctli collect cdr a$ 442 write"Choose a function to be substituted from "$ 443 listprint(fl)$terpri()$ 444 change_prompt_to ""$ 445 repeat h:=termread() until not freeof(fl,h); 446 restore_interactive_prompt()$ 447 while h neq cdar fctli do fctli:=cdr fctli; 448 return car fctli 449end$ 450 451symbolic procedure do_one_subst(ex,f,a,ftem,vl,level,eqn,pdes)$ 452% substitute ex for f in a 453% ex is in {'!*sq,SQ,t} form, f is a kernel, a is in SQ form 454begin scalar l,l1,p,oldstarde$ 455 % l is the list of factors of a 456 l:=get(a,'fac)$ 457 if not pairp l then l:={get(a,'sqval)}$ 458 oldstarde:=get(a,'starde)$ 459 while l do << % for each factor 460 if smember(f,car l) then << 461 % p:=subsq(car l,{(f . ex)})$ 462 p:=simp!* {'!*sq,subsq(car l,{(f . ex)}),nil}$ 463 % - simp!* {'!*sq,..,nil} to simplify poly using identities, like i^2 -> -1 464 if sqzerop p then <<l:=list nil$l1:=list 0>> 465 else << 466 l1:=cons(p,l1) 467 >> 468 >> else l1:=cons(car l,l1)$ 469 l:=cdr l 470 >>$ 471 l:=nil$ 472 while l1 do << 473 if not member(car l1,cdr l1) then 474 l:=union(simplifySQ(car l1,ftem,t,nil,nil),l)$ 475 l1:=cdr l1 476 >>$ 477 l:=delete((1 . 1),l)$ % delete all non-vanishing factors 478 if null l then << 479 if print_ then << 480 terpri()$ % new 481 write"Substitution of "$ 482 fctprint list f$ 483 if cdr get(eqn,'fcts) then << 484 write " by an expression in "$terpri()$ 485 fctprint delete(f,get(eqn,'fcts)) 486 >>$ 487 write " found in ",eqn," : "$ 488 eqprint(list('equal,f,ex)) 489 >>$ 490 raise_contradiction({'!*sq,get(a,'sqval),t}, 491 "leads to a contradiction in : ")$ 492 a:=nil 493 >> else 494 if member((nil . 1),l) then << % one factor is zero 495 drop_pde(a,nil,0)$ 496 a:=nil 497 >> else << 498% if pairp cdr l then l:=cons('times,l) 499% else l:=car l$ 500 if get(a,'level) neq level then << 501 pdes:=drop_pde(a,pdes,0)$ % pdes updated as it is used in mkeqSQ: 502 a:=mkeqSQ(nil,l,nil,ftem,vl,allflags_,nil,list(0),nil,pdes) 503 >> else << 504 p:=get(a,'derivs); % keep the leading derivative to check 505 if p then p:=caar p; % whether it changed in the substitution 506 for each b in allflags_ do flag(list a,b)$ 507 if null updateSQ(a,nil,l,nil,ftem,vl,nil,list(0),pdes) then << 508 % l is a list of sq-form factors 509 drop_pde(a,pdes,0)$ 510 a:=nil 511 >> else << 512 % If the leading derivative has changed then drop 513 % the 'dec_with and the 'dec_with_rl list. 514 l1:=get(a,'derivs); 515 if l1 then l1:=caar l1; 516 if l1 neq p then << 517 put(a,'dec_with,nil); 518 put(a,'dec_with_rl,nil); 519 for each l1 in pdes do if l1 neq a then << 520 drop_dec_with(a,l1,t)$ 521 drop_dec_with(a,l1,nil) 522 >> 523 >>; 524 drop_pde_from_idties(a,pdes,nil)$ 525 % nil as second argument for safety, for not knowing better 526 put(a,'rl_with,nil); 527 for each l1 in pdes do if l1 neq a then drop_rl_with(a,l1) 528 >> 529 >>$ 530 put(a,'level,level) 531 >>$ 532 if oldstarde and not get(a,'starde) then put(a,'dec_with,nil); 533 return a$ 534end$ 535 536symbolic procedure sub_in_forg(f,ex,forg); 537% f is the function to be substituted in forg 538% ex has form {'!*sq,..,t} 539begin scalar fl,h,was_subst,dnr$ 540 % fl:=delete(f,smemberl(ftem_,ex))$ % functs occur. in ex, delete should not be necess. 541 fl:=smemberl(append(ftem_,for each h in fsub_ collect car h),ex)$ % functions occuring in ex 542 forg:=for each h in forg collect 543 if atom h then 544 if f=h then <<put(f,'fcts,fl)$was_subst:=t$ 545 list('equal,f,if pairp ex and (car ex='!*sq) then cadr ex 546 else simp!* ex )>> 547 else h 548 else 549 if (car h='equal) and member(f,get(cadr h,'fcts)) then << 550 was_subst:=t$ 551 dnr:=simp!* {'!*sq,subf(denr caddr h,{(f . ex)}),nil}; 552 if sqzerop dnr then <<contradiction_:=t$ 553 terpri()$write"##### ERROR: When substituting ",f," by ",ex, 554 " in the denominator of the forg entry: ",h, 555 " then the denominator becomes zero!! #####"$ 556 terpri()$ 557 nil 558 >> else << 559 h:=list('equal,cadr h, 560 simp!* {'!*sq,quotsq(subf(numr caddr h,{(f . ex)}),dnr),nil}); 561 % h:=list('equal,cadr h,simp!* {'!*sq,subsq(caddr h,{(f . ex)}),nil}); 562 % - simp!* {'!*sq,..,nil} to simplify poly using identities, like i^2 -> -1 563 put(cadr h,'fcts, 564 smemberl(union(fl,delete(f,get(cadr h,'fcts))),caddr h)); 565 h 566 >> 567 >> else h$ 568 return (forg . was_subst) 569end$ 570 571symbolic procedure do_subst(md,p,l,pde,forg,vl,plim,keep_eqn)$ 572% md is the mode of substitution (1='fcteval_lin,..., 6='fct_nli_nus) 573% (md is needed in case of an ISE) 574% Use equation p for substitution. 575% l incodes the substitution itself = (cf . f) 576% Substitute a function in all pdes. 577begin scalar f,fl,h,ex,res,slim,too_large,was_subst,new_user_rules, 578 ruli,ise,cf,vl,nof,stde,partial_subs,ineq_bak,ineq_or_bak$ 579% l:=get(p,'fcteval_lin)$ 580% if null l then l:=get(p,'fcteval_nca)$ 581% if null l then l:=get(p,'fcteval_nli)$ 582% if l then << % l:=car l$ 583 f:=cdr l$ 584 cf:=car l$ 585 586 % at first check whether f can be added to flin_ 587 if flin_ and vl_ 588 % - We can not do a test 'if get(p,'fcteval_n2l) then ' 589 % because fcteval_n2l does only know about allvarfct functions in flin_ 590 % - We test 'and vl_' because this is much cheaper than freeof and 591 % currently (Dec 2010) non-linear substitutions have a lower priority 592 % than any flin_ substitutions, even if they are case generating. 593 and freeof(flin_,f) and null add2flin(pde,f) then << 594 % f is not occuring linearly wrt. flin_ therefore remove all functions 595 % of this equation from flin_ 596 for each h in get(p,'fcts) do 597 if not freeof(flin_,h) then flin_:=delete(h,flin_) 598 >>$ 599 600 if get(p,'starde) then ise:=t; 601 slim:=get(p,'length)$ 602 ruli:=start_let_rules()$ 603 if modular_comp then on modular$ 604 ex:={'!*sq,subs2 quotsq(subtrsq(multsq(cf,simp f),get(p,'sqval)),cf),t}; 605 % ex:={'!*sq,simp!* {'!*sq,quotsq(subtrsq(multsq(cf,simp f),get(p,'sqval)),cf),nil},t}; 606 % - simp!* {'!*sq,..,nil} to simplify poly using identities, like i^2 -> -1 607 608 % When the right hand side of a substitution becomes identical to the left 609 % hand side due to LET rules then there is no point in performing the 610 % substitution. But the equation used to formulate the substitution 611 % is deleted because it is used for the substitution so all that is left is 612 % the set of LET rules. :-( 613 614%####################### 615% CHECK WHETHER f==ex and then drop all LET RULES and add them as equations. 616%####################### 617 618 % This substitution might give a contradiction when being done to an 619 % equation. But if this equation is used as a LET rule then that 620 % contradicting substitution may turn a denominator of in forg to zero when 621 % the contradicted LET rule is used in the substitution in forg. To avoid 622 % the error message that a denominator in forg is zero, in the following 623 % it is checked whether a LET rule is contradicted by the substitution. 624 625 for each h in cdr userrules_ do % h = {replaceby,cadr h,caddr h}, cadr h => caddr h 626 if not freeof(cdr h,f) and 627 ({(1 . 1)} = simplifySQ(subsq(subtrsq(simp cadr h,simp caddr h), 628 {(f . ex)}), ftem_, t, nil, nil)) then << 629 contradiction_:=t$ 630 if print_ then << 631 write"The current substitution"$terpri()$ 632 algebraic(write lisp f," = ", lisp ex)$ 633 write"results in a contradiction in the LET rule:"$ 634 algebraic (write lisp {'list,h})$ 635 >> 636 >>$ 637 638 %---- specification of substitution in case of expert_mode (user guided) 639 if not contradiction_ then 640 if expert_mode then << 641 terpri()$ 642 write"Enter a list of equations in which substitution should take place."$ 643 terpri()$ 644 write"Substitution into the expressions for the original functions and"$ 645 terpri()$ 646 write"the inequalities is only done if you select all equations with `;' ."$ 647 l:=select_from_list(pde,nil)$ 648 if l then << 649 if not_included(pde,l) then partial_subs:=t 650 else partial_subs:=nil; 651 l:=delete(p,l) 652 >>; 653 if partial_subs then << 654 change_prompt_to ""$ 655 terpri()$ 656 write"Should substitutions be done in the inequalities? (y/nil)"$ 657 repeat h:=reval termread() until (h=t) or (h=nil)$ 658 restore_interactive_prompt() 659 >> 660 >> else l:=delete(p,pde)$ 661 662 %---- substitution in LET-rules 663 if not contradiction_ then << 664 for each h in cdr userrules_ do % h = {replaceby,cadr h,caddr h}, cadr h => caddr h 665 if freeof(cadr h,f) then 666 if freeof(caddr h,f) then new_user_rules:=cons(h,new_user_rules) 667 else << 668 if print_ then << 669 terpri()$ 670 write"The current substitution"$terpri()$ 671 algebraic(write lisp f," = ", lisp ex)$ 672 write"affects the following LET rule:"$terpri()$ 673 algebraic (write lisp {'list,h}) 674 >>$ 675 algebraic(clearrules lisp {'list,h}); 676 h:={car h, cadr h, mk!*sq subsq(simp caddr h,{(f . ex)})}; 677 if print_ then << 678 write"which therefore is modified to:"$ 679 algebraic (write lisp {'list,h})$ 680 >>$ 681 new_user_rules:=cons(h,new_user_rules); 682 algebraic(let lisp {'list,h}); 683 >> else << 684 if print_ then << 685 terpri()$ 686 write"The current substitution affects the following LET rule"$terpri()$ 687 write"which therefore has to be deleted:"$ 688 algebraic (write lisp {'list,h}) 689 >>$ 690 pde:=moverule2eqn(h,pde)$ 691 >>$ 692 userrules_:=cons('list,reverse new_user_rules)$ 693 >>$ 694 695 %---- substitution in inequalities 696 if not contradiction_ then << 697 ineq_bak:=ineq_$ % for the case that no substitution is done 698 ineq_or_bak:=ineq_or$ % " 699 if (not ise) and ((not partial_subs) or h) then %ineqsubst(ex,f,ftem_,pde)$ 700 simp_all_ineq_with_subst_SQ(ex,f,pde)$ 701 >>$ 702 703 if not contradiction_ then << 704 705 %--- substitution in forg 706 if (not ise) and (not partial_subs) then 707 if lazy_eval and null lin_problem then << 708 was_subst:=t$ 709 fsub_:=cons((f . ex),fsub_)$ 710 % alternatively: 711 %h:=cons(0,forg)$ 712 %forg:=h$ 713 %while cdr h and f neq cadr forg do h:=cdr h$ 714 %if cdr h then <<% f is found and still unevaluated 715 % replacd(h,cons((f . ex),cddr forg))$ 716 % put(f,'fcts,smemberl(ftem_,ex)) 717 %>> else fsub_:=cons((f . ex),fsub_)$ 718 %forg:=cdr forg$ 719 >> else << 720 was_subst:=sub_in_forg(f,ex,forg); 721 forg:=car was_subst; 722 was_subst:=cdr was_subst 723 >>$ 724 725 % The following test depends on the global structure, taken out 726 % for the time being: 727 %% no substitution in equations which do not include all functions 728 %% of all variables in ex 729 %h:=nil; 730 %fl:=get(p,'allvarfcts); 731 %while l do << 732 % if not_included(fl,get(car l,'fcts)) then too_large:=t 733 % else h:=cons(car l,h); 734 % l:=cdr l 735 %>>; 736 %l:=h; 737 % Do the substitution in all suitable equations 738 739 if ise and not contradiction_ then << % if the equation to be used, ie. p, 740 % is a star-equation then collecting suitable equations for substitution 741 h:=nil; 742 vl:=get(p,'vars)$ 743 fl:=get(p,'fcts)$ 744 nof:=caar get(p,'starde)$ 745 while l do << 746 if (stde:=get(car l,'starde)) and 747 (nof<=caar stde) and 748 (not not_included(vl,get(car l,'vars))) and 749 (not not_included(fl,get(car l,'fcts))) then h:=cons(car l,h); 750 l:=cdr l 751 >>$ 752 l:=h; 753 >>$ 754 755 while l and not contradiction_ do << 756 if member(f,get(car l,'fcts)) then 757 if not expert_mode and plim and (slim*get(car l,'length)>plim) 758 then too_large:=t 759 else << 760 pde:=eqinsert(do_one_subst(ex,f,car l,ftem_,union(vl,get(car l,'vars)), 761 get(p,'level),p,pde),delete(car l,pde))$ 762 for each h in pde do drop_rl_with(car l,h); 763 put(car l,'rl_with,nil); 764 for each h in pde do drop_dec_with(car l,h,'dec_with_rl); 765 put(car l,'dec_with_rl,nil); 766 flag(list car l,'to_int); 767 was_subst:=t 768 >>$ 769 l:=cdr l 770 >>$ 771 if print_ and (not contradiction_) and was_subst then << 772 terpri()$write "Substitution of "$ 773 fctprint list f$ 774 if cdr get(p,'fcts) then << 775 write " by an "$ 776 if ise then write"(separable) "$ 777 write "expression in "$terpri()$ 778 fctprint delete(f,get(p,'fcts)) 779 >>$ 780 write " found in ",p," : "$ 781 eqprint(list('equal,f,ex)) 782 >>$ 783 % To avoid using p repeatedly for substitutions of different 784 % functions in the same equations: 785 if ise and not contradiction_ then << 786 put(p,'fcteval_lin,nil); 787 put(p,'fcteval_nca,nil); 788 put(p,'fcteval_nli,nil); 789 remflag1(p,'to_eval)$ % otherwise 'fcteval_??? would be computed again 790 md:=md; % only in order to do something with md if the next 791 % statement is commented out 792 % if too_large then 793 % if md=1 then put(p,'fcteval_lin,list((cf . f))) else 794 % if md=2 then put(p,'fcteval_nca,list((cf . f))) else 795 % put(p,'fcteval_nli,list((cf . f)))$ 796 % could probably unnecessarily be repeated 797 >>; 798 % delete f and p if not anymore needed 799 if (not ise) and 800 (not keep_eqn) and 801 (not too_large) and 802 (not partial_subs) and 803 (not contradiction_) then << 804 %if not assoc(f,depl_copy_) then << 805 if (null lazy_eval) or lin_problem then << 806 h:=t; 807 for each l in forg do 808 if pairp l then if cadr l=f then h:=nil else 809 else if l=f then h:=nil; 810 if h then drop_fct(f)$ 811 >>$ 812 %>>$ 813 was_subst:=t$ % in the sense that pdes have been updated 814 ftem_:=delete(f,ftem_)$ 815 flin_:=delete(f,flin_)$ 816 pde:=drop_pde(p,pde,0)$ 817 >>$ 818% if was_subst then 819 res:=list(pde,forg,p) 820 % also if not used to delete the pde if the function to be 821 % substituted does not appear anymore 822 >>$ 823 if modular_comp then off modular$ 824 stop_let_rules(ruli)$ 825% >>$ 826 if null was_subst then <<ineq_:=ineq_bak; ineq_or:=ineq_or_bak>>$ 827 if not contradiction_ then return cons(was_subst,res)$ 828end$ 829 830 831symbolic procedure make_subst(pdes,forg,vl,l1,length_limit,pdelimit, 832 less_vars,no_df,no_cases,lin_subst, 833 min_growth,cost_limit,keep_eqn,sub_fc)$ 834% make a subst. 835% l1 is the list of possible "candidates" 836begin scalar p,q,r,l,h,hh,h3,cases_,w,md,tempchng,plim,mod_switched$ % ,ineq,cop,newfdep 837 838 for each h in l1 do 839 if not freeof(pdes,h) then l:=cons(h,l); 840 if l1 and null l then return nil 841 else l1:=reverse l$ 842 l:=nil; 843 if expert_mode then << 844 write"Which PDE should be used for substitution?"$ terpri()$ 845 l1:=selectpdes(pdes,1)$ 846 >>; 847 % a fully specified substitution from to_do_list 848 if sub_fc and % a specific function sub_fc is to be substituted using a 849 % specific equation car l1 850 l1 and null cdr l1 then << 851 p:=car l1$ 852 853 if null fcteval(p) then 854 if print_ then << 855 write"##### Strange: equation ",p," was to be solved for ",sub_fc, 856 " but facteval says that is not possible."$ 857 terpri() 858 >>$ 859 860 h:=get(p,'fcteval_lin); 861 while h and (sub_fc neq cdar h) do h:=cdr h; 862 if h then hh:=1 863 else << 864 h:=get(p,'fcteval_nca); 865 while h and (sub_fc neq cdar h) do h:=cdr h; 866 if h then hh:=2 867 else << 868 h:=get(p,'fcteval_nli); 869 while h and (sub_fc neq cdar h) do h:=cdr h; 870 if h then hh:=3 871 >> 872 >>; 873 if h then w:={hh,car l1,car h} 874 >>; 875 876 if sub_fc and null w then return nil; 877 878again: 879 if (sub_fc and w) or 880 (min_growth and (w:=err_catch_minsub(pdes,l1,cost_limit,no_cases)) ) or 881 ((null min_growth ) and 882 (w:=get_subst(pdes,l1,length_limit,less_vars,no_df,no_cases)) ) then 883 % w has form {mdu,p,(cf . f)} , cf is in SQ form 884 if null !*batch_mode and null expert_mode and confirm_subst and << 885 886 if print_ then << 887 terpri()$ 888 write"Proposal: Substitution of ",cdaddr w$terpri()$ 889 write" using equation ",cadr w,": "$ 890 >>; 891 if print_ and (get(cadr w,'printlength)<=print_) then print_stars(cadr w)$ 892 if print_ then <<typeeq(cadr w)$terpri()>>$ 893 894 if car w<=2 then write"No case distinctions will be necessary." 895 else write"Case distinctions will be necessary."$ 896 terpri()$ 897 write"The coefficient is:"$ mathprint factorize {'!*sq,caaddr w,t}$ 898 write"Accept? (Enter y or n or p for stopping substitution) "$ 899 change_prompt_to ""$ 900 repeat h:=termread() until (h='y) or (h='n) or (h='p); 901 restore_interactive_prompt()$ 902 if h='n then << 903 tempchng:=cons(w,tempchng); 904 if car w=1 then <<hh:=get(cadr w,'fcteval_lin); 905 hh:=delete(caddr w,hh); 906 put(cadr w,'fcteval_lin,hh)>> else 907 if car w=2 then <<hh:=get(cadr w,'fcteval_nca); 908 hh:=delete(caddr w,hh); 909 put(cadr w,'fcteval_nca,hh)>> else 910 <<hh:=get(cadr w,'fcteval_nli); 911 hh:=delete(caddr w,hh); 912 put(cadr w,'fcteval_nli,hh); 913 if car w=3 then <<hh:=get(cadr w,'fct_nli_lin); 914 hh:=delete(caddr w,hh); 915 put(cadr w,'fct_nli_lin,hh) 916 >> else 917 if car w=4 then <<hh:=get(cadr w,'fct_nli_nca); 918 hh:=delete(caddr w,hh); 919 put(cadr w,'fct_nli_nca,hh) 920 >> else 921 if car w=5 then <<hh:=get(cadr w,'fct_nli_nli); 922 hh:=delete(caddr w,hh); 923 put(cadr w,'fct_nli_nli,hh) 924 >> else 925 if car w=6 then <<hh:=get(cadr w,'fct_nli_nus); 926 hh:=delete(caddr w,hh); 927 put(cadr w,'fct_nli_nus,hh) 928 >> 929 >>; 930 if null hh and 931 null get(cadr w,'fcteval_lin) and 932 null get(cadr w,'fcteval_nca) and 933 null get(cadr w,'fcteval_nli) then remflag1(cadr w,'to_eval) 934 % otherwise 'fcteval_lin,... will be reassigned 935 >>; 936 if (h='p) then l1:=nil; 937 if (h='n) or (h='p) then t else nil 938 >> then goto again 939 else 940 if ( car w = 1) or 941 ((lin_subst=nil) and 942 ( (car w = 2) or 943 ((car w > 2) and 944 member(caaddr w,ineq_) ) ) ) then << 945 946 if pdelimit and in_cycle({'subst,stepcounter_,cdaddr w, 947 get(cadr w,'printlength),pdelimit}) 948 % function, printlength of equation 949 then plim:=nil 950 else plim:=pdelimit; 951 952 l:=do_subst(car w,cadr w,caddr w,pdes,forg,vl,plim,keep_eqn)$ 953 if l and null car l then << % not contradiction but not used 954 l1:=delete(cadr w,l1); 955 if l1 then << 956 pdes:=cadr l; 957 forg:=caddr l; 958 l:=nil; 959 goto again 960 >> else l:=nil 961 >>; 962 if l then << 963 l:=cdr l; 964 add_to_last_steps({'subst,stepcounter_,cdaddr w, 965 get(cadr w,'printlength),pdelimit}) 966 >> 967 >> else 968 if (null lin_subst) and (null no_cases) then << 969 md:=car w; % md = type of substitution, needed in case of ISE 970 p:=cadr w; % p = the equation 971 w:=caddr w; % w = (coeff . function) 972 if pdelimit and in_cycle({'subst,stepcounter_,w, 973 get(p,'printlength),pdelimit}) % (eqn,function) 974 then pdelimit:=nil; 975 backup_to_file(pdes,forg,nil)$ 976 977 % make an equation from the coefficient 978 q:=mkeqSQ(car w,nil,nil,get(p,'fcts),get(p,'vars),allflags_, 979 t,list(0),nil,pdes)$ 980 % and an equation from the remainder 981 r:=nil; 982 if not contradiction_ then << 983 hh:=subtrsq(get(p,'sqval),multsq(car w,simp cdr w))$ 984 if not sqzerop hh then 985 r:=mkeqSQ(hh,nil,nil,get(p,'fcts),get(p,'vars), 986 allflags_,t,list(0),nil,pdes) 987 >>; 988 if contradiction_ then << 989 if print_ then << 990 write"Therefore no special investigation whether the "$ 991 terpri()$ 992 write"coefficient of a function to be substituted is zero."$ 993 >>$ 994 contradiction_:=nil$ 995 h:=restore_backup_from_file(pdes,forg,nil)$ 996 pdes:= car h; 997 forg:=cadr h; 998 delete_backup()$ 999 if null q then h:={car w} 1000 else << 1001 h:=get(q,'fac)$ 1002 if not pairp h then h:={get(q,'sqval)} 1003 >>$ 1004 for each l in h do % (was: .. in cdr h) % pdes:= % <=<=<=<= 1005 addSQineq(pdes,l,if q then nil else t)$ % nil if l already simplified 1006 drop_pde(q,nil,nil)$ 1007 if r then drop_pde(r,nil,nil)$ 1008 l:=do_subst(md,p,w,pdes,forg,vl,pdelimit,keep_eqn)$ 1009 1010 if l and null car l then << % not contradiction but not used 1011 l1:=delete(p,l1); 1012 if l1 then << 1013 pdes:=cadr l; 1014 forg:=caddr l; 1015 l:=nil; 1016 goto again 1017 >> 1018 >>; 1019 if l then << 1020 l:=cdr l; 1021 add_to_last_steps({'subst,stepcounter_,cdr w, 1022 get(p,'printlength),pdelimit}) 1023 >> 1024 1025 >> else << 1026 % Generation of two equation does not instantly generate a contradiction 1027 remflag1(p,'to_eval)$ 1028 if print_ then << 1029 terpri()$ 1030 write "for the substitution of ",cdr w," by ",p$ 1031 write " we have to consider the case 0=",q,": "$ 1032 eqprint list('equal,0,{'!*sq,car w,t}) 1033 >>$ 1034 pdes:=eqinsert(q,drop_pde(p,pdes,nil))$ 1035 if freeof(pdes,q) then << 1036 % The coefficient must be zero --> only one case, no substitution 1037 if print_ then << 1038 terpri()$ 1039 write "It turns out that the coefficient of ",cdr w," in ", 1040 p," is zero due"$ 1041 terpri()$ 1042 write "to other equations. Therefore no substitution is made and"$ 1043 terpri()$ 1044 write "equation ",p," will be updated instead."$ 1045 terpri() 1046 >>$ 1047 h:=restore_backup_from_file(pdes,forg,nil)$ 1048 pdes:= car h; 1049 forg:=cadr h; 1050 delete_backup()$ 1051 if modular_comp and null !*modular then <<on modular$ mod_switched:=t>>$ 1052 hh:=subtrsq(get(p,'sqval),multsq(car w,simp cdr w))$ 1053 if mod_switched then off modular$ 1054 if sqzerop hh then <<drop_pde(p,pdes,nil); pdes:=delete(p,pdes)>> 1055 else << 1056 updateSQ(p,hh,nil,nil,get(p,'fcts),get(p,'vars),t,list(0),pdes)$ 1057 drop_pde_from_idties(p,pdes,nil)$ % new history is nil as r has no history 1058 drop_pde_from_properties(p,pdes) 1059 >>$ 1060 drop_pde(q,pdes,nil); % q is not in pdes but nevertheless 1061 if r then drop_pde(r,pdes,nil); % r is not in pdes but nevertheless 1062 l:=list(pdes,forg,p) 1063 >> else << 1064 % Setting the coefficient to zero does not give a contradiction 1065 % and is not an obvious consequence of other equations 1066 if r then pdes:=eqinsert(r,pdes)$ 1067 if print_ then << 1068 write"The coefficient to be set = 0 in the first subcase is:"$ 1069 %h:=print_all; print_all:=t; 1070 hh:=print_; print_:=30; 1071 typeeqlist(list q); 1072 %print_all:=h; 1073 print_:=hh 1074 >>$ 1075 % To try to use q=0 for a substitution if possible: 1076 to_do_list:=cons(list('subst_level_35,{q}),to_do_list)$ 1077 cp_sq2p_val(q)$ 1078 h3:=get(q,'pval)$ 1079 start_level(1,list {'equal,0,h3})$ 1080 1081 h:=get(q,'fac)$ % to add it to ineq_ afterwards 1082 if not pairp h then h:={get(q,'sqval)}$ 1083 1084 recycle_fcts:=nil$ 1085 l:=crackmain_if_possible_remote(pdes,forg)$ 1086 hh:=restore_and_merge(l,pdes,forg)$ 1087 1088 pdes:= car hh; 1089 forg:=cadr hh; % was not assigned above as it has not changed probably 1090 delete_backup()$ 1091 start_level(2,list {'ineq,0,h3})$ 1092 if print_ then << 1093 terpri()$ 1094 write "now back to the substitution of ",cdr w," by ",p$ 1095 >>$ 1096 for each h3 in h do % pdes:= % <=<=<=<= 1097 addSQineq(pdes,h3,nil); 1098 fcteval p$ % fcteval_... properties were deleted in addSQineq 1099 if contradiction_ then <<cases_:=t$contradiction_:=nil>> % but no 1100 % further investigation of this case 1101 else << 1102 % If one of the factors of q was an overall factor of the 1103 % equation, i.e. if car w is not anymore the coefficient of f 1104 % in p then new determination of car w (coeff. of f) is needed 1105 hh:=coeffn({'!*sq,get(p,'sqval),t},cdr w,1); 1106 w:=(if pairp hh and (car hh='!*sq) then cadr hh 1107 else simp hh) . cdr w; 1108 1109% The following lines are commented out (16.5.2008) because: 1110% - if contradiction_ then this must come from addSQineq (because it was 1111% set to nil in restore_and_merge()) and then there is no solution in 1112% this case, 1113% - if the substitution crashes then there is no backup file for this 1114% second case, so the substitution should be done in the crackmain() call 1115% - there should be no jump to again: and no misslabeling of cases occur 1116% instead this substitution is written into the to_do list to be done first 1117% in the crackmain call 1118% 1119% if contradiction_ or null l then << 1120% contradiction_:=nil$ 1121% l:=do_subst(md,p,w,pdes,forg,vl,pdelimit,keep_eqn)$ 1122%% if l and null car l then << % not contradiction but not used 1123% l1:=delete(p,l1); 1124% if l1 then << 1125% pdes:=cadr l; 1126% forg:=caddr l; 1127% l:=nil; 1128% goto again 1129% >> 1130% >>; 1131% if l then << 1132% l:=cdr l; 1133% add_to_last_steps({'subst,stepcounter_,cdr w, 1134% get(p,'printlength),pdelimit}) 1135% >> 1136% 1137% >> else << 1138 1139% New, instead of substitution: 1140 to_do_list:=cons(cons('subst_level_35, 1141 if sub_fc then {{p},sub_fc} 1142 else {{p}}), 1143 to_do_list)$ 1144 1145 % To avoid a loop the picked w='fcteval_nli is now stored as 1146 % w='fcteval_nca 1147 if md>2 then << 1148 h:=get(p,'fcteval_nli)$ 1149 if member(w,h) then << % otherwise p had just one term 1150 % where the non-zero coefficient was a factor which 1151 % is dropped by now, i.e. no further fix needed. 1152 % More generally, in addineq() and updateSQfcteval() 1153 % the following should be unnecessary by now 1154 h:=delete(w,h)$ 1155 put(p,'fcteval_nli,h)$ 1156 put(p,'fcteval_nca,cons(w,get(p,'fcteval_nca)))$ 1157 if md=3 then << 1158 h:=get(p,'fct_nli_lin)$ 1159 h:=delete(w,h)$ 1160 put(p,'fct_nli_lin,h)$ 1161 >> else 1162 if md=4 then << 1163 h:=get(p,'fct_nli_nca)$ 1164 h:=delete(w,h)$ 1165 put(p,'fct_nli_nca,h)$ 1166 >> else 1167 if md=5 then << 1168 h:=get(p,'fct_nli_nli)$ 1169 h:=delete(w,h)$ 1170 put(p,'fct_nli_nli,h)$ 1171 >> else 1172 if md=6 then << 1173 h:=get(p,'fct_nli_nus)$ 1174 h:=delete(w,h)$ 1175 put(p,'fct_nli_nus,h)$ 1176 >> 1177 >> 1178 >>$ 1179 cases_:=t$ 1180 % no backup of global data 1181 h:=crackmain_if_possible_remote(pdes,forg)$ 1182 % No recovery of global data as this crackmain will end now too. 1183 1184 % Because no data are changed, computation could just continue 1185 % without crackmain() sub-call but then combining the 1186 % different results would be difficult. 1187 % No delete_backup() as this has already been done. 1188 if contradiction_ then contradiction_:=nil 1189 else l:=merge_crack_returns(h,l) 1190% >> 1191 >> 1192 >> 1193 >>$ 1194 >>$ 1195 1196 if null !*batch_mode and null expert_mode and confirm_subst then 1197 while tempchng do << 1198 w:=car tempchng; tempchng:=cdr tempchng; 1199 if car w=1 then <<hh:=get(cadr w,'fcteval_lin); 1200 hh:=cons(caddr w,hh); 1201 put(cadr w,'fcteval_lin,hh)>> else 1202 if car w=2 then <<hh:=get(cadr w,'fcteval_nca); 1203 hh:=cons(caddr w,hh); 1204 put(cadr w,'fcteval_nca,hh)>> else 1205 <<hh:=get(cadr w,'fcteval_nli); 1206 hh:=cons(caddr w,hh); 1207 put(cadr w,'fcteval_nli,hh); 1208 if car w=3 then <<hh:=get(cadr w,'fct_nli_lin); 1209 hh:=cons(caddr w,hh); 1210 put(cadr w,'fct_nli_lin,hh)>> else 1211 if car w=4 then <<hh:=get(cadr w,'fct_nli_nca); 1212 hh:=cons(caddr w,hh); 1213 put(cadr w,'fct_nli_nca,hh)>> else 1214 if car w=5 then <<hh:=get(cadr w,'fct_nli_nli); 1215 hh:=cons(caddr w,hh); 1216 put(cadr w,'fct_nli_nli,hh)>> else 1217 if car w=6 then <<hh:=get(cadr w,'fct_nli_nus); 1218 hh:=cons(caddr w,hh); 1219 put(cadr w,'fct_nli_nus,hh)>> 1220 >>; 1221 flag1(cadr w,'to_eval) 1222 >>$ 1223 1224 return if contradiction_ then nil % list(nil,nil) 1225 else if cases_ then list l 1226 else l$ 1227end$ 1228 1229 1230symbolic procedure subst_level_0(arglist)$ % module 3 1231make_subst( car arglist, % all pdes 1232 cadr arglist,caddr arglist, % forg,vl 1233 cadddr arglist, % pdes to choose from 1234 subst_0, % length_limit for pde to use 1235 target_limit_3, % pdelimit for pdes to be changed 1236 t, % less_vars 1237 nil, % no_df 1238 t, % no_cases 1239 nil, % lin_subst 1240 nil, % min_growth 1241 nil, % cost_limit 1242 nil, % keep_eqn 1243 if length arglist > 4 then nth(arglist,5) 1244 else nil % sub_fc 1245 )$ 1246 1247symbolic procedure subst_level_03(arglist)$ % module 4 1248make_subst( car arglist, % all pdes 1249 cadr arglist,caddr arglist, % forg,vl 1250 cadddr arglist, % pdes to choose from 1251 subst_0, % length_limit for pde to use 1252 target_limit_3, % pdelimit for pdes to be changed 1253 nil, % less_vars 1254 t, % no_df 1255 t, % no_cases 1256 nil, % lin_subst 1257 nil, % min_growth 1258 nil, % cost_limit 1259 nil, % keep_eqn 1260 if length arglist > 4 then nth(arglist,5) 1261 else nil % sub_fc 1262 )$ 1263 1264symbolic procedure subst_level_04(arglist)$ % module 45 1265make_subst( car arglist, % all pdes 1266 cadr arglist,caddr arglist, % forg,vl 1267 cadddr arglist, % pdes to choose from 1268 subst_1, % length_limit for pde to use 1269 target_limit_1, % pdelimit for pdes to be changed 1270 nil, % less_vars 1271 t, % no_df 1272 t, % no_cases 1273 nil, % lin_subst 1274 nil, % min_growth 1275 nil, % cost_limit 1276 nil, % keep_eqn 1277 if length arglist > 4 then nth(arglist,5) 1278 else nil % sub_fc 1279 )$ 1280 1281symbolic procedure subst_level_05(arglist)$ % module 5 1282make_subst( car arglist, % all pdes 1283 cadr arglist,caddr arglist, % forg,vl 1284 cadddr arglist, % pdes to choose from 1285 subst_3, % length_limit for pde to use 1286 target_limit_3, % pdelimit for pdes to be changed 1287 nil, % less_vars 1288 t, % no_df 1289 t, % no_cases 1290 nil, % lin_subst 1291 nil, % min_growth 1292 nil, % cost_limit 1293 nil, % keep_eqn 1294 if length arglist > 4 then nth(arglist,5) 1295 else nil % sub_fc 1296 )$ 1297 1298symbolic procedure subst_level_1(arglist)$ % module 15 1299make_subst( car arglist, % all pdes 1300 cadr arglist,caddr arglist, % forg,vl 1301 cadddr arglist, % pdes to choose from 1302 subst_1, % length_limit for pde to use 1303 target_limit_2, % pdelimit for pdes to be changed 1304 t, % less_vars 1305 nil, % no_df 1306 nil, % no_cases 1307 nil, % lin_subst 1308 nil, % min_growth 1309 nil, % cost_limit 1310 nil, % keep_eqn 1311 if length arglist > 4 then nth(arglist,5) 1312 else nil % sub_fc 1313 )$ 1314 1315symbolic procedure subst_level_2(arglist)$ % module 18 1316make_subst( car arglist, % all pdes 1317 cadr arglist,caddr arglist, % forg,vl 1318 cadddr arglist, % pdes to choose from 1319 subst_3, % length_limit for pde to use 1320 target_limit_3, % pdelimit for pdes to be changed 1321 t, % less_vars 1322 nil, % no_df 1323 t, % no_cases 1324 nil, % lin_subst 1325 nil, % min_growth 1326 nil, % cost_limit 1327 nil, % keep_eqn 1328 if length arglist > 4 then nth(arglist,5) 1329 else nil % sub_fc 1330 )$ 1331 1332symbolic procedure subst_level_3(arglist)$ % module 16 1333make_subst( car arglist, % all pdes 1334 cadr arglist,caddr arglist, % forg,vl 1335 cadddr arglist, % pdes to choose from 1336 subst_2, % length_limit for pde to use 1337 target_limit_1, % pdelimit for pdes to be changed 1338 nil, % less_vars 1339 nil, % no_df 1340 nil, % no_cases 1341 nil, % lin_subst 1342 nil, % min_growth 1343 nil, % cost_limit 1344 nil, % keep_eqn 1345 if length arglist > 4 then nth(arglist,5) 1346 else nil % sub_fc 1347 )$ 1348 1349symbolic procedure subst_level_33(arglist)$ % module 19 1350make_subst( car arglist, % all pdes 1351 cadr arglist,caddr arglist, % forg,vl 1352 cadddr arglist, % pdes to choose from 1353 subst_3, % length_limit for pde to use 1354 target_limit_3, % pdelimit for pdes to be changed 1355 nil, % less_vars 1356 nil, % no_df 1357 t, % no_cases 1358 t, % lin_subst 1359 nil, % min_growth 1360 nil, % cost_limit 1361 nil, % keep_eqn 1362 if length arglist > 4 then nth(arglist,5) 1363 else nil % sub_fc 1364 )$ 1365 1366symbolic procedure subst_level_35(arglist)$ % module 20 1367make_subst( car arglist, % all pdes 1368 cadr arglist,caddr arglist, % forg,vl 1369 cadddr arglist, % pdes to choose from 1370 subst_3, % length_limit for pde to use 1371 target_limit_3, % pdelimit for pdes to be changed 1372 nil, % less_vars 1373 nil, % no_df 1374 t, % no_cases 1375 nil, % lin_subst 1376 nil, % min_growth 1377 nil, % cost_limit 1378 nil, % keep_eqn 1379 if length arglist > 4 then nth(arglist,5) 1380 else nil % sub_fc 1381 )$ 1382 1383symbolic procedure subst_level_4(arglist)$ % module 21 1384make_subst( car arglist, % all pdes 1385 cadr arglist,caddr arglist, % forg,vl 1386 cadddr arglist, % pdes to choose from 1387 subst_3, % length_limit for pde to use 1388 target_limit_3, % pdelimit for pdes to be changed 1389 nil, % less_vars 1390 nil, % no_df 1391 nil, % no_cases 1392 nil, % lin_subst 1393 nil, % min_growth 1394 nil, % cost_limit 1395 nil, % keep_eqn 1396 if length arglist > 4 then nth(arglist,5) 1397 else nil % sub_fc 1398 )$ 1399 1400symbolic procedure subst_level_45(arglist)$ % module 6 1401make_subst( car arglist, % all pdes 1402 cadr arglist,caddr arglist, % forg,vl 1403 cadddr arglist, % pdes to choose from 1404 subst_3, % length_limit for pde to use 1405 target_limit_3, % pdelimit for pdes to be changed 1406 nil, % less_vars 1407 nil, % no_df 1408 t, % no_cases 1409 nil, % lin_subst 1410 t, % min_growth 1411 cost_limit5, % cost_limit 1412 nil, % keep_eqn 1413 if length arglist > 4 then nth(arglist,5) 1414 else nil % sub_fc 1415 )$ 1416 1417symbolic procedure subst_level_5(arglist)$ % module 17 1418make_subst( car arglist, % all pdes 1419 cadr arglist,caddr arglist, % forg,vl 1420 cadddr arglist, % pdes to choose from 1421 subst_3, % length_limit for pde to use 1422 target_limit_3, % pdelimit for pdes to be changed 1423 nil, % less_vars 1424 nil, % no_df 1425 nil, % no_cases 1426 nil, % lin_subst 1427 t, % min_growth 1428 nil, % cost_limit 1429 nil, % keep_eqn 1430 if length arglist > 4 then nth(arglist,5) 1431 else nil % sub_fc 1432 )$ 1433 1434 1435symbolic procedure best_fac_pde(pdes)$ 1436% pdes must be pdes for which their 'fac property is pairp 1437begin scalar p,md,mdgr,mtm,f,dgr,tm,bestp; 1438 md:=1000; mtm:=100000; 1439 for each p in pdes do << 1440 % compute the max degree of any factor 1441 mdgr:=0$ 1442 for each f in get(p,'fac) do << 1443 dgr:=pde_degree_SQ(f,smemberl(get(p,'rational),f))$ 1444 if dgr>mdgr then mdgr:=dgr 1445 >>$ 1446 tm:=get(p,'length)$ 1447 if (mdgr<md) or ((mdgr=md) and (tm<mtm)) then 1448 <<bestp:=p; md:=mdgr; mtm:=tm>>; 1449 >>; 1450 return {bestp,md,mtm} 1451end$ 1452 1453algebraic procedure start_let_rules$ 1454begin scalar ruli; 1455 lisp (oldrules!*:=nil)$ % to fix a REDUCE bug 1456 ruli:={}; 1457 if cot(!%x) neq (1/tan(!%x)) then let explog_$ 1458 if lisp(userrules_) neq {} then let lisp userrules_$ 1459 if sin(!%x)**2+cos(!%x)**2 neq 1 then <<ruli:=cons(1,ruli);let trig1_>> else ruli:=cons(0,ruli)$ 1460 if cosh(!%x)**2 neq (sinh(!%x)**2 + 1) then <<ruli:=cons(1,ruli);let trig2_>> else ruli:=cons(0,ruli)$ 1461 if sin(!%x)*tan(!%x/2)+cos(!%x) neq 1 then <<ruli:=cons(1,ruli);let trig3_>> else ruli:=cons(0,ruli)$ 1462 if sin(!%x)*cot(!%x/2)-cos(!%x) neq 1 then <<ruli:=cons(1,ruli);let trig4_>> else ruli:=cons(0,ruli)$ 1463 if cos(2*!%x) + 2*sin(!%x)**2 neq 1 then <<ruli:=cons(1,ruli);let trig5_>> else ruli:=cons(0,ruli)$ 1464 if sin(2*!%x) neq 2*cos(!%x)*sin(!%x) then <<ruli:=cons(1,ruli);let trig6_>> else ruli:=cons(0,ruli)$ 1465 if sinh(2*!%x) neq 2*sinh(!%x)*cosh(!%x) then <<ruli:=cons(1,ruli);let trig7_>> else ruli:=cons(0,ruli)$ 1466 if cosh(2*!%x) neq 2*cosh(!%x)**2-1 then <<ruli:=cons(1,ruli);let trig8_>> else ruli:=cons(0,ruli)$ 1467 if sqrt(!%x*!%y) neq sqrt(!%x)*sqrt(!%y) then <<ruli:=cons(1,ruli);let sqrt1_>> else ruli:=cons(0,ruli)$ 1468 if sqrt(!%x/!%y) neq sqrt(!%x)/sqrt(!%y) then <<ruli:=cons(1,ruli);let sqrt2_>> else ruli:=cons(0,ruli)$ 1469 return ruli; 1470end$ 1471 1472algebraic procedure stop_let_rules(ruli)$ 1473begin 1474 clearrules explog_$ 1475 if (lisp(userrules_) neq {}) and 1476 (lisp (zerop reval {'difference, 1477 car cdadr userrules_, 1478 cadr cdadr userrules_})) 1479 then clearrules lisp userrules_$ 1480 if first ruli = 1 then clearrules sqrt2_$ ruli:=rest ruli$ 1481 if first ruli = 1 then clearrules sqrt1_$ ruli:=rest ruli$ 1482 if first ruli = 1 then clearrules trig8_$ ruli:=rest ruli$ 1483 if first ruli = 1 then clearrules trig7_$ ruli:=rest ruli$ 1484 if first ruli = 1 then clearrules trig6_$ ruli:=rest ruli$ 1485 if first ruli = 1 then clearrules trig5_$ ruli:=rest ruli$ 1486 if first ruli = 1 then clearrules trig4_$ ruli:=rest ruli$ 1487 if first ruli = 1 then clearrules trig3_$ ruli:=rest ruli$ 1488 if first ruli = 1 then clearrules trig2_$ ruli:=rest ruli$ 1489 if first ruli = 1 then clearrules trig1_$ ruli:=rest ruli$ 1490end$ 1491 1492%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1493% procedures for finding an optimal substitution % 1494%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1495 1496symbolic procedure fbts(a,b)$ 1497% fbts ... first better than second 1498(cadr a <= cadr b) and 1499(caddr a <= caddr b) and 1500(cadddr a <= cadddr b)$ 1501 1502 1503symbolic procedure list_subs(p,fevl,fli,mdu)$ 1504% p is an equation, fevl a substitution list of p, 1505% fli is a list of lists (f,p1,p2,..) where 1506% f is a function, 1507% pi are lists (eqn,nco,nte,mdu) where 1508% eqn is an equation that can be used for substituting f 1509% nco is the number of terms of the coefficient of f in the eqn 1510% nte is the number of terms without f in the eqn 1511% mdu is the kind of substitution (1:lin, 2:nca, 3:nli) 1512% returns an extended fli 1513% 1514begin 1515 scalar a,f,nco,nte,cpy,cc,ntry; 1516 for each a in fevl do << 1517 f:=cdr a; 1518% nco:=no_of_terms(car a); 1519 nco:=no_of_tm_sf numr car a; 1520 nte:=get(p,'terms); 1521 nte:=if nte=1 then 0 1522 else nte-nco$ 1523 ntry:={p,nco,nte,mdu}$ 1524 1525 % Is there already any substitution list for f? 1526 cpy:=fli; 1527 while cpy and (f neq caar cpy) do cpy:=cdr cpy$ 1528 if null cpy then fli:=cons({f,ntry},fli) % no, there was not 1529 else << % yes, there was one 1530 cc:=cdar cpy$ 1531 while cc and (null fbts(car cc,ntry)) do cc:=cdr cc$ 1532 if null cc then << % ntry is at least in one criterium better 1533 % than a known one 1534 rplaca(cpy,cons(f,cons(ntry,cdar cpy))); 1535 cc:=cdar cpy$ % just the list of derivatives with ntry as the first 1536 while cdr cc do 1537 if fbts(ntry,cadr cc) then rplacd(cc,cddr cc) 1538 else cc:=cdr cc$ 1539 >> 1540 >> 1541 >>; 1542 return fli 1543end$ 1544 1545symbolic procedure cwrno(n,r)$ 1546% number of terms of (a1+a2+..+an)**r if ai are pairwise prime 1547% number of combinations of r factors out of n possible factors 1548% with repititions and without order = (n+r-1 over r) 1549<<n:=n+r-1; 1550 % The rest of the procedure computes binomial(n,r). 1551 if 2*r>n then r:=n-r; 1552 (for i:=1:r product (n+1-i))/ 1553 (for i:=1:r product i ) 1554>>$ 1555 1556symbolic procedure besu(ic1,mdu1,ic2,mdu2)$ 1557% Is the first substitution better than the second? 1558((mdu1<mdu2) and (ic1<=ic2)) or 1559((mdu1=mdu2) and (ic1< ic2)) or 1560% ########## difficult + room for improvement as the decision is 1561% actually dependent on how precious memory is 1562% (more memory --> less cases and less time): 1563((mdu1=2) and (ic1<(ic2+ 4))) or 1564((mdu1=3) and (ic1<(ic2+25)))$ 1565 1566symbolic procedure search_subs(pdes,sbpdes,cost_limit,no_cases)$ 1567% returns a list {mdu ,p ,car get(p,'fcteval_???)}, i.e. the 3rd argument is SQ 1568begin 1569 scalar fli,p,el,f,fpl,dv,drf,d,ffl,hp,ff,nco,be,s,nte,ic,fp, 1570 rm,mc,subli,mdu,tr_search,h$ 1571 1572 % at first find the list of all functions that could be substituted 1573 % using one of the equations sbpdes together with 1574 % a list of such sbpdes, the number of terms in the coeff and 1575 % the type of substitution 1576 1577% tr_search:=t$ 1578 1579 for each p in sbpdes do fcteval(p)$ 1580 1581 % assuming that equations in sbpdes are sorted by size beginning with shortest 1582 fp:=sbpdes; 1583 h:=nil; 1584 while fp and null h do 1585 if get(car fp,'terms)>2 then fp:=nil else 1586 if null (h:=get(car fp,'fcteval_lin)) then fp:=cdr fp; 1587 if fp then return {1,car fp,car h}$ 1588 1589 while sbpdes do << 1590 p:=car sbpdes; sbpdes:=cdr sbpdes; 1591 fli:=list_subs(p,get(p,'fcteval_lin),fli,1)$ 1592 fli:=list_subs(p,get(p,'fcteval_nca),fli,2)$ 1593 if null no_cases then fli:=list_subs(p,get(p,'fcteval_nli),fli,3)$ 1594 1595 if s then if get(p,'length)>(3*s) then sbpdes:=nil else else 1596 if fli then s:=get(p,'length) 1597 >>$ 1598 1599 if tr_search then << 1600 write"equations substitution: (eqn, no of coeff. t., no of other t., mdu)"$ 1601 terpri()$ 1602 for each el in fli do <<write el;terpri()>>$ 1603 >>$ 1604 1605 if fli then 1606 if (null cdr fli) and % one function 1607 (null cddar fli) then % one equation, i.e. no choice 1608 return << 1609 fli:=cadar fli; % fli is now (eqn,nco,nte,mdu) 1610 mdu:=cadddr fli; 1611 {mdu,car fli,car get(car fli,if mdu = 1 then 'fcteval_lin else 1612 if mdu = 2 then 'fcteval_nca else 1613 'fcteval_nli) } 1614 >> else 1615 % (more than 1 fct.) or (only 1 function and more than 1 eqn.) 1616 for each el in fli do << % for any function to be substituted 1617 % (for the format of fli see proc list_subs) 1618 1619 f:=car el$ el:=cdr el$ 1620 % el is now a list of possible eqn.s to use for subst. of f 1621 1622 fpl:=nil$ % fpl will be a list of lists (p,hp,a1,a2,..) where 1623 % p is an equation that involves f, 1624 % hp the highest power of f in p 1625 % ai are lists {ff,cdr d,nco} where ff is a derivative of f, 1626 % cdr d its power and nco the number of coefficients 1627 for each p in pdes do << % for each equation in which f could be subst. 1628 dv:=get(p,'derivs)$ % ((fct var1 n1 ...).pow) 1629 drf:=nil$ 1630 for each d in dv do 1631 if caar d = f then drf:=cons(d,drf)$ 1632 % drf is now the list of powers of derivatives of f in p 1633 1634 ffl:=nil$ % ffl will be a list of derivatives of f in p 1635 % together with the power of f and number of 1636 % terms in the coeff. 1637 if drf then << % f occurs in this equation and we estimate the increase 1638 hp:=0$ 1639 for each d in drf do << 1640 if cdar d then ff:=cons('df,car d) 1641 else ff:=caar d; 1642 nco:=coeffn({'!*sq,get(p,'sqval),t},ff,cdr d); 1643 nco:=if pairp nco and (car nco='!*sq) then no_of_tm_sf numr cadr nco 1644 else 1; 1645 if cdr d > hp then hp:=cdr d$ 1646 ffl:=cons({ff,cdr d,nco},ffl); 1647 >> 1648 >>; 1649 1650 if drf then fpl:=cons(cons(p,cons(hp,ffl)),fpl); 1651 >>$ 1652 1653 % now all information about all occurences of f is collected and for 1654 % all possible substitutions of f the cost will be estimated and the 1655 % cheapest substitution for f will be determined 1656 1657 be:=nil; % be will be the best equation with an associated min. cost mc 1658 for each s in el do << 1659 % for each possible equation that can be used to subst. for f 1660 1661 % number of terms of (a1+a2+..+an)**r = n+r-1 over r 1662 % f = (a1+a2+..+a_nte) / (b1+b2+..+b_nco) 1663 nco:=cadr s; 1664 nte:=caddr s; 1665 ic:= - get(car s,'terms); % ic will be the cost associated with 1666 % substituting f by car s and car s 1667 % will be dropped after the substitution 1668 for each fp in fpl do 1669 if (car s) neq (car fp) then << 1670 rm:=get(car fp,'terms); % to become the number of terms without f 1671 hp:=cadr fp; 1672 ic:=ic - rm; % as the old eqn. car fp will be replaced 1673 for each ff in cddr fp do << % for each power of each deriv. of f 1674 ic:=ic + (caddr ff)* % number of terms of coefficient of ff 1675 cwrno(nte,cadr ff)* % (numerator of f)**(power of ff) 1676 cwrno(nco,hp - cadr ff); % (denom. of f)**(hp - power of ff) 1677 rm:=rm - caddr ff; % caddr ff is the number of terms with ff 1678 >>; 1679 % Now all terms containing f in car fp have been considered. The 1680 % remaining terms are multiplied with (denom. of f)**hp 1681 ic:=ic + rm*cwrno(nco,hp) 1682 >>; 1683 1684 % Is this substitution better than the best previous one? 1685 if (null be) or besu(ic,cadddr s,mc,mdu) then 1686 <<be:=car s; mc:=ic; mdu:=cadddr s>>; 1687 1688 >>; 1689 1690 % It has been estimated that the substitution of f using the 1691 % best eqn be has an additional cost of ic terms 1692 1693 if tr_search and (length el > 1) then << 1694 terpri()$ 1695 write"Best substitution for ",f," : ",{ic,f,be,mdu}$ 1696 >>$ 1697 1698 if (null cost_limit) or (ic<cost_limit) then 1699 subli:=cons({ic,mdu,f,be},subli)$ 1700 >>$ 1701 1702 % Now pick the best substitution 1703 if subli then << 1704 s:=car subli; 1705 subli:=cdr subli; 1706 for each el in subli do 1707 if besu(car el,cadr el,car s,cadr s) then s:=el$ 1708 1709 if tr_search then << 1710 terpri()$ 1711 write"Optimal substitution:"$terpri()$ 1712 write" replace ",caddr s," with the help of ",cadddr s,","$terpri()$ 1713 if car s < 0 then write" saving ", - car s," terms, " 1714 else write" with a cost of ",car s," additional terms, "$ 1715 terpri()$ 1716 write if cadr s = 1 then " linear substitution" else 1717 if cadr s = 2 then " nonlinearity inceasing substitution" else 1718 " with case distinction" $ 1719 >>$ 1720 1721 el:=get(cadddr s,if (cadr s) = 1 then 'fcteval_lin else 1722 if (cadr s) = 2 then 'fcteval_nca else 1723 'fcteval_nli); 1724 while (caddr s) neq (cdar el) do el:=cdr el; 1725 1726 return {cadr s,cadddr s,car el} 1727 % = {mdu ,p ,car get(p,'fcteval_???)} 1728 >>$ 1729 1730end$ 1731 1732%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1733% procedures for doing a substitution `bottom up' % 1734%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1735 1736symbolic procedure bottom_up_subst(arglist)$ 1737% This procedure uses the gobal variable largest_fully_substituted_in 1738% which is the latest/largest equation in the list pdes of 1739% equations (when starting to count from the shortest) in which 1740% all possible substitutions of shorter equations have been done 1741% that do not lead to case distinctions. 1742% Only one substitution will be made because the result could be 1743% big and then it would be better to continue substituting in a 1744% different shorter equation. 1745 1746if currently_to_be_substituted_in = '!*SUBST_DONE!* then nil else 1747 1748if (null car arglist) or (null cdar arglist) then 1749<<currently_to_be_substituted_in:='!*SUBST_DONE!*;nil>> else 1750 1751begin scalar pcp,found,pdes,fns,p,a,h; 1752 currently_to_be_substituted_in:=nil; % <-- for now as not in list 1753 % not_passed_back in crinit.red 1754 if null currently_to_be_substituted_in then 1755 currently_to_be_substituted_in:=cadar arglist; 1756 fns:=get(currently_to_be_substituted_in,'fcts); 1757 1758 pcp:=car arglist; 1759 while (currently_to_be_substituted_in neq '!*SUBST_DONE!*) and 1760 null found do 1761 1762 if car pcp=currently_to_be_substituted_in then 1763 if null cdr pcp then currently_to_be_substituted_in:='!*SUBST_DONE!* 1764 else << 1765 currently_to_be_substituted_in:=cadr pcp; 1766 fns:=get(currently_to_be_substituted_in,'fcts); 1767 pcp:=car arglist 1768 >> else 1769 if << % Is substitution possible? 1770 p:=car pcp; 1771 if get(p,'starde) then nil 1772 else << 1773 a:=get(p,'fcteval_lin); 1774 if null a or freeof(fns,cdar a) then << 1775 a:=get(p,'fcteval_nca); 1776 if a and freeof(fns,cdar a) then a:=nil 1777 >>; 1778 a 1779 >> 1780 >> then << % do substitution 1781 pdes:=car arglist$ 1782 pdes:=eqinsert(do_one_subst(car a,cdr a,currently_to_be_substituted_in, 1783 ftem_,get(currently_to_be_substituted_in,'vars), 1784 get(p,'level),p,pdes), 1785 delete(currently_to_be_substituted_in,pdes))$ 1786 for each h in pdes do drop_rl_with(currently_to_be_substituted_in,h); 1787 put(currently_to_be_substituted_in,'rl_with,nil); 1788 for each h in pdes do drop_dec_with(currently_to_be_substituted_in,h,'dec_with_rl); 1789 put(currently_to_be_substituted_in,'dec_with_rl,nil); 1790 flag(list currently_to_be_substituted_in,'to_int); 1791 found:=t 1792 >> else pcp:=cdr pcp; % Check next equation for substitution 1793 1794 return 1795 if currently_to_be_substituted_in='!*SUBST_DONE!* then nil 1796 else {pdes,cadr arglist} 1797end$ 1798 1799 1800%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1801% procedures for substitution of a derivative by a new function % 1802%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1803 1804symbolic procedure check_subst_df(pdes,forg)$ 1805% yields a list of derivatives which occur in all 1806% pdes and in forg 1807begin scalar l,l1,l2,n,cp,not_to_substdf$ 1808 if pdes then << 1809 for each s in pdes do l:=union(for each a in get(s,'derivs) 1810 collect car a,l)$ % all derivs 1811 for each s in forg do 1812 if pairp s then 1813 l:=union(for each a in all_deriv_search_SF(numr caddr s,ftem_) collect car a, 1814 union(for each a in all_deriv_search_SF(denr caddr s,ftem_) collect car a, 1815 l))$ 1816 1817 for each s in fsub_ do 1818 l:=union(for each a in all_deriv_search_SF(numr cadr cdr s,ftem_) collect car a, 1819 union(for each a in all_deriv_search_SF(denr cadr cdr s,ftem_) collect car a, 1820 l))$ 1821 1822 l1:=df_min_list(l)$ 1823 l:=nil$ 1824 for each s in l1 do 1825 if pairp s and not member(car s,not_to_substdf) then << 1826 l:=cons(cons('df,s),l)$ 1827 not_to_substdf:=cons(car s,not_to_substdf) 1828 >> $ 1829 1830 % Derivatives of functions should only be substituted if the 1831 % function occurs in at least 2 equations or forg functions 1832 while l do << 1833 n:=0; % counter 1834 cp:=pdes; 1835 while cp and (n<2) do << 1836 if member(cadar l,get(car cp,'fcts)) then n:=add1 n; 1837 cp:=cdr cp 1838 >>; 1839 cp:=forg; 1840 while cp and (n<2) do << 1841 if (pairp car cp) and (caar cp = 'equal) and 1842 member(cadar l,get(cadr car cp,'fcts)) then n:=add1 n; 1843 cp:=cdr cp 1844 >>; 1845 if n=2 then l2:=cons(car l,l2); 1846 l:=cdr l 1847 >> 1848 >>$ 1849 return l2$ 1850end$ 1851 1852symbolic procedure df_min_list(dflist)$ 1853% yields the intersection of the derivatives of each function 1854% in the list of deriv. dflist. 1855% e.g. dflist='((f x z) (g x) (g) (f y) (h x y) (h x z)) 1856% ==> result='(f g (h x)) 1857if dflist then 1858begin scalar l,d,m,lmax$ 1859while dflist do 1860 <<m:=car dflist$ 1861 dflist:=cdr dflist$ 1862 l:=nil$ 1863 while dflist do 1864 <<if (d:=df_min(car dflist,m)) then m:=d 1865 else l:=cons(car dflist,l)$ 1866 dflist:=cdr dflist$ 1867 >>$ 1868 if pairp m and null cdr m then lmax:=cons(car m,lmax) 1869 else lmax:=cons(m,lmax)$ 1870 dflist:=l$ 1871 >>$ 1872return lmax$ 1873end$ 1874 1875symbolic procedure df_min(df1,df2)$ 1876% yields the intersection of derivatives df1,df2 1877% e.g. df_min('(f x y),'(f x z))='(f x), df_min('(f x z),'(g x))=nil 1878<<if not pairp df1 then df1:=list df1$ 1879 if not pairp df2 then df2:=list df2$ 1880 if car df1=car df2 then 1881 if (df1:=df_min1(cdr df1,cdr df2)) then cons(car df2,df1) 1882 else car df2 1883>>$ 1884 1885symbolic procedure df_min1(df1,df2)$ 1886begin scalar l,a$ 1887while df1 do 1888 <<a:=car df1$ 1889 if not zerop (a:=min(dfdeg(df1,car df1),dfdeg(df2,car df1))) then 1890 l:=cons(car df1,l)$ 1891 if a>1 then l:=cons(a,l)$ 1892 df1:=cdr df1$ 1893 if df1 and numberp car df1 then df1:=cdr df1>>$ 1894return reverse l$ 1895end$ 1896 1897symbolic procedure dfsubst_forg(p,g,d,forg)$ 1898% substitute the function d in forg by an integral g 1899% of the function p 1900for each h in forg collect 1901if pairp h and member(d,get(cadr h,'fcts)) then << 1902 put(cadr h,'fcts,fctinsert(p,delete(d,get(cadr h,'fcts))))$ 1903 % reval subst(g,d,h) 1904 {'equal,cadr h,simp {'!*sq,subsq(caddr h,{(d . {'!*sq,g,t})}),nil}} 1905 % {'equal,cadr h,simp!* {'!*sq,subsq(caddr h,{(d . {'!*sq,g,t})}),nil}} 1906 % - simp!* {'!*sq,..,nil} to simplify poly using identities, like i^2 -> -1 1907>> else h$ 1908 1909symbolic procedure expand_INT(p,varlist)$ 1910if null varlist then p 1911else begin scalar v,n$ 1912 v:=car varlist$ 1913 varlist:=cdr varlist$ 1914 if pairp(varlist) and numberp(car varlist) then 1915 <<n:=car varlist$ 1916 varlist:=cdr varlist>> 1917 else n:=1$ 1918 for i:=1:n do p:=list('int,p,v)$ 1919 return expand_INT(p,varlist) 1920end$ 1921 1922symbolic procedure rational_less(a,b)$ 1923% a and b are two reval'ued rational numbers in prefix form 1924% It returns the boolean value of a<b 1925if (pairp a) and 1926 (car a='quotient) then rational_less(cadr a,reval{'times,caddr a,b}) else 1927if (pairp b) and 1928 (car b='quotient) then rational_less(reval{'times,caddr b,a},cadr b) else 1929if (pairp a) and (car a='minus) then 1930if (pairp b) and (car b='minus) then cadr a > cadr b 1931 else not rational_less(cadr a,reval{'minus,b}) 1932 else 1933if (pairp b) and (car b='minus) then 1934if a<0 then not rational_less(reval{'minus,a},cadr b) 1935 else nil 1936 else a<b$ 1937 1938symbolic procedure substitution_weight(k,l,m,n,p,q)$ 1939 % This function computes a weight for an equation to 1940 % be used for a substitution 1941 % k .. number of occurences as factor, 1942 % l .. total degree of factor as homogeneous polynomial, 1943 % m .. number of appearances in eqns, 1944 % n .. number of terms 1945 % p .. a factor taking care of the number of lists in ineq_or 1946 % and their length of which f is an element 1947 % q .. =1 if contains flin_ unknowns else =0 1948reval {'quotient,{'times,l,n,q},{'times,p,{'expt,{'plus,{'times,2,k},m},2}}}$ 1949 1950symbolic procedure test_factors_for_substitution(p,pv)$ 1951% p is the equation, pv the list of factors 1952begin scalar h,h1,h4,h3$ 1953% pv:=get(p,'fac)$ 1954% if not pairp pv then return nil; 1955 h:=get(p,'split_test); 1956 if null h then << % check all factors whether they can be used for subst. 1957 h1:=pv$ % the list of factors in p 1958 h4:=t$ 1959 % make an equation from the coefficient 1960 while h1 and h4 do << 1961 h3:=mkeqSQ(car h1,nil,nil,get(p,'fcts),get(p,'vars),allflags_, 1962 t,list(0),nil,nil)$ 1963 contradiction_:=nil$ 1964 % the last argument is nil to avoid having a lasting effect on pdes 1965 h1:=cdr h1$ 1966 fcteval(h3)$ 1967 if not(get(h3,'fcteval_lin) or get(h3,'fcteval_nca)) then h4:=nil; 1968 drop_pde(h3,nil,nil) 1969 >>$ 1970 h:=if h4 then 1 % p can be splited into substitutable equations 1971 else 0; % p can not be splited into only " equations 1972 put(p,'split_test,h) 1973 >>; 1974 return h 1975end$ 1976 1977symbolic procedure get_fact_pde(pdes,aim_at_subst)$ 1978% look for pde in pdes which can be factorized 1979begin scalar p,pv,f,fcl,fcc,h,h1,h2,h3,h4,h5,h6,h7,h8,eql,tr_gf$ 1980 % tr_gf:=t$ 1981 1982 h1:=pdes; 1983 if null aim_at_subst then 1984 1985 % Highest priority has an equation with a case2sep entry, i.e. which 1986 % allows to conclude either the independence of a function on a 1987 % variable or allows to do a direct separation. 1988 1989 while h1 and null h2 do << 1990 h3:=get(car h1,'case2sep)$ 1991 if h3 then if not member(h3,ineq_) then h2:=car h1 1992 else << 1993 put(car h1,'case2sep,nil)$ 1994 h4:=stardep3(get(car h1,'vars),get(car h1,'kern),get(car h1,'derivs))$ 1995 if h4 then <<put(car h1,'starde,{(0 . car h4)})$ flag1(car h1,'to_sep)>>$ 1996 h1:=cdr h1 1997 >> else h1:=cdr h1 1998 >>$ 1999 if h2 then return h2$ 2000 2001 % now choose an equation that minimizes a weight computed from 2002 % the weights of its factors, 2003 % the weight of a factor = 2004 % (if an atom then number of all equations it occurs 2005 % else the number of equations it occurs as a factor)/ 2006 % the total degree of this factor/ 2007 % the number of factors of the equation 2008 % The factor with the highest weight is to be set to 0 first. 2009 2010 % 1) collecting a list of all suitable equations eql and a list 2011 % of all factors of any equation, listing for each factor 2012 % in how many equations it appears 2013 for each p in pdes do << 2014 pv:=get(p,'fac)$ 2015 if pairp pv then << % otherwise not (simply) factorizable 2016 % increment the counter of appearances of each factor 2017 h1:=pv$ 2018 while h1 do << % for each factor 2019 f:=car h1; h1:=cdr h1; % f is in SQ form 2020 2021 fcc:=fcl$ 2022 2023 % fcl is list of lists 2024 % (factor itself, 2025 % no of occurences as factor, 2026 % total degree of factor as homogeneous polynomial, 2027 % number of appearances in eqns, 2028 % number of terms in factor 2029 % a factor taking care of the number of lists in ineq_or 2030 % and their length of which f is an element 2031 % does the factor contain flin_ (=1) or not (=0) ) 2032 2033 while fcc and (caar fcc neq f) do fcc:=cdr fcc$ 2034 2035 if fcc then << % factor had already appeared 2036 h:=cons(f,cons(add1 cadar fcc,cddar fcc)); 2037 rplaca(fcc,h); 2038 >> else << % factor is new 2039 % Computing the total degree of the factor 2040 if fhom_ then << 2041 h2:=find_hom_deg_SF(numr f)$ 2042 h2:=(car h2) + (cadr h2) 2043 >> else h2:=1; 2044 % If f is a function then counting in how many equations it appears 2045 if no_number_atom_SQ f then << % count in how many equations f does occur 2046 h5:=mvar numr f; % the function 2047 h3:=0; % the counter 2048 h4:=pdes; 2049 while h4 do << 2050 if not freeof(get(car h4,'fcts),h5) then h3:=add1 h3; 2051 h4:=cdr h4 2052 >> 2053 >> else h3:=1$ 2054 2055 % The number of terms of f: 2056 h4:=no_of_tm_sf numr f$ 2057 2058 % Is f element of a list-element (with at most 8 elements) of ineq_or? 2059 h5:=1$ 2060 h6:=ineq_or$ 2061 while h6 do << 2062 h7:=length h6; 2063 if h7 < 9 then 2064 if member({f},car h6) then h5:={'times,h5,{'expt,2,{'difference,9,h7}}}$ 2065 h6:=cdr h6 2066 >>$ 2067 h5:=reval h5; 2068 2069 if flin_ and not freeoflist(f,flin_) then h6:=1 else h6:=0$ 2070 2071 fcl:=cons({f,1,h2,h3,h4,h5,h6},fcl) 2072 >> 2073 >>$ % done for all factors 2074 2075 % check whether each factor can be used for subst., i.e. whether 2076 % this equation should be factorized 2077 if null aim_at_subst then h:=1 2078 else h:=test_factors_for_substitution(p,pv)$ 2079 2080 % adding the equation to the ones suited for factorization 2081 if not zerop h then eql:=cons(p,eql) 2082 2083 >> 2084 >>$ % looked at all factorizable equations 2085 2086 % Any equation worth to be case-split? 2087 if null eql then return nil; 2088 2089 % Now that it is known how often each factor appears in all equations, 2090 % each factor can be given a weight and each equation be given a weight 2091 2092 fcl:=for each h in fcl collect cons(car h,cons( 2093 substitution_weight(cadr h,caddr h,cadddr h,car cddddr h, 2094 cadr cddddr h,caddr cddddr h), 2095 cdr h ))$ 2096 2097 % fcl is now list of lists 2098 % (factor itself, 2099 % substitution_weight, 2100 % no of occurences as factor, 2101 % total degree of factor as homogeneous polynomial, 2102 % number of appearances in eqns, 2103 % number of terms in factor 2104 % a factor taking care of the number of lists in ineq_or 2105 % and their length of which f is an element 2106 % does the factor contain flin_ (=1) or not (=0) ) 2107 2108 h2:=nil; % h2 is the best equation, its weight will be h3 and the 2109 % factors of the best equation sorted by weight will be h4 2110 % in the new order they will be set to zero 2111 2112 for each p in eql do << 2113 pv:=get(p,'fac)$ % list of factors in SQ-form 2114 h8:=length pv$ % number of factors of p 2115 h5:=nil; % the list of (weight of factor . factor) 2116 h6:=1; % the weight of equation p, the less the better 2117 while pv do << % for each factor 2118 h:=assoc(car pv,fcl); % the properties of the factor 2119 if tr_gf then << write "h assoc= ",h$terpri()>>$ 2120 h5:=cons(cons(reval cadr h,car h),h5); % cadr h is substitution_weight 2121 % h6:={'plus,h6,cadr h}; % if adding up the weights for the factors and 2122 % then taking the equation with the minimal SUM of weights will 2123 % favour equations where the weight of the highest weight factor in 2124 % the equation is lowest. So, how good the best of the factors is 2125 % does little matter. We therefore change to product. 2126 h6:={'times,h6,cadr h}; 2127 %if not zerop nth(h,8) then h6:={'times,30,h6}; 2128 % The following change was made to use an equation less 2129 % if more flin_ functions appear. 2130 if (null lin_problem) and 2131 (not zerop nth(h,8)) then h6:={'times, 2132 {'expt,2,{'plus,4, 2133 get(p,'nfct_lin)}}, 2134 h6}; 2135 % evaluating flin_ functions 2136 % has lower priority as quicker progress is made by evaluating the 2137 % fewer non-flin_ unknowns (at least in bi-lin alg problems) 2138 pv:=cdr pv 2139 >>$ 2140 h6:=reval {'times,{'expt,2,h8},h6}; % punishment of many factors 2141 if null h2 or rational_less(h6,h3) then << 2142 h2:=p; 2143 h3:=h6; 2144 h4:=h5 2145 >> 2146 >>$ 2147 2148 % sort factors of the equation (h2) by their substitution_weight, 2149 h4:=rat_idx_sort h4$ 2150 if tr_gf then << write "h4(rat_idx_sort'ed)= ",h4$terpri()>>$ 2151 2152 % Putting the flin_ factor last is bad if this factor comes up in 2153 % many equations, like in the case of bi-linear systems when at the 2154 % end only one flin_ function is left being a factor of all equations 2155 % 2156% if flin_ then << 2157% h5:=h4; % car h5 will be the factor involving flin_ functions 2158% while h5 and freeoflist(cdar h5,flin_) do h5:=cdr h5; 2159% if h5 then h4:=append(delete(car h5,h4),list car h5) 2160% >>$ 2161 2162 put(h2,'fac,for each a in h4 collect cdr a)$ 2163 return h2 2164end$ 2165 2166endmodule$ 2167 2168end$ 2169 2170------------------------------------------------------------------------------------------------------ 2171Substitutions with: 2172 2173procedure # length_li pdelimit less_vars no_df no_cases lin_subst min_growth cost_limit keep_eqn 2174 2175subst_level_0 3 subst_0 target_limit_3 t nil t nil nil nil nil 2176subst_level_03 4 subst_0 target_limit_3 nil t t nil nil nil nil 2177subst_level_04 45 subst_1 target_limit_1 nil t t nil nil nil nil 2178subst_level_05 5 subst_3 target_limit_3 nil t t nil nil nil nil 2179subst_level_1 15 subst_1 target_limit_2 t nil nil nil nil nil nil 2180subst_level_2 18 subst_3 target_limit_3 t nil t nil nil nil nil 2181subst_level_3 16 subst_2 target_limit_1 nil nil nil nil nil nil nil 2182subst_level_33 19 subst_3 target_limit_3 nil nil t t nil nil nil 2183subst_level_35 20 subst_3 target_limit_3 nil nil t nil nil nil nil 2184subst_level_4 21 subst_3 target_limit_3 nil nil nil nil nil nil nil 2185subst_level_45 6 subst_3 target_limit_3 nil nil t nil t cost_limit5 nil 2186subst_level_5 17 subst_3 target_limit_3 nil nil nil nil t nil nil 2187 2188subst_0:=2$ maximal printlength of an expression to be substituted 2189subst_1:=8$ 2190subst_2:=20$ 2191subst_3:=10^3$ 2192subst_4:=nil$ 2193cost_limit5:=100$ % maximal number of extra terms generated by a subst. 2194target_limit_0:=10^2$ % maximal product length(pde)*length(subst_expr) 2195target_limit_1:=10^3$ % nil=no length limit 2196target_limit_2:=10^5$ 2197target_limit_3:=nil$ 2198 2199make_subst 2200 get_subst 2201 do_subst 2202 simp_all_ineq_with_subst_SQ 2203 do_one_subst 2204 mkeqSQ 2205 updatesq 2206 mkeqSQ 2207 updatesq 2208 2209 2210tr make_subst 2211tr get_subst 2212tr do_subst 2213tr do_one_subst 2214tr mkeqSQ 2215tr updatesq 2216tr fcteval 2217tr simp_all_ineq_with_subst_SQ 2218