On Wed, Dec 20, 2023 at 07:40:00AM +0800, Qian Yun wrote:
> integrate(atan(x*(x^2+1)^(1/2)),x)
> integrate(tanh(x)/(exp(x)+exp(2*x))^(1/2),x)
> integrate((2+2*tan(x)+tan(x)^2)^(1/2),x)
> ...
>
> These integrals are different between trunk and this patch.
> (same as 1.3.9.)
Newer patch, fixes bunch of problems with previous patches.
Previous patch did not handle correclty the following:
integrate((14*x^2-18*x-4)/(x^4+(-7)*x^2+6*x+1), x)
-- poly x^4 - 7*x^2 + 6*x + 1, 4 real roots
integrate((-8*x^2-24*x-4)/(x^4+4*x^2+8*x+1), x)
-- poly x^4 + 4*x^2 + 8*x + 1, 2 complex roots and 2 real roots
--
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/ZYMxdnX7ym5Lb6Kj%40fricas.org.
diff --git a/src/algebra/rsimp.spad b/src/algebra/rsimp.spad
index 04322e8..22f00f5 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) -> Union(r_rec, "failed")
++ quartic2(p) should be local but conditional
+ my_sqrt : F -> F
+ ++ my_sqrt(a) should be local but conditional
my_root3 : F -> F
++ my_root3(a) should be local but conditional
@@ -24,7 +26,7 @@ RootUtilities(R, F) : Exports == Implementation where
P ==> SparseMultivariatePolynomial(R, K)
dummy := create()$SingletonAsOrderedSet
-
+
root_pair(a : F, b : F) : C_rec == [a, b]
root4(a : F) : F ==
@@ -54,58 +56,78 @@ RootUtilities(R, F) : Exports == Implementation where
degree(fr1) ~= 1 => error "impossible"
-coefficient(fr1, 0)/coefficient(fr1, 1)
- -- Q = (-4)*u*v^2+(4*u^3+2*a*u+b)
- quartic3(r1 : F, a0 : F, a : F, b : F) : List(C_rec) ==
+ do_rp(r : F, q : F, res1 : r__rec) : Union(r_rec, "failed") ==
+ (su := sign(q)) case "failed" => "failed"
+ (0 < su@Integer) =>
+ qr := my_sqrt(q)
+ [res1.reals, cons(root_pair(r, qr), res1.complexes)]
+ qr := my_sqrt(-q)
+ [cons(r + qr, cons(r - qr, res1.reals)), res1.complexes]
+
+ -- Square of real part u and square of imaginary part v
+ -- satisfy (-4)*u*v^2+(4*u^3+2*a*u+b) = 0
+ -- so given real part we compute imaginary parts.
+ -- In fact, above u is square of average of two roots and
+ -- -v^2 is square od half of difference between roots, so
+ -- this also works for real roots.
+ quartic3(r1 : F, a0 : F, a : F, b : F) : Union(r_rec, "failed") ==
rr := my_sqrt(r1)
qq := r1 + (1/2::F)*a
q1 := qq + b*rr/(4*r1)
+ res1 := do_rp(rr + a0, q1, [[], []])
+ res1 case "failed" => "failed"
q2 := qq - b*rr/(4*r1)
- [root_pair(rr - a0/(4::F), my_sqrt(q1)),
- root_pair(-rr - a0/(4::F), my_sqrt(q2))]
+ do_rp(-rr + a0, q2, res1)
- do_f2(del : F, a2 : F, b2 : F, a0 : F, a : F, b : F) : List(C_rec) ==
+ do_f2(del : F, a2 : F, b2 : F, a0 : F, a : F, b : F
+ ) : Union(r_rec, "failed") ==
d2 := my_sqrt(del)
- (su := sign(d2)) case "failed" => []
+ (su := sign(d2)) case "failed" => "failed"
if su@Integer < 0 then
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)
+
+ -- Squares of averages of pairs of roots of quartic satisfy
+ -- equation 64u^3 + 32au^2 + (-16c+4a^2)*u - b^2 = 0.
+ -- If the qubic factors, and we can determine needed
+ -- signs, then this works, otherwise fails.
+ quartic2(p : UP, a0 : F) : Union(r_rec, "failed") ==
a := coefficient(p, 2)
b := coefficient(p, 1)
c := coefficient(p, 0)
xx := monomial(1, 1)$UP
- r := 64*xx^3 + 32*a*xx^2
+ r := 64*xx^3 + 32*a*xx^2
r := r + (-16*c+4*a^2)*xx
r := r - (b^2)::UP
fr := factor(r)
fl := factorList(fr)
- #fl = 1 => []
- #fl = 3 =>
+ #fl = 1 => "failed"
+ flf : List(UP) := []
+ for fac in fl repeat
+ f1 := fac.factor
+ fac.exponent = 2 =>
+ flf := cons(f1, cons(f1, flf))
+ flf := cons(f1, flf)
+ #flf = 3 =>
r1 : F
cnt : Integer := 0
- for fac in fl repeat
- f1 := fac.factor
+ for f1 in flf repeat
cc := coefficient(f1, 0)/coefficient(f1, 1)
- (su := sign(cc)) case "failed" => return []
+ (su := sign(cc)) case "failed" => return "failed"
if su@Integer < 0 then
- cnt > 0 => return []
- cnt := 1
+ cnt := cnt + 1
+ cnt > 1 => iterate
r1 := cc
cnt = 0 => error "impossible"
quartic3(-r1, a0, a, b)
- #fl = 2 =>
+ #flf = 2 =>
r1 : F
- f1 := first(fl).factor
- f2 := second(fl).factor
- if (degree(f2) = 1) then
+ f1 := first(flf)
+ f2 := second(flf)
+ if degree(f2) = 1 then
(f1, f2) := (f2, f1)
r1 := -coefficient(f1, 0)/coefficient(f1, 1)
- (su := sign(r1)) case "failed" => return []
+ (su := sign(r1)) case "failed" => return "failed"
s1 := su@Integer
a2 := coefficient(f2, 2)
b2 := coefficient(f2, 1)
@@ -122,17 +144,22 @@ RootUtilities(R, F) : Exports == Implementation where
return do_f2(del, a2, b2, a0, a, b)
return quartic3(r1, a0, a, b)
su := sign(b2)
- su case "failed" or not(0 < su@Integer) => return []
+ su case "failed" or not(0 < su@Integer) => return "failed"
su := sign(a2*coefficient(f2, 0))
if su case "failed" or not(0 < su@Integer) then
- return []
+ return "failed"
quartic3(r1, a0, a, b)
error "impossible"
else
- quartic2(p : UP) : List(C_rec) == []
-
+ quartic2(p : UP, a0 : F) : Union(r_rec, "failed") == "failed"
+
+ 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)
rec.exponent = 3 => a^(1/3)
@@ -156,22 +183,61 @@ 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
- si = 1 =>
+ a = 0 => [[a0, a0, a0, a0], []]
+ (su := sign(a)) case "failed" => "failed"
+ (si := su@Integer) = 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)
- not(empty?(res1)) => [[], res1]
- "failed"
+ 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], []]
+ quartic2(p, a0)
+ quartic2(p, a0)
complex_roots(p : UP) : Union(r_rec, "failed") ==
d := degree(p)
