On Tue, Dec 19, 2023 at 06:14:16PM +0800, Qian Yun wrote:
>
>
> On 12/19/23 01:37, Waldek Hebisch wrote:
> > 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.
>
> Do you think this part is problematic?
> "alpha := rootSimp zeroOf p"
Yes, we should have something better. We probably should limit
work done in 'irexpand.spad'. Namely, by necessity routines there
work locally without view of the whole expression. I think
that outside we should do cleanup removing dependencies between
roots.
> for p is (x^5+1), alpha is (-1)^(1/5) instead of -1.
> (unlike radicalSovle)
>
> With this patch I can get better result.
>
> =====
> diff --git a/src/algebra/irexpand.spad b/src/algebra/irexpand.spad
> index 4195f3b7..760f14bd 100644
> --- a/src/algebra/irexpand.spad
> +++ b/src/algebra/irexpand.spad
> @@ -185,7 +185,12 @@ IntegrationResultToFunction(R, F) : Exports ==
> Implementation where
> 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)],
> + sn := sign(r)
> + if sn case Integer and (sn::Integer) > 0 then
> + alpha := - rootSimp nthRoot(r, d)
> + else
> + alpha := rootSimp zeroOf p
> + pairsum([cmplex(alpha, lg.logand)],
> lg2func([lg.scalar,
> (p exquo (monomial(1, 1)$UP - alpha::UP))::UP,
> lg.logand], x))
> =====
For now it looks good.
> > > 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,
>
> Just to clarify, you revert some class of biquadratic to old behavior,
> right?
'quartic2' only handles case of complex roots, if is did anything to
equations with reals roots, then this was a bug. In the patch added
case for biquadratic should handle all cases that where correctly
handled by version in trunk and handle also cases with real roots.
The only unhandled biquadratics are ones where we can not determine
configuration of roots. I think that version in trunk is unable
to handle them too.
--
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/ZYHa8_w5Anv6ay1b%40fricas.org.
