SF.net SVN: fricas:[2630] trunk

[wheb...@users.sourceforge.net]

Revision: 2630
http://sourceforge.net/p/fricas/code/2630
Author: whebisch
Date: 2020-03-04 17:00:29 +0000 (Wed, 04 Mar 2020)
Log Message:
-----------
Handle recursion for algebraic integrals

Modified Paths:
--------------
trunk/ChangeLog
trunk/src/algebra/intaf.spad
trunk/src/algebra/intalg.spad
trunk/src/algebra/intef.spad
trunk/src/input/integ.input

Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2020-03-03 19:38:13 UTC (rev 2629)
+++ trunk/ChangeLog 2020-03-04 17:00:29 UTC (rev 2630)
@@ -1,3 +1,9 @@
+2020-03-04 Waldek Hebisch <heb...@math.uni.wroc.pl>
+
+ * src/algebra/intaf.spad, src/algebra/intalg.spad,
+ src/algebra/intef.spad, src/input/integ.input:
+ Handle recursion for algebraic integrals
+
2020-03-03 Waldek Hebisch <heb...@math.uni.wroc.pl>

* src/interp/i-eval.boot: Initialize variable

Modified: trunk/src/algebra/intaf.spad
===================================================================
--- trunk/src/algebra/intaf.spad 2020-03-03 19:38:13 UTC (rev 2629)
+++ trunk/src/algebra/intaf.spad 2020-03-04 17:00:29 UTC (rev 2630)
@@ -725,7 +725,7 @@
FAIL==> error "failed - cannot handle that integrand"

Exports ==> with
- algint : (F, K, K, UP -> UP) -> IR
+ algint : (F, K, K, UP -> UP, F -> IR) -> 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]}.
@@ -746,8 +746,8 @@
import from PolynomialCategoryQuotientFunctions(IndexedExponents K,
K, R, P, F)

- rootintegrate : (F, K, K, UP -> UP) -> IR
- algintegrate : (F, K, K, UP -> UP) -> IR
+ rootintegrate : (F, K, K, UP -> UP, F -> IR) -> IR
+ algintegrate : (F, K, K, UP -> UP, F -> IR) -> IR
UPUP2F : (UPUP, RF, K, K) -> F
F2UPUP : (F, K, K, UP) -> UPUP
UP2UPUP : (UP, K) -> UPUP
@@ -762,7 +762,7 @@

F2UPUP(f, kx, k, p) == UP2UPUP(univariate(f, k, p), kx)

- rootintegrate(f, t, k, derivation) ==
+ rootintegrate(f, t, k, derivation, rec_int) ==
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)
@@ -769,15 +769,17 @@
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(reduce f1, derivation, rec_int
+ )$ALG)$IR2

- algintegrate(f, t, k, derivation) ==
+ algintegrate(f, t, k, derivation, rec_int) ==
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
+ algintegrate(reduce f1, derivation, rec_int
+ )$ALG)$IR2

