errors in rsimp.spad

[Qian Yun]

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)

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)

- Qian

[Waldek Hebisch]

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

[Waldek Hebisch]

On Mon, Dec 18, 2023 at 06:37:28PM +0100, 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...

Dependent roots are introduced by 'hackroot' in src/algebra/algfunc.spad
AFAICS current code is from open source Axiom era. NAG version
was doing no splitting, but when fraction had nominator 1 or -1 it
replaced root by inverse of root of inverse, which caused undesirable
effect for things like sqrt(-1/x), so the code was replaced by
version doing splitting.

I think now that we should disable this splitting. ATM I am not
sure if we should just drop this part. Or maybe restore restricted
variant of old code.

Anyway, current

(5) -> (-1/108)^(1/5)

5+---+
\|- 1
(5) ------
5+---+
\|108
Type: AlgebraicNumber
(6) -> sqrt(3/1501)

+-+
\|3
(6) -------
+----+
\|1501
Type: AlgebraicNumber

is arguably wrong: root of degree 5 satisfies equation of degree 5, while
above we produce two root suggesting extention of degree 25. If they are
independent we merely are doing more work than needed. With dependent
roots various things no longer work. In particular, eqation for
(-1)^(1/5) is reducible, so we should avoid producing such things.

Similarly in the second case, we get two roots giving extention of
degree 4. This time roots are independent, but we if we already
have sqrt(2) and sqrt(6), then we will get dependence. And for
things like

(7) -> sqrt(3^(2/3)/11)

+-----+
|3+-+2
\|\|3
(7) --------
+--+
\|11
Type: AlgebraicNumber

we get reducible exuation for the nominator, while expression for
the whole root is irreducible.

--
Waldek Hebisch

[Qian Yun]

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"

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))
=====

>> 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?

Otherwise, your patch works.

- Qian

[Qian Yun]

If revert to old behavior, I expect there will be lots of test
breakage?

I have not tested it -- just think this will have extensive
impacts in various places...

(I'm not saying we shouldn't do this -- I think we should
do mathematically correct changes, and this is one of them.)

- Qian

[Waldek Hebisch]

On Tue, Dec 19, 2023 at 06:18:09PM +0800, Qian Yun wrote:
> If revert to old behavior, I expect there will be lots of test
> breakage?
>
> I have not tested it -- just think this will have extensive
> impacts in various places...
>
> (I'm not saying we shouldn't do this -- I think we should
> do mathematically correct changes, and this is one of them.)

Extensive impact, yes. AFAICS no real breakage. One test
in 'integ.input' fails but can be changed to pass. Results in
fixed, mapleok, r20abugs, tutchap3 have large changes.
tutchap3 in a sense is an improvement: after change comments
agree with computed results. Other are neutral. Changed
results are slightly bigger, but not much.

Most changes are to results from integrator, here we could
improve results by using extra simplification pass in
final processing (after irexpand).

--
Waldek Hebisch

[Waldek Hebisch]

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

[Qian Yun]

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.)

- Qian

[Qian Yun]

Oh, they are fixed by the most recent commit
"Disable splitting of roots".

But here is the list that broke by it:
(error: impossible)

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)
...

- Qian

[Qian Yun]

Ignore this mail, I forgot to apply sum6a.diff.

- Qian

[Waldek Hebisch]

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.)

Yes, thanks. In trunk case of biquadratic with two complex
roots is unhandled. We can send it to 'quartic2'. But
I am looking at 'quartic2', currently it assumes that
all roots are complex and was intended to fail otherwise.
But it seems that I missed some cases where 'quartic2'
goes on but roots are real. And in fact, it should be
possible to handle real roots in 'quartic2'.

--
Waldek Hebisch

[Waldek Hebisch]

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

[Qian Yun]

Yes, with this patch, there's much improvement.

And only 2 cases are different, are they expected?

integrate(1/(a^5+x^5),x)
integrate((-x^4+p*x^2+1)^(1/2)/(x^4+1),x)

typo: "od" -> "of"

Minor problem: double underscore in "r__rec". BTW, why not capitalize
to "R_rec"?

- Qian

[Waldek Hebisch]

On Thu, Dec 21, 2023 at 06:16:16PM +0800, Qian Yun wrote:
> Yes, with this patch, there's much improvement.
>
> And only 2 cases are different, are they expected?
>
> integrate(1/(a^5+x^5),x)
> integrate((-x^4+p*x^2+1)^(1/2)/(x^4+1),x)

They both are parametric integrals and AFAICS we can not determine
configuration of roots, so transformations do no apply. There
is possible change because now 'real' is gone from rational
function integrator.

--
Waldek Hebisch

[Qian Yun]

Off topic a little bit:

For parametric integrals, do you think these are good ideas?

1. Integration by substitution.

For 1/(a^5+x^5), recognize and transform it to 1/a^5*1/(1+(x/a)^5).
Substitute x/a with y.

2. partialFraction instead of current resultant method.

e.g. integrate(1/((x-a)*(x-b)*(x-c)*(x-d)),x)

- Qian

[Waldek Hebisch]

On Thu, Dec 21, 2023 at 08:32:46PM +0800, Qian Yun wrote:
> Off topic a little bit:
>
> For parametric integrals, do you think these are good ideas?
>
> 1. Integration by substitution.
>
> For 1/(a^5+x^5), recognize and transform it to 1/a^5*1/(1+(x/a)^5).
> Substitute x/a with y.

If you literally mean cases like 1/(a^n + x^n), then this is probably
too special. If you want more general cases, than is becomes
more tricky. For example a lot of integrals in Rubi testsute
contains 'a*x + b' as a factors. One could try substitution
'y = a*x + b', but this may complicate other parts of the integral.

Currently FriCAS useses substitution when they give "substantial"
simplification:
- allow completely getting rid of a transcendental kernel
- reduce degree of algebriac kernels

> 2. partialFraction instead of current resultant method.
>
> e.g. integrate(1/((x-a)*(x-b)*(x-c)*(x-d)),x)

That is promising idea. In principle partialFraction should
always be faster than resultant and AFAICS we will not
loose elementary integrability by splitting into partial fractions.
But it needs testing.

--
Waldek Hebisch