Artifact bcfc9fe5df6cf56a1491dd362e4cf1f8c92e42530280d5293a2b35bc456e9404:
- Executable file
r38/packages/support/patches.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 34590) [annotate] [blame] [check-ins using] [more...]
module patches; % Patches to correct problems in REDUCE 3.8. % Author: Anthony C. Hearn. % Copyright (c) 2004, 2005, 2006, 2007 Anthony C. Hearn. All Rights Reserved. global '(patch!-date!* patch!-url!-list!*); patch!-date!* := "11-Jan-2007"; patch!-url!-list!* := '("http://reduce-algebra.com/support/patches/patches.fsl"); % Bugs fixed by these patches. % 26 Jun 04. With rounded arithmetic, solving some linear equation % problems could lead to a catastrophic error. % 8 Jul 04. Some non-zero integrals (e.g., int(e^(a^(1/3)*x)*sin x,x)) % returned zero. % 5 Aug 04. Using RLFI with latex on could lead to invalid operator errors. % 2 Sep 04. In rare circumstances, floating point conversion could give % an extraneous error. % 6 Sep 04. With rational on, some non-zero factorizations could produce % a zero coefficient (e.g., on rational; % factorize(r^((1/4*n^2 - 1/4*n + 1)/(n - 1)));). % 28 Sep 04. Some integrals would not return a closed form solution % with algint on that would with algint off % (e.g., int(sqrt(x-1)/(sqrt x*(x-1)),x)). % 10 Dec 04. With dfprin on, some products and sums printed incorrectly. % 31 Jan 05. Some integrals involving square roots could run forever. % 12 Feb 05. SOLVE could produce a spurious recursive loop (e.g., % solve((4*e^(y^3/3)*cte+2x^2+y^3+3)/e^(y^3/3),y)). % 20 Apr 05. int(e^(-a^(1/4)*(-1)^(1/4)*x),x); terminated with an error. % 2 May 05. int(e^(-a^(1/4)*(-1)^(1/4)*i*x)*b+(1/4)*e^(-a^(1/4) % *(-1)^(1/4)*i*x)*x,x); terminated with an error. % 22 May 05. Some integrals, e.g., int(e^((3sqrt 5+1)*x)*(sqrt 5+1) % +e^((3sqrt 5-1)*x)*(sqrt 5+1),x), never completed. % 30 May 05. SOLVE could produce a spurious "Zero divisor" error % (e.g., solve({log tan(y/2),y+1/x},{x,y})). % 4 Oct 05. DEG did not work with rational coefficients (e.g., % deg(x**3/a-x/5+1/4,x)). % 5 Oct 05. Some SOLVE calculations could give a spurious "Zero Divisor" % error (e.g., ex0:= sqrt(a^2-y^2); solve((-log(( - x + a + y) % /ex0) + log((x + a + y)/ex0) + x - (a^2 - y^2)/ex0),y); % 16 Nov 05. System errors could occur with rounded and combineexpt on. % (E.g., on rounded,combineexpt; 0.183*e^x*t^4.39;). % 22 Nov 05. Some definite integrals with variables other than x could % give a wrong answer, e.g., int(e^(-y),y,0,x). % 9 Dec 05. With combineexpt on, expressions could be dropped (e.g., % on combineexpt; 4*e^(-3*h/2) - 3*h*e^(-h) + 2*e^(-h)). % 4 Feb 06. Setcrackflags() was not set in crack, but needed to be. % 20 Feb 06. The rule for df(Jacobidn(~u,~m),~u) was wrong. % 21 Feb 06. Evaluating some integrals could suppress the printing of % the results. % 22 Feb 06. Some sub evaluations could include superfluous terms like % x = (x^(1/7))^7. % 23 May 06. Derivatives and integrals of matrices were not computed. % 18 Aug 06. After nospur, some traces were still evaluated. % 29 Sep 06. With dfprint on, derivatives of integrals would print in a % truncated form. % 11 Jan 07. With rounded arithmetic and factor on, a non-numeric % argument error could occur. % Alg declarations. fluid '(sublist!*); patch alg; % 16 Nov 05, 9 Dec 05. symbolic procedure exptunwind(u,v); begin scalar x,x1,x2,y,z,z2; a: if null v then return u; x := caar v; x1 := cadr x; x2 := caddr x; y := cdar v; v := cdr v; if !*combineexpt and null domainp u and null red u and (z2 := kernels u) and null cdr z2 then u := {(({'expt,car z2,ldeg u} . 1) . lc u)}; while (z := assocp1(x1,v)) and (z2 := simp {'plus,{'times,x2,y},{'times,caddar z,cdr z}}) and (!*combineexpt or (fixp numr z2 and fixp denr z2)) do <<if fixp numr z2 and fixp denr z2 then <<x2 := divide(numr z2,denr z2); if car x2>0 then <<if fixp x1 then u := multf(x1**car x2,u) else u := multpf(mksp(x1,car x2),u); z2 := cdr x2 ./ denr z2>>; y := numr z2>> else if domainp numr z2 then y := 1 else <<y := lcoeffgcd cdr comfac numr z2; if not fixp y then y := 1>>; x2 := prepsq(quotf(numr z2,y) ./ denr z2); v := delete(z,v)>>; if !*combineexpt and y=1 and fixp x1 then <<while (z := assocp2(x2,v)) and cdr z=1 and fixp cadar z do <<x1 := cadar z * x1; v := delete(z,v)>>; if eqcar(x2,'quotient) and fixp cadr x2 and fixp caddr x2 and cadr x2<caddr x2 then <<z := nrootn(x1**cadr x2,caddr x2); if cdr z = 1 then u := multd(car z,u) else if car z = 1 then u := multf(formsf(x1,x2,1),u) else <<u := multd(car z,u); v := (list('expt,cdr z,x2) . 1) . v>>>> else u := multf(formsf(x1,x2,y),u)>> else u := multf(formsf(x1,x2,y),u); go to a end; % 22 Feb 06. symbolic procedure subeval0 u; begin scalar x,y,z,ns; while cdr u do <<if not eqcar(car u,'equal) then x := car u . x else if not(cadar u = (y := reval caddar u)) then x := {caar u,cadar u,y} . x; u := cdr u>>; if null x then return car u else u := nconc(reversip x,u); if u member sublist!* then return mk!*sq !*p2q mksp('sub . u,1) else sublist!* := u . sublist!*; if null(u and cdr u) then rederr "SUB requires at least 2 arguments"; (while cdr u do <<x := reval car u; if getrtype x eq 'list then u := append(cdr x,cdr u) else <<if not eqexpr x then errpri2(car u,t); y := cadr x; if null getrtype y then y := !*a2kwoweight y; if getrtype caddr x then ns := (y . caddr x) . ns else z := (y . caddr x) . z; u := cdr u>>>>) where !*evallhseqp=nil; x := aeval car u; return subeval1(append(ns,z),x) end; symbolic procedure subsubf(l,expn); begin scalar x,y; for each j in l do if car j neq (y := prepsq!* simp!* cdr j) then x := (car j . y) . x; l := reversip x; if null l then return expn; y := nil; for each j in cddr expn do if (x := assoc(j,l)) then <<y := x . y; l := delete(x,l)>>; expn := sublis(l,car expn) . for each j in cdr expn collect subsublis(l,j); if null y then return expn; expn := aconc!*(for each j in reversip!* y collect list('equal,car j,aeval cdr j),expn); return if l then subeval expn else mk!*sq !*p2q mksp('sub . expn,1) end; % 23 May 06. symbolic procedure reval1(u,v); (begin scalar x,y; if null u then return nil else if stringp u then return u else if fixp u then return if flagp(dmode!*,'convert) then reval2(u,v) else u else if atom u then if null subfg!* then return u else if idp u and (x := get(u,'avalue)) then if u memq varstack!* then recursiveerror u else <<varstack!* := u . varstack!*; return if y := get(car x,'evfn) then apply2(y,u,v) else reval1(cadr x,v)>> else nil else if not idp car u then errpri2(u,t) else if car u eq '!*sq then return if caddr u and null !*resimp then if null v then u else prepsqxx cadr u else reval2(u,v) else if flagp(car u,'remember) then return rmmbreval(u,v) else if flagp(car u,'opfn) then return reval1(opfneval u,v) else if x := get(car u,'psopfn) then <<u := apply1(x,cdr u); if x := get(x,'cleanupfn) then u := apply2(x,u,v); return u>> else if arrayp car u then return reval1(getelv u,v); return if x := getrtype u then if y := get(x,'evfn) then apply2(y,u,v) else rerror(alg,101, list("Missing evaluation for type",x)) else if not atom u and not atom cdr u and (y := getrtype cadr u) and null(y eq 'list and cddr u) and (x := get(y,'aggregatefn)) and (not(x eq 'matrixmap) or flagp(car u,'matmapfn)) and not flagp(car u,'boolean) and not !*listargs and not flagp(car u,'listargp) then apply2(x,u,v) else reval2(u,v) end) where varstack!* := varstack!*; symbolic procedure getrtype2 u; begin scalar x; return if (x := get(car u,'rtype)) and (x := get(x,'rtypefn)) then apply1(x,cdr u) else if x := get(car u,'rtypefn) then apply1(x,cdr u) else if flagp(car u,'matmapfn) and cdr u and getrtype cadr u eq 'matrix then 'matrix else nil end; endpatch; patch arith; % 2 Sep 04. symbolic procedure read!:num(n); if fixp n then make!:ibf(n, 0) else if not(numberp n or stringp n) then bflerrmsg 'read!:num else begin integer j,m,sign; scalar ch,u,v,l,appear!.,appear!/; j := m := 0; sign := 1; u := v := appear!. := appear!/ := nil; l := explode n; loop: ch := car l; if digit ch then << u := ch . u; j := j + 1 >> else if ch eq '!. then << appear!. := t; j := 0 >> else if ch eq '!/ then << appear!/ := t; v := u; u := nil >> else if ch eq '!- then sign := -1 else if ch memq '(!E !D !B !e !d !b) then go to jump; if l := cdr l then goto loop else goto make; jump: while l := cdr l do <<if digit(ch := car l) or ch eq '!- then v := ch . v >>; l := reverse v; if car l eq '!- then m := - compress cdr l else m:= compress l; make: u := reverse u; v := reverse v; if appear!/ then return conv!:r2bf(make!:ratnum(sign*compress v,compress u), if !:bprec!: then !:bprec!: else 170); if appear!. then j := - j else j := 0; if sign = 1 then u := compress u else u := - compress u; return round!:mt (decimal2internal (u, j + m), !:bprec!:) where !:bprec!: := if !:bprec!: then !:bprec!: else msd!: abs u end; endpatch; patch crack; setcrackflags(); endpatch; % Defint declarations. symbolic smacro procedure listsq(u); for each uu in u collect simp!* uu; patch defint; symbolic procedure new_meijer(u); begin scalar f,y,mellin,new_mellin,m,n,p,q,old_num,old_denom,temp,a1, b1,a2,b2,alpha,num,denom,n1,temp1,temp2,coeff,v,var,new_var,new_y, new_v,k; f := prepsq simp car u; y := caddr u; mellin := bastab(car f,cddr f); temp := car cddddr mellin; var := cadr f; if not idp VAR then RETURN error(99,'FAIL); temp := reval algebraic(sub(x=var,temp)); mellin := {car mellin,cadr mellin,caddr mellin,cadddr mellin,temp}; temp := reduce_var(cadr u,mellin,var); alpha := simp!* car temp; new_mellin := cdr temp; if car cddddr new_mellin neq car cddddr mellin then << k := car cddddr mellin; y := reval algebraic(sub(var=y,k)); new_y := simp y>> else << new_var := car cddddr new_mellin; new_y := simp reval algebraic(sub(x=y,new_var))>>; n1 := addsq(alpha,'(1 . 1)); temp1 := {'expt,y,prepsq n1}; temp2 := cadddr new_mellin; coeff := simp!* reval algebraic(temp1*temp2); m := caar new_mellin; n := cadar new_mellin; p := caddar new_mellin; q := car cdddar new_mellin; old_num := cadr new_mellin; old_denom := caddr new_mellin; for i:=1 :n do << if old_num = nil then a1 := append(a1,{simp!* old_num }) else << a1 := append(a1,{simp!* car old_num}); old_num := cdr old_num>>; >>; for j:=1 :m do << if old_denom = nil then b1 := append(b1,{simp!* old_denom }) else << b1 := append(b1,{simp!* car old_denom}); old_denom := cdr old_denom>>; >>; a2 := listsq old_num; b2 := listsq old_denom; if a1 = nil and a2 = nil then num := list({negsq alpha}) else if a2 = nil then num := list(append(a1,{negsq alpha})) else << num := append(a1,{negsq alpha}); num := append({num},a2)>>; if b1 = nil and b2 = nil then denom := list({subtrsq(negsq alpha,'(1 . 1))}) else if b2 = nil then denom := list(b1,subtrsq(negsq alpha,'(1 . 1))) else << denom := list(b1,subtrsq(negsq alpha,'(1 . 1))); denom := append(denom,b2)>>; v := gfmsq(num,denom,new_y); if v = 'fail then return simp 'fail else v := prepsq subsq(v,list(prepsq new_y . y)); if eqcar(v,'meijerg) then new_v := v else new_v := simp v; return multsq(new_v,coeff); end; endpatch; patch hephys; symbolic procedure nospur u; <<rmsubs(); !*nospurp := t; flag(u,'nospur)>>; endpatch; % Int declarations. fluid '(!*purerisch !*trdint gaussiani indexlist intvar sqrt!-places!-alist loglist !*intflag!* listofnewsqrts listofallsqrts sqrt!-intvar basic!-listofallsqrts basic!-listofnewsqrts !*precise dmode!* !*exp !*gcd !*keepsqrts !*limitedfactors !*mcd !*rationalize !*structure !*uncached kord!*); smacro procedure argof u; cadr u; patch int; % 8 Jul 04, 31 Jan 05, 20 Apr 05, 2 May 05. symbolic procedure df2q p; begin scalar n,d,w,x,y,z; if null p then return nil ./ 1; d:=denr lc p; w:=red p; while w do <<d := multf(d,quotf(denr lc w,gcdf(d,denr lc w))); w := red w>>; while p do begin w := sqrt2top lc p; x := multf(xl2f(lpow p,zlist,indexlist),multf(numr w,d)); if null x then return (p := red p); y := denr w; z := quotf(x,y); if null z then <<z := rationalizesq(x ./ y); if denr z neq 1 then <<d := multf(denr z,d); n := multf(denr z,n)>>; z := numr z>>; n := addf(n,z); p := red p end; return tidy!-powersq (n ./ d) end; % 8 Jul 04, 22 May 05. symbolic procedure tidy!-powersq x; begin scalar expts,!*precise,!*keepsqrts; !*keepsqrts := t; x := subs2q x; expts := find!-expts(numr x,find!-expts(denr x,nil)); if null expts then return x; x := subsq(x,for each v in expts collect (car v . list('expt,cadr v,cddr v))); x := subsq(x,for each v in expts collect (cadr v . list('expt,car v,list('quotient,1,cddr v)))); return x end; symbolic procedure find!-expts(ff,l); begin scalar w; if domainp ff then return l; l := find!-expts(lc ff,find!-expts(red ff, l)); ff := mvar ff; if eqcar(ff,'sqrt) then ff := list('expt, cadr ff,'(quotient 1 2)) else if eqcar(ff,'expt) and eqcar(caddr ff,'quotient) and numberp caddr caddr ff then <<w := assoc(cadr ff,l); if null w then <<w := cadr ff . gensym() . 1; l := w . l >>; rplacd(cdr w,lcm(cddr w,caddr caddr ff))>>; return l end; % 28 Sep 04. symbolic procedure look_for_quad(integrand, var, zz); begin if (car zz = 'sqrt and listp cadr zz and caadr zz = 'plus) or (car zz = 'expt and listp cadr zz and caadr zz = 'plus and listp caddr zz and car caddr zz = 'quotient and fixp caddr caddr zz) then << zz := simp cadr zz; if (cdr zz = 1) then << zz := cdr coeff1(prepsq zz, var, nil); if length zz = 2 then return begin scalar a, b; scalar nvar, res, ss; a := car zz; b := cadr zz; if (depends(a,var) or depends(b,var)) then return nil; nvar := gensym(); if !*trint then << prin2 "Linear shift suggested "; prin2 a; prin2 " "; prin2 b; terpri(); >>; integrand := subsq(integrand, list(var . list('quotient, list('difference, list('expt,nvar,2),a), b))); integrand := multsq(integrand, simp list('quotient,list('times,nvar,2), b)); if !*trint then << prin2 "Integrand is transformed by substitution to "; printsq integrand; prin2 "using substitution "; prin2 var; prin2 " -> "; printsq simp list('quotient, list('difference,list('expt,nvar,2),a), b); >>; res := integratesq(integrand, nvar, nil, nil); ss := list(nvar . list('sqrt,list('plus,list('times,var,b), a))); res := subsq(car res, ss) . subsq(multsq(cdr res, simp list('quotient,b, list('times,nvar,2))), ss); return res; end else if length zz = 3 then return begin scalar a, b, c; a := car zz; b := cadr zz; c:= caddr zz; if (depends(a,var) or depends(b,var) or depends(c,var)) then return nil; a := simp list('difference, a, list('times,b,b, list('quotient,1,list('times,4,c)))); if null numr a then return nil; b := simp list('quotient, b, list('times, 2, c)); c := simp c; return if minusf numr c then << if minusf numr a then begin scalar !*hyperbolic; !*hyperbolic := t; return look_for_invhyp(integrand,nil,var,a,b,c) end else look_for_asin(integrand,var,a,b,c)>> else << if minusf numr a then look_for_invhyp(integrand,t,var,a,b,c) else look_for_invhyp(integrand,nil,var,a,b,c) >> end else if length zz = 5 then return begin scalar a, b, c, d, e, nn, dd, mm; a := car zz; b := cadr zz; c:= caddr zz; d := cadddr zz; e := car cddddr zz; if not(b = 0) or not(d = 0) then return nil; if (depends(a,var) or depends(c,var)) or depends(e,var) then return nil; nn := numr integrand; dd := denr integrand; if denr(mm :=quotsq(nn ./ 1, !*kk2q var)) = 1 and even_power(numr mm, var) and even_power(dd, var) then << return sqrt_substitute(numr mm, dd, var); >>; if denr(mm :=quotsq(dd ./ 1, !*kk2q var)) = 1 and even_power(nn, var) and even_power(numr mm, var) then << return sqrt_substitute(nn, multf(dd,!*kk2f var), var); >>; return nil; end; >>>>; return nil end; % 21 Feb 06. symbolic procedure simpint u; if atom u or null cdr u or cddr u and (null cdddr u or cddddr u) then rerror(int,1,"Improper number of arguments to INT") else if cddr u then simpdint u else begin scalar ans,dmod,expression,variable,loglist,oldvarstack, !*intflag!*,!*purerisch,cflag,intvar,listofnewsqrts, listofallsqrts,sqrtfn,sqrt!-intvar,sqrt!-places!-alist, basic!-listofallsqrts,basic!-listofnewsqrts,coefft, varchange,w,!*precise; !*intflag!* := t; variable := !*a2k cadr u; if not(idp variable or pairp variable and numlistp cdr variable) then <<varchange := variable . intern gensym(); if !*trint then printc {"Integration kernel", variable, "replaced by simple variable", cdr varchange}; variable := cdr varchange>>; intvar := variable; w := cddr u; if w then rerror(int,3,"Too many arguments to INT"); listofnewsqrts:= list mvar gaussiani; listofallsqrts:= list (argof mvar gaussiani . gaussiani); sqrtfn := get('sqrt,'simpfn); put('sqrt,'simpfn,'proper!-simpsqrt); if dmode!* then << if (cflag:=get(dmode!*, 'cmpxfn)) then onoff('complex, nil); if (dmod := get(dmode!*,'dname)) then onoff(dmod,nil)>> where !*msg := nil; begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*limitedfactors,!*mcd, !*rationalize,!*structure,!*uncached,kord!*, ans1,badbit,denexp,erfg,nexp,oneterm; !*keepsqrts := !*limitedfactors := t; !*exp := !*gcd := !*mcd := !*structure := !*uncached := t; dmode!* := nil; if !*algint then << sqrt!-intvar:=!*q2f simpsqrti variable; if (red sqrt!-intvar) or (lc sqrt!-intvar neq 1) or (ldeg sqrt!-intvar neq 1) then interr "Sqrt(x) not properly formed" else sqrt!-intvar:=mvar sqrt!-intvar; basic!-listofallsqrts:=listofallsqrts; basic!-listofnewsqrts:=listofnewsqrts; sqrtsave(basic!-listofallsqrts,basic!-listofnewsqrts, list(variable . variable))>>; coefft := (1 ./ 1); expression := int!-simp car u; if varchange then <<depend1(car varchange,cdr varchange,t); expression := int!-subsq(expression,{varchange})>>; denexp := 1 ./ denr expression; nexp := numr expression; while not atom nexp and null cdr nexp and not depends(mvar nexp,variable) do <<coefft := multsq(coefft,(((caar nexp) . 1) . nil) ./ 1); nexp := lc nexp>>; ans1 := nil; while nexp do begin scalar x,zv,tmp; if atom nexp then <<x := !*f2q nexp; nexp := nil>> else <<x := !*t2q car nexp; nexp := cdr nexp>>; x := multsq(x,denexp); zv := zvars(getvariables x,zv,variable,t); tmp := ans1; while tmp do <<if zv=caar tmp then <<rplacd(car tmp,addsq(cdar tmp,x)); tmp := nil; zv := nil>> else tmp := cdr tmp>>; if zv then ans1 := (zv . x) . ans1 end; if length ans1 = 1 then oneterm := t; nexp := ans1; ans := nil ./ 1; badbit:=nil ./ 1; while nexp do <<u := cdar nexp; if !*trdint then <<princ "Integrate"; printsq u; princ "with Zvars "; print caar nexp>>; erfg := erfg!*; ans1 := errorset!*(list('integratesq,mkquote u, mkquote variable,mkquote loglist, mkquote caar nexp), !*backtrace); erfg!* := erfg; nexp := cdr nexp; if errorp ans1 then badbit := addsq(badbit,u) else <<ans := addsq(caar ans1, ans); badbit:=addsq(cdar ans1,badbit)>>>>; if !*trdint then <<prin2 "Partial answer="; printsq ans; prin2 "To do="; printsq badbit>>; if badbit neq '(nil . 1) then <<setkorder nil; badbit := reordsq badbit; ans := reordsq ans; coefft := reordsq coefft; if !*trdint then <<princ "Retrying..."; printsq badbit>>; if oneterm and ans = '(nil . 1) then ans1 := nil else ans1 := errorset!*(list('integratesq,mkquote badbit, mkquote variable,mkquote loglist,nil), !*backtrace); if null ans1 or errorp ans1 then ans := addsq(ans,simpint1(badbit . variable . w)) else <<ans := addsq(ans,caar ans1); if not smemq(variable, ans) then ans := nil ./ 1; if cdar ans1 neq '(nil . 1) then ans := addsq(ans, simpint1(cdar ans1 . variable . w)) >>>>; end; ans := multsq(coefft,ans); if !*trdint then << printc "Resimp and all that"; printsq ans >>; put('int,'simpfn,'simpiden); put('sqrt,'simpfn,sqrtfn); << if dmod then onoff(dmod,t); if cflag then onoff('complex,t)>> where !*msg := nil; oldvarstack := varstack!*; varstack!* := nil; ans := errorset!*(list('int!-resub,mkquote ans,mkquote varchange),t); put('int,'simpfn,'simpint); varstack!* := oldvarstack; return if errorp ans then error1() else car ans end; endpatch; patch mathpr; % 10 Dec 04, 29 Sep 06. symbolic procedure dflayout u; (begin scalar op, args, w; w := car (u := cdr u); u := cdr u; if smemq('int,w) then !*noarg := nil; if !*noarg and (atom w or not get(car w, 'op)) then << if atom w then << op := w; args := assoc(op, depl!*); if args then args := cdr args >> else << op := car w; args := cdr w >>; remember!-args(op, args); w := op >>; maprin w; if u then << u := layout!-formula('!!dfsub!! . u, 0, nil); if null u then return 'failed; w := 1 + cddr u; putpline((update!-pline(0, -w, caar u) . cdar u) . ((cadr u - w) . (cddr u - w))) >> end) where !*noarg = !*noarg; endpatch; patch matrix; % 26 Jun 04. symbolic procedure sparse_backsub(exlis,varlis); begin scalar d,z,c; if null exlis then return nil; d := lc car exlis; foreach x in exlis do begin scalar s,p,v,r; p := lc x; v := mvar x; x := red x; while not domainp x and mvar x member varlis do <<if (c := atsoc(mvar x,z)) then s := addf(multf(lc x,cdr c),s) else r := addf(!*t2f lt x,r); x := red x>>; s := negf quotff(addf(multf(addf(r,x),d),s),p); z := (v . s) . z; end; for each p in z do cdr p := cancel(cdr p ./ d); return z end; symbolic procedure quotff(u,v); if null u then nil else (if x then x else (if denr y = 1 then numr y else rederr "Invalid division in backsub") where y=rationalizesq(u ./ v)) where x=quotf(u,v); % 23 May 06. symbolic procedure matsm1 u; begin scalar x,y,z; integer n; a: if null u then return z else if eqcar(car u,'!*div) then go to d else if atom car u then go to er else if caar u eq 'mat then go to c1 else if flagp(caar u,'matmapfn) and cdar u and getrtype cadar u eq 'matrix then x := matsm matrixmap(car u,nil) else <<x := lispapply(caar u,cdar u); if eqcar(x,'mat) then x := matsm x>>; b: z := if null z then x else if null cdr z and null cdar z then multsm(caar z,x) else multm(x,z); c: u := cdr u; go to a; c1: if not lchk cdar u then rerror(matrix,3,"Matrix mismatch"); x := for each j in cdar u collect for each k in j collect xsimp k; go to b; d: y := matsm cadar u; if (n := length car y) neq length y then rerror(matrix,4,"Non square matrix") else if (z and n neq length z) then rerror(matrix,5,"Matrix mismatch") else if cddar u then go to h else if null cdr y and null cdar y then go to e; x := subfg!*; subfg!* := nil; if null z then z := apply1(get('mat,'inversefn),y) else if null(x := get('mat,'lnrsolvefn)) then z := multm(apply1(get('mat,'inversefn),y),z) else z := apply2(get('mat,'lnrsolvefn),y,z); subfg!* := x; z := for each j in z collect for each k in j collect <<!*sub2 := t; subs2 k>>; go to c; e: if null caaar y then rerror(matrix,6,"Zero divisor"); y := revpr caar y; z := if null z then list list y else multsm(y,z); go to c; h: if null z then z := generateident n; go to c; er: rerror(matrix,7,list("Matrix",car u,"not set")) end; symbolic procedure matrixmap(u,v); if flagp(car u,'matmapfn) then matsm!*1 for each j in matsm cadr u collect for each k in j collect simp!*(car u . mk!*sq k . cddr u) else if flagp(car u,'matfn) then reval2(u,v) else typerr(car u,"matrix operator"); put('matrix,'aggregatefn,'matrixmap); flag('(int df taylor),'matmapfn); flag('(det trace),'matfn); endpatch; patch poly; % 6 Sep 04. symbolic procedure fctrf u; (begin scalar !*ezgcd,!*gcd,denom,x,y; if domainp u then return list u else if ncmp!* and not noncomfp u then ncmp!* := nil; !*gcd := t; if null !*limitedfactors and null dmode!* then !*ezgcd := t; if null !*mcd then rerror(poly,15,"Factorization invalid with MCD off") else if null !*exp then <<!*exp := t; u := !*q2f resimp !*f2q u>>; if dmode!* eq '!:rn!: then <<dmode!* := nil; alglist!* := nil . nil; x := simp prepf u; if atom denr x then <<denom := denr x; u := numr x>> else denom := 1>>; if null ncmp!* then <<x := sf2ss u; if homogp x then <<if !*trfac then prin2t "This polynomial is homogeneous - variables scaled"; y := caaar x . listsum caaadr x; x := fctrf1 ss2sf(car(x) . (reverse subs0 cadr x . 1)); x := rconst(y,x); return car x . sort!-factors cdr x>>>>; u := fctrf1 u; if denom then <<alglist!* := nil . nil; dmode!* := '!:rn!:; car u := quotf!*(car u,denom)>>; return car u . sort!-factors cdr u end) where !*exp = !*exp, ncmp!* = ncmp!*; % 4 Oct 05. symbolic procedure deg(u,kern); <<u := simp!* u; tstpolyarg2(u,kern); numrdeg(numr u,kern)>> where dmode!* = gdmode!*; symbolic procedure tstpolyarg2(u,kern); <<for each j in kernels numr u do if j=kern then nil else if depends(j,kern) then typerr(prepsq u,"polynomial"); for each j in kernels denr u do if depends(j,kern) then typerr(prepsq u,"polynomial")>>; % 11 Jan 07. symbolic procedure rnfactor!: u; begin scalar x,y,dmode!*; integer m,n; x := subf(u,nil); if not domainp denr x then return {1,(u . 1)}; y := factorf numr x; n := car y; dmode!* := '!:rn!:; y := for each j in cdr y collect <<n := n*(m := (lnc ckrn car j)**cdr j); quotfd(car j,m) . cdr j>>; return int!-equiv!-chk mkrn(n,denr x) . y end; endpatch; patch rlfi; put('tex,'simpfn,'simpcar); endpatch; % Solve declarations. fluid '(!*cramer bareiss!-step!-size!*); global '(assumptions); patch solve; % 26 Jun 04. symbolic procedure solvelnrsys(exlis,varlis); begin scalar w,x; if w := solvesparsecheck(exlis,varlis) then exlis := w else exlis := exlis . varlis; if null !*cramer and null errorp(x := errorset2{'solvebareiss,mkquote car exlis, mkquote cdr exlis} where bareiss!-step!-size!* = if w then 4 else 2) then exlis := car x else exlis := solvecramer(car exlis,cdr exlis); return solvesyspost(exlis,varlis) end; % 12 Feb 05, 5 Oct 05. symbolic procedure solvesq (ex,var,mul); begin scalar r,x; r:= for each w in solvesq1(ex,var,mul) join if null cadr w or eqcar(x := prepsq caar w,'root_of) or numr subfx(denr ex,{caadr w . x}) then {w}; if r and not domainp denr ex then assumptions:=append(assumptions,{prepf denr ex}); return r end; % 5 Oct 05. symbolic procedure subfx(u,v); begin scalar x; x := errorset2 {'subf,mkquote u,mkquote v}; return if errorp x then 1 ./ 1 else car x end; % 12 Feb 05 symbolic procedure polypeval u; begin scalar bool,v; v := cadr u; u := simpcar u; if cdr u neq 1 then return nil else u := kernels car u; while u and null bool do <<if v neq car u and smember(v,car u) then bool := t; u := cdr u>>; return null bool end; put('polyp,'psopfn,'polypeval); (algebraic << depend(!~p,!~x); clearrules {root_of(~p,~x,~tg)^~n => sub(x=root_of(p,x,tg), -reduct(p,x)/coeffn(p,x,deg(p,x)))^(n-deg(p,x)+1) when fixp n and deg(p,x)>=1 and n>=deg(p,x)}; let root_of(~p,~x,~tg)^~n => sub(x=root_of(p,x,tg), -reduct(p,x)/coeffn(p,x,deg(p,x))) ^ (n-deg(p,x)+1) when polyp(p,x) and fixp n and deg(p,x)>=1 and n>=deg(p,x); nodepend(!~p,!~x); >>) where dmode!*=nil,!*modular=nil,!*rounded=nil,!*complex=nil; % 30 May 05. symbolic procedure solvenonlnrtansolve(u,x,w); begin scalar v,s,z,r,y; integer ar; ar:=!!arbint; v:=caar u;u:=prepf numr simp cdr u; s:=solveeval{u,'tg!-}; !!arbint:=ar; for each q in cdr s do <<z:=reval caddr q; z:=reval sublis(solvenonlnrtansolve1 z,z); !!arbint:=ar; y:=solve0({'equal,{'tan,{'quotient,V,2}},z},x); r:=union(y,r)>>; y := errorset2 {'subf,mkquote w,mkquote{x . 'pi}}; if null errorp y and null numr y then <<!!arbint:=ar; r:=union(solve0({'equal,{'cos,x},-1},x),r)>>; return t.r end; % 5 Oct 05. symbolic procedure check!-solns(z,ex,var); begin scalar x,y; if not errorp (x := errorset2 {'check!-solns1,mkquote z,mkquote ex,mkquote var}) then return car x else if ex = (y := (numr simp!* prepf ex where !*reduced=t)) or errorp (x := errorset2 {'check!-solns1,mkquote z,mkquote y,mkquote var}) then return 'unsolved else return car x end; symbolic procedure check!-solns1(z,ex,var); begin scalar x,y,fv,sx,vs; fv := freevarl(ex,var); for each z1 in z do fv := union(fv,union(freevarl(numr caar z1,var), freevarl(denr caar z1,var))); fv := delete('i,fv); if fv then for each v in fv do if not flagp(v,'constant) then vs := (v . list('quotient,1+random 999,1000)) . vs; sx := if vs then numr subf(ex,vs) else ex; while z do if null cadar z or errorp(y := errorset2 {'check!-solns2,mkquote ex,mkquote z}) then <<z := nil; x := 'unsolved>> else if null(y := car y) or fv and null(y := numr subf(sx,list(caadar z . mk!*sq subsq(caaar z,vs)))) or null numvalue y then <<x := car z . x; z := cdr z>> else z := cdr z; return if null x then 'unsolved else x end; symbolic procedure check!-solns2(ex,z); if smemq('root_of,z) then rederr 'check!-solns else numr subf(ex,{caadar z . mk!*sq caaar z}); endpatch; patch specfn; % 20 Feb 06. algebraic (for all u,m let df(Jacobidn(u,m),u) = -m^2 *Jacobisn(u,m)*Jacobicn(u,m)); endpatch; patch rlisp; !#if (member 'psl lispsystem!*) symbolic procedure global idlist; fluid idlist; symbolic procedure global1 id1; if not get(id1,'vartype) then fluid1 id1; !#endif endpatch; endmodule; end;