Curv_Rec ==> Record(funs : List(UPUP),
irec : Record(coef : Fraction(UP), poly : UPUP),
@@ -837,9 +839,9 @@
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)
+ algint(f, t, y, derivation, rec_int) ==
+ is?(y, 'nthRoot) => rootintegrate(f, t, y, derivation, rec_int)
+ is?(y, 'rootOf) => algintegrate(f, t, y, derivation, rec_int)
FAIL

--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.

Modified: trunk/src/algebra/intalg.spad
===================================================================
--- trunk/src/algebra/intalg.spad 2020-03-03 19:38:13 UTC (rev 2629)
+++ trunk/src/algebra/intalg.spad 2020-03-04 17:00:29 UTC (rev 2630)
@@ -153,6 +153,7 @@
GP ==> LaurentPolynomial(F, UP)
K ==> Kernel F
IR ==> IntegrationResult R
+ IRF ==> IntegrationResult F
UPR ==> SparseUnivariatePolynomial R
FD ==> FiniteDivisor(F, UP, UPUP, R)
LOG ==> Record(scalar : Q, coeff : UPR, logand : UPR)
@@ -166,8 +167,8 @@
NOTI ==> error "Not integrable (provided residues have no relations)"

Exports ==> with
- algintegrate : (R, UP -> UP) -> IR
- ++ algintegrate(f, d) integrates f with respect to the derivation d.
+ algintegrate : (R, UP -> UP, F -> IRF) -> IR
+ ++ algintegrate(f, d, rec) integrates f with respect to the derivation d.
palgintegrate : (R, UP -> UP) -> IR
++ palgintegrate(f, d) integrates f with respect to the derivation d.
++ Argument f must be a pure algebraic function.
@@ -196,8 +197,8 @@
nonQ : (DIV, UP) -> Union(List LOG, "failed")
rlift : (F, K, K) -> R
varRoot? : (UP, F -> F) -> Boolean
- algintexp : (R, UP -> UP) -> IR
- algintprim : (R, UP -> UP) -> IR
+ algintexp : (R, UP -> UP, F -> IRF) -> IR
+ algintprim : (R, UP -> UP, F -> IRF) -> IR

dummy : R := 0

@@ -207,7 +208,19 @@
F2UPR f == F2R(f)::UPR
F2R f == f::UP::QF::R

- algintexp(f, derivation) ==
+ IRF_to_IR(irf : IRF) : IR ==
+ map(F2R, irf)$IntegrationResultFunctions2(F, R)
+
+ RF ==> Fraction UP
+
+ R_to_F(f : R) : Union(F, "failed") ==
+ (u1 := retractIfCan(f)@Union(RF, "failed")) case "failed" => "failed"
+ f1 := u1::RF
+ (u2 := retractIfCan(f1)@Union(UP, "failed")) case "failed" => "failed"
+ f2 := u2::UP
+ retractIfCan(f2)@Union(F, "failed")
+
+ algintexp(f, derivation, rec_int) ==
d := (c := integralCoordinates f).den
v := c.num
vp : Vector(GP) := new(n := #v, 0)
@@ -218,16 +231,23 @@
qsetelt!(vp, i, r.polyPart)
ff := represents(vf, w := integralBasis())
h := HermiteIntegrate(ff, derivation)
- p := represents(map((x1 : GP) : QF+->convert(x1)@QF, vp)$VectorFunctions2(GP, QF), w)
- zero?(h.logpart) and zero? p => h.answer::IR
+ p := represents(map((x1 : GP) : QF+->convert(x1)@QF,
+ vp)$VectorFunctions2(GP, QF), w)
+ not(zero?(p)) => FAIL3
+ zero?(h.logpart) => h.answer::IR
+ (p1u := R_to_F(h.logpart)) case F =>
+ res1 := rec_int(p1u::F)
+ mkAnswer(h.answer, empty(), empty()) + IRF_to_IR(res1)
(u := alglogint(h.logpart, derivation)) case "failed" =>
mkAnswer(h.answer, empty(), [[p + h.logpart, dummy]])
- zero? p => mkAnswer(h.answer, u::List(LOG), empty())
- FAIL3
+ mkAnswer(h.answer, u::List(LOG), empty())

- algintprim(f, derivation) ==
+ algintprim(f, derivation, rec_int) ==
h := HermiteIntegrate(f, derivation)
zero?(h.logpart) => h.answer::IR
+ (p1u := R_to_F(h.logpart)) case F =>
+ res1 := rec_int(p1u::F)
+ mkAnswer(h.answer, empty(), empty()) + IRF_to_IR(res1)
(u := alglogint(h.logpart, derivation)) case "failed" =>
mkAnswer(h.answer, empty(), [[h.logpart, dummy]])
mkAnswer(h.answer, u::List(LOG), empty())
@@ -572,13 +592,13 @@
mkAnswer(h.answer, u::List(LOG), [[difFirstKind, dummy]])

-- for mixed functions. f dx not assumed locally integral at infinity
- algintegrate(f, derivation) ==
+ algintegrate(f, derivation, rec_int) ==
x' := derivation(x := monomial(1, 1)$UP)
zero? degree(x') =>
- algintprim(f, derivation)
+ algintprim(f, derivation, rec_int)
((xx := x' exquo x) case UP) and
(retractIfCan(xx::UP)@Union(F, "failed") case F) =>
- algintexp(f, derivation)
+ algintexp(f, derivation, rec_int)
error "should not happen"

alglogint(f, derivation) ==
@@ -585,6 +605,7 @@
r := primitivePart(doubleResultant(f, derivation))
varRoot?(r, x1+->retract(derivation(x1::UP))@F) => "failed"
FAIL0
+
)boot $tryRecompileArguments := true

--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.

Modified: trunk/src/algebra/intef.spad
===================================================================
--- trunk/src/algebra/intef.spad 2020-03-03 19:38:13 UTC (rev 2629)
+++ trunk/src/algebra/intef.spad 2020-03-04 17:00:29 UTC (rev 2630)
@@ -211,7 +211,8 @@
algint(f, t, y,
(x1 : UP) : UP+->differentiate(x1,
(x2 : F) : F+->differentiate(x2, x),
- monomial(differentiate(first argument t, x), 1)))
+ monomial(differentiate(first argument t, x), 1)),
+ f1 +-> lfintegrate(f1, x))

algprimint(f, t, y, x) ==
(u := tryChangeVar(f, t, x)) case IR => u::IR
@@ -218,7 +219,8 @@
algint(f, t, y,
(x1 : UP) : UP+->differentiate(x1,
(x2 : F) : F+->differentiate(x2, x),
- differentiate(t::F, x)::UP))
+ differentiate(t::F, x)::UP),
+ f1 +-> lfintegrate(f1, x))

lfextendedint(f, x, g) ==
res1 := extendedint(f, x, [g]).particular

Modified: trunk/src/input/integ.input
===================================================================
--- trunk/src/input/integ.input 2020-03-03 19:38:13 UTC (rev 2629)
+++ trunk/src/input/integ.input 2020-03-04 17:00:29 UTC (rev 2630)
@@ -85,6 +85,21 @@
fun4 := D(log(log(log(x)))*exp(log(x) + x)*(x + 1/(log(x)^2 + x)), x)
testIntegrate("fun4", "x", "rde")

+-- Just algebraic Hermite
+testIntegrate(
+ "(4*(log(x)-1)*(log(x)^2+x+1)^(1/2)+(2*log(x)^2-2*log(x)+x+2))/"
+ "((10*log(x)^2+2*x+2)*(log(x)^2+x+1)^(1/2)+8*log(x)^3+(8*x+8)*log(x))",
+ "x", "alg")
+-- Needs recursive call to integrator
+testIntegrate(
+ "(4*(log(x)-1)*(log(x)^2+x+1)^(1/2)+(2*log(x)^2-2*log(x)+x+2))/"
+ "((10*log(x)^2+2*x+2)*(log(x)^2+x+1)^(1/2)+8*log(x)^3+(8*x+8)*log(x))+"
+ "1/x", "x", "alg")
+testIntegrate(
+ "(4*(log(x)-1)*(log(x)^2+x+1)^(1/2)+(2*log(x)^2-2*log(x)+x+2))/"
+ "((10*log(x)^2+2*x+2)*(log(x)^2+x+1)^(1/2)+8*log(x)^3+(8*x+8)*log(x))+"
+ "1/x^2", "x", "alg")
+
-- Needs logarithmic derivative in algebraic extension
testIntegrate("sqrt(z)^I^2", "x", "rde")
-- D(exp(ellipticF(x, m)), x)

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.