login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit
 Topics FrontPage SandBox SanBoxIntPar <-- You are here. ---------- Forwarded message ----------:  From: Waldek Hebisch Date: 31 August 2014 20:19 Subject: [fricas-devel] Reworking extended integration To: fricas-devel@googlegroups.com  As I wrote I am reworking extented integration routines. Now there is some code that you can try. At http://www.math.uni.wroc.pl/~hebisch/fricas/intpar.tar.bz2 there is tarball containing current version. To build run current FriCAS (tarball depend on changes and fixes present in trunk) in the directory obtained by unpacking the tarball and do:  )read compi.input  to compile. After compilation you can load the routines using:  )read load.input  Two examples that fails in trunk, but runs whith new version are:  (1) -> integrate((polylog(3, x) - x*polylog(2, x))/(1 - x)^2, x) 2 - 2x polylog(3,x) + (x - 1)log(- x + 1) + 2x dilog(- x + 1) (1) ------------------------------------------------------------ 2x - 2 Type: Union(Expression(Integer),...) (2) -> integrate(log(x+exp(x))/(x^2 + 3*x + 2) + (log(x*(1+x)) - log(x*(x+2)))*(exp(x) + 1)/(exp(x) + x), x) 2 2 x (2) (- log(x + 2x) + log(x + x))log(%e + x) Type: Union(Expression(Integer),...)  In general you will see difference only for relatively complicated integrals (when there are multiple levels of recursion in Risch algorithm). New code has still significant gaps. One is that parametric logarithmic derivative solver (needed to get bounds in third case of RDE) uses factorization. Currently factorizer I use will fail when base ring is different than Integer or if coefficients of factored polynomials contains nontrivial algebraics. Also I left unimplemented one degenerate case in RDE solver. It seems that this are only gaps in transcendental Liouvillean case. In algebraic case RDE is completely unhandled and extendend integration only works in pure root extensions. Compared to current code new one will probably loose in algebraic case: old code can solve RDE and do integration in general pure algebraic extensions (when there is no nontrivial transcendentals) and some other cases via change of variables. I have some idea to extend parametric RDE routines to algebraic case. This should allow eventually replacing current RDE solver by parametic one (non-parametric RDE can be reduced to parametric RDE). -- Waldek Hebisch )compile LINCOMB.spad spad)abbrev package LINCOMB LinearCombinationUtilities LinearCombinationUtilities(F, UP) : Exports == Implementation where F : Field UP : UnivariatePolynomialCategory F N ==> NonNegativeInteger RF ==> Fraction UP GP ==> LaurentPolynomial(F, UP) Param_Rec_F ==> Record(ratpart : F, coeffs : Vector F) L_Param_F ==> List Param_Rec_F Partial_F ==> Union(Param_Rec_F, "failed") Both_F ==> Record(particular : Partial_F, basis : L_Param_F) Q ==> Fraction(Integer) Param_Rec_Q2 ==> Record(ratpart : F, coeffs : Vector Q) L_Param_Q2 ==> List Param_Rec_Q2 Partial_Q2 ==> Union(Param_Rec_Q2, "failed") Both_Q2 ==> Record(particular : Partial_Q2, basis : L_Param_Q2) Exports ==> with dehomogenize : L_Param_F -> Both_F ++ dehomogenize(ls) converts list of solutions ++ (a, [c0, c1, ..., cn]) to homogeneous equation ++ L(a) + c0 f + c1 g1 + ... + cn gn = 0 into list of ++ solutions of inhomogeneous equation ++ L(a) + f + c1 g1 + ... + cn gn = 0. This transformation ++ works the same for all equations, so we only ++ need list of solutions as argument )if false dehomogenize : L_Param_Q2 -> Both_Q2 ++ dehomogenize(ls) converts list of solutions ++ (a, [c0, c1, ..., cn]) to homogeneous equation ++ L(a) + c0 f + c1 g1 + ... + cn gn = 0 into list of ++ solutions of inhomogeneous equation ++ L(a) + f + c1 g1 + ... + cn gn = 0. This transformation ++ works the same for all equations, so we only ++ need list of solutions as argument )endif lin_comb : (Vector F, List(F)) -> F ++ lin_comb(v, [f1, ..., fn]) computes linear combination ++ v(1) f1 + ... v(n) fn. Vector v and list [f1, ..., fn] ++ must be of equal length. lin_comb : (Vector Q, List(F)) -> F ++ lin_comb(v, [f1, ..., fn]) computes linear combination ++ v(1) f1 + ... v(n) fn. Vector v and list [f1, ..., fn] ++ must be of equal length. lin_comb! : (Vector F, Vector F, List(Vector F)) -> Vector F ++ lin_comb(v, w, [f1, ..., fn]) computes linear combination ++ w + v(1) f1 + ... v(n) fn by modifying w in place. ++ Vector v and list [f1, ..., fn] must be of equal length. lin_comb : (Vector F, List(Vector F)) -> Vector F ++ lin_comb(v, [f1, ..., fn]) computes linear combination ++ v(1) f1 + ... v(n) fn. Vector v and list [f1, ..., fn] ++ must be of equal positive length. lin_comb : (Vector F, List(RF)) -> RF ++ lin_comb(v, [f1, ..., fn]) computes linear combination ++ v(1) f1 + ... v(n) fn. Vector v and list [f1, ..., fn] ++ must be of equal length. lin_comb : (Vector Q, List(RF)) -> RF ++ lin_comb(v, [f1, ..., fn]) computes linear combination ++ v(1) f1 + ... v(n) fn. Vector v and list [f1, ..., fn] ++ must be of equal length. lin_comb : (Vector F, List(GP)) -> GP ++ lin_comb(v, [f1, ..., fn]) computes linear combination ++ v(1) f1 + ... v(n) fn. Vector v and list [f1, ..., fn] ++ must be of equal length. lin_comb : (Vector F, List(UP)) -> UP ++ lin_comb(v, [f1, ..., fn]) computes linear combination ++ v(1) f1 + ... v(n) fn. Vector v and list [f1, ..., fn] ++ must be of equal length. Implementation ==> add dehomogenize_body(Lpar, Par) ==> empty?(ls) => ["failed", []] nn := #(ls(1).coeffs) be1 : Par found : Boolean := false res2 : Lpar := [] for be in ls repeat c := (be.coeffs)(1) c = 0 or found => res2 := cons(be, res2) be1 := be found := true not(found) => ["failed", []] c1inv := 1/(be1.coeffs)(1) ppa := c1inv*be1.ratpart ppv := (c1inv*be1.coeffs)(2..nn) res3 : Lpar := [] for be in res2 repeat c := (be.coeffs)(1) bv := (be.coeffs)(2..nn) rp := be.ratpart if c ~= 0 then rp := rp - c*ppa bv := bv - c*ppv res3 := cons([rp, bv], res3) [[ppa, ppv], res3] dehomogenize(ls : L_Param_F) : Both_F == dehomogenize_body(L_Param_F, Param_Rec_F) )if false dehomogenize(ls : L_Param_Q2) : Both_Q2 == dehomogenize_body(L_Param_Q2, Param_Rec_Q2) )endif lin_comb(v : Vector F, lf : List(F)) : F == res : F := 0 for i in 1..#v for f in lf repeat res := res + v(i)*f res lin_comb(v : Vector Q, lf : List(F)) : F == res : F := 0 for i in 1..#v for f in lf repeat res := res + v(i)::F*f res lin_comb!(u : Vector F, v : Vector F, lw : List(Vector F)) : Vector F == res := v n := #res for i in 1..#u for w in lw repeat c := u(i) for j in 1..n repeat res(j) := res(j) + c*w(j) res lin_comb(u : Vector F, lw : List(Vector F)) : Vector F == n := #first(lw) lin_comb!(u, new(n, 0)$Vector(F), lw) lin_comb(v : Vector F, lf : List(RF)) : RF == res : RF := 0 for i in 1..#v for f in lf repeat res := res + v(i)::UP*f res lin_comb(v : Vector Q, lg : List RF) : RF == res : RF := 0 for i in 1..#v for g in lg repeat res := res + v(i)::F::UP::RF*g res lin_comb(v : Vector F, lf : List(GP)) : GP == res : GP := 0 for i in 1..#v for f in lf repeat res := res + v(i)::GP*f res lin_comb(v : Vector F, lf : List(UP)) : UP == res : UP := 0 for i in 1..#v for f in lf repeat res := res + v(i)*f res spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/8542904851967799624-25px001.spad using old system compiler. LINCOMB abbreviates package LinearCombinationUtilities ------------------------------------------------------------------------ initializing NRLIB LINCOMB for LinearCombinationUtilities compiling into NRLIB LINCOMB processing macro definition dehomogenize_body(Lpar,Par) ==> SEQ(LET(G653: Boolean,empty? ls),exit(1,IF(G653,construct(failed,construct),SEQ(LET(nn,# (ls (One)) coeffs),be1: Par,LET(found: Boolean,false),LET(res2: Lpar,construct),REPEAT(IN(be,ls),SEQ(LET(c,(be coeffs) One),IF(=(c,Zero),exit(1,LET(res2,cons(be,res2))),IF(found,exit(1,LET(res2,cons(be,res2))),noBranch)),LET(be1,be),exit(1,LET(found,true)))),exit(1,IF(found,SEQ(LET(c1inv,/(One,(be1 coeffs) One)),LET(ppa,*(c1inv,be1 ratpart)),LET(ppv,(* c1inv (be1 coeffs)) SEGMENT(2,nn)),LET(res3: Lpar,construct),REPEAT(IN(be,res2),SEQ(LET(c,(be coeffs) One),LET(bv,(be coeffs) SEGMENT(2,nn)),LET(rp,be ratpart),SEQ(LET(G654: Boolean,~=(c,Zero)),exit(1,IF(G654,SEQ(LET(rp,-(rp,*(c,ppa))),exit(1,LET(bv,-(bv,*(c,ppv))))),noBranch))),exit(1,LET(res3,cons(construct(rp,bv),res3))))),exit(1,construct(construct(ppa,ppv),res3))),construct(failed,construct))))))) compiling exported dehomogenize : List Record(ratpart: F,coeffs: Vector F) -> Record(particular: Union(Record(ratpart: F,coeffs: Vector F),failed),basis: List Record(ratpart: F,coeffs: Vector F)) Time: 0.10 SEC. compiling exported lin_comb : (Vector F,List F) -> F Time: 0.01 SEC. compiling exported lin_comb : (Vector Fraction Integer,List F) -> F Time: 0.03 SEC. compiling exported lin_comb! : (Vector F,Vector F,List Vector F) -> Vector F Time: 0.01 SEC. compiling exported lin_comb : (Vector F,List Vector F) -> Vector F Time: 0 SEC. compiling exported lin_comb : (Vector F,List Fraction UP) -> Fraction UP Time: 0.02 SEC. compiling exported lin_comb : (Vector Fraction Integer,List Fraction UP) -> Fraction UP Time: 0.19 SEC. compiling exported lin_comb : (Vector F,List LaurentPolynomial(F,UP)) -> LaurentPolynomial(F,UP) Time: 0.01 SEC. compiling exported lin_comb : (Vector F,List UP) -> UP Time: 0.02 SEC. (time taken in buildFunctor: 0) ;;; *** |LinearCombinationUtilities| REDEFINED ;;; *** |LinearCombinationUtilities| REDEFINED Time: 0 SEC. Warnings: [1] dehomogenize: coeffs has no value [2] dehomogenize: found has no value [3] dehomogenize: be1 has no value [4] dehomogenize: ratpart has no value [5] dehomogenize: res2 has no value Cumulative Statistics for Constructor LinearCombinationUtilities Time: 0.39 seconds finalizing NRLIB LINCOMB Processing LinearCombinationUtilities for Browser database: --->-->LinearCombinationUtilities(constructor): Not documented!!!! --------(dehomogenize ((Record (: particular (Union (Record (: ratpart F) (: coeffs (Vector F))) failed)) (: basis (List (Record (: ratpart F) (: coeffs (Vector F)))))) (List (Record (: ratpart F) (: coeffs (Vector F))))))--------- --------(lin_comb (F (Vector F) (List F)))--------- --------(lin_comb (F (Vector (Fraction (Integer))) (List F)))--------- --------(lin_comb! ((Vector F) (Vector F) (Vector F) (List (Vector F))))--------- --->-->LinearCombinationUtilities((lin_comb! ((Vector F) (Vector F) (Vector F) (List (Vector F))))): Improper first word in comments: lin_comb "lin_comb(\\spad{v},{} \\spad{w},{} [\\spad{f1},{} ...,{} \\spad{fn}]) computes linear combination \\spad{w} + \\spad{v}(1) \\spad{f1} + ... \\spad{v}(\\spad{n}) \\spad{fn} by modifying \\spad{w} in place. Vector \\spad{v} and list [\\spad{f1},{} ...,{} \\spad{fn}] must be of equal length." --------(lin_comb ((Vector F) (Vector F) (List (Vector F))))--------- --------(lin_comb ((Fraction UP) (Vector F) (List (Fraction UP))))--------- --------(lin_comb ((Fraction UP) (Vector (Fraction (Integer))) (List (Fraction UP))))--------- --------(lin_comb ((LaurentPolynomial F UP) (Vector F) (List (LaurentPolynomial F UP))))--------- --------(lin_comb (UP (Vector F) (List UP)))--------- --->-->LinearCombinationUtilities(): Missing Description ; compiling file "/var/aw/var/LatexWiki/LINCOMB.NRLIB/LINCOMB.lsp" (written 05 SEP 2014 03:02:55 AM): ; /var/aw/var/LatexWiki/LINCOMB.NRLIB/LINCOMB.fasl written ; compilation finished in 0:00:00.095 ------------------------------------------------------------------------ LinearCombinationUtilities is now explicitly exposed in frame initial LinearCombinationUtilities will be automatically loaded when needed from /var/aw/var/LatexWiki/LINCOMB.NRLIB/LINCOMB )compile EFACTOR.spad spad)abbrev package EFACTOR ExpressionFactorPolynomial ExpressionFactorPolynomial(R, F) : Exports == Implementation where R : Join(GcdDomain, Comparable, CharacteristicZero, RetractableTo Integer, LinearlyExplicitRingOver Integer) F : Join(AlgebraicallyClosedField, TranscendentalFunctionCategory, FunctionSpace R) UP ==> SparseUnivariatePolynomial F Exports ==> with factorPolynomial : UP -> Factored(UP) Implementation ==> add K ==> Kernel F Exp ==> IndexedExponents(K) MP ==> SparseMultivariatePolynomial(R, K) if R is Integer and F is Expression(Integer) then MULTFACT ==> MultivariateFactorize(K, Exp, R, MP) dummy := create()$SingletonAsOrderedSet factorPolynomial(p) == n := degree(p) cn := leadingCoefficient(p) cn = 0 => nilFactor(p, 1) p := p/cn cnp := cn::UP n < 1 => makeFR(cnp, []) n = 1 => makeFR(cnp, [["prime", p, 1]]) kl := kernels(coefficients(p)) for k in kl repeat not(is?(k, 'nthRoot) or is?(k, 'rootOf)) => "iterate" error "factorPolynomial: algebraic coefficients" xs := new()$Symbol xk := kernel(xs) pf := retract(eval(p, dummy, xk::F))@F fres1 := factor(numer(pf))$MULTFACT res : Factored(UP) := makeFR(cnp, []) for frec in factorList(fres1) repeat fr1 := frec.fctr degree(fr1, xk) < 1 => "iterate" frec.flg ~= "prime" => error "impossible" fru := numer(univariate(fr1::F, xk)) fr2 := fru/leadingCoefficient(fru) res := res*primeFactor(fr2, frec.xpnt) res else factorPolynomial(p) == error "factorPolynomial unimplemented" spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/4099479118789942951-25px002.spad using old system compiler. EFACTOR abbreviates package ExpressionFactorPolynomial ------------------------------------------------------------------------ initializing NRLIB EFACTOR for ExpressionFactorPolynomial compiling into NRLIB EFACTOR ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: F already in scope ****** Domain: F already in scope processing macro definition K ==> Kernel F processing macro definition Exp ==> IndexedExponents Kernel F processing macro definition MP ==> SparseMultivariatePolynomial(R,Kernel F) processing macro definition MULTFACT ==> MultivariateFactorize(Kernel F,IndexedExponents Kernel F,R,SparseMultivariatePolynomial(R,Kernel F)) compiling exported factorPolynomial : SparseUnivariatePolynomial F -> Factored SparseUnivariatePolynomial F Time: 0.09 SEC. compiling exported factorPolynomial : SparseUnivariatePolynomial F -> Factored SparseUnivariatePolynomial F EFACTOR;factorPolynomial;SupF;2 is replaced by errorfactorPolynomial unimplemented Time: 0 SEC. compiling exported factorPolynomial : SparseUnivariatePolynomial F -> Factored SparseUnivariatePolynomial F EFACTOR;factorPolynomial;SupF;3 is replaced by errorfactorPolynomial unimplemented Time: 0.01 SEC. (time taken in buildFunctor: 0) ;;; *** |ExpressionFactorPolynomial| REDEFINED ;;; *** |ExpressionFactorPolynomial| REDEFINED Time: 0 SEC. Warnings: [1] factorPolynomial: fctr has no value [2] factorPolynomial: flg has no value [3] factorPolynomial: xpnt has no value [4] factorPolynomial: res has no value Cumulative Statistics for Constructor ExpressionFactorPolynomial Time: 0.10 seconds finalizing NRLIB EFACTOR Processing ExpressionFactorPolynomial for Browser database: --->-->ExpressionFactorPolynomial(constructor): Not documented!!!! --->-->ExpressionFactorPolynomial((factorPolynomial ((Factored (SparseUnivariatePolynomial F)) (SparseUnivariatePolynomial F)))): Not documented!!!! --->-->ExpressionFactorPolynomial(): Missing Description ; compiling file "/var/aw/var/LatexWiki/EFACTOR.NRLIB/EFACTOR.lsp" (written 05 SEP 2014 03:02:57 AM): ; /var/aw/var/LatexWiki/EFACTOR.NRLIB/EFACTOR.fasl written ; compilation finished in 0:00:00.053 ------------------------------------------------------------------------ ExpressionFactorPolynomial is now explicitly exposed in frame initial ExpressionFactorPolynomial will be automatically loaded when needed from /var/aw/var/LatexWiki/EFACTOR.NRLIB/EFACTOR )compile PFUTIL.spad spad)abbrev package PFUTIL PartialFractionUtilities PartialFractionUtilities(F, UP) : Exports == Implementation where F : Field UP : UnivariatePolynomialCategory F Exports ==> with decompose : (UP, List(UP)) -> List(UP) ++ decompse(n, [p1, ..., pn]) returns numerators of partial ++ fraction decomposition of n divided by product of pi-s ++ Note: decomposition is assumed to have no whole part. Implementation ==> add FUP ==> Factored(UP) PFUP ==> PartialFraction(UP) decompose(nn : UP, dens : List(UP)) : List(UP) == fdens := [nilFactor(nden, 1)$FUP for nden in dens] nd := reduce(_*, fdens, 1) pfr := partialFraction(nn, nd)$PFUP wholePart(pfr) ~= 0 => error "decompose: wholePart(pfr) ~= 0" nfacs := numberOfFractionalTerms(pfr) res : List(UP) := [] for fden in fdens repeat for i in 1..nfacs repeat ct := nthFractionalTerm(pfr, i) cnum := firstNumer(ct) cden := firstDenom(ct) if cden = fden then res := cons(cnum, res) reverse!(res) spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/3035627895780083912-25px003.spad using old system compiler. PFUTIL abbreviates package PartialFractionUtilities ------------------------------------------------------------------------ initializing NRLIB PFUTIL for PartialFractionUtilities compiling into NRLIB PFUTIL processing macro definition FUP ==> Factored UP processing macro definition PFUP ==> PartialFraction UP compiling exported decompose : (UP,List UP) -> List UP Time: 0.02 SEC. (time taken in buildFunctor: 0) ;;; *** |PartialFractionUtilities| REDEFINED ;;; *** |PartialFractionUtilities| REDEFINED Time: 0 SEC. Warnings: [1] decompose: res has no value Cumulative Statistics for Constructor PartialFractionUtilities Time: 0.02 seconds finalizing NRLIB PFUTIL Processing PartialFractionUtilities for Browser database: --->-->PartialFractionUtilities(constructor): Not documented!!!! --------(decompose ((List UP) UP (List UP)))--------- --->-->PartialFractionUtilities((decompose ((List UP) UP (List UP)))): Improper first word in comments: decompse "decompse(\\spad{n},{} [\\spad{p1},{} ...,{} \\spad{pn}]) returns numerators of partial fraction decomposition of \\spad{n} divided by product of \\spad{pi}-\\spad{s} Note: decomposition is assumed to have no whole part." --->-->PartialFractionUtilities(): Missing Description ; compiling file "/var/aw/var/LatexWiki/PFUTIL.NRLIB/PFUTIL.lsp" (written 05 SEP 2014 03:02:59 AM): ; /var/aw/var/LatexWiki/PFUTIL.NRLIB/PFUTIL.fasl written ; compilation finished in 0:00:00.026 ------------------------------------------------------------------------ PartialFractionUtilities is now explicitly exposed in frame initial PartialFractionUtilities will be automatically loaded when needed from /var/aw/var/LatexWiki/PFUTIL.NRLIB/PFUTIL )compile INTPAR1.spad spad)abbrev package INTPAR1 ParametricTranscendentalIntegration ParametricTranscendentalIntegration(F, UP) : Exports == Implementation where F : Field UP : UnivariatePolynomialCategory F Z ==> Integer Q ==> Fraction Z RF ==> Fraction UP Param_Rec_F ==> Record(ratpart : F, coeffs : Vector F) L_Param_F ==> List Param_Rec_F Param_Rec_Q ==> Record(ratpart : RF, coeffs : Vector F) L_Param_Q ==> List Param_Rec_Q Param_Rec_QF ==> Record(logands : List F, basis : List Vector Q) Param_Rec_Q2 ==> Record(logands : List RF, basis : List Vector Q) -- L_Param_Q2 ==> List Param_Rec_Q2 Exports ==> with primextint : (UP -> UP, List F -> L_Param_F, Matrix F -> List Vector F, List RF) -> L_Param_Q ++ primext(', ext, csolve, [g1, ..., gn]) returns ++ a basis of solutions of the homogeneous system ++ \spad{h' + c1*g1 + ... + cn*gn = 0}. ++ Argument ext is an extended integration function on F. ++ csolve is solver over constants. expextint : (UP -> UP, (Z, List F) -> L_Param_F, Matrix F -> List Vector F, List RF) -> L_Param_Q ++ expextint(', rde, csolve, [g1, ..., gn]) returns ++ a basis of solution of the system homogeneous ++ \spad{h' + c1*g1 + ... + cn*gn = 0} ++ Argument foo is an parametric rde solver on F. ++ csolve is solver over constants. diffextint : (List UP -> L_Param_F, Matrix F -> List Vector F, List RF) -> L_Param_F ++ diffext(ext, csolve, [g1, ..., gn]) unkextint : (List F -> L_Param_F, Matrix F -> List Vector F, List RF) -> L_Param_F ++ unkext(ext, csolve, [g1, ..., gn]) monologint : (List UP, UP, UP -> UP, UP -> UP, Matrix F -> List Vector F) -> List Vector F ++ monologint(la, b, der) see Raab logextint : (UP -> UP, UP -> Factored(UP), Matrix F -> List Vector Q, List UP -> Param_Rec_Q2, List(RF)) -> Param_Rec_Q2 ++ logextint(der, ufactor, csolve, rec, [g1, ..., gn]) ++ returns [[u1, ..., um], bas] giving basis of solution of ++ the homogeneous systym ++ \spad{c1*g1 + ... + cn*gn + c_{n+1}u1'/u1 + ... c_{n+m}um'/um = 0} monologextint : (List UP, Matrix F -> List Vector Q, List F -> Param_Rec_QF) -> Param_Rec_Q2 ++ monologextint(lup, csolve, rec) is a helper for logextint Implementation ==> add N ==> NonNegativeInteger GP ==> LaurentPolynomial(F, UP) import from LinearCombinationUtilities(F, UP) import from LinearCombinationUtilities(Q, SparseUnivariatePolynomial(Q)) import from TranscendentalHermiteIntegration(F, UP) monologextint(lup, csolve, rec1) == n0 := #lup lc0 := [coefficient(p, 0) for p in lup] lup2 := [p - c0::UP for p in lup for c0 in lc0] m1 := matrix([lup2])$Matrix(UP) rs1 : Matrix F := reducedSystem(m1) cb := csolve(rs1) nlc0 := [lin_comb(bv, lc0) for bv in cb] (ll, bl) := rec1(nlc0) empty?(bl) => [[], []] n1 := #cb n2 := #ll n3 := n0 + n2 rbl := [new(n3, 0)$Vector(Q) for bv in bl] nl := [le::UP::RF for le in ll] for rbv in rbl for bv in bl repeat pv := lin_comb(bv(1..n1), cb) for i in 1..n0 repeat rbv(i) := pv(i) for i in n0+1..n3 for j in n1+1.. repeat rbv(i) := bv(j) [nl, rbl] logextint(der, ufactor, csolve, rec, lg) == empty?(lg) => [[], []] n0 := #lg lghr := [HermiteIntegrate(g, der) for g in lg] lg1 := [ghr.answer for ghr in lghr] m1 := matrix([lg1])$Matrix(RF) rs1 : Matrix UP := reducedSystem(m1) rs2 : Matrix F := reducedSystem(rs1) cb := csolve(rs2) lg2 := [ghr.logpart for ghr in lghr] lg3 := [lin_comb(bv, lg2) for bv in cb] lden1 := [denom(g) for g in lg3] (mbas, m2) := gcdDecomposition(vector(lden1))$GcdBasis(UP) n1 := #lg3 n2 := #mbas mbasl := entries(mbas) mbasfl := [[frr.factor for frr in factors(ufactor(mbasp))] for mbasp in mbasl] basl1 := reduce(concat, mbasfl, []) -- print(basl1::OutputForm) sl := [#fl1 for fl1 in mbasfl] n3 := reduce(_+, sl, 0) m3 := new(n3, n1 + n3, 0)$Matrix(UP) for i in 1..n1 for g in lg3 repeat fl : List(UP) := [] jl : List(Z) := [] sl1 := sl fl1 : List(UP) := [] j0 : Z := 1 j1 : Z := 0 for j in 1..n3 for fj in basl1 repeat j1 := j1 + 1 if j1 > first(sl1) then sl1 := rest(sl1) j1 := 0 j0 := j0 + 1 if m2(j0, i) = 1 then fl := cons(fj, fl) jl := cons(j, jl) nl := decompose(numer(g), fl)$PartialFractionUtilities(F, UP) for num in nl for j in jl repeat m3(j, i) := num -- print("m3 part 1 done"::OutputForm) lpc : List(RF) := [] for j in 1..n3 for bj in basl1 repeat dbj := der(bj) (q, r) := divide(dbj, bj) m3(j, n1 + j) := r lpc := cons(q::RF, lpc) lpc := reverse!(lpc) -- print(lpc::OutputForm) -- print(m3::OutputForm) rs3 := reducedSystem(m3) cb2 := csolve(rs3) cb3 := [bv(1..n1) for bv in cb2] cb3e := [bv(n1+1..n1+n3) for bv in cb2] -- print(cb3e::OutputForm) ncb := [lin_comb(bv, cb) for bv in cb3] nlpc := [lin_comb(bv, lpc) for bv in cb3e] lg4 := [ghr.polypart::RF + ghr.specpart for ghr in lghr] lg5 := [lin_comb(bv, lg4) + pc for bv in ncb for pc in nlpc] -- print(lg5::OutputForm) lrf : List(RF) := [] lg6 : List(RF) := [] for g in lg5 repeat den := denom(g) (q, r) := divide(numer(g), den) lrf := cons(r/den, lrf) lg6 := cons(q::RF, lg6) lrf := reverse!(lrf) lg6 := reverse!(lg6) m1 := matrix([lrf])$Matrix(RF) rs1 : Matrix UP := reducedSystem(m1) rs2 : Matrix F := reducedSystem(rs1) cb4 := csolve(rs2) ncb2 := [lin_comb(bv, ncb) for bv in cb4] cb4e := [lin_comb(bv, cb3e) for bv in cb4] -- print(cb4e::OutputForm) -- lg7 := [lin_comb(bv, lg6) + lin_comb(bv, nlpc) for bv in cb4] lg7 := [lin_comb(bv, lg6) for bv in cb4] n4 := #lg7 (flog, fbas) := rec([retract(g)@UP for g in lg7]) empty?(fbas) => [[], []] rbasl1 := [up::RF for up in basl1] nlog := concat(rbasl1, flog) n5 := #flog n6 := (n0 + n3 + n5)::N rbas := [new(n6, 0)$Vector(Q) for bv in fbas] -- print("[n0, n3, n4, n5, n6], cb4e"::OutputForm) -- print([n0, n3, n4, n5, n6]::OutputForm) -- print(cb4e::OutputForm) for rbv in rbas for bv in fbas repeat bv1 := bv(1..n4) pv := lin_comb(bv1, ncb2) for i in 1..n0 repeat rbv(i) := pv(i) pv := lin_comb(bv1, cb4e) for i in n0 + 1..n0 + n3 for j in 1..n3 repeat rbv(i) := pv(j) for i in n0 + n3 + 1..n6 for j in n4 + 1..n4 + n5 repeat rbv(i) := bv(j) [nlog, rbas] monologint(la : List UP, b : UP, der1 : UP -> UP, der2 : UP -> UP, csolve : Matrix F -> List Vector F ) : List Vector F == ee1 := extendedEuclidean(b, der1(b)) ee1.generator ~= 1 => error "impossible" c1 := ee1.coef2 lpi := [(a*c1) rem b for a in la] d2b := der2(b) ee2 := extendedEuclidean(b, differentiate(b)) ee2.generator ~= 1 => error "impossible" c2 := ee2.coef2 lp2i := [(differentiate(pi)*d2b*c2) rem b for pi in lpi] lqi := [der2(pi) - p2i for pi in lpi for p2i in lp2i] rs1 := reducedSystem(matrix([lqi])$Matrix(UP)) csolve(rs1) RF_to_GP(f : RF) : GP == (numer(f)::GP exquo denom(f)::GP)::GP primextint(der : UP -> UP, ext : List F -> L_Param_F, csolve : Matrix F -> List Vector F, lg : List RF) : L_Param_Q == empty?(lg) => [] n := #lg lghr := [HermiteIntegrate(g, der) for g in lg] -- print("after HermiteIntegrate"::OutputForm) -- print(lghr::OutputForm) lg1 := [ghr.logpart for ghr in lghr] m1 := matrix([lg1])$Matrix(RF) rs1 : Matrix UP := reducedSystem(m1) rs2 : Matrix F := reducedSystem(rs1) cb := csolve(rs2) -- print("after csolve"::OutputForm) -- print(cb::OutputForm) a1l := [ghr.answer for ghr in lghr] lba : List(RF) := [0 for bv in cb] lg2 := [ghr.polypart for ghr in lghr] vg2 := vector([lg2])$Vector(UP) ldg : List(N) := [degree(g2) for g2 in lg2] d := reduce(max, ldg) dk := retract(der(monomial(1, 1)$UP))@F cba := [0$F for bv in cb] nlba : List(RF) -- print("before for"::OutputForm) for j in d..0 by -1 repeat n1 := #cb lgj : List(F) := [] for i in 1..n repeat gi := vg2(i) gij : F := 0 if degree(gi) = j then gij := leadingCoefficient(gi) vg2(i) := reductum(gi) lgj := cons(gij, lgj) lgj := reverse!(lgj) -- print("lgj ="::OutputForm) -- print(lgj::OutputForm) lgj1 : List(F) := [] for bv in cb for aa in cba repeat ff := lin_comb(bv, lgj) + (j + 1)::F*dk*aa lgj1 := cons(ff, lgj1) lgj1 := reverse!(lgj1) -- print("before recursive call"::OutputForm) b2 := ext(cons(dk, lgj1)) -- print("after recursive call"::OutputForm) n1p := n1 + 1 -- print(n1::OutputForm) cb0 := [(be.coeffs)(2..n1p) for be in b2] -- print("got subparts"::OutputForm) ncb := [lin_comb(bv, cb) for bv in cb0] -- print("transformed subparts"::OutputForm) cba := [be.ratpart for be in b2] nlba := [lin_comb(bv, lba) + monomial(be.ratpart, j)$UP::RF + monomial((be.coeffs)(1)/(j + 1)::F, j + 1)$UP::RF for be in b2 for bv in cb0] -- print("j, nlba"::OutputForm) -- print(j::OutputForm) -- print(nlba::OutputForm) cb := ncb lba := nlba nlba := [ba - lin_comb(bv, a1l) for bv in cb for ba in lba] [[ba, bv] for bv in cb for ba in nlba] expextint(der : UP -> UP, rde : (Z, List F) -> L_Param_F, csolve : Matrix F -> List Vector F, lg : List RF) : L_Param_Q == empty?(lg) => [] lghr := [HermiteIntegrate(g, der) for g in lg] lg1 := [ghr.logpart for ghr in lghr] -- print(lg1::OutputForm) m1 := matrix([lg1])$Matrix(RF) rs1 : Matrix UP := reducedSystem(m1) rs2 : Matrix F := reducedSystem(rs1) cb := csolve(rs2) a1l := [ghr.answer for ghr in lghr] lba : List(GP) := [0 for bv in cb] lg2 := [ghr.polypart::GP + RF_to_GP(ghr.specpart) for ghr in lghr] vg2 := vector([lg2])$Vector(GP) ldg : List(Z) := [degree(g2) for g2 in lg2] d := reduce(max, ldg) cba := [0$F for bv in cb] j := d repeat last_iter : Boolean := true n1 := #cb lgj : List(F) := [] for i in 1..#vg2 repeat gi := vg2(i) gij : F := 0 if gi ~= 0 then last_iter := false if degree(gi) = j then gij := leadingCoefficient(gi) vg2(i) := reductum(gi) lgj := cons(gij, lgj) last_iter => break lgj := reverse!(lgj) lgj1 : List(F) := [] for bv in cb repeat ff := lin_comb(bv, lgj) lgj1 := cons(ff, lgj1) lgj1 := reverse!(lgj1) -- print("before rde"::OutputForm) -- print(j::OutputForm) -- print(lgj1::OutputForm) b2 := rde(j, lgj1) -- print("after rde"::OutputForm) -- print(b2::OutputForm) empty?(b2) => return [] ncb := [lin_comb(be.coeffs, cb) for be in b2] nlba := [lin_comb(be.coeffs, lba) + monomial(be.ratpart, j)$GP for be in b2] cb := ncb lba := nlba -- print("end of interation"::OutputForm) -- print(j::OutputForm) j := j - 1 lbar := [convert(ba)@RF - lin_comb(bv, a1l) for bv in cb for ba in lba] [[bar, bv] for bv in cb for bar in lbar] diffextint1(trim : RF -> UP, ext : List UP -> L_Param_F, csolve : Matrix(F) -> List Vector(F), lg : List RF) : L_Param_F == lup := [trim(g) for g in lg] lg1 := [g - up::RF for g in lg for up in lup] m1 := matrix([lg1])$Matrix(RF) rs1 : Matrix(UP) := reducedSystem(m1) rs2 : Matrix F := reducedSystem(rs1) cb := csolve(rs2) lup1 := [lin_comb(bv, lup) for bv in cb] res1 := ext(lup1) [[re.ratpart, lin_comb(re.coeffs, cb)] for re in res1] lin_part(f : RF) : UP == p := numer(f) quo denom(f) monomial(coefficient(p, 1), 1)$UP + coefficient(p, 0)::UP diffextint(ext, csolve, lg) == diffextint1((x : RF) : UP +-> lin_part(x), ext, csolve, lg) coeff0(f : RF) : UP == p := numer(f) quo denom(f) coefficient(p, 0)::UP unkextint(ext, csolve, lg) == ext1 := (lup : List UP) : L_Param_F +-> lf := [retract(p)@F for p in lup] ext(lf) diffextint1((x : RF) : UP +-> coeff0(x), ext1, csolve, lg) spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/1208729076270301804-25px004.spad using old system compiler. INTPAR1 abbreviates package ParametricTranscendentalIntegration ------------------------------------------------------------------------ initializing NRLIB INTPAR1 for ParametricTranscendentalIntegration compiling into NRLIB INTPAR1 processing macro definition N ==> NonNegativeInteger processing macro definition GP ==> LaurentPolynomial(F,UP) importing LinearCombinationUtilities(F,UP) importing LinearCombinationUtilities(Fraction Integer,SparseUnivariatePolynomial Fraction Integer) importing TranscendentalHermiteIntegration(F,UP) compiling exported monologextint : (List UP,Matrix F -> List Vector Fraction Integer,List F -> Record(logands: List F,basis: List Vector Fraction Integer)) -> Record(logands: List Fraction UP,basis: List Vector Fraction Integer) Time: 0.08 SEC. compiling exported logextint : (UP -> UP,UP -> Factored UP,Matrix F -> List Vector Fraction Integer,List UP -> Record(logands: List Fraction UP,basis: List Vector Fraction Integer),List Fraction UP) -> Record(logands: List Fraction UP,basis: List Vector Fraction Integer) Time: 0.35 SEC. compiling exported monologint : (List UP,UP,UP -> UP,UP -> UP,Matrix F -> List Vector F) -> List Vector F Time: 0.12 SEC. compiling local RF_to_GP : Fraction UP -> LaurentPolynomial(F,UP) Time: 0.01 SEC. compiling exported primextint : (UP -> UP,List F -> List Record(ratpart: F,coeffs: Vector F),Matrix F -> List Vector F,List Fraction UP) -> List Record(ratpart: Fraction UP,coeffs: Vector F) Time: 0.28 SEC. compiling exported expextint : (UP -> UP,(Integer,List F) -> List Record(ratpart: F,coeffs: Vector F),Matrix F -> List Vector F,List Fraction UP) -> List Record(ratpart: Fraction UP,coeffs: Vector F) Time: 0.13 SEC. compiling local diffextint1 : (Fraction UP -> UP,List UP -> List Record(ratpart: F,coeffs: Vector F),Matrix F -> List Vector F,List Fraction UP) -> List Record(ratpart: F,coeffs: Vector F) Time: 0.01 SEC. compiling local lin_part : Fraction UP -> UP Time: 0.01 SEC. compiling exported diffextint : (List UP -> List Record(ratpart: F,coeffs: Vector F),Matrix F -> List Vector F,List Fraction UP) -> List Record(ratpart: F,coeffs: Vector F) Time: 0 SEC. compiling local coeff0 : Fraction UP -> UP Time: 0.01 SEC. compiling exported unkextint : (List F -> List Record(ratpart: F,coeffs: Vector F),Matrix F -> List Vector F,List Fraction UP) -> List Record(ratpart: F,coeffs: Vector F) Time: 0.01 SEC. (time taken in buildFunctor: 0) ;;; *** |ParametricTranscendentalIntegration| REDEFINED ;;; *** |ParametricTranscendentalIntegration| REDEFINED Time: 0 SEC. Warnings: [1] monologextint: logands has no value [2] monologextint: basis has no value [3] logextint: answer has no value [4] logextint: logpart has no value [5] logextint: basis has no value [6] logextint: transform has no value [7] logextint: j0 has no value [8] logextint: fl has no value [9] logextint: jl has no value [10] logextint: quotient has no value [11] logextint: remainder has no value [12] logextint: polypart has no value [13] logextint: specpart has no value [14] logextint: logands has no value [15] monologint: generator has no value [16] monologint: coef2 has no value [17] primextint: logpart has no value [18] primextint: answer has no value [19] primextint: polypart has no value [20] primextint: gij has no value [21] primextint: coeffs has no value [22] primextint: ratpart has no value [23] expextint: logpart has no value [24] expextint: answer has no value [25] expextint: polypart has no value [26] expextint: specpart has no value [27] expextint: gij has no value [28] expextint: last_iter has no value [29] expextint: coeffs has no value [30] expextint: ratpart has no value [31] diffextint1: ratpart has no value [32] diffextint1: coeffs has no value Cumulative Statistics for Constructor ParametricTranscendentalIntegration Time: 1.01 seconds finalizing NRLIB INTPAR1 Processing ParametricTranscendentalIntegration for Browser database: --->-->ParametricTranscendentalIntegration(constructor): Not documented!!!! --------(primextint ((List (Record (: ratpart (Fraction UP)) (: coeffs (Vector F)))) (Mapping UP UP) (Mapping (List (Record (: ratpart F) (: coeffs (Vector F)))) (List F)) (Mapping (List (Vector F)) (Matrix F)) (List (Fraction UP))))--------- --->-->ParametricTranscendentalIntegration((primextint ((List (Record (: ratpart (Fraction UP)) (: coeffs (Vector F)))) (Mapping UP UP) (Mapping (List (Record (: ratpart F) (: coeffs (Vector F)))) (List F)) (Mapping (List (Vector F)) (Matrix F)) (List (Fraction UP))))): Improper first word in comments: primext "primext(',{} ext,{} csolve,{} [\\spad{g1},{} ...,{} \\spad{gn}]) returns a basis of solutions of the homogeneous system \\spad{h' + c1*g1 + ... + cn*gn = 0}. Argument ext is an extended integration function on \\spad{F}. csolve is solver over constants." --------(expextint ((List (Record (: ratpart (Fraction UP)) (: coeffs (Vector F)))) (Mapping UP UP) (Mapping (List (Record (: ratpart F) (: coeffs (Vector F)))) (Integer) (List F)) (Mapping (List (Vector F)) (Matrix F)) (List (Fraction UP))))--------- --------(diffextint ((List (Record (: ratpart F) (: coeffs (Vector F)))) (Mapping (List (Record (: ratpart F) (: coeffs (Vector F)))) (List UP)) (Mapping (List (Vector F)) (Matrix F)) (List (Fraction UP))))--------- --->-->ParametricTranscendentalIntegration((diffextint ((List (Record (: ratpart F) (: coeffs (Vector F)))) (Mapping (List (Record (: ratpart F) (: coeffs (Vector F)))) (List UP)) (Mapping (List (Vector F)) (Matrix F)) (List (Fraction UP))))): Improper first word in comments: diffext "diffext(\\spad{ext},{} \\spad{csolve},{} [\\spad{g1},{} ...,{} \\spad{gn}])" --------(unkextint ((List (Record (: ratpart F) (: coeffs (Vector F)))) (Mapping (List (Record (: ratpart F) (: coeffs (Vector F)))) (List F)) (Mapping (List (Vector F)) (Matrix F)) (List (Fraction UP))))--------- --->-->ParametricTranscendentalIntegration((unkextint ((List (Record (: ratpart F) (: coeffs (Vector F)))) (Mapping (List (Record (: ratpart F) (: coeffs (Vector F)))) (List F)) (Mapping (List (Vector F)) (Matrix F)) (List (Fraction UP))))): Improper first word in comments: unkext "unkext(\\spad{ext},{} \\spad{csolve},{} [\\spad{g1},{} ...,{} \\spad{gn}])" --------(monologint ((List (Vector F)) (List UP) UP (Mapping UP UP) (Mapping UP UP) (Mapping (List (Vector F)) (Matrix F))))--------- --------(logextint ((Record (: logands (List (Fraction UP))) (: basis (List (Vector (Fraction (Integer)))))) (Mapping UP UP) (Mapping (Factored UP) UP) (Mapping (List (Vector (Fraction (Integer)))) (Matrix F)) (Mapping (Record (: logands (List (Fraction UP))) (: basis (List (Vector (Fraction (Integer)))))) (List UP)) (List (Fraction UP))))--------- --------(monologextint ((Record (: logands (List (Fraction UP))) (: basis (List (Vector (Fraction (Integer)))))) (List UP) (Mapping (List (Vector (Fraction (Integer)))) (Matrix F)) (Mapping (Record (: logands (List F)) (: basis (List (Vector (Fraction (Integer)))))) (List F))))--------- --->-->ParametricTranscendentalIntegration(): Missing Description ; compiling file "/var/aw/var/LatexWiki/INTPAR1.NRLIB/INTPAR1.lsp" (written 05 SEP 2014 03:03:03 AM): ; /var/aw/var/LatexWiki/INTPAR1.NRLIB/INTPAR1.fasl written ; compilation finished in 0:00:00.354 ------------------------------------------------------------------------ ParametricTranscendentalIntegration is now explicitly exposed in frame initial ParametricTranscendentalIntegration will be automatically loaded when needed from /var/aw/var/LatexWiki/INTPAR1.NRLIB/INTPAR1 )compile RDEPAR.spad spad)abbrev package RDEPAR ParametricRischDE ParametricRischDE(R, F) : Exports == Implementation where R : Join(GcdDomain, Comparable, CharacteristicZero, RetractableTo Integer, LinearlyExplicitRingOver Integer) F : Join(TranscendentalFunctionCategory, AlgebraicallyClosedField, FunctionSpace R) N ==> NonNegativeInteger Z ==> Integer Q ==> Fraction(Integer) SE ==> Symbol LF ==> List F K ==> Kernel F LK ==> List K Param_Rec_F ==> Record(ratpart : F, coeffs : Vector F) L_Param_F ==> List Param_Rec_F Partial_F ==> Union(Param_Rec_F, "failed") Both_F ==> Record(particular : Partial_F, basis : L_Param_F) UP ==> SparseUnivariatePolynomial F Param_Rec_UP ==> Record(ratpart : UP, coeffs : Vector F) L_Param_UP ==> List Param_Rec_UP Param_Rec_QF ==> Record(logands : List F, basis : List Vector Q) Exports ==> with param_rde : (Z, F, LF, SE, LK, (LK, LF) -> L_Param_F, (LK, LF) -> Param_Rec_QF) -> L_Param_F ++ param_rde(n, f, lg, x, lk, ext, logi) finds basis of ++ solution to the equation ++ dy/dx + n df/dx y + c1 g1 + ... cn gn = 0 where ++ y is in field generated by lk and ci are ++ constants. param_rde : (Z, F, F, LF, SE, LK, (LK, LF) -> L_Param_F, (LK, LF) -> Param_Rec_QF) -> Both_F ++ param_rde(n, f, h, lg, x, lk, ext, logi) finds a particular ++ solution and basis of solutions to homogeneous ++ equation for equation ++ dy/dx + n df/dx y + c1 g1 + ... cn gn = h where ++ y is in field generated by lk and ci are ++ constants. param_rde2 : (F, LF, SE, LK, (LK, LF) -> L_Param_F, (LK, LF) -> Param_Rec_QF) -> L_Param_F do_spde1 : (UP, List(UP), Z, UP -> UP, Matrix(F) -> Matrix(F)) -> L_Param_UP Implementation ==> add RDEEF ==> ElementaryRischDE(R, F) MET ==> MonomialExtensionTools(F, UP) GP ==> LaurentPolynomial(F, UP) RF ==> Fraction UP Param_Rec_Q ==> Record(ratpart : RF, coeffs : Vector F) L_Param_Q ==> List Param_Rec_Q RSOL ==> Record(ans : UP, remainder : UP) DSOL ==> Record(ans : List(UP), acoeff : UP, eegen : UP, bpar : UP, lcpar : List UP, dpar : Z) P ==> SparseMultivariatePolynomial(R, K) import from Integer import from List(Integer) import from IntegrationTools(R, F) import from PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, P, F) import from LinearCombinationUtilities(F, UP) ALGOP := '%alg PRIM := 'prim do_param_rde : (F, LF, SE, K, LK, (LK, LF) -> L_Param_F, (LK, LF) -> Param_Rec_QF) -> L_Param_F -- param_rde2 : (F, LF, SE, LK, (LK, LF) -> L_Param_F, -- (LK, LF) -> Param_Rec_QF) -> L_Param_F get_denom(f : RF, lg : List RF, der : UP -> UP) : List(UP) == d := normalDenom(f, der)$MET e0 := lcm [denom(g) for g in lg]$List(UP) (e, s) := split(e0, der)$MET gg := gcd(d, e) h := (gcd(e, differentiate e) exquo gcd(gg, differentiate gg))::UP [d, h] Frec ==> Record(fctr : UP, xpnt : Z) normalize(f : RF, der : UP -> UP) : List(Frec) == d := normalDenom(f, der)$MET g := gcd(d, differentiate(d)) dd := gcd(d, g) d1 := (d exquo dd)::UP d2 := (denom(f) exquo d1)::UP eeu := extendedEuclidean(d2, d1, numer(f)) (a, b) := eeu::Record(coef1 : UP, coef2 : UP) zk := kernel(new()$Symbol)$K dd1 := der(d1) r := resultant(a - zk::F*dd1, d1) rql := get_rational_roots(r, zk )$FunctionSpaceRationalRoots(R, F) rl : List Frec := [] for rq in rql repeat if (mu := retractIfCan(rq)@Union(Z, "failed")) case Z then m := mu::Z if m > 0 then pi := gcd(a - m::F*dd1, d1) rl := cons([pi, m], rl) rl RF_to_GP(f : RF) : GP == (numer(f)::GP exquo denom(f)::GP)::GP do_spde1(b : UP, lc : List(UP), d : Z, der : UP -> UP, get_rs : Matrix(F) -> Matrix(F)) : L_Param_UP == import TranscendentalRischDE(F, UP) la := [SPDEnocancel1(b, c, d, der).ans for c in lc] lrem := [der(a1) + b*a1 - c for a1 in la for c in lc] rs1 : Matrix(F) := reducedSystem(matrix([lrem])) rs2 := get_rs(rs1) lkv := nullSpace(rs2) [[lin_comb(kv, la), kv] for kv in lkv] param_SPDE(a : UP, b : UP, lc : List(UP), d : Z, der : UP -> UP, get_rs : Matrix(F) -> Matrix(F), x : SE) : L_Param_UP == -- print("param_SPDE"::OutputForm) -- print(a::OutputForm) -- print(b::OutputForm) -- print(lc::OutputForm) -- print(d::OutputForm) dt := der monomial(1, 1) degree(a) = 0 => a ~= 1 => error "param_SPDE: degree(a) = 0 but a ~= 1" degt := degree(dt) - 1 base_case := dt = 1 b = 0 => error "param_SPDE: b = 0" base_case or degree(b) > max(0, degt) => do_spde1(b, lc, d, der, get_rs) error "param_SPDE: degree(a) = 0 unimplemented" n1 := #lc s1 := multi_SPDE(a, b, lc, d, der)$RDEaux(F) s1 case List(RSOL) => lrs := s1::List(RSOL) m1 := matrix([[rsol.remainder for rsol in lrs]])$Matrix(UP) rs1 : Matrix(F) := reducedSystem(m1) rs2 := get_rs(rs1) lkv := nullSpace(rs2) a1l := [rsol.ans for rsol in lrs] [[lin_comb(kv, a1l), kv] for kv in lkv] dres := s1::DSOL g := dres.eegen a := (a exquo g)::UP b := (dres.bpar exquo g)::UP aa := dres.acoeff oans := dres.ans lq := [divide(c, g) for c in dres.lcpar] lr := [q.remainder for q in lq] rs1 : Matrix(F) := reducedSystem(matrix([lr])$Matrix(UP)) rs2 := get_rs(rs1) lkv := nullSpace(rs2) empty?(lkv) => [] lc := [q.quotient for q in lq] nlc : List(UP) := [] for kv in lkv repeat nlc := cons(lin_comb(kv, lc), nlc) nlc := reverse!(nlc) n2 := #lkv s2 := param_SPDE(a, b, nlc, dres.dpar, der, get_rs, x) nres : L_Param_UP := [] for be in s2 repeat bv := lin_comb(be.coeffs, lkv) ans1 := lin_comb(bv, oans) + aa*be.ratpart nres := cons([ans1, bv], nres) reverse!(nres) integer_vector(v : Vector(Q)) : Union(Vector(Integer), "failed") == (nv, d) := splitDenominator(v)$CommonDenominator(Integer, Q, Vector(Q)) d ~= 1 => "failed" nv(1) ~= 1 => "failed" vector([retract(nv(i))@Integer for i in 1..#nv]) do_SPDE_prim0(b : F, lc : List(UP), lk : LK, ext : (LK, LF) -> L_Param_F, logi : (LK, LF) -> Param_Rec_QF, der : UP -> UP, get_rs : Matrix(F) -> Matrix(F), x : SE ) : L_Param_Q == -- print("do_SPDE_prim0"::OutputForm) n := #lc vg2 := vector([lc])$Vector(UP) cb := [new(n, 0)$Vector(F) for i in 1..n] for i in 1..n for bv in cb repeat bv(i) := 1 d := reduce(max, [degree(c) for c in lc], 0) dk := retract(der(monomial(1, 1)$UP))@F lba : List(RF) := [0 for bv in cb] cba := [0$F for bv in cb] for j in d..0 by -1 repeat -- print("for j, j ="::OutputForm) -- print(j::OutputForm) n1 := #cb lgj : List(F) := [] for i in 1..n repeat gi := vg2(i) gij : F := 0 if degree(gi) = j then gij := leadingCoefficient(gi) vg2(i) := reductum(gi) lgj := cons(gij, lgj) lgj := reverse!(lgj) lgj1 : List(F) := [] for bv in cb for aa in cba repeat ff := lin_comb(bv, lgj) + (j + 1)::F*dk*aa lgj1 := cons(ff, lgj1) lgj1 := reverse!(lgj1) s2f := param_rde2(b, lgj1, x, lk, ext, logi) ncb := [lin_comb(be.coeffs, cb) for be in s2f] cba := [be.ratpart for be in s2f] nlba := [lin_comb(be.coeffs, lba) + monomial(be.ratpart, j::N)$UP::RF for be in s2f] cb := ncb lba := nlba [[ba, bv] for ba in lba for bv in cb] do_SPDE_prim(a : UP, bbr : RF, lcr : List RF, lk : LK, ext : (LK, LF) -> L_Param_F, logi : (LK, LF) -> Param_Rec_QF, der : UP -> UP, get_rs : Matrix(F) -> Matrix(F), x : SE ) : L_Param_Q == -- print("do_SPDE_prim"::OutputForm) fp := retract(der(monomial(1, 1)))@F base_case : Boolean := fp = 1 b := retract(bbr)@UP lc : List(UP) := [retract(cr)@UP for cr in lcr] da := degree(a) db := degree(b) dc := reduce(max, [degree(c) for c in lc]$List(Z)) not(base_case) and da = 0 and db = 0 => b1 := retract(b)@F/retract(a)@F b1 = 0 => error "need call ext" (ll, bl) := logi(lk, [b1]) if not(empty?(bl)) then bv := first(bl) bvu := integer_vector(bv) if bvu case Vector(Integer) then error "need transform an call ext" do_SPDE_prim0(b1, lc, lk, ext, logi, der, get_rs, x) n := db > da => max(0, dc - db) max(0, dc - da + 1) if da = db + 1 then f0 := -leadingCoefficient(b)/leadingCoefficient(a) base_case => if (mu := retractIfCan(f0)@Union(Z, "failed")) case Z then n := max(n, mu::Z) r0 := dehomogenize(ext(lk, [-f0, fp])).particular if not(r0 case "failed") then mf : F := ((r0::Param_Rec_F).coeffs)(1) if (mu := retractIfCan(mf)@Union(Z, "failed")) case Z then n := max(n, mu::Z) if not base_case and da = db then f0 := -leadingCoefficient(b)/leadingCoefficient(a) b1 := f0*a + b if degree(b1) + 1 = da then f1 := -leadingCoefficient(b1)/leadingCoefficient(a) r0 := dehomogenize(ext(lk, [-f1, fp])).particular if not(r0 case "failed") then mf : F := ((r0::Param_Rec_F).coeffs)(1) if (mu := retractIfCan(mf)@Union(Z, "failed")) case Z then n := max(n, mu::Z) -- print("calling param_SPDE"::OutputForm) res1 := param_SPDE(a, b, lc, n, der, get_rs, x) [[re.ratpart::RF, re.coeffs] for re in res1] do_SPDE_exp0(a : F, b : F, lcr : List(GP), lk : LK, eta : F, ext : (LK, LF) -> L_Param_F, logi : (LK, LF) -> Param_Rec_QF, x : SE ) : L_Param_Q == n := #lcr vg2 := vector([lcr])$Vector(GP) cb := [new(n, 0)$Vector(F) for cr in lcr] for i in 1..n for bv in cb repeat bv(i) := 1 lba : List(GP) := [0 for cr in lcr] d := reduce(max, [degree(cr) for cr in lcr]) j := d f0 := b/a repeat last_iter : Boolean := true n1 := #cb lgj : List(F) := [] for i in 1..n repeat gi := vg2(i) gij : F := 0 if gi ~= 0 then last_iter := false if degree(gi) = j then gij := leadingCoefficient(gi) vg2(i) := reductum(gi) lgj := cons(gij, lgj) last_iter => break lgj := reverse!(lgj) lgj1 : List(F) := [] for bv in cb repeat ff := lin_comb(bv, lgj) lgj1 := cons(ff, lgj1) lgj1 := reverse!(lgj1) s2f := param_rde2(f0 + j::F*eta, lgj1, x, lk, ext, logi) ncb := [lin_comb(be.coeffs, cb) for be in s2f] nlba := [lin_comb(be.coeffs, lba) + monomial(be.ratpart, j)$GP for be in s2f] cb := ncb lba := nlba j := j - 1 [[convert(re)@RF, bv] for re in lba for bv in cb] exp_lower_bound(a : UP, b : GP, ob : Z, nc0 : Z, lk : LK, eta : F, logi : (LK, LF) -> Param_Rec_QF) : Z == ob < 0 => min(0, nc0 - ob) n0 := min(0, nc0) 0 < ob => n0 c0 := coefficient(b, 0)/coefficient(a, 0) (ll, bl) := logi(lk, [c0, eta]) empty?(bl) => n0 bv := first(bl) bvu := integer_vector(bv) bvu case "failed" => n0 min((bvu::Vector(Integer))(2), n0) exp_upper_bound(a : UP, b : UP, nc1 : Z, lk : LK, eta : F, logi : (LK, LF) -> Param_Rec_QF) : Z == -- print("exp_upper_bound"::OutputForm) -- print(a::OutputForm) -- print(b::OutputForm) da := degree(a) db := degree(b) da < db => nc1 - db n0 := max(0, nc1 - da) db < da => n0 c1 := leadingCoefficient(b)/leadingCoefficient(a) (ll, bl) := logi(lk, [c1, eta]) empty?(bl) => n0 bv := first(bl) bvu := integer_vector(bv) bvu case "failed" => n0 max((bvu::Vector(Integer))(2), n0) do_SPDE_exp(a : UP, bbr : RF, lcr : List RF, lk : LK, ext : (LK, LF) -> L_Param_F, logi : (LK, LF) -> Param_Rec_QF, der : UP -> UP, get_rs : Matrix(F) -> Matrix(F), x : SE) : L_Param_Q == -- print("do_SPDE_exp"::OutputForm) -- print(a::OutputForm) b := RF_to_GP(bbr) -- print(b::OutputForm) lc := [RF_to_GP(cr) for cr in lcr] -- print(lc::OutputForm) nb0 := order(b) nc0 := reduce(min, [order(c) for c in lc]) eta := retract((der(monomial(1, 1)) exquo monomial(1, 1))::UP)@F degree(a) = 0 and degree(b) = 0 and nb0 = 0 => do_SPDE_exp0(retract(a), retract(b), lc, lk, eta, ext, logi, x) n0 := exp_lower_bound(a, b, nb0, nc0, lk, eta, logi) if n0 < 0 then -- print("n0 < 0"::OutputForm) -- print(n0::OutputForm) b := b + (n0::F*eta)::GP*a::GP if nb0 < 0 then -- print("nb0 < 0"::OutputForm) -- print(nb0::OutputForm) t1 := monomial(1, (-nb0)::N)$UP b := t1::GP*b a := t1*a bu := retract(b)@UP m0 := min(0, n0) + min(0, nb0) if m0 < 0 then t1 := monomial(1, (-m0)::N)$GP lc := [t1::GP*c for c in lc] lcu := [retract(c)@UP for c in lc] nc1 := reduce(max, [degree(cu) for cu in lcu]) n1 := exp_upper_bound(a, bu, nc1, lk, eta, logi) res1 := param_SPDE(a, bu, lcu, n1, der, get_rs, x) tt : RF := n0 < 0 => (monomial(1, (-n0)::N)$UP)::RF 1 [[re.ratpart::RF/tt, re.coeffs] for re in res1] param_rde(m, f, g0, lg, x, lk, ext, logi) == lg1 := cons(-g0, lg) -- nn := #lg1 res1 := param_rde(m, f, lg1, x, lk, ext, logi) dehomogenize(res1) param_rde2(fp, lg, x, lk, ext, logi) == -- print("param_rde2"::OutputForm) -- print(fp::OutputForm) -- print(lg::OutputForm) -- print(x::OutputForm) -- print(lk::OutputForm) k := kmax(lk) -- print(k::OutputForm) lk := [k1 for k1 in lk | k1 ~= k] dku := univariate(differentiate(k::F, x), k) -- print(dku::OutputForm) denom(dku) ~= 1 => -- print("denom(dku) ~= 1"::OutputForm) [] dk := numer(dku) fpu := univariate(fp, k) denfp := denom(fpu) der1 : UP -> UP := z1 +-> differentiate(z1, (z2 : F) : F +-> differentiate(z2, x), dk) rl := normalize(fpu, der1) p : UP := 1 for re in rl repeat (pi, ni) := re fpu := fpu - ni::F*der1(pi)/pi p := pi^ni*p fp := multivariate(fpu, k) pf := multivariate(p::RF, k) lg1 := [pf*g for g in lg] -- print("calling do_param_rde") do_param_rde(fp, lg, x, k, lk, ext, logi) param_rde(m, f, lg, x, lk, ext, logi) == (fp := D(m*f, x)) = 0 => ext(lk, lg) k := kmax(lk) lk := [k1 for k1 in lk | k1 ~= k] has?(operator k, ALGOP)$BasicOperator => error "param_rde: ALGOP unimplemented" -- (fp, f, p) := normalize(m, f, fp, x, k)$RDEEF -- lg1 := [p*g for g in lg] do_param_rde(fp, lg, x, k, lk, ext, logi) -- assumes fp is weakly normalized and main kernel k -- is a monomial do_param_rde(fp, lg, x, k, lk, ext, logi) == -- print("do_param_rde"::OutputForm) -- print(lg::OutputForm) dku := univariate(differentiate(k::F, x), k) denom(dku) ~= 1 => [] dk := numer(dku) fpu := univariate(fp, k) denfp := denom(fpu) nfp := numer(fpu) lgu := [univariate(g, k) for g in lg] der1 : UP -> UP := z1 +-> differentiate(z1, (z2 : F) : F +-> differentiate(z2, x), dk) (d, h) := get_denom(fpu, lgu, der1) aa := d*h bbr := aa*fpu - (d*der1(h))::RF aa1 := aa*h lgu := [aa1*gu for gu in lgu] lgd := [decompose(gu, der1)$MET for gu in lgu] -- print(lgd::OutputForm) lnor1 := [dr.normal for dr in lgd] -- print(lnor1::OutputForm) rs1 : Matrix(UP) := reducedSystem(matrix([lnor1])$Matrix(RF)) rs2 : Matrix(F) := reducedSystem(rs1) get_rs := (m : Matrix(F)) : Matrix(F) +-> reducedSystem(m, [(ff : F) : F +-> differentiate(ff, x)] )$ConstantLinearDependence(R, F) rs3 := get_rs(rs2) lker := nullSpace(rs3) empty?(lker) => [] lgu := [dr.poly::RF + dr.special for dr in lgd] n1 := #lgu n2 := #lker lgu1 : List(RF) := [] for kv in lker repeat lgu1 := cons(lin_comb(kv, lgu), lgu1) lgu1 := reverse!(lgu1) -- print("lgu1 = "::OutputForm) -- print(lgu1::OutputForm) res1 := symbolIfCan(k) case SE or is?(k, 'log) or has?(operator k, PRIM)$BasicOperator => -- print("calling do_SPDE_prim"::OutputForm) do_SPDE_prim(aa, bbr, lgu1, lk, ext, logi, der1, get_rs, x) is?(k, 'exp) => do_SPDE_exp(aa, bbr, lgu1, lk, ext, logi, der1, get_rs, x) error "param_rde: unimplemented kernel" res2 : L_Param_F := [] for re in res1 repeat bv := lin_comb(re.coeffs, lker) ans1 : RF := -re.ratpart/h::RF anf : F := multivariate(ans1, k) res2 := cons([anf, bv], res2) reverse!(res2) spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/1790683143857815936-25px005.spad using old system compiler. RDEPAR abbreviates package ParametricRischDE ------------------------------------------------------------------------ initializing NRLIB RDEPAR for ParametricRischDE compiling into NRLIB RDEPAR ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: F already in scope ****** Domain: F already in scope processing macro definition RDEEF ==> ElementaryRischDE(R,F) processing macro definition MET ==> MonomialExtensionTools(F,SparseUnivariatePolynomial F) processing macro definition GP ==> LaurentPolynomial(F,SparseUnivariatePolynomial F) processing macro definition RF ==> Fraction SparseUnivariatePolynomial F processing macro definition Param_Rec_Q ==> Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F) processing macro definition L_Param_Q ==> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F) processing macro definition RSOL ==> Record(ans: SparseUnivariatePolynomial F,remainder: SparseUnivariatePolynomial F) processing macro definition DSOL ==> Record(ans: List SparseUnivariatePolynomial F,acoeff: SparseUnivariatePolynomial F,eegen: SparseUnivariatePolynomial F,bpar: SparseUnivariatePolynomial F,lcpar: List SparseUnivariatePolynomial F,dpar: Integer) processing macro definition P ==> SparseMultivariatePolynomial(R,Kernel F) importing Integer importing List Integer importing IntegrationTools(R,F) importing PolynomialCategoryQuotientFunctions(IndexedExponents Kernel F,Kernel F,R,SparseMultivariatePolynomial(R,Kernel F),F) importing LinearCombinationUtilities(F,SparseUnivariatePolynomial F) compiling local get_denom : (Fraction SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> List SparseUnivariatePolynomial F Time: 0.05 SEC. processing macro definition Frec ==> Record(fctr: SparseUnivariatePolynomial F,xpnt: Integer) compiling local normalize : (Fraction SparseUnivariatePolynomial F,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> List Record(fctr: SparseUnivariatePolynomial F,xpnt: Integer) Semantic Errors: [1] normalize: FunctionSpaceRationalRoots is not a known type Warnings: [1] get_denom: IN has no value [2] get_denom: g has no value [3] get_denom: normal has no value [4] get_denom: special has no value [5] normalize: coef1 has no value [6] normalize: coef2 has no value ****** comp fails at level 3 with expression: ****** error in function normalize (SEQ (LET |d| ((|Sel| (|MonomialExtensionTools| F (|SparseUnivariatePolynomial| F)) |normalDenom|) |f| |der|)) (LET |g| (|gcd| |d| (|differentiate| |d|))) (LET |dd| (|gcd| |d| |g|)) (LET |d1| (|::| (|exquo| |d| |dd|) (|SparseUnivariatePolynomial| F))) (LET |d2| (|::| (|exquo| (|denom| |f|) |d1|) (|SparseUnivariatePolynomial| F))) (LET |eeu| (|extendedEuclidean| |d2| |d1| (|numer| |f|))) (LET (|@Tuple| |a| |b|) (|::| |eeu| (|Record| (|:| |coef1| (|SparseUnivariatePolynomial| F)) (|:| |coef2| (|SparseUnivariatePolynomial| F))))) (LET |zk| ((|Sel| (|Kernel| F) |kernel|) ((|Sel| (|Symbol|) |new|)))) (LET |dd1| (|der| |d1|)) (LET |r| (|resultant| (- |a| (* (|::| |zk| F) |dd1|)) |d1|)) (LET |rql| | << | ((|Sel| (|FunctionSpaceRationalRoots| R F) |get_rational_roots|) |r| |zk|) | >> |) (LET (|:| |rl| (|List| (|Record| (|:| |fctr| (|SparseUnivariatePolynomial| F)) (|:| |xpnt| (|Integer|))))) (|construct|)) (REPEAT (IN |rq| |rql|) (SEQ (LET |mu| (@ (|retractIfCan| |rq|) (|Union| (|Integer|) "failed"))) (|exit| 1 (IF (|case| |mu| (|Integer|)) (SEQ (LET |m| (|::| |mu| (|Integer|))) (|exit| 1 (IF (> |m| 0) (SEQ (LET |pi| (|gcd| (- |a| (* (|::| |m| F) |dd1|)) |d1|)) (|exit| 1 (LET |rl| (|cons| (|construct| |pi| |m|) |rl|)))) |noBranch|))) |noBranch|)))) (|exit| 1 |rl|)) ****** level 3 ****** $x:= ((Sel (FunctionSpaceRationalRoots R F) get_rational_roots) r zk)$m:= $EmptyMode$f:= ((((|r| #) (|dd1| #) (|zk| #) (|b| #) ...) ((|get_denom| #) (|#| #) (< #) (<= #) ...))) >> Apparent user error: NoValueMode is an unknown mode )compile INTALG2.spad spad)abbrev package INTALG2 AlgebraicIntegrate2 AlgebraicIntegrate2(R0, F, R) : Exports == Implementation where R0 : Join(Comparable, IntegralDomain, RetractableTo Integer) F : Join(AlgebraicallyClosedField, FunctionSpace R0) R : FunctionFieldCategory(F, UP, UPUP) UP ==> SparseUnivariatePolynomial F UPUP ==> SparseUnivariatePolynomial QF QF ==> Fraction UP Param_Rec_F ==> Record(ratpart : F, coeffs : Vector(F)) L_Param_F ==> List Param_Rec_F Param_Rec_Q ==> Record(ratpart : QF, coeffs : Vector(F)) L_Param_Q ==> List Param_Rec_Q Param_Rec_R ==> Record(ratpart : R, coeffs : Vector(F)) L_Param_R ==> List Param_Rec_R Exports ==> with algextint : (UP -> UP, List(QF) -> L_Param_Q, (QF, List QF) -> L_Param_Q, Matrix F -> List Vector F, List(R)) -> L_Param_R ++ algextint(der, ext, rde, csolve, lg) algprimextint : (UP -> UP, List(QF) -> L_Param_Q, (QF, List QF) -> L_Param_Q, Matrix F -> List Vector F, List(R)) -> L_Param_R algexpextint : (UP -> UP, List(QF) -> L_Param_Q, (QF, List QF) -> L_Param_Q, Matrix F -> List Vector F, List(R)) -> L_Param_R Implementation ==> add GP ==> LaurentPolynomial(F, UP) AHR ==> AlgebraicHermiteIntegration(F, UP, UPUP, R) HER_Rec ==> Record(answer : R, logpart : R, polypart : R) Param_Rec_UP ==> Record(ratpart : UP, coeffs : Vector(F)) L_Param_UP ==> List Param_Rec_UP RSOL ==> Record(ans : UP, remainder : UP) exp_hermite1(f : R, der : UP -> UP) : HER_Rec == d := (c := integralCoordinates f).den v := c.num vp : Vector(QF) := new(n := #v, 0) vf : Vector(QF) := new(n, 0) for i in minIndex v .. maxIndex v repeat r := separate(qelt(v, i) / d)$GP qsetelt!(vf, i, r.fracPart) qsetelt!(vp, i, convert(r.polyPart)@QF) ff := represents(vf, w := integralBasis()) fp := represents(vp, w) h := HermiteIntegrate(ff, der)$AHR [h.answer, h.logpart, fp] prim_hermite1(f : R, der : UP -> UP) : HER_Rec == h := HermiteIntegrate(f, der)$AHR zero?(hh := h.logpart) => [h.answer, 0, 0] d := (c := integralCoordinates(hh)).den v := c.num vp : Vector(QF) := new(n := #v, 0) vf : Vector(QF) := new(n, 0) for i in minIndex v .. maxIndex v repeat r := divide(qelt(v, i), d)$UP qsetelt!(vf, i, (r.remainder)/d) qsetelt!(vp, i, (r.quotient)::QF) ff := represents(vf, w := integralBasis()) fp := represents(vp, w) [h.answer, ff, fp] list_hermite(lf : List(R), hermite1 : R -> HER_Rec) : List(HER_Rec) == [hermite1(f) for f in lf] import from LinearCombinationUtilities(F, UP) lin_comb2(v : Vector F, lr : List(R)) : R == res : R := 0 for i in 1..#v for r in lr repeat res := res + v(i)::UP::QF*r res )if false diagextint(ndw : UP, ddw : UP, lup : List UP, der : UP -> UP, rde : (QF, List QF) -> L_Param_Q, csolve : Matrix F -> List Vector F) : L_Param_UP == g := gcd(ndw, ddw) b := (ndw exquo g)::UP a := (ddw exquo g)::UP degree(a) < 1 => error "diagextint: degree(a) < 1" d := reduce(max, [degree(up) for up in lup], 0)$List(Integer) lup1 := [a*up for up in lup] s1 := multi_SPDE(a, b, lup1, d + 1, der)$RDEaux(F) s1 case List(RSOL) => lrs := s1::List(RSOL) m1 := matrix([[rsol.remainder for rsol in lrs]])$Matrix(UP) rs1 : Matrix(F) := reducedSystem(m1) lkv := csolve(rs1) a1l : List(UP) := [rsol.ans for rsol in lrs] [[-lin_comb(kv, a1l), kv] for kv in lkv] error "impossible" )endif diagextint(dden : UP, dm : Matrix(UP), w : Vector(R), lpv : List Vector QF, cb : List Vector F, ca0 : List R, ext : List QF -> L_Param_Q, rde : (QF, List QF) -> L_Param_Q) : L_Param_R == lup1 := [pv(1) for pv in lpv] lrf := [lin_comb(bv, lup1) for bv in cb] res0 := ext(lrf) ncb := [lin_comb(be0.coeffs, cb) for be0 in res0] wi := w(1) ca : List(R) := [be0.ratpart*wi for be0 in res0] for i in 2..#w repeat lup1 := [pv(i) for pv in lpv] lup2 := [lin_comb(bv, lup1) for bv in cb] res1 := rde(dm(i, i)/dden, lup2) ncb := [lin_comb(be.coeffs, cb) for be in res1] wi := w(i) nca := [be.ratpart*wi + lin_comb2(be.coeffs, ca) for be in res1] cb := ncb ca := nca [[ai + lin_comb2(bv, ca0), bv] for ai in ca for bv in cb] algprimextint(der, ext, rde, csolve, lg) == hermite1 : R -> HER_Rec := g +-> prim_hermite1(g, der) lh := list_hermite(lg, hermite1) lg1 := [h.logpart for h in lh] m1 : Matrix(QF) := reducedSystem(matrix([lg1])) rs1 : Matrix UP := reducedSystem(m1)$QF rs2 : Matrix F := reducedSystem(rs1) cb := csolve(rs2) empty?(cb) => [] lp1 := [h.polypart for h in lh] lpv : List Vector(UP) := [] for p1 in lp1 repeat ci := integralCoordinates(p1) ci.den ~= 1 => error "ci.den ~= 1" lpv := cons(ci.num, lpv) lpv := reverse!(lpv) w := integralBasis()$R dm := integralDerivationMatrix(der)$R dden := dm.den ca0 := [h.answer for h in lh] w(1) = 1 and diagonal?(dm.num) => lpv2 : List(Vector(QF)) := [] for pv in lpv repeat n := #pv pv2 := new(n, 0)$Vector(QF) for i in 1..n repeat pv2(i) := pv(i)::QF lpv2 := cons(pv2, lpv2) lpv2 := reverse!(lpv2) )if false -- Why this fails to compile ???? for pv in lpv repeat c1 := (xp : UP) : QF +-> xp::QF pv2 : Vector(QF):= map(c1, pv )$VectorFunctions2(Vector(UP), Vector(QF)) lpv2 := cons(pv2, lpv2) lpv2 := reverse!(lpv2) )endif )if false -- Why this fails to compile ???? lpv2 := [(map((xp : UP) : QF +-> xp::QF, pv )$VectorFunctions2(Vector(UP), Vector(QF)) )@Vector(QF) for pv in lpv] )endif diagextint(dden, dm.num, w, lpv2, cb, ca0, ext, rde) )if false diagextint := dm.den lup1 := [pv(1) for pv in lpv] lrf := [lin_comb(bv, lup1)::QF for bv in cb] res0 := ext(lrf) ncb := [lin_comb(be0.coeffs, cb) for be0 in res0] wi := w(1) ca : List(R) := [be0.ratpart*wi for be0 in res0] for i in 2..#w repeat lup1 := [pv(i) for pv in lpv] lup2 := [lin_comb(bv, lup1) for bv in cb] -- res1 := diagextint((dm.num)(i, i), dden, -- lup2, der, rde, csolve) res1 := rde((dm.num)(i, i)/dden, ncb := [lin_comb(be.coeffs, cb) for be in res1] wi := w(i) nca := [be.ratpart::QF*wi + lin_comb2(be.coeffs, ca) for be in res1] cb := ncb ca := nca [[ai, bv] for ai in ca for bv in cb] )endif -- for i in 1..#w for wi in w repeat error "unimplemented" algexpextint(der, ext, rde, csolve, lg) == hermite1 : R -> HER_Rec := g +-> exp_hermite1(g, der) lh := list_hermite(lg, hermite1) lg1 := [h.logpart for h in lh] m1 : Matrix(QF) := reducedSystem(matrix([lg1])) rs1 : Matrix UP := reducedSystem(m1)$QF rs2 : Matrix F := reducedSystem(rs1) cb := csolve(rs2) empty?(cb) => [] lp1 := [h.polypart for h in lh] lpv : List Vector(QF) := [] w := integralBasis()$R n := #w for p1 in lp1 repeat ci := integralCoordinates(p1) deni := ci.den pv1 := ci.num pv2 := new(n, 0)$Vector(QF) for i in 1..n repeat pv2(i) := pv1(i)/deni lpv := cons(pv2, lpv) lpv := reverse!(lpv) dm := integralDerivationMatrix(der)$R dden := dm.den ca0 := [h.answer for h in lh] w(1) = 1 and diagonal?(dm.num) => diagextint(dden, dm.num, w, lpv, cb, ca0, ext, rde) error "unimplemented" algextint(der, ext, rde, csolve, lg) == x' := der(x := monomial(1, 1)$UP) zero? degree(x') => algprimextint(der, ext, rde, csolve, lg) ((xx := x' exquo x) case UP) and (retractIfCan(xx::UP)@Union(F, "failed") case F) => algexpextint(der, ext, rde, csolve, lg) error "algextint: unhandled case" spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/586361080427230223-25px006.spad using old system compiler. INTALG2 abbreviates package AlgebraicIntegrate2 ------------------------------------------------------------------------ initializing NRLIB INTALG2 for AlgebraicIntegrate2 compiling into NRLIB INTALG2 ****** Domain: R0 already in scope ****** Domain: R0 already in scope ****** Domain: F already in scope processing macro definition GP ==> LaurentPolynomial(F,SparseUnivariatePolynomial F) processing macro definition AHR ==> AlgebraicHermiteIntegration(F,SparseUnivariatePolynomial F,SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F,R) processing macro definition HER_Rec ==> Record(answer: R,logpart: R,polypart: R) processing macro definition Param_Rec_UP ==> Record(ratpart: SparseUnivariatePolynomial F,coeffs: Vector F) processing macro definition L_Param_UP ==> List Record(ratpart: SparseUnivariatePolynomial F,coeffs: Vector F) processing macro definition RSOL ==> Record(ans: SparseUnivariatePolynomial F,remainder: SparseUnivariatePolynomial F) compiling local exp_hermite1 : (R,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> Record(answer: R,logpart: R,polypart: R) Time: 0.03 SEC. compiling local prim_hermite1 : (R,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> Record(answer: R,logpart: R,polypart: R) Time: 0.01 SEC. compiling local list_hermite : (List R,R -> Record(answer: R,logpart: R,polypart: R)) -> List Record(answer: R,logpart: R,polypart: R) Time: 0 SEC. importing LinearCombinationUtilities(F,SparseUnivariatePolynomial F) compiling local lin_comb2 : (Vector F,List R) -> R Time: 0.09 SEC. compiling local diagextint : (SparseUnivariatePolynomial F,Matrix SparseUnivariatePolynomial F,Vector R,List Vector Fraction SparseUnivariatePolynomial F,List Vector F,List R,List Fraction SparseUnivariatePolynomial F -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),(Fraction SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F) -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F)) -> List Record(ratpart: R,coeffs: Vector F) Time: 0.32 SEC. compiling exported algprimextint : (SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),(Fraction SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F) -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),Matrix F -> List Vector F,List R) -> List Record(ratpart: R,coeffs: Vector F) Time: 0.14 SEC. compiling exported algexpextint : (SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),(Fraction SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F) -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),Matrix F -> List Vector F,List R) -> List Record(ratpart: R,coeffs: Vector F) Time: 0.14 SEC. compiling exported algextint : (SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),(Fraction SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F) -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),Matrix F -> List Vector F,List R) -> List Record(ratpart: R,coeffs: Vector F) Time: 0.02 SEC. (time taken in buildFunctor: 0) ;;; *** |AlgebraicIntegrate2| REDEFINED ;;; *** |AlgebraicIntegrate2| REDEFINED Time: 0 SEC. Warnings: [1] exp_hermite1: num has no value [2] exp_hermite1: fracPart has no value [3] exp_hermite1: polyPart has no value [4] exp_hermite1: answer has no value [5] exp_hermite1: logpart has no value [6] prim_hermite1: logpart has no value [7] prim_hermite1: answer has no value [8] prim_hermite1: num has no value [9] prim_hermite1: remainder has no value [10] prim_hermite1: quotient has no value [11] diagextint: coeffs has no value [12] diagextint: ratpart has no value [13] algprimextint: logpart has no value [14] algprimextint: polypart has no value [15] algprimextint: den has no value [16] algprimextint: num has no value [17] algprimextint: lpv has no value [18] algprimextint: answer has no value [19] algexpextint: logpart has no value [20] algexpextint: polypart has no value [21] algexpextint: den has no value [22] algexpextint: num has no value [23] algexpextint: answer has no value Cumulative Statistics for Constructor AlgebraicIntegrate2 Time: 0.75 seconds finalizing NRLIB INTALG2 Processing AlgebraicIntegrate2 for Browser database: --->-->AlgebraicIntegrate2(constructor): Not documented!!!! --------(algextint ((List (Record (: ratpart R) (: coeffs (Vector F)))) (Mapping (SparseUnivariatePolynomial F) (SparseUnivariatePolynomial F)) (Mapping (List (Record (: ratpart (Fraction (SparseUnivariatePolynomial F))) (: coeffs (Vector F)))) (List (Fraction (SparseUnivariatePolynomial F)))) (Mapping (List (Record (: ratpart (Fraction (SparseUnivariatePolynomial F))) (: coeffs (Vector F)))) (Fraction (SparseUnivariatePolynomial F)) (List (Fraction (SparseUnivariatePolynomial F)))) (Mapping (List (Vector F)) (Matrix F)) (List R)))--------- --->-->AlgebraicIntegrate2((algprimextint ((List (Record (: ratpart R) (: coeffs (Vector F)))) (Mapping (SparseUnivariatePolynomial F) (SparseUnivariatePolynomial F)) (Mapping (List (Record (: ratpart (Fraction (SparseUnivariatePolynomial F))) (: coeffs (Vector F)))) (List (Fraction (SparseUnivariatePolynomial F)))) (Mapping (List (Record (: ratpart (Fraction (SparseUnivariatePolynomial F))) (: coeffs (Vector F)))) (Fraction (SparseUnivariatePolynomial F)) (List (Fraction (SparseUnivariatePolynomial F)))) (Mapping (List (Vector F)) (Matrix F)) (List R)))): Not documented!!!! --->-->AlgebraicIntegrate2((algexpextint ((List (Record (: ratpart R) (: coeffs (Vector F)))) (Mapping (SparseUnivariatePolynomial F) (SparseUnivariatePolynomial F)) (Mapping (List (Record (: ratpart (Fraction (SparseUnivariatePolynomial F))) (: coeffs (Vector F)))) (List (Fraction (SparseUnivariatePolynomial F)))) (Mapping (List (Record (: ratpart (Fraction (SparseUnivariatePolynomial F))) (: coeffs (Vector F)))) (Fraction (SparseUnivariatePolynomial F)) (List (Fraction (SparseUnivariatePolynomial F)))) (Mapping (List (Vector F)) (Matrix F)) (List R)))): Not documented!!!! --->-->AlgebraicIntegrate2(): Missing Description ; compiling file "/var/aw/var/LatexWiki/INTALG2.NRLIB/INTALG2.lsp" (written 05 SEP 2014 03:03:06 AM): ; /var/aw/var/LatexWiki/INTALG2.NRLIB/INTALG2.fasl written ; compilation finished in 0:00:00.122 ------------------------------------------------------------------------ AlgebraicIntegrate2 is now explicitly exposed in frame initial AlgebraicIntegrate2 will be automatically loaded when needed from /var/aw/var/LatexWiki/INTALG2.NRLIB/INTALG2 )compile INTAF.spad spad)abbrev package INTAF AlgebraicIntegration ++ Mixed algebraic integration; ++ Author: Manuel Bronstein ++ Date Created: 12 October 1988 ++ Date Last Updated: 4 June 1988 ++ Description: ++ This package provides functions for the integration of ++ algebraic integrands over transcendental functions; AlgebraicIntegration(R, F) : Exports == Implementation where R : Join(Comparable, IntegralDomain) F : Join(AlgebraicallyClosedField, FunctionSpace R) SY ==> Symbol N ==> NonNegativeInteger K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) UP ==> SparseUnivariatePolynomial F RF ==> Fraction UP UPUP==> SparseUnivariatePolynomial RF IR ==> IntegrationResult F IR2 ==> IntegrationResultFunctions2(curve, F) Param_Rec_F ==> Record(ratpart : F, coeffs : Vector F) L_Param_F ==> List Param_Rec_F Param_Rec_Q ==> Record(ratpart : RF, coeffs : Vector(F)) L_Param_Q ==> List Param_Rec_Q ALG ==> AlgebraicIntegrate(R, F, UP, UPUP, curve) FAIL==> error "failed - cannot handle that integrand" Exports ==> with algint : (F, K, K, UP -> UP) -> IR ++ algint(f, x, y, d) returns the integral of \spad{f(x, y)dx} ++ where y is an algebraic function of x; ++ d is the derivation to use on \spad{k[x]}. algextint : (K, K, UP -> UP, List(RF) -> L_Param_Q, (RF, List(RF)) -> L_Param_Q, Matrix F -> List Vector F, List(F)) -> L_Param_F ++ algextint(x, y, d, ext, rde, csolve, [g1, ..., gn]) returns ++ \spad{[h, [c1, ..., cn]]} such that \spad= dh/dx + sum(ci gi) ++ and dci/dx = 0, if such \spad{[h, [c1, ..., cn]]} exist, ++ "failed" otherwise. Implementation ==> add import from ChangeOfVariable(F, UP, UPUP) import from PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, P, F) rootintegrate : (F, K, K, UP -> UP) -> IR algintegrate : (F, K, K, UP -> UP) -> IR UPUP2F : (UPUP, RF, K, K) -> F F2UPUP : (F, K, K, UP) -> UPUP UP2UPUP : (UP, K) -> UPUP F2UPUP(f, kx, k, p) == UP2UPUP(univariate(f, k, p), kx) rootintegrate(f, t, k, derivation) == r1 := mkIntegral(modulus := UP2UPUP(p := minPoly k, t)) f1 := F2UPUP(f, t, k, p) monomial(inv(r1.coef), 1) r := radPoly(r1.poly)::Record(radicand : RF, deg : N) q := retract(r.radicand) curve := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg) map(x1+->UPUP2F(lift x1, r1.coef, t, k), algintegrate(reduce f1, derivation)$ALG)$IR2 algintegrate(f, t, k, derivation) == r1 := mkIntegral(modulus := UP2UPUP(p := minPoly k, t)) f1 := F2UPUP(f, t, k, p) monomial(inv(r1.coef), 1) modulus := UP2UPUP(p := minPoly k, t) curve := AlgebraicFunctionField(F, UP, UPUP, r1.poly) map(x1+->UPUP2F(lift x1, r1.coef, t, k), algintegrate(reduce f1, derivation)$ALG)$IR2 Curv_Rec ==> Record(funs : List(UPUP), irec : Record(coef : Fraction(UP), poly : UPUP), curve_dom : FunctionFieldCategory(F, UP, UPUP)) rootcurve(lf : List(F), t : K, k : K, derivation : UP -> UP ) : Curv_Rec == r1 := mkIntegral(modulus := UP2UPUP(p := minPoly k, t)) mon1 := monomial(inv(r1.coef), 1)$UPUP lf1 := [F2UPUP(f, t, k, p) mon1 for f in lf]$List(UPUP) r := radPoly(r1.poly)::Record(radicand : RF, deg : N) q := retract(r.radicand) curve := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg) [lf1, r1, curve]$Curv_Rec algcurve(lf : List(F), t : K, k : K, derivation : UP -> UP ) : Curv_Rec == r1 := mkIntegral(modulus := UP2UPUP(p := minPoly k, t)) mon1 := monomial(inv(r1.coef), 1)$UPUP lf1 := [F2UPUP(f, t, k, p) mon1 for f in lf]$List(UPUP) curve := AlgebraicFunctionField(F, UP, UPUP, r1.poly) [lf1, r1, curve]$Curv_Rec algextint(t, y, der, ext, rde, csolve, lg) == c_rec := is?(y, 'nthRoot) => rootcurve(lg, t, y, der) is?(y, 'rootOf) => algcurve(lg, t, y, der) FAIL curve := c_rec.curve_dom red : UPUP -> curve:= reduce$curve cc := c_rec.irec.coef res1 := algextint(der, ext, rde, csolve, [red(f) for f in c_rec.funs] )$AlgebraicIntegrate2(R, F, curve) [[UPUP2F(lift(be.ratpart), cc, t, y), be.coeffs] for be in res1] UP2UPUP(p, k) == map(x1+->univariate(x1, k), p)$SparseUnivariatePolynomialFunctions2(F, RF) UPUP2F(p, cf, t, k) == map((x1 : RF) : F+->multivariate(x1, t), p)$SparseUnivariatePolynomialFunctions2(RF, F) (multivariate(cf, t) * k::F) algint(f, t, y, derivation) == is?(y, 'nthRoot) => rootintegrate(f, t, y, derivation) is?(y, 'rootOf) => algintegrate(f, t, y, derivation) FAIL spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/8454275636552109700-25px007.spad using old system compiler. INTAF abbreviates package AlgebraicIntegration ------------------------------------------------------------------------ initializing NRLIB INTAF for AlgebraicIntegration compiling into NRLIB INTAF ****** Domain: R already in scope ****** Domain: F already in scope importing ChangeOfVariable(F,SparseUnivariatePolynomial F,SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F) importing PolynomialCategoryQuotientFunctions(IndexedExponents Kernel F,Kernel F,R,SparseMultivariatePolynomial(R,Kernel F),F) compiling local F2UPUP : (F,Kernel F,Kernel F,SparseUnivariatePolynomial F) -> SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F Time: 0.01 SEC. compiling local rootintegrate : (F,Kernel F,Kernel F,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> IntegrationResult F Time: 0.04 SEC. compiling local algintegrate : (F,Kernel F,Kernel F,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> IntegrationResult F Time: 0.02 SEC. processing macro definition Curv_Rec ==> Record(funs: List SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F,irec: Record(coef: Fraction SparseUnivariatePolynomial F,poly: SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F),curve_dom: FunctionFieldCategory(F,SparseUnivariatePolynomial F,SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F)) compiling local rootcurve : (List F,Kernel F,Kernel F,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> Record(funs: List SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F,irec: Record(coef: Fraction SparseUnivariatePolynomial F,poly: SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F),curve_dom: FunctionFieldCategory(F,SparseUnivariatePolynomial F,SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F)) Time: 0.04 SEC. compiling local algcurve : (List F,Kernel F,Kernel F,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> Record(funs: List SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F,irec: Record(coef: Fraction SparseUnivariatePolynomial F,poly: SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F),curve_dom: FunctionFieldCategory(F,SparseUnivariatePolynomial F,SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F)) Time: 0.01 SEC. compiling exported algextint : (Kernel F,Kernel F,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),(Fraction SparseUnivariatePolynomial F,List Fraction SparseUnivariatePolynomial F) -> List Record(ratpart: Fraction SparseUnivariatePolynomial F,coeffs: Vector F),Matrix F -> List Vector F,List F) -> List Record(ratpart: F,coeffs: Vector F) Time: 0.03 SEC. compiling local UP2UPUP : (SparseUnivariatePolynomial F,Kernel F) -> SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F Time: 0 SEC. compiling local UPUP2F : (SparseUnivariatePolynomial Fraction SparseUnivariatePolynomial F,Fraction SparseUnivariatePolynomial F,Kernel F,Kernel F) -> F Time: 0.01 SEC. compiling exported algint : (F,Kernel F,Kernel F,SparseUnivariatePolynomial F -> SparseUnivariatePolynomial F) -> IntegrationResult F Time: 0 SEC. (time taken in buildFunctor: 0) ;;; *** |AlgebraicIntegration| REDEFINED ;;; *** |AlgebraicIntegration| REDEFINED Time: 0 SEC. Warnings: [1] rootintegrate: coef has no value [2] rootintegrate: poly has no value [3] rootintegrate: radicand has no value [4] rootintegrate: deg has no value [5] algintegrate: coef has no value [6] algintegrate: poly has no value [7] rootcurve: coef has no value [8] rootcurve: poly has no value [9] rootcurve: radicand has no value [10] rootcurve: deg has no value [11] algcurve: coef has no value [12] algcurve: poly has no value [13] algextint: curve_dom has no value [14] algextint: irec has no value [15] algextint: funs has no value [16] algextint: ratpart has no value [17] algextint: coeffs has no value Cumulative Statistics for Constructor AlgebraicIntegration Time: 0.16 seconds finalizing NRLIB INTAF Processing AlgebraicIntegration for Browser database: --------constructor--------- --------(algint ((IntegrationResult F) F (Kernel F) (Kernel F) (Mapping (SparseUnivariatePolynomial F) (SparseUnivariatePolynomial F))))--------- --------(algextint ((List (Record (: ratpart F) (: coeffs (Vector F)))) (Kernel F) (Kernel F) (Mapping (SparseUnivariatePolynomial F) (SparseUnivariatePolynomial F)) (Mapping (List (Record (: ratpart (Fraction (SparseUnivariatePolynomial F))) (: coeffs (Vector F)))) (List (Fraction (SparseUnivariatePolynomial F)))) (Mapping (List (Record (: ratpart (Fraction (SparseUnivariatePolynomial F))) (: coeffs (Vector F)))) (Fraction (SparseUnivariatePolynomial F)) (List (Fraction (SparseUnivariatePolynomial F)))) (Mapping (List (Vector F)) (Matrix F)) (List F)))--------- ; compiling file "/var/aw/var/LatexWiki/INTAF.NRLIB/INTAF.lsp" (written 05 SEP 2014 03:03:09 AM): ; /var/aw/var/LatexWiki/INTAF.NRLIB/INTAF.fasl written ; compilation finished in 0:00:00.105 ------------------------------------------------------------------------ AlgebraicIntegration is now explicitly exposed in frame initial AlgebraicIntegration will be automatically loaded when needed from /var/aw/var/LatexWiki/INTAF.NRLIB/INTAF )compile INTPAR2.spad spad)abbrev package INTPAR2 ParametricIntegration ParametricIntegration(R, F) : Exports == Implementation where R : Join(GcdDomain, Comparable, CharacteristicZero, RetractableTo Integer, LinearlyExplicitRingOver Integer) F : Join(AlgebraicallyClosedField, TranscendentalFunctionCategory, FunctionSpace R) Q ==> Fraction(Integer) SE ==> Symbol K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) UP ==> SparseUnivariatePolynomial F RF ==> Fraction UP OPDIFF ==> '%diff Partial_C ==> Union(Vector F, "failed") Both_C ==> Record(particular : Partial_C, basis : List Vector F) Param_Rec_F ==> Record(ratpart : F, coeffs : Vector F) L_Param_F ==> List Param_Rec_F Partial_F ==> Union(Param_Rec_F, "failed") Both_F ==> Record(particular : Partial_F, basis : L_Param_F) Param_Rec_Q ==> Record(ratpart : RF, coeffs : Vector F) L_Param_Q ==> List Param_Rec_Q Param_Rec_QF ==> Record(logands : List F, basis : List Vector Q) Param_Rec_Q2 ==> Record(logands : List RF, basis : List Vector Q) Exports ==> with extendedint : (SE, List K, List F) -> L_Param_F ++ extendedint(x, [k1, ..., kn], [g1, ..., gn]) returns ++ a basis of the homogeneous system ++ \spad{dh/dx + c1*g1 + ... + cn*gn = 0}. Solutions are ++ in the field generated by k1, ..., kn. extendedint : (F, SE, List K, List F) -> Both_F extendedint : (F, SE, List F) -> Both_F ++ extendedint(f, x, [g1, ..., gn]) returns ++ solution of the system \spad= dh/dx + c1*g1 + ... + cn*gn and ++ and a basis of the associated homogeneous system ++ \spad{dh/dx + c1*g1 + ... + cn*gn = 0}. Solutions are ++ in the field generated by kernels of f and g1, ..., gn. logextint : (SE, List K, List F) -> Param_Rec_QF ++ logextint(x, lk, lg) returns [[u1, ..., um], bas] giving ++ basis of solution of ++ the homogeneous systym ++ \spad{c1*g1 + ... + cn*gn + c_{n+1}u1'/u1 + ... c_{n+m}um'/um = 0} algextint : (SE, K, List K, List F) -> L_Param_F primextint : (SE, K, List K, List F) -> L_Param_F expextint : (SE, K, List K, List F) -> L_Param_F diffextint : (SE, K, List K, List F) -> L_Param_F unkextint : (SE, K, List K, List F) -> L_Param_F Implementation ==> add import from IntegrationTools(R, F) import from AlgebraicManipulations(R, F) import from FunctionSpacePrimitiveElement(R, F) import from PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, P, F) PRIM := 'prim ALGOP := '%alg alglfint : (F, K, List K, SE) -> IR algprimint : (F, K, K, SE) -> IR algexpint : (F, K, K, SE) -> IR primint : (F, SE, K) -> IR expint : (F, SE, K) -> IR lambint : (F, SE, K) -> IR tanint : (F, SE, K) -> IR prim? : (K, SE) -> Boolean isx? : (F, SE) -> Boolean addx : (IR, F) -> IR cfind : (F, LLG) -> F lfintegrate0 : (F, SE) -> IR unknownint : (F, SE) -> IR diffint : (F, SE, K) -> IR tryChangeVar : (F, K, SE) -> Union(IR, "failed") prim?(k, x) == is?(k, 'log) or has?(operator k, PRIM) csolve2(m : Matrix(F)) : List Vector Q == n := nrows(m) v := new(n, 0)$Vector(F) (solveLinearlyOverQ(m, v)$IntegerLinearDependence(F)).basis EFACT ==> ExpressionFactorPolynomial(R, F) primlogextint(x : SE, k : K, l : List(K), lg : List F) : Param_Rec_QF == rec1 := (lg1 : List F) : Param_Rec_QF +-> logextint(x, l, lg1) rec := (lg2 : List UP) : Param_Rec_Q2 +-> monologextint(lg2, csolve2, rec1 )$ParametricTranscendentalIntegration(F, UP) d1 := (x1 : F) : F +-> differentiate(x1, x) der := (x1 : UP) : UP+->differentiate(x1, d1, differentiate(k::F, x)::UP) uf : UP -> Factored(UP) := factorPolynomial$EFACT (ll, bl) := logextint(der, uf, csolve2, rec, [univariate(g, k) for g in lg] )$ParametricTranscendentalIntegration(F, UP) [[multivariate(le, k) for le in ll], bl] explogextint1(lg : List(UP), eta : F, rec1 : List F -> Param_Rec_QF ) : Param_Rec_Q2 == lg1 := concat(lg, eta::UP) (ll, bl) := monologextint(lg1, csolve2, rec1 )$ParametricTranscendentalIntegration(F, UP) ll1 := cons(monomial(1, 1)$UP::RF, ll) [ll1, bl] explogextint(x : SE, k : K, l : List(K), lg : List F) : Param_Rec_QF == eta := differentiate(first(argument(k)), x) d1 := (x1 : F) : F +-> differentiate(x1, x) der := (x1 : UP) : UP +-> differentiate(x1, d1, monomial(eta, 1)$UP) rec1 := (lg1 : List F) : Param_Rec_QF +-> logextint(x, l, lg1) rec := (lg2 : List UP) : Param_Rec_Q2 +-> explogextint1(lg2, eta, rec1) uf : UP -> Factored(UP) := factorPolynomial$EFACT (ll, bl) := logextint(der, uf, csolve2, rec, [univariate(g, k) for g in lg] )$ParametricTranscendentalIntegration(F, UP) [[multivariate(le, k) for le in ll], bl] logextint(x, l, lg) == empty?(l) => cb := csolve2(matrix([lg])$Matrix(F)) [[], cb] k := kmax(l) l := [k1 for k1 in l | k1 ~= k] symbolIfCan(k) case SE or prim?(k, x) => primlogextint(x, k, l, lg) is?(k, 'exp) => explogextint(x, k, l, lg) error "logextint: unhandled kernel" extendedint(f : F, x : SE, lg : List F) == l := varselect(tower(cons(x::F, cons(f, lg))), x) extendedint(f, x, l, lg) extendedint(x : SE, l : List K, lg : List F) == empty?(l) => cb := nullSpace(matrix([lg])$Matrix(F)) [[0, kv] for kv in cb] k := kmax(l) l := [k1 for k1 in l | k1 ~= k] symbolIfCan(k) case SE or prim?(k, x) => primextint(x, k, l, lg) is?(k, 'exp) => expextint(x, k, l, lg) has?(operator k, ALGOP) => algextint(x, k, l, lg) is?(k, OPDIFF) => diffextint(x, k, l, lg) unkextint(x, k, l, lg) extendedint(f, x, l, lg) == dehomogenize(extendedint(x, l, cons(-f, lg)) )$LinearCombinationUtilities(F, UP) csolve1(m : Matrix F, d1 : F -> F) : List Vector(F) == (solveLinearOverConstants(m, new(nrows(m), 0)$Vector(F), [d1] )$ConstantLinearDependence(R, F)).basis wrapfn(fn : List F -> L_Param_F, k : K) : List RF -> L_Param_Q == (lrf : List RF) : L_Param_Q +-> lf := [multivariate(rf, k) for rf in lrf] r1 := fn(lf) la := [univariate(be.ratpart, k) for be in r1] [[a, be.coeffs] for a in la for be in r1] algextint(x, k, l, lg) == k1 := kmax(l) l := [k2 for k2 in l | k2 ~= k1] symbolIfCan(k1) case SE or prim?(k1, x) or is?(k1, 'exp) => d1 := (x1 : F) : F +-> differentiate(x1, x) cs1 := (x2 : Matrix F) : List Vector(F) +-> csolve1(x2, d1) dk : UP := is?(k1, 'exp) => monomial(differentiate(first argument k1, x), 1) differentiate(k1::F, x)::UP der1 := (x1 : UP) : UP+->differentiate(x1, d1, dk) ext1 := (x3 : List F) : L_Param_F +-> extendedint(x, cons(k1, l), x3) ext2 := (x4 : List K, x3 : List F) : L_Param_F +-> extendedint(x, x4, x3) logi := (x2 : List K, x3 : List F) : Param_Rec_QF +-> logextint(x, x2, x3) rde1 := (x6 : F, x3 : List F) : L_Param_F +-> param_rde2(x6, x3, x, cons(k1, l), ext2, logi )$ParametricRischDE(R, F) -- ext := (x4 : List RF) : L_Param_Q +-> -- primextint(der1, ext1, cs1, x4 -- )$ParametricTranscendentalIntegration(F, UP) rde2 := (x5 : RF, x4 : List RF) : L_Param_Q +-> wrapfn((x3 : List F) : L_Param_F +-> rde1(multivariate(x5, k1), x3), k1)(x4) algextint(k1, k, der1, wrapfn(ext1, k1), rde2, cs1, lg )$AlgebraicIntegration(R, F) has?(operator k1, ALGOP) => error "algextint: nested algop" -- algalgextint(f, x, k, k1, lg) error "algextint unimplemented kernel" primextint(x, k, l, lu) == d1 := (x1 : F) : F +-> differentiate(x1, x) der := (x1 : UP) : UP+->differentiate(x1, d1, differentiate(k::F, x)::UP) ext := (x3 : List F) : L_Param_F +-> extendedint(x, l, x3) cs1 := (x4 : Matrix F) : List Vector F +-> csolve1(x4, d1) res1 := primextint(der, ext, cs1, [univariate(u, k) for u in lu] )$ParametricTranscendentalIntegration(F, UP) [[multivariate(si.ratpart, k), si.coeffs] for si in res1] expextint(x, k, l, lu) == eta := first argument k d1 := (x1 : F) : F +-> differentiate(x1, x) der := (x1 : UP) : UP+->differentiate(x1, d1, monomial(differentiate(eta, x), 1)) ext := (x2 : List K, x3 : List F) : L_Param_F +-> extendedint(x, x2, x3) logi := (x2 : List K, x3 : List F) : Param_Rec_QF +-> logextint(x, x2, x3) rde := (x1 : Integer, x3 : List F) : L_Param_F +-> param_rde(x1, eta, x3, x, l, ext, logi )$ParametricRischDE(R, F) cs1 := (x4 : Matrix F) : List Vector F +-> csolve1(x4, d1) res1 := expextint(der, rde, cs1, [univariate(u, k) for u in lu] )$ParametricTranscendentalIntegration(F, UP) [[multivariate(si.ratpart, k), si.coeffs] for si in res1] constant_subspace(b : List Vector(F), d1 : (F -> F) ) : Record(transform : Matrix(F), basis : List Vector(F)) == empty?(b) => [new(0, 0, 1), []] nc := #first(b) nr := #b m0b := matrix([entries(bv) for bv in b])$Matrix(F) m0 := horizConcat(m0b, scalarMatrix(nr, 1)) m1 := rowEchelon(m0) k : Integer := 1 lpiv : List Integer := [] lri : List Integer := [] for i in 1..nrows(m1) repeat k1 := k while m1(i, k1) = 0 and k1 <= nc repeat k1 := k1 + 1 nc < k1 => "iterate" k := k1 m1(i, k) ~= 1 => error "impossible" lri := cons(i, lri) lpiv := cons(k, lpiv) lpiv := cons(nc + 1, lpiv) lpiv := reverse!(lpiv) lr : List Vector F := [] ll : List List F := [] for i in lri repeat lpiv1 := lpiv kk := first(lpiv1) lpiv1 := rest(lpiv1) ll1 : List F := [] for j in 1..nc repeat j = kk => kk := first(lpiv1) lpiv1 := rest(lpiv1) ll1 := cons(d1(m1(i, j)), ll1) ll := cons(ll1, ll) lr := cons(row(m1, i), lr) m := transpose(matrix(ll)) v0 := new(ncols(m), 0)$Vector(F) s1 := (solveLinearOverConstants(m, v0, [d1] )$ConstantLinearDependence(R, F)).basis l3 := [lin_comb(bv, lr)$LinearCombinationUtilities(F, UP) for bv in s1] m2 := matrix([[bv(i) for i in nc+1..ncols(m0)] for bv in l3]) [m2, [bv(1..nc) for bv in l3]] -- [[bv(nc + 1), bv(1..nc)] for bv in l3] diffextint1(lup : List UP, x : SE, k : K, lk : List K, csolve : Matrix(F) -> List Vector(F)) : L_Param_F == import from LinearCombinationUtilities(F, UP) args := argument(k) #args ~= 3 => error "internal error, k is not a diff" arg3 := args(3) (da3 := differentiate(arg3, x)) = 0 => m1 : Matrix(UP) := matrix([lup]) rs1 : Matrix(F) := reducedSystem(m1) b0 := csolve(rs1) [[0, bv] for bv in b0] lg0 := [coefficient(up, 0) for up in lup] lg1 := [coefficient(up, 1) for up in lup] k1 := eval(args(1), retract(args(2))@K, arg3) dv := new()$Symbol dvf := dv::F lg2 := [eval(g/da3, k1, dvf) for g in lg1] lek := [eval(ki::F, k1, dvf) for ki in lk] lk1 := varselect(tower(cons(dvf, append(lg2, lek))), dv) res1 := extendedint(dv, lk1, lg2) empty?(res1) => [] lca := [be.ratpart for be in res1] cb := [be.coeffs for be in res1] lg3 := [lin_comb(bv, lg1) for bv in cb] nlg0 := [lin_comb(bv, lg0) for bv in cb] lca := map((x1 : F) : F+->eval(x1, kernel(dv), k1::F), lca) nlg : List(F) := [] for ca in lca for g1 in lg3 for g0 in nlg0 repeat du := univariate(differentiate(ca, x), k) nu := numer(du) denom(du) ~= 1 or degree(nu) > 1 => return [] g + coefficient(nu, 1) ~= 0 => return [] ng := g0 + coefficient(nu, 0) nlg := cons(ng, nlg) nlg := reverse!(nlg) res2 := extendedint(x, lk, nlg) cb0 := [be.coeffs for be in res2] d1 := (x1 : F) : F +-> differentiate(x1, x) (m2, cb1) := constant_subspace(cb0, d1) empty?(cb1) => [] v2 := vector([be.ratpart for be in res2])$Vector(F) v3 := m2*v2 ncb := [lin_comb(bv, cb) for bv in cb1] nlca := [lin_comb(bv, lca) for bv in cb1] [[v3(i) + ba, bv] for i in 1..#v3 for ba in nlca for bv in ncb] diffextint(x, k, l, lg) == d1 := (x1 : F) : F +-> differentiate(x1, x) cs1 := (x4 : Matrix F) : List Vector F +-> csolve1(x4, d1) diffi1 := (x1 : List UP) : L_Param_F +-> diffextint1(x1, x, k, l, cs1) diffextint(diffi1, cs1, [univariate(u, k) for u in lg] )$ParametricTranscendentalIntegration(F, UP) unkextint(x, k, l, lg) == d1 := (x1 : F) : F +-> differentiate(x1, x) cs1 := (x4 : Matrix F) : List Vector F +-> csolve1(x4, d1) ext := (x3 : List F) : L_Param_F +-> extendedint(x, l, x3) unkextint(ext, cs1, [univariate(u, k) for u in lg] )$ParametricTranscendentalIntegration(F, UP) spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/954272669615211323-25px008.spad using old system compiler. INTPAR2 abbreviates package ParametricIntegration ------------------------------------------------------------------------ initializing NRLIB INTPAR2 for ParametricIntegration compiling into NRLIB INTPAR2 ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: F already in scope ****** Domain: F already in scope importing IntegrationTools(R,F) importing AlgebraicManipulations(R,F) importing FunctionSpacePrimitiveElement(R,F) importing PolynomialCategoryQuotientFunctions(IndexedExponents Kernel F,Kernel F,R,SparseMultivariatePolynomial(R,Kernel F),F) compiling local prim? : (Kernel F,Symbol) -> Boolean Time: 0 SEC. compiling local csolve2 : Matrix F -> List Vector Fraction Integer Time: 0.01 SEC. processing macro definition EFACT ==> ExpressionFactorPolynomial(R,F) compiling local primlogextint : (Symbol,Kernel F,List Kernel F,List F) -> Record(logands: List F,basis: List Vector Fraction Integer) Time: 0.03 SEC. compiling local explogextint1 : (List SparseUnivariatePolynomial F,F,List F -> Record(logands: List F,basis: List Vector Fraction Integer)) -> Record(logands: List Fraction SparseUnivariatePolynomial F,basis: List Vector Fraction Integer) Time: 0.02 SEC. compiling local explogextint : (Symbol,Kernel F,List Kernel F,List F) -> Record(logands: List F,basis: List Vector Fraction Integer) Time: 0.02 SEC. compiling exported logextint : (Symbol,List Kernel F,List F) -> Record(logands: List F,basis: List Vector Fraction Integer) Time: 0 SEC. compiling exported extendedint : (F,Symbol,List F) -> Record(particular: Union(Record(ratpart: F,coeffs: Vector F),failed),basis: List Record(ratpart: F,coeffs: Vector F)) Time: 0.01 SEC. compiling exported extendedint : (Symbol,List Kernel F,List F) -> List Record(ratpart: F,coeffs: Vector F) Time: 0.01 SEC. compiling exported extendedint : (F,Symbol,List Kernel F,List F) -> Record(particular: Union(Record(ratpart: F,coeffs: Vector F),failed),basis: List Record(ratpart: F,coeffs: Vector F)) Time: 0.01 SEC. compiling local csolve1 : (Matrix F,F -> F) -> List Vector F ****** comp fails at level 4 with expression: ****** error in function csolve1 (((|Sel| (|ConstantLinearDependence| R F) |solveLinearOverConstants|) |m| ((|Sel| (|Vector| F) |new|) (|nrows| |m|) 0) (|construct| |d1|)) |basis|) ****** level 4 ******$x:= d1 $m:= (List (Fraction (Integer)))$f:= ((((|d1| # #) (|m| # . #1=#) (|csolve1| #) (|m| . #1#) ...) ((|explogextint| #) (|explogextint1| #)) ((|explogextint1| #) (|copy| #) (|setelt!| #) (|basis| #) ...) ((|primlogextint| #) (|copy| #) (|setelt!| #) (|basis| #) ...) ...)) >> Apparent user error: Cannot coerce d1 of mode (Mapping F F) to mode (Fraction (Integer)) )compile INTEF.spad spad)abbrev package INTEF ElementaryIntegration ++ Integration of elementary functions ++ Author: Manuel Bronstein ++ Date Created: 1 February 1988 ++ Date Last Updated: 24 October 1995 ++ Description: ++ This package provides functions for integration, limited integration, ++ extended integration and the risch differential equation for ++ elemntary functions. ++ Keywords: elementary, function, integration. ++ Examples: )r INTEF INPUT ElementaryIntegration(R, F) : Exports == Implementation where R : Join(GcdDomain, Comparable, CharacteristicZero, RetractableTo Integer, LinearlyExplicitRingOver Integer) F : Join(AlgebraicallyClosedField, TranscendentalFunctionCategory, FunctionSpace R) SE ==> Symbol K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) UP ==> SparseUnivariatePolynomial F RF ==> Fraction UP IR ==> IntegrationResult F IRRF ==> IntegrationResult RF FF ==> Record(ratpart : RF, coeff : RF) LLG ==> List Record(coeff : F, logand : F) Ext_Rec ==> Record(ratpart:F, coeff:F) U2 ==> Union(Ext_Rec, "failed") U3 ==> Union(Record(mainpart:F, limitedlogs:LLG), "failed") ANS ==> Record(special : F, integrand : F) PSOL==> Record(ans : F, right : F, sol? : Boolean) PSOL2 ==> Record(ans : F, right : F, primpart : F, sol? : Boolean) FAIL==> error "failed - cannot handle that integrand" OPDIFF ==> '%diff Param_Rec_F ==> Record(ratpart : F, coeffs : Vector F) Exports ==> with lfextendedint : (F, SE, F) -> U2 ++ lfextendedint(f, x, g) returns functions \spad{[h, c]} such that ++ \spad{dh/dx = f - cg}, if (h, c) exist, "failed" otherwise. lflogint : (F, SE, List F) -> U3 ++ lflimitedint(f, x, [g1, ..., gn]) returns functions \spad{[h, [[ci, gi]]]} ++ such that the gi's are among \spad{[g1, ..., gn]}, and ++ \spad{d(h+sum(ci log(gi)))/dx = f}, if possible, "failed" otherwise. lfinfieldint : (F, SE) -> Union(F, "failed") ++ lfinfieldint(f, x) returns a function g such that \spad{dg/dx = f} ++ if g exists, "failed" otherwise. lfintegrate : (F, SE) -> IR ++ lfintegrate(f, x) = g such that \spad{dg/dx = f}. Implementation ==> add import from IntegrationTools(R, F) import from ElementaryRischDE(R, F) import from RationalIntegration(F, UP) import from AlgebraicIntegration(R, F) import from AlgebraicManipulations(R, F) import from ElementaryRischDESystem(R, F) import from TranscendentalIntegration(F, UP) import from PureAlgebraicIntegration(R, F, F) import from IntegrationResultFunctions2(F, F) import from IntegrationResultFunctions2(RF, F) import from FunctionSpacePrimitiveElement(R, F) import from PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, P, F) import from ParametricIntegration(R, F) PRIM := 'prim ALGOP := '%alg alglfint : (F, K, List K, SE) -> IR algprimint : (F, K, K, SE) -> IR algexpint : (F, K, K, SE) -> IR primint : (F, SE, K) -> IR expint : (F, SE, K) -> IR lambint : (F, SE, K) -> IR tanint : (F, SE, K) -> IR prim? : (K, SE) -> Boolean isx? : (F, SE) -> Boolean addx : (IR, F) -> IR cfind : (F, LLG) -> F lfintegrate0 : (F, SE) -> IR unknownint : (F, SE) -> IR diffint : (F, SE, K) -> IR tryChangeVar : (F, K, SE) -> Union(IR, "failed") prim?(k, x) == is?(k, 'log) or has?(operator k, PRIM) lambint(f, x, k) == eta' := differentiate(eta := first argument k, x) dfac := (monomial(1, 1)$UP + 1$UP) derivative : UP -> UP := (x1 : UP) : UP +-> dfac*differentiate(x1, (x2 : F) : F +-> differentiate(x2, x), 0$UP) + differentiate(x1, (x2 : F) : F +-> 0, monomial(eta'/eta, 1)) -- XXX FIXME need proper extended integration function extint : F -> U2 := f +-> lfextendedint(f, x, eta'/eta) r1 := lambintegrate((monomial(1, 1)$UP + 1$UP)*univariate(f, k), eta'/eta, (x1 : F) : F +-> differentiate(x1, x), derivative, extint, (x2 : F) : IR +-> lfintegrate(x2, x)) map((x1 : RF) : F +-> multivariate(x1, k), r1.answer) + r1.a0 tanint(f, x, k) == eta' := differentiate(eta := first argument k, x) r1 := tanintegrate(univariate(f, k), (x1 : UP) : UP +-> differentiate(x1, (x2 : F) : F +-> differentiate(x2, x), monomial(eta', 2) + eta'::UP), (x6 : Integer, x2 : F, x3 : F) : Union(List F, "failed") +-> rischDEsys(x6, 2 * eta, x2, x3, x, (x4 : F, x5 : List F) : U3 +-> lflogint(x4, x, x5), (x4 : F, x5 : F) : U2 +-> lfextendedint(x4, x, x5))) map((x1 : RF) : F +-> multivariate(x1, k), r1.answer) + _ lfintegrate(r1.a0, x) -- tries various tricks since the integrand contains something not elementary unknownint(f, x) == (da := differentiate(a := denom(f)::F, x)) ~= 0 and zero? differentiate(c := numer(f)::F / da, x) => (c * log a)::IR mkAnswer(0, empty(), [[f, x::F]]) LOG ==> Record(scalar : Fraction(Integer), coeff : UP, logand : UP) dummy := create()$SingletonAsOrderedSet diffint1(f : F, x : SE, k : K) : Union(IR, "failed") == fu := univariate(f, k) denom(fu) ~= 1 => "failed" nfu := numer(fu) degree(nfu) ~= 1 => "failed" f1 := leadingCoefficient(nfu) args := argument(k) #args ~= 3 => error "internal error, k is not a diff" arg3 := args(3) da3 := differentiate(arg3, x) da3 = 0 => "failed" k1 := eval(args(1), retract(args(2))@K, arg3) dv := new()$Symbol f2 := eval(f1/da3, k1, dv::F) nres := lfintegrate(f2, dv) not(empty?(notelem(nres))) => "failed" logs := logpart(nres) nlogs : List(LOG) := [] nrat : F := 0 for lg in logs repeat cpol := lg.coeff every?((x1 : F) : Boolean +-> D(x1, x) = 0, coefficients(cpol)) => nlogs := cons(lg, nlogs) degree(cpol) ~= 1 => -- print(lg::OutputForm) return "failed" alpha := -coefficient(cpol, 0)/leadingCoefficient(cpol) nrat := nrat + lg.scalar*alpha* log(retract(eval(lg.logand, dummy, alpha))@F) nres2 := mkAnswer(ratpart(nres) + nrat, reverse!(nlogs), []) nres3 := map((x1 : F) : F+->eval(x1, kernel(dv), k1::F), nres2) f3 := f - differentiate(nres3, x) -- should not happen in differential transcendental -- case, but in principle possible member?(k, kernels(f3)) => "failed" nres3 + lfintegrate(f3, x) diffint(f, x, k) == (r1 := diffint1(f, x, k)) case IR => r1 mkAnswer(0, empty(), [[f, x::F]]) isx?(f, x) == (k := retractIfCan(f)@Union(K, "failed")) case "failed" => false (r := symbolIfCan(k::K)) case "failed" => false r::SE = x Alg_Rec ==> Record(fun : F, nroot : F, npow1 : Integer, npow2 : Integer) alg_split_root0(f : F, r : K, n : Integer) : List F == n = 2 => ef := eval(f, r, -(r::F)) f0 := (f + ef)/(2::F) f1 := (f - ef)/(2::F) member?(r, kernels(f0)) => print(f0::OutputForm) error "alg_split_root0: roots did not cancel 1" f1 := f1/(r::F) member?(r, kernels(f1)) => print(f1::OutputForm) error "alg_split_root0: roots did not cancel 2" [f0, f1] q := univariate(f, r, minPoly r )$PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, P, F) [coefficient(q, i) for i in 0..(n-1)] alg_split_roots(f : F, r1 : K, r2 : K) : List(Alg_Rec) == a1 := argument(r1) a2 := argument(r2) n1 : Integer := retract(a1(2))@Integer b1 := a1(1) n2 : Integer := retract(a2(2))@Integer b2 := a2(1) l1 : List F := alg_split_root0(f, r1, n1) res := empty()$List(Alg_Rec) pow1 : Integer rop := operator(r1) for f1 in l1 for pow1 in (0$Integer)..(n1 - 1) repeat f1 = 0 => "iterate" l2 := alg_split_root0(f1, r2, n2) g1 := gcd(n1, pow1)$Integer nn1 := (n1 exquo g1)::Integer np1 := (pow1 exquo g1)::Integer for f2 in l2 for pow2 in 0..(n2 - 1) repeat f2 = 0 => "iterate" g2 := gcd(n2, pow2) nn2 := (n2 exquo g2)::Integer np2 := (pow2 exquo g2)::Integer nn := lcm(nn1, nn2) bb1 := b1^(np1*(nn exquo nn1)::Integer) bb2 := b2^(np2*(nn exquo nn2)::Integer) nr := kernel(rop, [bb1*bb2, nn::F]) nrr := pow1 + pow2 > 0 => nr::F 1$F res := cons([f2*nrr, nrr, pow1, pow2]$Alg_Rec, res) res alglfint(f, k, l, x) == xf := x::F )if false if is?(k, 'nthRoot) then args := argument(k) -- a1 := args(1) n := retract(args(2))@Integer (m, c, r) := froot(args(1), n::NonNegativeInteger )$PolynomialRoots(IndexedExponents K, K, R, P, F) if c ~= 1 then rop := operator(k) k1 := kernel(rop, [r, m::F]) f1 := eval(f, k, c*k1) ir1 := lfintegrate(f1, x) return map((x1 : F) : F +-> eval(x1, k1, k::F/r), ir1) )endif symbolIfCan(kx := ksec(k, l, x)) case SE => addx(palgint(f, kx, k), xf) is?(kx, 'exp) => addx(algexpint(f, kx, k, x), xf) prim?(kx, x) => addx(algprimint(f, kx, k, x), xf) has?(operator kx, ALGOP) => is?(operator k, 'nthRoot) and is?(operator kx, 'nthRoot) and not(member?(kx, tower(k::F))) => al := alg_split_roots(f, k, kx) -- print(al::OutputForm) res := 0$IR for rec in al repeat ir1 := lfintegrate(rec.fun, x) if rec.npow1 + rec.npow2 > 0 then oroot := (k::F)^rec.npow1*(kx::F)^rec.npow2 ir1 := map((x1 : F) : F +-> eval(x1, retract(rec.nroot)@K, oroot), ir1) res := res + ir1 -- print(res::OutputForm) res rec := primitiveElement(kx::F, k::F) y := rootOf(rec.prim) map((x1 : F) : F+->eval(x1, retract(y)@K, rec.primelt), lfintegrate(eval(f, [kx, k], [(rec.pol1) y, (rec.pol2) y]), x)) unknownint(f, x) if R has Join(ConvertibleTo Pattern Integer, PatternMatchable Integer) and F has Join(LiouvillianFunctionCategory, RetractableTo SE) then import from PatternMatchIntegration(R, F) lfintegrate(f, x) == intPatternMatch(f, x, lfintegrate0, pmintegrate) else lfintegrate(f, x) == lfintegrate0(f, x) lfintegrate0(f, x) == zero? f => 0 xf := x::F empty?(l := varselect(kernels f, x)) => (xf * f)::IR symbolIfCan(k := kmax l) case SE => map((x1 : RF) : F+->multivariate(x1, k), integrate univariate(f, k)) is?(k, 'tan) => addx(tanint(f, x, k), xf) is?(k, 'exp) => addx(expint(f, x, k), xf) is?(k, 'lambertW) => addx(lambint(f, x, k), xf) prim?(k, x) => addx(primint(f, x, k), xf) has?(operator k, ALGOP) => alglfint(f, k, l, x) is?(k, OPDIFF) => diffint(f, x, k) unknownint(f, x) addx(i, x) == elem? i => i mkAnswer(ratpart i, logpart i, [[ne.integrand, x] for ne in notelem i]) tryChangeVar(f, t, x) == z := new()$Symbol g := subst(f / differentiate(t::F, x), [t], [z::F]) freeOf?(g, x) => -- can we do change of variables? map((x1 : F) : F+->eval(x1, kernel z, t::F), lfintegrate(g, z)) "failed" algexpint(f, t, y, x) == (u := tryChangeVar(f, t, x)) case IR => u::IR algint(f, t, y, (x1 : UP) : UP+->differentiate(x1, (x2 : F) : F+->differentiate(x2, x), monomial(differentiate(first argument t, x), 1))) algprimint(f, t, y, x) == (u := tryChangeVar(f, t, x)) case IR => u::IR algint(f, t, y, (x1 : UP) : UP+->differentiate(x1, (x2 : F) : F+->differentiate(x2, x), differentiate(t::F, x)::UP)) lfextendedint(f, x, g) == res1 := extendedint(f, x, [g]).particular res1 case "failed" => "failed" [res1.ratpart, (res1.coeffs)(1)] lflogint(f, x, lu) == lg := [differentiate(u, x)/u for u in lu] res1 := extendedint(f, x, lg).particular res1 case "failed" => "failed" s1 := res1::Param_Rec_F llg : LLG := [] for c in entries(s1.coeffs) for u in lu repeat if c ~= 0 then llg := cons([c, u], llg) [s1.ratpart, llg] lfinfieldint(f, x) == (u := lfextendedint(f, x, 0)) case "failed" => "failed" u.ratpart DREC ==> Record(answer : RF, logpart : RF, ir : IRRF) denint_dummy(f : RF) : DREC == [0, f, 0] DMap ==> (RF, K, SE) -> DREC denint_li : DMap := R is Integer and F is Expression(Integer) => (li_int$DenominatorIntegration(R, F))@DMap R is Complex(Integer) and F is Expression(Complex(Integer)) => (li_int$DenominatorIntegration(R, F))@DMap (rf : RF, k : K, x : SE) : DREC +-> denint_dummy(rf) denint_poly : DMap := R is Integer and F is Expression(Integer) => (poly_int$DenominatorIntegration(R, F))@DMap R is Complex(Integer) and F is Expression(Complex(Integer)) => (poly_int$DenominatorIntegration(R, F))@DMap (rf : RF, k : K, x : SE) : DREC +-> denint_dummy(rf) primint(f, x, k) == lk := varselect([a for a in tower f | a ~= k], x) denint : RF -> DREC := is?(k, 'log) => (rf : RF) : DREC +-> denint_li(rf, k, x) denint_dummy dk := differentiate(k::F, x) ext := (f1 : F) : U2 +-> r1 := extendedint(f1, x, lk, [dk]).particular r1 case "failed" => "failed" rf := r1::Param_Rec_F [rf.ratpart, (rf.coeffs)(1)] r1 := primintegrate(univariate(f, k), (x1 : UP) : UP+->differentiate(x1, (x2 : F) : F+->differentiate(x2, x), dk::UP), denint, ext) map((x1 : RF) : F+->multivariate(x1, k), r1.answer) + lfintegrate(r1.a0, x) cfind(f, l) == for u in l repeat f = u.logand => return u.coeff 0 risch_de_solver(x3 : Integer, x4 : F, eta : F, x : SE) : PSOL2 == risch_de_ext(x3, eta, x4, x, (x5 : F, x6 : List F) : U3+-> lflogint(x5, x, x6), (x7 : F, x8 : F) : U2+->lfextendedint(x7, x, x8) )$ElementaryRischDEX2(R, F) expint(f, x, k) == eta := first argument k r1 := expintegrate(univariate(f, k), (x1 : UP) : UP+->differentiate(x1, (x2 : F) : F+->differentiate(x2, x), monomial(differentiate(eta, x), 1)), (rf : RF) : DREC +-> denint_poly(rf, k, x), (x3 : Integer, x4 : F) : PSOL2 +-> risch_de_solver(x3, x4, eta, x)) map((x1 : RF) : F+->multivariate(x1, k), r1.answer) + lfintegrate(r1.a0, x) --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- SPAD files for the integration world should be compiled in the -- following order: -- -- intaux rderf intrf curve curvepkg divisor pfo -- intalg intaf efstruc rdeef intpm INTEF irexpand integrat spad Compiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/7929300622890310071-25px009.spad using old system compiler. INTEF abbreviates package ElementaryIntegration ------------------------------------------------------------------------ initializing NRLIB INTEF for ElementaryIntegration compiling into NRLIB INTEF ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: R already in scope ****** Domain: F already in scope ****** Domain: F already in scope importing IntegrationTools(R,F) importing ElementaryRischDE(R,F) importing RationalIntegration(F,SparseUnivariatePolynomial F) importing AlgebraicIntegration(R,F) importing AlgebraicManipulations(R,F) importing ElementaryRischDESystem(R,F) importing TranscendentalIntegration(F,SparseUnivariatePolynomial F) importing PureAlgebraicIntegration(R,F,F) importing IntegrationResultFunctions2(F,F) importing IntegrationResultFunctions2(Fraction SparseUnivariatePolynomial F,F) importing FunctionSpacePrimitiveElement(R,F) importing PolynomialCategoryQuotientFunctions(IndexedExponents Kernel F,Kernel F,R,SparseMultivariatePolynomial(R,Kernel F),F) importing ParametricIntegration(R,F) compiling local prim? : (Kernel F,Symbol) -> Boolean Time: 0.01 SEC. compiling local lambint : (F,Symbol,Kernel F) -> IntegrationResult F Time: 0.05 SEC. compiling local tanint : (F,Symbol,Kernel F) -> IntegrationResult F Time: 0.01 SEC. compiling local unknownint : (F,Symbol) -> IntegrationResult F Time: 0.03 SEC. processing macro definition LOG ==> Record(scalar: Fraction Integer,coeff: SparseUnivariatePolynomial F,logand: SparseUnivariatePolynomial F) compiling local diffint1 : (F,Symbol,Kernel F) -> Union(IntegrationResult F,failed) Time: 0.19 SEC. compiling local diffint : (F,Symbol,Kernel F) -> IntegrationResult F Time: 0.01 SEC. compiling local isx? : (F,Symbol) -> Boolean Time: 0 SEC. processing macro definition Alg_Rec ==> Record(fun: F,nroot: F,npow1: Integer,npow2: Integer) compiling local alg_split_root0 : (F,Kernel F,Integer) -> List F Time: 0.02 SEC. compiling local alg_split_roots : (F,Kernel F,Kernel F) -> List Record(fun: F,nroot: F,npow1: Integer,npow2: Integer) Time: 0.05 SEC. compiling local alglfint : (F,Kernel F,List Kernel F,Symbol) -> IntegrationResult F Time: 0.09 SEC. ****** Domain: R already in scope augmenting R: (ConvertibleTo (Pattern (Integer))) ****** Domain: R already in scope augmenting R: (PatternMatchable (Integer)) ****** Domain: F already in scope augmenting F: (LiouvillianFunctionCategory) importing PatternMatchIntegration(R,F) compiling exported lfintegrate : (F,Symbol) -> IntegrationResult F Time: 0.02 SEC. compiling exported lfintegrate : (F,Symbol) -> IntegrationResult F Time: 0 SEC. compiling exported lfintegrate : (F,Symbol) -> IntegrationResult F Time: 0 SEC. compiling exported lfintegrate : (F,Symbol) -> IntegrationResult F Time: 0 SEC. compiling local lfintegrate0 : (F,Symbol) -> IntegrationResult F Time: 0.02 SEC. compiling local addx : (IntegrationResult F,F) -> IntegrationResult F Time: 0 SEC. compiling local tryChangeVar : (F,Kernel F,Symbol) -> Union(IntegrationResult F,failed) Time: 0.01 SEC. compiling local algexpint : (F,Kernel F,Kernel F,Symbol) -> IntegrationResult F Time: 0.01 SEC. compiling local algprimint : (F,Kernel F,Kernel F,Symbol) -> IntegrationResult F Time: 0 SEC. compiling exported lfextendedint : (F,Symbol,F) -> Union(Record(ratpart: F,coeff: F),failed) ****** comp fails at level 4 with expression: ****** error in function lfextendedint (SEQ (LET |res1| ((|extendedint| |f| |x| (|construct| |g|)) | << particular >> |)) (|exit| 1 (IF (|case| |res1| "failed") "failed" (|construct| (|res1| |ratpart|) ((|res1| |coeffs|) 1))))) ****** level 4 ******$x:= particular $m:=$EmptyMode \$f:= ((((|coeff| #) (|ratpart| # . #1=#) (|coeff| #) (|ratpart| # . #1#) ...) ((|alg_split_roots| #) (|#| #) (< #) (<= #) ...) ((|alg_split_root0| #) (|#| #) (< #) (<= #) ...) ((|diffint1| #) (|case| #) (|autoCoerce| #) (|coerce| #) ...))) >> Apparent user error: no modemap for extendedint with 3 arguments )read load.input fricas)lib LINCOMB LinearCombinationUtilities is already explicitly exposed in frame initial LinearCombinationUtilities will be automatically loaded when needed from /var/aw/var/LatexWiki/LINCOMB.NRLIB/LINCOMB fricas)lib EFACTOR ExpressionFactorPolynomial is already explicitly exposed in frame initial ExpressionFactorPolynomial will be automatically loaded when needed from /var/aw/var/LatexWiki/EFACTOR.NRLIB/EFACTOR fricas)lib PFUTIL PartialFractionUtilities is already explicitly exposed in frame initial PartialFractionUtilities will be automatically loaded when needed from /var/aw/var/LatexWiki/PFUTIL.NRLIB/PFUTIL fricas)lib INTPAR1 ParametricTranscendentalIntegration is already explicitly exposed in frame initial ParametricTranscendentalIntegration will be automatically loaded when needed from /var/aw/var/LatexWiki/INTPAR1.NRLIB/INTPAR1 fricas)lib RDEPAR )library cannot find the file RDEPAR. fricas)lib INTALG2 AlgebraicIntegrate2 is already explicitly exposed in frame initial AlgebraicIntegrate2 will be automatically loaded when needed from /var/aw/var/LatexWiki/INTALG2.NRLIB/INTALG2 fricas)lib INTAF AlgebraicIntegration is already explicitly exposed in frame initial AlgebraicIntegration will be automatically loaded when needed from /var/aw/var/LatexWiki/INTAF.NRLIB/INTAF fricas)lib INTPAR2 )library cannot find the file INTPAR2. fricas)lib INTEF )library cannot find the file INTEF. fricasintegrate((polylog(3, x) - x*polylog(2, x))/(1 - x)^2, x) >> Error detected within library code: Sorry - cannot handle that integrand yet

 Subject:   Be Bold !! ( 14 subscribers )