On Mon, Dec 18, 2023 at 06:12:01PM +0800, Qian Yun wrote:
> I'm opening a new thread for this newly added rsimp.spad.
>
> I have tested with first 10k of Nasser's integrals, the recent
> changes to integrator is positive:
>
> It affects the result of a few percent integrals, almost all are
> positive improvement.
>
> Only a few regressions, in 3 classes:
>
> 1. not invertiable:
> integrate(x^6/(2*x^5+3)^3,x)
This is caused by troubles with other routines. Namely 'complex_roots'
gets polynomial witgh dependent roots which causes failure in 'factor'.
To reproduce do:
a := (-1)^(1/5)/(108)^(1/5)
p1 := x^4 + a*x^3/250 + a^2*x^2/62500 + a^3*x/15625000 + a^4/3906250000
factor(p1::UP('x, EXPR(INT))::SUP(EXPR(INT)))
If one replaces a by
a := subst(b^(1/5), b = -1/108)
then there is error later, again coused by dependent roots. This
time 'factor' returns empty list of factors for polynomial of
degree 2...
I have a workaround, but ATM I am not sure how to fix this properly.
> 2. third "impossible" in rsimp.spad:
> integrate(1/(x^4+4*x^2+1),x)
>
> 3. fourth "impossible" in rsimp.spad:
>
> This error happens when resolvent contains zero factor, it seems.
>
> integrate(log(x^2+(-x^2+1)^(1/2)),x)
> integrate(log(1+x*(x^2+1)^(1/2)),x)
> integrate(atan(x*(-x^2+1)^(1/2)),x)
> integrate(1/(x^4+3*x^2-1),x)
> integrate(1/(x^4-3*x^2-1),x)
> integrate(1/(x^2-1)^(1/2)/(x^(1/2)+(x^2-1)^(1/2))^2,x)
> integrate((x^(1/2)-(x^2-1)^(1/2))^2/(-x^2+x+1)^2/(x^2-1)^(1/2),x)
All of the above are fixed by adding special case for biquadratic.
I missed possibility that irreducible quartic may have 4 purely imaginary
roots, so real part may be 0. Case of general quartic needs rechecking,
while in our case real part may be 0 implies biquadratic, so will
no longer happen but it is possible that I also missed something
else.
Patch for biquadratic in the attachment.
--
Waldek Hebisch
--
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/ZYCDWLJ_9DyVq23V%40fricas.org.
diff --git a/src/algebra/rsimp.spad b/src/algebra/rsimp.spad
index 04322e8..20f0bdb 100644
--- a/src/algebra/rsimp.spad
+++ b/src/algebra/rsimp.spad
@@ -10,8 +10,10 @@ RootUtilities(R, F) : Exports == Implementation where
UP ==> SUP(F)
Exports ==> with
complex_roots : UP -> Union(r_rec, "failed")
- quartic2 : UP -> List(C_rec)
+ quartic2 : (UP, F) -> List(C_rec)
++ quartic2(p) should be local but conditional
+ my_sqrt : F -> F
+ ++ my_sqrt(a) should be local but conditiona
my_root3 : F -> F
++ my_root3(a) should be local but conditional
@@ -60,8 +62,8 @@ RootUtilities(R, F) : Exports == Implementation where
qq := r1 + (1/2::F)*a
q1 := qq + b*rr/(4*r1)
q2 := qq - b*rr/(4*r1)
- [root_pair(rr - a0/(4::F), my_sqrt(q1)),
- root_pair(-rr - a0/(4::F), my_sqrt(q2))]
+ [root_pair(rr + a0, my_sqrt(q1)),
+ root_pair(-rr + a0, my_sqrt(q2))]
do_f2(del : F, a2 : F, b2 : F, a0 : F, a : F, b : F) : List(C_rec) ==
d2 := my_sqrt(del)
@@ -70,11 +72,7 @@ RootUtilities(R, F) : Exports == Implementation where
d2 := -d2
quartic3((-b2 + d2)/(2*a2), a0, a, b)
- quartic2(p : UP) : List(C_rec) ==
- lc := leadingCoefficient(p)
- p := (1/lc)*p
- a0 := coefficient(p, 3)
- p := eval(p, dummy, monomial(1, 1)$UP - (a0/4::F)::UP)
+ quartic2(p : UP, a0 : F) : List(C_rec) ==
a := coefficient(p, 2)
b := coefficient(p, 1)
c := coefficient(p, 0)
@@ -131,7 +129,12 @@ RootUtilities(R, F) : Exports == Implementation where
else
- quartic2(p : UP) : List(C_rec) == []
+ quartic2(p : UP, a0 : F) : List(C_rec) == []
+
+ my_sqrt(a : F) ==
+ rec := froot(a, 2)$PolynomialRoots(IndexedExponents(K), K, R, P, F)
+ rec.exponent = 2 => sqrt(a)
+ rec.coef*rec.radicand
my_root3(a : F) : F ==
rec := froot(a, 3)$PolynomialRoots(IndexedExponents(K), K, R, P, F)
@@ -156,20 +159,62 @@ RootUtilities(R, F) : Exports == Implementation where
[[-r1 + a0], [root_pair(r1/(2::F) + a0, sqrt(3::F)*r1/(2::F))]]
"failed"
+ quartic0(p : UP, a0 : F) : Union(r_rec, "failed") ==
+ pu := p exquo monomial(1, 1)$UP
+ pu case "failed" => error "impossible"
+ p := pu@UP
+ ru := qubic(p)
+ ru case "failed" => "failed"
+ r1 := ru@r_rec
+ rl : List(F) := [a0]
+ for r in r1.reals repeat
+ rl := cons(r + a0, rl)
+ cl : List(C_rec) :=
+ empty?(r1.complexes) => []
+ cp := first(r1.complexes)
+ [root_pair(cp.real + a0, cp.imag)]
+ [rl, cl]
+
quartic(p : UP) : Union(r_rec, "failed") ==
+ coefficient(p, 0) = 0 => quartic0(p, 0)
+ lc := leadingCoefficient(p)
+ p := (1/lc)*p
+ a0 := -coefficient(p, 3)/(4::F)
+ p := eval(p, dummy, monomial(1, 1)$UP + a0::UP)
ground?(rp := reductum(p)) =>
a := ground(rp)
- s := sign(a)
- s case "failed" => "failed"
- si := s@Integer
+ a = 0 => [[a0, a0, a0, a0], []]
+ (su := sign(a)) case "failed" => "failed"
+ si := su@Integer
si = 1 =>
r1 := root4(a/(4::F*leadingCoefficient(p)))
- [[], [root_pair(r1, r1), root_pair(-r1, r1)]]
+ [[], [root_pair(r1 + a0, r1), root_pair(-r1 + a0, r1)]]
si = -1 =>
r1 := rootSimp(zeroOf(p))
- [[r1, -r1], [root_pair(0, r1)]]
+ [[r1 + a0, -r1 + a0], [root_pair(a0, r1)]]
error "impossible"
- res1 := quartic2(p)
+ coefficient(p, 0) = 0 => quartic0(p, a0)
+ coefficient(p, 1) = 0 =>
+ b := coefficient(p, 2)
+ c := coefficient(p, 0)
+ del := b*b - 4*c
+ (su := sign(del)) case "failed" => "failed"
+ not((si := su@Integer) < 0) =>
+ r1 := my_sqrt(del)
+ x1 := (-b + r1)/(2::F)
+ x2 := (-b - r1)/(2::F)
+ (su := sign(x1)) case "failed" => "failed"
+ (su@Integer < 0) =>
+ [[], [root_pair(a0, my_sqrt(-x1)),
+ root_pair(a0, my_sqrt(-x2))]]
+ (su := sign(x2)) case "failed" => "failed"
+ r2 := my_sqrt(x1)
+ (su@Integer < 0) =>
+ [[r2 + a0, -r2 + a0], [root_pair(a0, my_sqrt(-x2))]]
+ r3 := my_sqrt(x2)
+ [[r2 + a0, -r2 + a0, r3 + a0, -r3 + a0], []]
+ "failed"
+ res1 := quartic2(p, a0)
not(empty?(res1)) => [[], res1]
"failed"
