Reorganized and somewhat extended version. I have split purely
algebraic part to a separate file and added another case.
In principle we could add a few more cases, but unless somebody
want to do this or modify the patch I will commit it.
--
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 fricas-devel+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/fricas-devel/ZX4i3TCpfd1k3unI%40fricas.org.
diff --git a/src/algebra/Makefile.in b/src/algebra/Makefile.in
index 45b1298c..6fe52f39 100644
--- a/src/algebra/Makefile.in
+++ b/src/algebra/Makefile.in
@@ -75,7 +75,7 @@ SPAD_SRCS= \
ptranfn puiseux qalgset quat radeigen radix \
random ratfact rdeefx real0q \
realzero reclos rec regset rep1 rep2 \
- resring retract rf riccati rinterp \
+ resring retract rf riccati rinterp rsimp \
rule scene seg serexp setorder sets sex sf \
sgcf sign si skpol smith smith2 solvedio \
solvelin solverad sortpak space special special2 sregset \
@@ -266,7 +266,7 @@ SPADLIST7=\
RESULT RETFROM RETRACT RETSOL RFDIST RFFACT \
RFSSPLIT RF RGCHAIN RIDIST RING RINTERP RMATCAT \
RMATRIX RMCAT2 RMODULE RNG RNORM RNS ROIRC \
- ROMAN RPOLCAT RRCC RSDCMPK RSETCAT \
+ ROMAN ROOTUT RPOLCAT RRCC RSDCMPK RSETCAT \
RSETGCD RULECOLD RULESET RULE RURPK \
SAE SAOS SARGND SBOUND SCACHE SCANUTIL \
SCELL SCENE SCIFS SCNP SCONF SCPKG SCRT SDPOL \
diff --git a/src/algebra/irexpand.spad b/src/algebra/irexpand.spad
index ead11d7f..756d5463 100644
--- a/src/algebra/irexpand.spad
+++ b/src/algebra/irexpand.spad
@@ -150,10 +150,39 @@ 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
+
+ C_rec ==> Record(real : F, imag : F)
+ r_rec ==> Record(reals : List(F), complexes : List(C_rec))
+
+ do_rr(rr : r_rec, lg : UP, x : Symbol) : List(F) ==
+ rl := rr.reals
+ cl := rr.complexes
+ res := 0$F
+ for r1 in rl repeat
+ res := res + cmplex(r1, lg)
+ for cp1 in cl repeat
+ res := res + root_pair(lg, cp1.real, cp1.imag, x)
+ [res]
+
+ quartic(p : UP, lg : UP, x : Symbol) : List(F) ==
+ rru := complex_roots(p)$RootUtilities(R, F)
+ (rru case "failed") => [lg2cfunc2(p, lg)]
+ do_rr(rru, lg, x)
+
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 +191,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)
