Third version in the attachemtnt, this time handling rather
general quartics where we can find real parts of roots. It
uses precomputed resolvent. After normalizing so that sum
of roots is zero resolvent is an even polynomial, so we
get equation of degree 3 for square of real part. If
this has one posivite real root which can be obtained via
factoring, then we are in business, otherwise this fails.
The routine uses factorization, so works only when F has
PolynomialFactorizationExplicit.
This case is somewhat rare, example is:
integrate(tan(x)/(sqrt(tan(x)) - 1)^2, x)
(1)
+-+ +------+ +-+ +-+ +------+
((- \|2 + 2)\|tan(x) + \|2 - 2)log(\|2 \|tan(x) + tan(x) + 1)
+
+-+ +------+ +-+ +------+
(4 \|2 \|tan(x) - 4 \|2 )log(\|tan(x) - 1)
+
+-+ +------+ +-+ +-+ +------+
((- \|2 - 2)\|tan(x) + \|2 + 2)log(- \|2 \|tan(x) + tan(x) + 1)
+
+-+ +------+ +-+ +-+ +------+
((2 \|2 - 4)\|tan(x) - 2 \|2 + 4)atan(\|2 \|tan(x) + 1)
+
+-+ +------+ +-+ +-+ +------+ +-+
((- 2 \|2 - 4)\|tan(x) + 2 \|2 + 4)atan(\|2 \|tan(x) - 1) - 4 \|2
/
+-+ +------+ +-+
4 \|2 \|tan(x) - 4 \|2
Type: Union(Expression(Integer),...)
--
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/ZW0i10xUfWBXTmvg%40fricas.org.
diff --git a/src/algebra/irexpand.spad b/src/algebra/irexpand.spad
index ead11d7f..d7fb0dbf 100644
--- a/src/algebra/irexpand.spad
+++ b/src/algebra/irexpand.spad
@@ -34,6 +34,8 @@ IntegrationResultToFunction(R, F) : Exports == Implementation where
complexExpand : IR -> F
++ complexExpand(i) returns the expanded complex function
++ corresponding to i.
+ quartic2 : (UP, UP, Symbol) -> List(F)
+ ++ quartic2(p, lg, x) should be local but conditional
Implementation ==> add
import from AlgebraicManipulations(R, F)
@@ -150,10 +152,121 @@ IntegrationResultToFunction(R, F) : Exports == Implementation where
bb := -(r.coef1) k
tantrick(aa * a + bb * b, g) + ilog0(aa, bb, r.coef2, -r.coef1, k)
+ -- Computes (a + b*%i)log(lg(a + b*%i)) + (a - b*%i)log(lg(a - b*%i))
+ -- using Rioboo transformation.
+ root_pair(lg : UP, a : F, b : F, x : Symbol) : F ==
+ lge := quadeval(lg, a, b, -1)
+ f := lge.ans1
+ g := lge.ans2
+ a*log(f*f + g*g) + b*ilog(f, g, x)
+
+ lg2cfunc2 : (UP, UP) -> F
+
+ root4(a : F) : F ==
+ rec := froot(a, 4)$PolynomialRoots(IndexedExponents(K), K, R, P, F)
+ p1 := monomial(1, rec.exponent)$UP - rec.radicand::UP
+ rec.coef*zeroOf(p1)
+
+ SUP ==> SparseUnivariatePolynomial
+
+ if F has PolynomialFactorizationExplicit then
+
+ dummy := create()$SingletonAsOrderedSet
+
+ my_sqrt(s : F) : F ==
+ p := monomial(1, 2)$UP - s::UP
+ fr := factor(p)
+ fl := factorList(fr)
+ #fl = 1 => sqrt(s)
+ fr1 := fl(1).factor
+ 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, lg : UP, x : Symbol
+ ) : List(F) ==
+ rr := my_sqrt(r1)
+ qq := r1 + (1/2::F)*a
+ q1 := qq + b*rr/(4*r1)
+ q2 := qq - b*rr/(4*r1)
+ [root_pair(lg, rr - a0/(4::F), my_sqrt(q1), x) +
+ root_pair(lg, -rr - a0/(4::F), my_sqrt(q2), x)]
+
+
+ quartic2(p : UP, lg : UP, x : Symbol) : List(F) ==
+ lc := leadingCoefficient(p)
+ p := (1/lc)*p
+ a0 := coefficient(p, 3)
+ p := eval(p, dummy, monomial(1, 1)$UP - (a0/4::F)::UP)
+ 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 := r + (-16*c+4*a^2)*xx
+ r := r - (b^2)::UP
+ fr := factor(r)
+ fl := factorList(fr)
+ #fl = 1 => []
+ #fl = 3 =>
+ r1 : F
+ cnt : Integer := 0
+ for fac in fl repeat
+ f1 := fac.factor
+ cc := coefficient(f1, 0)
+ (su := sign(cc)) case "failed" => return []
+ if su@Integer < 0 then
+ cnt > 0 => return []
+ cnt := 1
+ r1 := cc
+ cnt = 0 => error "impossible"
+ quartic3(-r1, a0, a, b, lg, x)
+ #fl = 2 =>
+ r1 : F
+ for fac in fl repeat
+ f1 := fac.factor
+ (degree(f1) = 1) =>
+ r1 := coefficient(f1, 0)
+ (su := sign(r1)) case "failed" => return []
+ if not (su@Integer < 0) then return []
+ del := coefficient(f1, 1)^2 -
+ 4*coefficient(f1, 2)*coefficient(f1, 0)
+ if not((su := sign(del)) case "failed") then
+ if (su@Integer < 0) then iterate
+ su := sign(coefficient(f1, 1))
+ su case "failed" or not(0 < su@Integer) => return []
+ su := sign(coefficient(f1, 2)*coefficient(f1, 0))
+ if su case "failed" or not(0 < su@Integer) then
+ return []
+ quartic3(-r1, a0, a, b, lg, x)
+ error "impossible"
+
+ else
+
+ quartic2(p : UP, lg : UP, x : Symbol) : List(F) == []
+
+ quartic(p : UP, lg : UP, x : Symbol) : List(F) ==
+ ground?(rp := reductum(p)) =>
+ a := ground(rp)
+ s := sign(a)
+ s case "failed" => [lg2cfunc2(p, lg)]
+ si := s@Integer
+ si = 1 =>
+ r1 := root4(a/(4::F*leadingCoefficient(p)))
+ [root_pair(lg, r1, r1, x) + root_pair(lg, -r1, r1, x)]
+ si = -1 =>
+ r1 := rootSimp(zeroOf(p))
+ [cmplex(r1, lg) + cmplex(-r1, lg) + root_pair(lg, 0, r1, x)]
+ error "impossible"
+ res1 := quartic2(p, lg, x)
+ not(empty?(res1)) => res1
+ [lg2cfunc2(p, lg)]
+
lg2func(lg, x) ==
zero?(d := degree(p := lg.coeff)) => error "poly has degree 0"
(d = 1) => [linear(p, lg.logand)]
d = 2 => quadratic(p, lg.logand, x)
+ d = 4 => quartic(p, lg.logand, x)
odd? d and
((r := retractIfCan(reductum p)@Union(F, "failed")) case F) =>
pairsum([cmplex(alpha := rootSimp zeroOf p, lg.logand)],
@@ -162,8 +275,10 @@ IntegrationResultToFunction(R, F) : Exports == Implementation where
lg.logand], x))
[lg2cfunc lg]
- lg2cfunc lg ==
- +/[cmplex(alpha, lg.logand) for alpha in zerosOf(lg.coeff)]
+ lg2cfunc(lg) == lg2cfunc2(lg.coeff, lg.logand)
+
+ lg2cfunc2(p : UP, lg : UP) ==
+ +/[cmplex(alpha, lg) for alpha in zerosOf(p)]
mkRealFunc(l, x) ==
ans := empty()$List(F)
