On 11/19/23 20:28, Waldek Hebisch wrote:
Until recently I found that Rioboo's algorithm described in
Manuel Bronstein's book, at the end of section 2.8, covers
this situation.
Hmm, I am not sure what you want to say here. We use Rioboo's
algorithm, but in some cases it does not give us enough
improvement.
My mistake. We are using Rioboo's algorithm for rational integration.
I was talking about section 2.8 "Riboo's Algorithm for Real Rational
Functions", page 65 to 70 for second edition of "Symbolic
Integration I".
AFAICS most significant reason for troubles is "atan doubling"
or even tripling: trying to keep things real we get square
(or higer power) of argument to logaritm. And we use changes
of variables which require kind of "to real" transform when
coming back. This in principle could be solved by factoring
arguments of logs.
This algorithm prevents generation of (unnecessary) complex
expression from the beginning.
For even higher degree polynomials, I think it's wise to
return rootSum instead of give results containing algebraic
numbers, which causes trouble for sagemath.
I posted a patch to do this some time ago. But it causes
regression in definite integrator. AFAICS, currently we
pass complex algebraic numbers to limit code which is wrong.
But in some cases we get sensible results. Once we have
rootSum limit code sees that there is trouble and gives up.
General case of limit of rootSum looks tricky. There is
good chance that with moderate effort we could improve
limit to handle cases that currently work using algebraic
numbers (and fail on other cases).
Now my patch has no test failures, because I limit it to
work on polynomials degree >= 4. Essentially replacing "lg2cfunc".
(For degree < 4, my patch is not wrong -- just not as optimized
as "quadratic".)
I'll post patch after I reviewed the test failures and do more
testing.
Could you test interaction of your patch with attached patch.
Attached patch remove (I think) last call to 'real' from the
integrator. We should do this as call to 'real' causes several
wrong results. But removing call to 'real' breaks some definite
integrals, so I am waiting with this patch for some solution for
definite integrals.
Brief test shows no bad interaction -- my patch works at a place
before yours have effect.
Removing call to 'real' means that following transformations of
integrals will see complex expression. For example, with 'real'
removed we have:
(2) -> integrate(1/(z^4+1),z)
(2)
+---+ +-+ +---+ +-+
(\|- 1 + 1)\|2 log((\|- 1 + 1)\|2 + 2 z)
+
+---+ +-+ +---+ +-+
(\|- 1 - 1)\|2 log((\|- 1 - 1)\|2 + 2 z)
+
+---+ +-+ +---+ +-+
(- \|- 1 + 1)\|2 log((- \|- 1 + 1)\|2 + 2 z)
+
+---+ +-+ +---+ +-+
(- \|- 1 - 1)\|2 log((- \|- 1 - 1)\|2 + 2 z)
/
8
Isn't a real result even better than this?
My patch is in attachment, the variable names are bad, I know.
- Qian
--
You received this message because you are subscribed to the Google Groups "FriCAS -
computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/fricas-devel/906651f9-758e-482d-aa89-2806c21f5a75%40gmail.com.
diff --git a/src/algebra/irexpand.spad b/src/algebra/irexpand.spad
index ead11d7f..97ffc16a 100644
--- a/src/algebra/irexpand.spad
+++ b/src/algebra/irexpand.spad
@@ -31,6 +31,7 @@ IntegrationResultToFunction(R, F) : Exports == Implementation where
expand : (IR, Symbol) -> List F
++ expand(i, x) returns the list of possible real functions
++ of x corresponding to i.
+ lg2func2 : (LOG, Symbol) -> List F
complexExpand : IR -> F
++ complexExpand(i) returns the expanded complex function
++ corresponding to i.
@@ -38,6 +39,7 @@ IntegrationResultToFunction(R, F) : Exports == Implementation where
Implementation ==> add
import from AlgebraicManipulations(R, F)
import from ElementaryFunctionSign(R, F)
+ import from TrigonometricManipulations(R, F)
IR2F : IR -> F
insqrt : F -> Record(sqrt : REC, sgn : Z)
@@ -50,6 +52,7 @@ IntegrationResultToFunction(R, F) : Exports == Implementation where
ilog0 : (F, F, UP, UP, F) -> F
nlogs : LOG -> List LOG
lg2func : (LOG, Symbol) -> List F
+ --lg2func2 : (LOG, Symbol) -> List F
quadratic : (UP, UP, Symbol) -> List F
mkRealFunc : (List LOG, Symbol) -> List F
lg2cfunc : LOG -> F
@@ -160,8 +163,63 @@ IntegrationResultToFunction(R, F) : Exports == Implementation where
lg2func([lg.scalar,
(p exquo (monomial(1, 1)$UP - alpha::UP))::UP,
lg.logand], x))
+ l := lg2func2(lg, x)
+ not empty? l => l
[lg2cfunc lg]
+ if R is Integer and F is Expression Integer then
+ import from TransSolvePackage(R)
+
+ lg2func2(lg, x) ==
+ RR := lg.coeff
+ SS := lg.logand
+ u := new('u)$Symbol
+ v := new('v)$Symbol
+ i := new('i)$Symbol
+ Ruvi : F := RR(u::F + v::F * i::F)
+ not real? Ruvi => []
+ Suvi : F := SS(u::F + v::F * i::F)
+ not real? Suvi => []
+ Ruv := eval(Ruvi, i::F = sqrt(-1::F))
+ PP := real Ruv
+ QQ := imag Ruv
+ Suv := eval(Suvi, i::F = sqrt(-1::F))
+ AA := real Suv
+ BB := imag Suv
+ all_roots := solve([PP=0,QQ=0],[u,v])
+ real_roots := [s for s in all_roots | real? rhs first s and real? rhs second s]
+ degree(RR) ~= #real_roots => []
+ SOL ==> List Equation F
+ l2 : List List SOL := []
+ real_real_roots : List SOL := []
+ while not empty? real_roots repeat
+ s := first real_roots
+ zero? rhs second s =>
+ real_real_roots := cons(s, real_real_roots)
+ real_roots := rest real_roots
+ member?(s, concat l2) =>
+ real_roots := rest real_roots
+ lu := [t for t in rest real_roots | rhs first s = rhs first t]
+ ll := [t for t in lu | rhs second s + rhs second t = 0]
+ not size?(ll, 1) => error "not expected"
+ l2 := cons([s, first ll], l2)
+ real_roots := rest real_roots
+ g : F := 0
+ h : F := 0
+ atans : F := 0
+ for root_uv in real_real_roots repeat
+ s : F := elt(SS, rhs first root_uv)
+ g := g + rhs first root_uv * log s
+ for pair in l2 repeat
+ root_uv := first pair
+ a := eval(AA, root_uv)
+ b := eval(BB, root_uv)
+ h := h + rhs first root_uv * log(a^2+b^2)
+ atans := atans + rhs second root_uv * ilog(a,b,x)
+ [g + h + atans]
+ else
+ lg2func2(lg, x) == []
+
lg2cfunc lg ==
+/[cmplex(alpha, lg.logand) for alpha in zerosOf(lg.coeff)]
