better integration result -- avoid complicated roots

[Qian Yun]

Some of the integration results are pages long --
because IntegrationResult often has a 'logpart', aka
a rootSum form.

'lg2func' in IR2F uses 'zeroOf' which can introduce
large radicals into the result.

For example

internalIntegrate(1/(z^4+1),z)

internalIntegrate((1+z^2)^4/ (1+z^4)^4,z)

Until recently I found that Rioboo's algorithm described in
Manuel Bronstein's book, at the end of section 2.8, covers
this situation.

The core idea is that for polynomials with real coefficient,
its roots come in pairs (symmetric over the real axis).
So in the following computation of the rootSum of ilogs,
many terms can cancel with each other.

Sympy implements it, see 'log_to_real':

https://github.com/sympy/sympy/blob/67e4a5c07c390d648601e8fbd6937c877c7438bb/sympy/integrals/rationaltools.py#L327C15-L327C15

Of course, this algorithm is not limited to rational function
integration, it can improve other integrals as well.

I have worked out a version in FriCAS, it's not a replace for
'lg2func' yet -- so my new function 'lg2func2' will try first,
if it can't handle, then pass to 'lg2func'.

The very rough first version only has 3 failures. I'm still
working on them.

Now, some questions:

My code needs to use solve$TransSolvePackage. Because the
resultant method mentioned in the book can't handle
coefficients of EXPR. Are there alternatives?

And will solve$TransSolvePackage only introduce necessary extra
root kernels?

TransSolvePackage is fixed to work with EXPR, not a generalized
domain like the second parameter to IR2F.

This doesn't work very well with integrals with parameters,
because there will be branches -- is sqrt(-a) real or not?

For even higher degree polynomials, I think it's wise to
return rootSum instead of give results containing algebraic
numbers, which causes trouble for sagemath.

I'll post patch after I reviewed the test failures and do more
testing.

I think this will has positive effect on Nasser's CAS independent
integration tests, to what extent, I don't know.

- Qian

This is the result with patch:

(1) -> integrate(1/(z^4+1),z)

(1)
+-+ 2 +-+ 2 +-+
log(z\|2 + z + 1) - log(- z\|2 + z + 1) + 2 atan(z\|2 + 1)
+
+-+
2 atan(z\|2 - 1)
/
+-+
4 \|2

(2) -> integrate((1+z^2)^4/ (1+z^4)^4,z)

(2)
12 8 4 +-+ 2
(33 z + 99 z + 99 z + 33)log(z\|2 + z + 1)
+
12 8 4 +-+ 2
(- 33 z - 99 z - 99 z - 33)log(- z\|2 + z + 1)
+
12 8 4 +-+
(306 z + 918 z + 918 z + 306)atan(z\|2 + 1)
+
12 8 4 +-+
(306 z + 918 z + 918 z + 306)atan(z\|2 - 1)
+
11 9 7 5 3 +-+
(240 z + 124 z + 672 z + 264 z + 432 z + 12 z)\|2
/
12 8 4 +-+
(384 z + 1152 z + 1152 z + 384)\|2

[Waldek Hebisch]

On Sun, Nov 19, 2023 at 07:17:14PM +0800, Qian Yun wrote:
> Some of the integration results are pages long --
> because IntegrationResult often has a 'logpart', aka
> a rootSum form.
>
> 'lg2func' in IR2F uses 'zeroOf' which can introduce
> large radicals into the result.
>
> For example
>
> internalIntegrate(1/(z^4+1),z)
>
> internalIntegrate((1+z^2)^4/ (1+z^4)^4,z)
>
> Until recently I found that Rioboo's algorithm described in
> Manuel Bronstein's book, at the end of section 2.8, covers
> this situation.

Hmm, I am not sure what you want to say here. We use Rioboo's
algorithm, but in some cases it does not give us enough
improvement.

AFAICS most significant reason for troubles is "atan doubling"
or even tripling: trying to keep things real we get square
(or higer power) of argument to logaritm. And we use changes
of variables which require kind of "to real" transform when
coming back. This in principle could be solved by factoring
arguments of logs.

> The core idea is that for polynomials with real coefficient,
> its roots come in pairs (symmetric over the real axis).
> So in the following computation of the rootSum of ilogs,
> many terms can cancel with each other.
>
> Sympy implements it, see 'log_to_real':
>
> https://github.com/sympy/sympy/blob/67e4a5c07c390d648601e8fbd6937c877c7438bb/sympy/integrals/rationaltools.py#L327C15-L327C15
>
> Of course, this algorithm is not limited to rational function
> integration, it can improve other integrals as well.
>
> I have worked out a version in FriCAS, it's not a replace for
> 'lg2func' yet -- so my new function 'lg2func2' will try first,
> if it can't handle, then pass to 'lg2func'.
>
> The very rough first version only has 3 failures. I'm still
> working on them.
>
> Now, some questions:
>
> My code needs to use solve$TransSolvePackage. Because the
> resultant method mentioned in the book can't handle
> coefficients of EXPR. Are there alternatives?
>
> And will solve$TransSolvePackage only introduce necessary extra
> root kernels?

I would rather avoid TransSolvePackage. Main approach to Expression
is convert it to polynomial (with Expression coefficients), say using
'univariate'. Once you have polynomials you can use resultants
or Groebner bases to solve systems of equations.

> TransSolvePackage is fixed to work with EXPR, not a generalized
> domain like the second parameter to IR2F.
>
> This doesn't work very well with integrals with parameters,
> because there will be branches -- is sqrt(-a) real or not?

Yes, integrals with parameters cause troubles. IMO those troubles
can be mostly ignored for indefinite integration. We can not
avoid branch troubles for definite integration, but for definite
integration we have more info (bound of integration interval)
and some extra possiblities (for example split integral into
sum of integrals on smaller intervals).

Anyway, my opinion is that we should do generaly applicable
simplification on indefinte integral, but leave thing which
say require looking at sign to definite integrator.

> For even higher degree polynomials, I think it's wise to
> return rootSum instead of give results containing algebraic
> numbers, which causes trouble for sagemath.

I posted a patch to do this some time ago. But it causes
regression in definite integrator. AFAICS, currently we
pass complex algebraic numbers to limit code which is wrong.
But in some cases we get sensible results. Once we have
rootSum limit code sees that there is trouble and gives up.
General case of limit of rootSum looks tricky. There is
good chance that with moderate effort we could improve
limit to handle cases that currently work using algebraic
numbers (and fail on other cases).

> I'll post patch after I reviewed the test failures and do more
> testing.

Could you test interaction of your patch with attached patch.
Attached patch remove (I think) last call to 'real' from the
integrator. We should do this as call to 'real' causes several
wrong results. But removing call to 'real' breaks some definite
integrals, so I am waiting with this patch for some solution for
definite integrals.

Removing call to 'real' means that following transformations of
integrals will see complex expression. For example, with 'real'
removed we have:

(2) -> integrate(1/(z^4+1),z)

(2)
+---+ +-+ +---+ +-+
(\|- 1 + 1)\|2 log((\|- 1 + 1)\|2 + 2 z)
+
+---+ +-+ +---+ +-+
(\|- 1 - 1)\|2 log((\|- 1 - 1)\|2 + 2 z)
+
+---+ +-+ +---+ +-+
(- \|- 1 + 1)\|2 log((- \|- 1 + 1)\|2 + 2 z)
+
+---+ +-+ +---+ +-+
(- \|- 1 - 1)\|2 log((- \|- 1 - 1)\|2 + 2 z)
/
8


--
Waldek Hebisch

[Qian Yun]

On 11/19/23 20:28, Waldek Hebisch wrote:
>>
>> Until recently I found that Rioboo's algorithm described in
>> Manuel Bronstein's book, at the end of section 2.8, covers
>> this situation.
>
> Hmm, I am not sure what you want to say here. We use Rioboo's
> algorithm, but in some cases it does not give us enough
> improvement.

My mistake. We are using Rioboo's algorithm for rational integration.
I was talking about section 2.8 "Riboo's Algorithm for Real Rational
Functions", page 65 to 70 for second edition of "Symbolic
Integration I".

> AFAICS most significant reason for troubles is "atan doubling"
> or even tripling: trying to keep things real we get square
> (or higer power) of argument to logaritm. And we use changes
> of variables which require kind of "to real" transform when
> coming back. This in principle could be solved by factoring
> arguments of logs.

This algorithm prevents generation of (unnecessary) complex
expression from the beginning.

>
>> For even higher degree polynomials, I think it's wise to
>> return rootSum instead of give results containing algebraic
>> numbers, which causes trouble for sagemath.
>
> I posted a patch to do this some time ago. But it causes
> regression in definite integrator. AFAICS, currently we
> pass complex algebraic numbers to limit code which is wrong.
> But in some cases we get sensible results. Once we have
> rootSum limit code sees that there is trouble and gives up.
> General case of limit of rootSum looks tricky. There is
> good chance that with moderate effort we could improve
> limit to handle cases that currently work using algebraic
> numbers (and fail on other cases).

Now my patch has no test failures, because I limit it to
work on polynomials degree >= 4. Essentially replacing "lg2cfunc".
(For degree < 4, my patch is not wrong -- just not as optimized
as "quadratic".)

>> I'll post patch after I reviewed the test failures and do more
>> testing.
>
> Could you test interaction of your patch with attached patch.
> Attached patch remove (I think) last call to 'real' from the
> integrator. We should do this as call to 'real' causes several
> wrong results. But removing call to 'real' breaks some definite
> integrals, so I am waiting with this patch for some solution for
> definite integrals.

Brief test shows no bad interaction -- my patch works at a place
before yours have effect.

> Removing call to 'real' means that following transformations of
> integrals will see complex expression. For example, with 'real'
> removed we have:
>
> (2) -> integrate(1/(z^4+1),z)
>
> (2)
> +---+ +-+ +---+ +-+
> (\|- 1 + 1)\|2 log((\|- 1 + 1)\|2 + 2 z)
> +
> +---+ +-+ +---+ +-+
> (\|- 1 - 1)\|2 log((\|- 1 - 1)\|2 + 2 z)
> +
> +---+ +-+ +---+ +-+
> (- \|- 1 + 1)\|2 log((- \|- 1 + 1)\|2 + 2 z)
> +
> +---+ +-+ +---+ +-+
> (- \|- 1 - 1)\|2 log((- \|- 1 - 1)\|2 + 2 z)
> /
> 8
>
>

Isn't a real result even better than this?

My patch is in attachment, the variable names are bad, I know.

- Qian

[Waldek Hebisch]

On Sun, Nov 19, 2023 at 08:47:26PM +0800, Qian Yun wrote:
>
>
> On 11/19/23 20:28, Waldek Hebisch wrote:
> I was talking about section 2.8 "Riboo's Algorithm for Real Rational
> Functions", page 65 to 70 for second edition of "Symbolic
> Integration I".
>
> > AFAICS most significant reason for troubles is "atan doubling"
> > or even tripling: trying to keep things real we get square
> > (or higer power) of argument to logaritm. And we use changes
> > of variables which require kind of "to real" transform when
> > coming back. This in principle could be solved by factoring
> > arguments of logs.
>
> This algorithm prevents generation of (unnecessary) complex
> expression from the beginning.

"atan doubling" makes thing harder (increases degrees) and after
doubling Rioboo was failing. Doubling may happen both before
irexpand and after.

> > > For even higher degree polynomials, I think it's wise to
> > > return rootSum instead of give results containing algebraic
> > > numbers, which causes trouble for sagemath.
> >
> > I posted a patch to do this some time ago. But it causes
> > regression in definite integrator. AFAICS, currently we
> > pass complex algebraic numbers to limit code which is wrong.
> > But in some cases we get sensible results. Once we have
> > rootSum limit code sees that there is trouble and gives up.
> > General case of limit of rootSum looks tricky. There is
> > good chance that with moderate effort we could improve
> > limit to handle cases that currently work using algebraic
> > numbers (and fail on other cases).
>
> Now my patch has no test failures, because I limit it to
> work on polynomials degree >= 4. Essentially replacing "lg2cfunc".
> (For degree < 4, my patch is not wrong -- just not as optimized
> as "quadratic".)

There are several changes in 'mapleok.input'. Did you look
that they all are improvements?

> > > I'll post patch after I reviewed the test failures and do more
> > > testing.
> >
> > Could you test interaction of your patch with attached patch.
> > Attached patch remove (I think) last call to 'real' from the
> > integrator. We should do this as call to 'real' causes several
> > wrong results. But removing call to 'real' breaks some definite
> > integrals, so I am waiting with this patch for some solution for
> > definite integrals.
>
> Brief test shows no bad interaction -- my patch works at a place
> before yours have effect.

Well, there is interaction, in the following integral.

> > Removing call to 'real' means that following transformations of
> > integrals will see complex expression. For example, with 'real'
> > removed we have:
> >
> > (2) -> integrate(1/(z^4+1),z)
> >
> > (2)
> > +---+ +-+ +---+ +-+
> > (\|- 1 + 1)\|2 log((\|- 1 + 1)\|2 + 2 z)
> > +
> > +---+ +-+ +---+ +-+
> > (\|- 1 - 1)\|2 log((\|- 1 - 1)\|2 + 2 z)
> > +
> > +---+ +-+ +---+ +-+
> > (- \|- 1 + 1)\|2 log((- \|- 1 + 1)\|2 + 2 z)
> > +
> > +---+ +-+ +---+ +-+
> > (- \|- 1 - 1)\|2 log((- \|- 1 - 1)\|2 + 2 z)
> > /
> > 8
> >
> >
>
> Isn't a real result even better than this?

Real result is better. But applying 'real' is wrong, we need
honest simplification that preserves value. This is one example
which code in trunk was unable to simplify without 'real'.
AFAICS your patch handles this.

--
Waldek Hebisch

[Qian Yun]

On 11/20/23 02:19, Waldek Hebisch wrote:
>
> "atan doubling" makes thing harder (increases degrees) and after
> doubling Rioboo was failing. Doubling may happen both before
> irexpand and after.
>

What is "atan doubling"?

>
> There are several changes in 'mapleok.input'. Did you look
> that they all are improvements?
>

I'll take some time to do more testing and report back.

- Qian

[Waldek Hebisch]

On Sun, Nov 19, 2023 at 07:17:14PM +0800, Qian Yun wrote:
>

> The core idea is that for polynomials with real coefficient,
> its roots come in pairs (symmetric over the real axis).
> So in the following computation of the rootSum of ilogs,
> many terms can cancel with each other.

<snip<

> Now, some questions:
>
> My code needs to use solve$TransSolvePackage. Because the
> resultant method mentioned in the book can't handle
> coefficients of EXPR. Are there alternatives?
>
> And will solve$TransSolvePackage only introduce necessary extra
> root kernels?

AFAICS TransSolvePackage is potentially problematic. In general,
working simultanelously with all roots of polynomials is
problematic. Done rigorously it is equivalent to computing
in splitting field, which usually is quite expensive due to
high degree algberaic extentions and need for algebraic
factoring. AFAIK our solvers do not compute splitting field,
and "all roots" commands are fake: they are reasonable when
one works with single root, but may give wrong result for
expressions involving multiple roots. Is seems that
TransSolvePackage uses formulas for solution of degree 3
and degree 4 polynomial. Those formulas are problematic too,
IIRC one has to treat them in special way to avoid troubles.

Also, you (and IIRC TransSolvePackage) use eval. 'eval' is
doing various simplifications which may transform kernels
in undesired ways. In integrator I switched to using 'subst',
but such change requires care, so we can not just mass-edit
the algebra. We could use something like 'quadeval' to
avoid 'eval'.

There is question in which cases your apprach works. Examples
in mapleok.input which work have minimal polynomial for residues
which is variation of x^4 + a^4 where a is a rational constant
or x^4 + 2*x^2 + 2 or x^4 + x^2/64 + 1/8192. x^4 + a^4 factors
in extention by square root of 2:

factor(x^4 + a^4, [sqrt(2)])

2 +-+ 2 2 +-+ 2
(13) (x - \|2 a x + a )(x + \|2 a x + a )

The two other need extention by sqrt(sqrt(2) - 1).

In general, your approach builds equations so
that real u, v

p(u + i*v) = 0 <=> PP(u, v) = 0 & QQ(u, v) = 0

Computing Groebner basis from PP and QQ we can find polynomial
r such that

r(u) = 0 <=> exists v p(u + i*v) = 0

AFAICS your approach needs set of real roots of r (so splitting
of a factor of r) and for each real root of r pair of
real v.

However, in

integrate(((-3)*x^3+6)/(x^6+(-1)*x^3+1), x)

minimal polynomial of residues is cyclotomic(9) (minimal
polynomial for primitive root of 1 of degree 9). Equation
for real parts of roots is 8*u^3 - 6*u + 1. But solution
via radicals looks complex:

(44) -> solve(8*u^3 - 6*u + 1, 1.0e-12)

(44)
[u = - 0.9396926207_8599058441, u = 0.7660444431_1897168518,
u = 0.1736481776_6656415188]

(48) -> radicalSolve(8*u^3 - 6*u + 1)

(48)
+----------+2
+---+ 3| +---+ 3+-+2
(- \|- 3 + 1)\|\|- 3 - 1 - 2 \|2
[u = --------------------------------------,
+----------+
+---+ 3+-+3| +---+
(2 \|- 3 + 2)\|2 \|\|- 3 - 1
+----------+2 +----------+2
+---+ 3| +---+ 3+-+2 3| +---+ 3+-+2
(- \|- 3 - 1)\|\|- 3 - 1 + 2 \|2 \|\|- 3 - 1 + \|2
u = --------------------------------------, u = ----------------------]
+----------+ +----------+
+---+ 3+-+3| +---+ 3+-+3| +---+
(2 \|- 3 - 2)\|2 \|\|- 3 - 1 2 \|2 \|\|- 3 - 1

It is not clear to me if TransSolvePackage produces the same
solutions, but it looks that solutions from TransSolvePackage
contain complex quantities.

Groebner produces sensible thing:

(50) -> groebner([PP::POLY(INT), QQ::POLY(INT), 8*u^3 - 6*u + 1])

2 2 3
(50) [v + u - 1, 8 u - 6 u + 1]

and solution of this could be passed to Rioboo. But TransSolvePackage
can not usefully handle it.

--
Waldek Hebisch

[Waldek Hebisch]

On Mon, Nov 20, 2023 at 08:48:30AM +0800, Qian Yun wrote:
>
>
> On 11/20/23 02:19, Waldek Hebisch wrote:
> >
> > "atan doubling" makes thing harder (increases degrees) and after
> > doubling Rioboo was failing. Doubling may happen both before
> > irexpand and after.
> >
>
> What is "atan doubling"?

instead of

log(u)

one can take

(1/2)*log(u^2)

as integration result. This tends to happen for logs that will
be later converted to atans. It seems that it may happen that
u is irrational, u^2 is rational and atan produced from u
has rational argument. If u is rational then factoring
argument to log should simplify it. But if u is irrational,
then things are more tricky.

> > There are several changes in 'mapleok.input'. Did you look
> > that they all are improvements?
> >
>
> I'll take some time to do more testing and report back.

I looked at them and they seem good.

--
Waldek Hebisch

[Qian Yun]

On 11/20/23 10:08, Waldek Hebisch wrote:
>
> There is question in which cases your apprach works. Examples
> in mapleok.input which work have minimal polynomial for residues
> which is variation of x^4 + a^4 where a is a rational constant
> or x^4 + 2*x^2 + 2 or x^4 + x^2/64 + 1/8192. x^4 + a^4 factors
> in extention by square root of 2:

In theory this approach should work with higher degree of polynomial.
After all as you said, it is about finding a algebraic extension to
factor the polynomial into simpler ones.

Let's see if I can find out such an example.

>
> Groebner produces sensible thing:
>
> (50) -> groebner([PP::POLY(INT), QQ::POLY(INT), 8*u^3 - 6*u + 1])
>
> 2 2 3
> (50) [v + u - 1, 8 u - 6 u + 1]
>
> and solution of this could be passed to Rioboo. But TransSolvePackage
> can not usefully handle it.
>

How can this go further? Return rootSum in Rioboo?

- Qian

[Waldek Hebisch]

On Tue, Nov 21, 2023 at 08:33:50AM +0800, Qian Yun wrote:
>
>
> On 11/20/23 10:08, Waldek Hebisch wrote:
> >
> > There is question in which cases your apprach works. Examples
> > in mapleok.input which work have minimal polynomial for residues
> > which is variation of x^4 + a^4 where a is a rational constant
> > or x^4 + 2*x^2 + 2 or x^4 + x^2/64 + 1/8192. x^4 + a^4 factors
> > in extention by square root of 2:
>
> In theory this approach should work with higher degree of polynomial.
> After all as you said, it is about finding a algebraic extension to
> factor the polynomial into simpler ones.

Yes. However, explict formulas for root are possible only in
special cases. And even when possible are likely to use
complex expressions for real roots. AFAICS biquadratic case,
that is degree 4 polynomial containing only even powers
should work. I am not sure if out solver is smart enough,
but symmetric degree 4 polynomial should work. Basically
polynomials which one can solve repeatedly solving quadratic
equations. Once solution involves general formula for degree 3
or degree 4 polynomial we will get complex expressions.

> Let's see if I can find out such an example.
>
> >
> > Groebner produces sensible thing:
> >
> > (50) -> groebner([PP::POLY(INT), QQ::POLY(INT), 8*u^3 - 6*u + 1])
> >
> > 2 2 3
> > (50) [v + u - 1, 8 u - 6 u + 1]
> >
> > and solution of this could be passed to Rioboo. But TransSolvePackage
> > can not usefully handle it.
> >
>
> How can this go further? Return rootSum in Rioboo?

Possibly. Or do what we do now, that is use symbolic roots.
For degree 3 after you add symbolic root polynomial splits
into linear term and quadratic (or into linear terms).

rootSum is the future, but as long as limit can not handle
rootSum we have to use symbolic roots.

AFAICS we should first set up equations, then use Groebner,
then factor resulting univariate polynomial. Unfortunately,
beside needed roots this polynomial contains extra roots
and IMO main issue is to get rid of extra roots. In case above
this was moderately easy: only one factor had real roots.
As long as polynomial has constant coefficients and we are
able to resolve branch problems for coefficient you can
use semi-numeric methods to check which factors have real
roots.

There are theoretical questions, in examples I looked at
there was always a single factor giving good real roots.
Situation would be much worse if some factor gave say one
good root and a bunch of bad (complex) roots.

For "generic" degree 4 polynomial, that is

p := x^4 + a*x^3 + b*x^2 + c*x + d

I got the following univarite equation for real parts of
roots:

2 3 2 2
6 3 a 5 2 b + 3 a 4 4 a b + a 3 - 4 d + a c + b + 2 a b 2
u + --- u + ---------- u + ---------- u + ------------------------ u
2 4 8 16
+
2 2 2 2
- 4 a d + a c + a b - a d - c + a b c
-------------------- u + ------------------
32 64

This polynomial is irreducible. Writing p in factored form:

p := (x - a)*(x - b)*(x - c)*(x - d)

We get

- d - c - d - b - d - a - c - b - c - a
(u + -------)(u + -------)(u + -------)(u + -------)(u + -------)
2 2 2 2 2
*
- b - a
(u + -------)
2

Generic polynomial of degree 4 has Galois group S(4). So at purely
algebraic level it is impossible to say which roots are paired.
In nice case the above will have factor of degree 2 with
real roots.

--
Waldek Hebisch

[Qian Yun]

On 11/21/23 09:52, Waldek Hebisch wrote:
> On Tue, Nov 21, 2023 at 08:33:50AM +0800, Qian Yun wrote:
>>
>>
>> On 11/20/23 10:08, Waldek Hebisch wrote:
>>>
>>> There is question in which cases your apprach works. Examples
>>> in mapleok.input which work have minimal polynomial for residues
>>> which is variation of x^4 + a^4 where a is a rational constant
>>> or x^4 + 2*x^2 + 2 or x^4 + x^2/64 + 1/8192. x^4 + a^4 factors
>>> in extention by square root of 2:
>>
>> In theory this approach should work with higher degree of polynomial.
>> After all as you said, it is about finding a algebraic extension to
>> factor the polynomial into simpler ones.
>
> Yes. However, explict formulas for root are possible only in
> special cases. And even when possible are likely to use
> complex expressions for real roots. AFAICS biquadratic case,
> that is degree 4 polynomial containing only even powers
> should work. I am not sure if out solver is smart enough,
> but symmetric degree 4 polynomial should work. Basically
> polynomials which one can solve repeatedly solving quadratic
> equations. Once solution involves general formula for degree 3
> or degree 4 polynomial we will get complex expressions.

How do you suggest to proceed? Write a special case for biquadratic?

>> Let's see if I can find out such an example.
>>
>>>
>>> Groebner produces sensible thing:
>>>
>>> (50) -> groebner([PP::POLY(INT), QQ::POLY(INT), 8*u^3 - 6*u + 1])
>>>
>>> 2 2 3
>>> (50) [v + u - 1, 8 u - 6 u + 1]
>>>
>>> and solution of this could be passed to Rioboo. But TransSolvePackage
>>> can not usefully handle it.
>>>
>>
>> How can this go further? Return rootSum in Rioboo?
>
> Possibly. Or do what we do now, that is use symbolic roots.
> For degree 3 after you add symbolic root polynomial splits
> into linear term and quadratic (or into linear terms).
>
> rootSum is the future, but as long as limit can not handle
> rootSum we have to use symbolic roots.
>

If current limit can handle result from integrate, then the limit
can expand the rootSum just like integrate does. Not an improvement,
but also not a step back, right?

- Qian

[Waldek Hebisch]

On Tue, Nov 21, 2023 at 06:31:27PM +0800, Qian Yun wrote:
>
>
> On 11/21/23 09:52, Waldek Hebisch wrote:
> > On Tue, Nov 21, 2023 at 08:33:50AM +0800, Qian Yun wrote:
> > >
> > >
> > > On 11/20/23 10:08, Waldek Hebisch wrote:
> > > >
> > > > There is question in which cases your apprach works. Examples
> > > > in mapleok.input which work have minimal polynomial for residues
> > > > which is variation of x^4 + a^4 where a is a rational constant
> > > > or x^4 + 2*x^2 + 2 or x^4 + x^2/64 + 1/8192. x^4 + a^4 factors
> > > > in extention by square root of 2:
> > >
> > > In theory this approach should work with higher degree of polynomial.
> > > After all as you said, it is about finding a algebraic extension to
> > > factor the polynomial into simpler ones.
> >
> > Yes. However, explict formulas for root are possible only in
> > special cases. And even when possible are likely to use
> > complex expressions for real roots. AFAICS biquadratic case,
> > that is degree 4 polynomial containing only even powers
> > should work. I am not sure if out solver is smart enough,
> > but symmetric degree 4 polynomial should work. Basically
> > polynomials which one can solve repeatedly solving quadratic
> > equations. Once solution involves general formula for degree 3
> > or degree 4 polynomial we will get complex expressions.
>
> How do you suggest to proceed? Write a special case for biquadratic?

Well, as I wrote we can recognize favourable case by factoring.
This is probably most general method.

> > > Let's see if I can find out such an example.
> > >
> > > >
> > > > Groebner produces sensible thing:
> > > >
> > > > (50) -> groebner([PP::POLY(INT), QQ::POLY(INT), 8*u^3 - 6*u + 1])
> > > >
> > > > 2 2 3
> > > > (50) [v + u - 1, 8 u - 6 u + 1]
> > > >
> > > > and solution of this could be passed to Rioboo. But TransSolvePackage
> > > > can not usefully handle it.
> > > >
> > >
> > > How can this go further? Return rootSum in Rioboo?
> >
> > Possibly. Or do what we do now, that is use symbolic roots.
> > For degree 3 after you add symbolic root polynomial splits
> > into linear term and quadratic (or into linear terms).
> >
> > rootSum is the future, but as long as limit can not handle
> > rootSum we have to use symbolic roots.
> >
>
> If current limit can handle result from integrate, then the limit
> can expand the rootSum just like integrate does. Not an improvement,
> but also not a step back, right?

Well, we want to handle more cases. But there is also problem
of incorrect results. Old limit code has a few fundamental
problems, basically parts of it work on principle "this
usually produces correct result". But sometimes we get
wrong results. Handling of symbolic roots is in this category.
In MRVLIM I tried to avoid such things. There is still possiblity
that MRVLIM gets something which breaks assumptions from
other parts of FriCAS (or from the user), but I do not want
to deliberatly include unsound code in MRVLIM. To put it
differently, if we allow the same liberty for new code as
for old code, then we will never get rid of unsound parts.

You may ask why unsound parts were included at all? Well,
one thing is research: in research you concentrate on
one part and take other as given, that is the only sane
way of doing research. But when it comes to implementation
it may happen that part treated "as given" is in fact
unavailable. Then you need to substitute something
which works well enough to test the new part. There
is also competitive pressure: if competitor uses
unsound method there is strong pressure to "do something"
which frequently boils down to using the same unsound
method. But IMO there is no compeling reason _now_
to include unsound methods in MRVLIM.

--
Waldek Hebisch

[Qian Yun]

Can you give a concrete example that limit fails to do
rootSum? Also, is limit on rootSum generally solvable?

- Qian

[Waldek Hebisch]

ATM limit just refuses to do rootSum:

(1) -> ff2 := rootSum(l*log(x - l), univariate(x^2 + 1)::SUP(EXPR(INT)), 'l)

--+
(1) > l log(x - l)
--+
2
l + 1 = 0
Type: Expression(Integer)
(2) -> limit(ff2, x=%minusInfinity)

(2) "failed"
Type: Union("failed",...)

In my developement vesrion I pass rootSum to to MRVLIM and I get:

(12) -> limit(ff2, x=%minusInfinity)

(12) 0
Type: Union(OrderedCompletion(Expression(Integer)),...)
(13) -> limit(ff2, x=%plusInfinity)

(13) 0
Type: Union(OrderedCompletion(Expression(Integer)),...)
(14) -> D(ff2, x)

2
(14) - ------
2
x + 1
Type: Expression(Integer)

As you can see the result from limit is wrong.

In general, rootSum leads to tricky problems because theory underlying
MRVLIM is real, while complex roots introduce complex functions.
But cases like above should be doable with moderate effort. Namely,

log(x - c) = log(x)*log(1 - c/x)

the second term goes to 1 when x goes to infinity, and we can
write asymptitic expansion of log(x - c). For single term we
get complex coefficients, but rootSum has property that coeficients
of expansion are again rootSum-s and give real values. So
given

\sum_{p(c) = 0} c*log(x - c)

we can compute its asymptotic expansion when x goes to infinity.
In case of convergent integral real parts will give 0. So value
is decided by imaginary parts, that is argument of logs, but this
also goes to 0. A bit more tricky is case when x goes to -infinity.
We can write this as log(-x - c). As first approximation, we
rewrite it as log(-1) + log(x + c). The log(x + c) parts can be
handled as before. The log(-1) part requires extra care, it is
%i*%pi when imaginary part of c is positive and -%i*%pi when imaginary
part of c is negative. So we need to compute

\sum_{Im(c) > 0} c - sum_{Im(c) < 0} c

Multiplied by %i this is real algebraic number and in principle
we should be able to compute its minimal polynomial. If p
has degree 2*k, than degree of minimal polynomial may be
as high as binomial(2*k, k)/2. In first relevant case, that is k=2,
degree will be 3 or less, so quite managable. When degree is 3,
there is one real root and two complex ones. So we get reasonably
satisfactory answer.

When x has finite limit we can simply return rootSum with x replaced
by limiting value. Simplifying such rootSum may be quite tricky,
but is valid result.

In more general case we may have some rational function of c before
log, the computation should go on with little changes.

The above is most general integral of rational function over
real numbers, that is when we split polynomials into linear
factors over complex numbers. In practice we get logs of
polynomials of degree bigger than 1. Splitting into linear
factors theoreticaly should work, but is quite heavy. ATM I have
no significantly better idea.

--
Waldek Hebisch

[Waldek Hebisch]

On Wed, Nov 22, 2023 at 04:35:30AM +0100, Waldek Hebisch wrote:
> On Wed, Nov 22, 2023 at 08:20:30AM +0800, Qian Yun wrote:
> >
> > Can you give a concrete example that limit fails to do
> > rootSum? Also, is limit on rootSum generally solvable?

A little correction and addition to previous post.

> But cases like above should be doable with moderate effort. Namely,
>
> log(x - c) = log(x)*log(1 - c/x)
>
> the second term goes to 1 when x goes to infinity, and we can
> write asymptitic expansion of log(x - c). For single term we
> get complex coefficients, but rootSum has property that coeficients
> of expansion are again rootSum-s and give real values. So
> given
>
> \sum_{p(c) = 0} c*log(x - c)
>
> we can compute its asymptotic expansion when x goes to infinity.
> In case of convergent integral real parts will give 0. So value
> is decided by imaginary parts, that is argument of logs, but this
> also goes to 0. A bit more tricky is case when x goes to -infinity.
> We can write this as log(-x - c). As first approximation, we
> rewrite it as log(-1) + log(x + c). The log(x + c) parts can be
> handled as before. The log(-1) part requires extra care, it is
> %i*%pi when imaginary part of c is positive and -%i*%pi when imaginary
> part of c is negative. So we need to compute
>
> \sum_{Im(c) > 0} c - sum_{Im(c) < 0} c
>
> Multiplied by %i this is real algebraic number and in principle
> we should be able to compute its minimal polynomial. If p
> has degree 2*k, than degree of minimal polynomial may be
> as high as binomial(2*k, k)/2. In first relevant case, that is k=2,

Should be binomial(2*k, k)

> degree will be 3 or less, so quite managable. When degree is 3,

^ 6

> there is one real root and two complex ones.

2 real roots and 4 complex ones. Actually, in this case when
r is a root also -r is a root. We need to choose positive
real root.

> In more general case we may have some rational function of c before
> log, the computation should go on with little changes.

ATM most troubling part is how to choose correct root. I expect
two real roots but it is not clear which one is the roght one.

> The above is most general integral of rational function over
> real numbers, that is when we split polynomials into linear
> factors over complex numbers. In practice we get logs of
> polynomials of degree bigger than 1. Splitting into linear
> factors theoreticaly should work, but is quite heavy. ATM I have
> no significantly better idea.

AFAICS splitting is not needed. Just knowing degree d we know
that split will have form

P(x) = \prod((x - c[i]), i = 1..d)

which as before gives 0 in the limit. For limit when x goes to
minus infinity we may observe that it is minus the corresponding
definite integral over the real line. So we can differentiate
to recover the integrand and use residue method to compute
integral over the real line. AFAICS the residue method is
equivalent to computing sum that we considered above, so looks
doable.

There is information loss, for P(x) of even degree trying to
compute limit when x goes to minus infinity by computing
limit of expression built from -x when x goes to plus infinity
leads to wrong result. So it seems that we need special handling
in definite integration code. Also, for limit at finite point
we get rootSum which may be tricky to correctly evaluate
numerically: we need to use correct brach of log and once
polynomial has degree bigger than 1 principal branch may
be wrong.

--
Waldek Hebisch

[Qian Yun]

Here are my thoughts:

Couldn't we utilize the same trick used in section 2.8
"Riboo's Algorithm for Real Rational Functions"?

From my observations, most (if not all) rootSum has the
following structure:

rootSum(p(a), a*log(q(a)+f(x))), where p,q are polynomials
and f is a general function.

{
rootSum(p(a), a*log(q(a)*g(x)+f(x))) can be written as
rootSum(p(a), a*log(q(a)+f(x)/g(x))+a*log(g(x)))
}

So if p, q have real (this condition can be relaxed?) coefficients, then

rootSum(p(a), a*log(q(a)+f(x))) =
rootSum(p(a) | real? a, a*log(q(a)+f(x))) +
rootSum((u,v), 1/2*u*log((A(u,v) + f(x))^2 + B(u,v)^2)) +
rootSum((u,v), 1/2*v*atan((A(u,v) + f(x))/B(u,v)))

{where u+%i*v is complex root of p(a). q(u+%i*v)= A(u,v) + %i*B(u,v)}

Not that the summation is over all of v, not limited to positive v.

So, then take the limit on x:

The atan is finite, so whether the rootSum is finite depends on the
rootSum part.

If f(x) goes to +infinity, then the rootSum is going to infinity too,
the sign is determined by

rootSum(p(a), a * log(f(x))), which is coefficient(p(a), degree(p) - 1),
which is equal to 0 most of the time?

If f(x) goes to -infinity, well, isn't integrate(1/x,x) is log(abs(x))?
If the rootSum comes from real integral, then taking logs on negative
value could be avoided by abs.

- Qian

[Waldek Hebisch]

On Thu, Nov 23, 2023 at 08:35:02AM +0800, Qian Yun wrote:
> Here are my thoughts:
>
> Couldn't we utilize the same trick used in section 2.8
> "Riboo's Algorithm for Real Rational Functions"?

When it works, yes.

> From my observations, most (if not all) rootSum has the
> following structure:
>
> rootSum(p(a), a*log(q(a)+f(x))), where p,q are polynomials
> and f is a general function.
>
> {
> rootSum(p(a), a*log(q(a)*g(x)+f(x))) can be written as
> rootSum(p(a), a*log(q(a)+f(x)/g(x))+a*log(g(x)))
> }
>
> So if p, q have real (this condition can be relaxed?) coefficients, then
>
> rootSum(p(a), a*log(q(a)+f(x))) =
> rootSum(p(a) | real? a, a*log(q(a)+f(x))) +
> rootSum((u,v), 1/2*u*log((A(u,v) + f(x))^2 + B(u,v)^2)) +
> rootSum((u,v), 1/2*v*atan((A(u,v) + f(x))/B(u,v)))
>
> {where u+%i*v is complex root of p(a). q(u+%i*v)= A(u,v) + %i*B(u,v)}
>
> Not that the summation is over all of v, not limited to positive v.

Trouble here is with A(u,v) and B(u,v). Already when p is of degree
4 with Galois group S(4) we have trouble: beside right solutions
there are spurious ones and it is hard to separate good ones from
spurious one.

Consider

p := 10*(a^4 + 1) + a

and

rootSum(a*log(x - a), p::SUP(EXPR(INT)))

How you want to find right u, v? The approach I described
works around difficulty with u, v by working with coefficients
of p. We still have difficulty with final result, but it is
less problematic: we need to choose positive real root of
a polyniomial. Note that the polynomial is minimal polynomial
for the answer, so we need to choose this root directly or
indirectly. But we can avoid trouble with u and v.

--
Waldek Hebisch

[Qian Yun]

In this particular case, B(u, v) = - v, when v is positive, B is
negative. So for this case, limit when x is +infinity is
rootSum(-%pi*v, V(v) where v is real and positive)

V(v) can be achieved by resultant:

p := 10*(a^4 + 1) + a

puv := eval(p,a=u+%i*v)
P := real puv
Q := imag puv
V := resultant(univariate(P, second kernels P), univariate(Q, second
kernels P))::POLY INT

40960000000*v^16+20480000000*v^12+(-166400000)*v^10+(-17920000000)*v^8+(-86400000)*v^6+2559973000*v^4

countRealRoots V returns 5, which is a 0 (ignore) and 2 positive v and
2 negative v (come in pairs), which is the correct number of v.

This does not scale for more complex B(u, v) though.

- Qian

[Waldek Hebisch]

The trouble is that AFAICS we have _no_ reasonable expression for
real roots. When V factors we can try to use countRealRoots to
find factors responsible for real roots. However, above V has
irreducible factor of degree 12 and need to sum over its real roots.

If wee look at u we get polynomial of degree 6, which may be easier
to handle (but then we need square roots to get v). In approach
that I sketched we get polynomial of degree 6 for final result.

BTW: Q is always divisible by v and v=0 can not give complex root,
so we can divide by it in "interesting" cases (that is when there
is no real root).

--
Waldek Hebisch

[Qian Yun]

ZeroDimensionalSolvePackage has "realSolve" and "positiveSolve"
which gives RealClosure as result, and the numeric value is right
(I computed another one and compared with wolframalpha.)

- Qian

[Waldek Hebisch]

On Thu, Nov 23, 2023 at 07:15:55PM +0800, Qian Yun wrote:
> ZeroDimensionalSolvePackage has "realSolve" and "positiveSolve"
> which gives RealClosure as result, and the numeric value is right
> (I computed another one and compared with wolframalpha.)

Well, ZeroDimensionalSolvePackage and CAD may help. But AFAICS
at least ATM there is no way to convert from RealClosure back
to expression.

--
Waldek Hebisch

[Waldek Hebisch]

On Tue, Nov 21, 2023 at 06:31:27PM +0800, Qian Yun wrote:
>

> How do you suggest to proceed? Write a special case for biquadratic?

AFAICS we can quickly add case for x^4 + a, that when we can
determine sign of a. Case for biquadratic is more complicated
but should be doable with moderate effort. Both should be
rather efficient as they do not need more general factoring
nor complicated solve.

We could try to handle more cases with easy Galois group
by factoring resolvent (resultant of P and Q). This is
more complicated, but still moderately hard.

I think that handling more general cases requires more
research and will take some time. So it makes sense to
do easier cases first.

--
Waldek Hebisch

[Qian Yun]

I have examined integrals where rootSum is over degree 4, from a
subset of Nasser's test.

Other less common cases are:

%A^8+1 integrate(log(1/x^4+x^4),x)
%A^4-1/2048 integrate(1/(x^4-2),x)
%A^6-1/1492992 integrate(1/(x^6-2),x)
%A^6+1/1492992 integrate(1/(x^6+2),x)
%A^5-7776/3125 integrate(x^(1/2)/(-1/x^(1/3)+x^(1/2)),x)
(factors to %A^4+6/5*%A^3+36/25*%A^2+216/125*%A+1296/625)
%A^6+%A^3+1 integrate(((-3)*x^3+6)/(x^6+(-1)*x^3+1), x)

Of course, there are corresponding cases where coefficients are
POLY INT.

I'll do more testing to see if the generic resultant method can
handle them or special cases are required.

- Qian

[Waldek Hebisch]

On Tue, Nov 28, 2023 at 07:02:16PM +0800, Qian Yun wrote:
> I have examined integrals where rootSum is over degree 4, from a
> subset of Nasser's test.
>
> Other less common cases are:

<snip>
> %A^4-1/2048 integrate(1/(x^4-2),x)

%A^4-1 is actually an example of x^4 + a with negative a, so
easy to handle once we extend scalars by sqrt(2048) (equivalently
sqrt(2))

> %A^6-1/1492992 integrate(1/(x^6-2),x)

similar to the above, using sqrt(1492992)

> %A^6+1/1492992 integrate(1/(x^6+2),x)

This leads to degree 3 equation for real parts of the roots.

--
Waldek Hebisch

[Waldek Hebisch]

On Tue, Nov 28, 2023 at 07:02:16PM +0800, Qian Yun wrote:

> I have examined integrals where rootSum is over degree 4, from a
> subset of Nasser's test.
>
> Other less common cases are:
>
> %A^8+1 integrate(log(1/x^4+x^4),x)

Here equation for real part of roots (obtained via resolvent) is of
degree 4 and we can solve it using real square roots:

(18) -> fl(2).factor

4 2
(18) 8 u - 8 u + 1
Type: Polynomial(Integer)
(19) -> radicalSolve(fl(2).factor)

(19)
+----------+ +----------+ +--------+ +--------+
| +-+ | +-+ | +-+ | +-+
\|- \|2 + 2 \|- \|2 + 2 \|\|2 + 2 \|\|2 + 2
[u = -------------, u = - -------------, u = -----------, u = - -----------]
2 2 2 2
Type: List(Equation(Expression(Integer)))


--
Waldek Hebisch

[Waldek Hebisch]

On Tue, Nov 28, 2023 at 07:02:16PM +0800, Qian Yun wrote:

> I have examined integrals where rootSum is over degree 4, from a
> subset of Nasser's test.

There are also examples of degree 3. When there is 1 real root
solver is doing reasonable job and real roots can be recognized,
like

integrate(((-3)*x^2+(-6)*x+6)/(x^3+3*x^2+3*x+(-2)), x)
-- a^3+3*a^2+3*a-2

or

integrate(((-10)*x+(-3))/(x^3+5*x+1), x)
-- a^3 + 5*a + 1

The first one is actually easy one, shifted root of degree 3.
The second one almost works, but the result is bigger than result
with implict root. And differentiating the integral gives
expression with unsimplified roots.

Actually there is general trouble with using roots of polynomial
of degree 3: splitting field F of polynomial (that is field
generated by all its roots) in normal, that is any irreducible
polynomial which has a root in F splits into linear factors in F.
In particular, if x^3 - a is irreducible and one of the roots
belongs to F, then also other roots belong to F. So if Cardano
fomula works inside F, then either root is reducible or roots
of 1 of degree 3 belong to F. Reducible root means that corresponding
extention is of degree lower than 3, so 1 or 2. As other
roots in Cardano formula are squere roots the polynomial would
have solution in extention of degree which is power of 2. But root
of irreducible polynomial of degree 3 generates extention of degree 3
and consequently any field containing its root has degree divisible
by 3. So reducible root is posible only when polynomial is
reducible. Now, in case of real roots it is clear that
F does not contain primitive root of degree 3 which is complex.
But actually, even if two roots are complex normally F does
not contain primitive root of degree 3. So Cardano formuals
force work in field of degree 12. In case of a^3 + 5*a + 1
this means field generated by sqrt(527/3), sqrt(-3) and
root of degree 3. Our routines rewrite sqrt(527/3) as
sqrt(527)/sqrt(3), so we get extra spurious roots.

Integration process introduces extra spurious roots, sqrt(1581)
should be rewriten as 3*sqrt(527/3), other root apparently
generate extra extentions which however are unnecessary.
So, to get simple result in this case we need special handling.
I think that we should look at discriminant. When discriminant
is negative we can try Rioobo for two complex roots, but we
need to avoid extra roots introduced by general solver.

--
Waldek Hebisch

[Waldek Hebisch]

On Tue, Nov 28, 2023 at 07:02:16PM +0800, Qian Yun wrote:

> I have examined integrals where rootSum is over degree 4, from a
> subset of Nasser's test.

Martin Weltz (clicliclic) points out to example from Timofeev:

integrate(tan(x)/(sqrt(tan(x)) - 1)^2, x)

where root sum is over roots of

a^5 - (1/2)*a^3 - (5/8)*a^2 + (9/64)*a - 1/64

which factors as

(54) -> factor(a^5 - (1/2)*a^3 - (5/8)*a^2 + (9/64)*a - 1/64)

4 3 1 2 1 1
(54) (a - 1)(a + a + - a - - a + --)
2 8 64
Type: Factored(Polynomial(Fraction(Integer)))

The quartic gives complicated roots when passed to radicalSolve.
'rsimp2.input' which I posted some time ago when applied kernel-by-kernel
can express roots in terms of sqrt(-1) and sqrt(2), but ATM it is
not clear to me if this is good method. I would prefer to get
simple answer in more direct way.

Equation for real parts of the roots splits into 3 factors:

[[flag = "prime", factor = 8 u + 4 u + 1, exponent = 2],
2
[flag = "prime", factor = 16 u + 8 u - 1, exponent = 2],
2
[flag = "prime", factor = 16 u + 8 u + 3, exponent = 2]]

and 'countRealRoots' finds out that the second one is the
correct factor. This leads to simple expressions for real
parts:

(28) -> radicalSolve(fl(2).factor)

+-+ +-+
- \|2 - 1 \|2 - 1
(28) [u = ----------, u = --------]
4 4
Type: List(Equation(Expression(Integer)))


and complex parts. So here resolvent works nicely.

--
Waldek Hebisch

[Qian Yun]

On 12/1/23 05:51, Waldek Hebisch wrote:
> On Tue, Nov 28, 2023 at 07:02:16PM +0800, Qian Yun wrote:
>> I have examined integrals where rootSum is over degree 4, from a
>> subset of Nasser's test.
>
> Martin Weltz (clicliclic) points out to example from Timofeev:
>
> integrate(tan(x)/(sqrt(tan(x)) - 1)^2, x)
>
> where root sum is over roots of
>
> a^5 - (1/2)*a^3 - (5/8)*a^2 + (9/64)*a - 1/64

Is this integral example wrong? "internalIntegrate" shows a different
rootSum.

So it looks like that this "Riboo real functions" method is effective
for rootSum over polynomial with FRAC(INT) coefficients.

As for POLY(INT) coefficients, 'countRealRoots' fails and we have to
do case by case treatment, and for some cases, we have to assume
different signs (branches), giving a list of possible forms as result.

If you are working on the code, please continue, so our efforts will
not be duplicated.

- Qian

[Waldek Hebisch]

On Fri, Dec 01, 2023 at 07:55:45AM +0800, Qian Yun wrote:
>
>
> On 12/1/23 05:51, Waldek Hebisch wrote:
> > On Tue, Nov 28, 2023 at 07:02:16PM +0800, Qian Yun wrote:
> > > I have examined integrals where rootSum is over degree 4, from a
> > > subset of Nasser's test.
> >
> > Martin Weltz (clicliclic) points out to example from Timofeev:
> >
> > integrate(tan(x)/(sqrt(tan(x)) - 1)^2, x)
> >
> > where root sum is over roots of
> >
> > a^5 - (1/2)*a^3 - (5/8)*a^2 + (9/64)*a - 1/64
>
> Is this integral example wrong? "internalIntegrate" shows a different
> rootSum.

This is essentially what ')trace INTEF )math' shows me. The
only difference is that I use 'a' as variable name istead of
generated symbol name.

> So it looks like that this "Riboo real functions" method is effective
> for rootSum over polynomial with FRAC(INT) coefficients.
>
> As for POLY(INT) coefficients, 'countRealRoots' fails and we have to
> do case by case treatment, and for some cases, we have to assume
> different signs (branches), giving a list of possible forms as result.

Yes, parameters in coefficients make things more tricky. We can
still do something if our sign deternination works. For example
'sign' treats 'a^2' as nonnegative and 'exp(a)' as positive.

> If you are working on the code, please continue, so our efforts will
> not be duplicated.

I have put almost all that I did into messages. As I wrote,
I think we should implement a few easy cases:

- degree 3 with one real root
- a^4 + c with c of know sign
- factorizable resolvent for degree 4
- maybe few other

AFAICS a^4 + c is most important as it appears in many examples.
I would limit this to cases when 'sign' gives us answer or
when split between cases is simple and easy.

I am not sure if you want to do them. I could do first two and
third for rational coefficients. Or maybe you want more general
code.


--
Waldek Hebisch

[Qian Yun]

On 12/1/23 08:58, Waldek Hebisch wrote:
>
> I have put almost all that I did into messages. As I wrote,
> I think we should implement a few easy cases:
>
> - degree 3 with one real root
> - a^4 + c with c of know sign
> - factorizable resolvent for degree 4
> - maybe few other
>
> AFAICS a^4 + c is most important as it appears in many examples.
> I would limit this to cases when 'sign' gives us answer or
> when split between cases is simple and easy.
>
> I am not sure if you want to do them. I could do first two and
> third for rational coefficients. Or maybe you want more general
> code.
>
>

You should do them, and if possible, commit separately for
different cases -- it's easier for me to testing against Nasser's
list.

- Qian

[Waldek Hebisch]

I prefer commit code when reasonably complete and ready.
If you want to test attached is my first shot at x^4 + a case
with negative a. It seems to be rare, but it was easy and
will share most code with case of positive a.

Example:

(1) -> integrate(1/(x^4-2),x)

4+-+
4+-+ 4+-+ x\|8
- log(x\|8 + 2) + log(x\|8 - 2) - 2 atan(-----)
2
(1) -------------------------------------------------
4+-+
4 \|8
Type: Union(Expression(Integer),...)

(before there were two complex logs).

--
Waldek Hebisch

[Qian Yun]

Thanks. I would not have thought of using "quadeval".
It is much more suitable than my "eval with new symbols and sqrt(-1)".

A minor comment: the right parenthesis at the end the second line
should move to the end of first line.

- Qian

[Waldek Hebisch]

Attached is at x^4 + a case handling both negative and positive a.
I added separate function to compute root of degree 4, but this
still needs improvement. Positive a catches several examples
in mapleok.input.

Example:

(1) -> integrate(1/(x^4 + a^4), x)

(1)
+-+
+-+ 2 2 +-+ 2 2 x\|2 + a
log(a x\|2 + x + a ) - log(- a x\|2 + x + a ) + 2 atan(---------)
a
+
+-+
x\|2 - a
2 atan(---------)
a
/
3 +-+
4 a \|2
Type: Union(Expression(Integer),...)

Unfortunately,

integrate(1/(x^4 + a^2), x)

produces (4*a^2)^(1/4) and consequently answer is more complicated
then it should be.

--
Waldek Hebisch

[Qian Yun]

Should "froot(4*a^2,4)" return
[exponent = 2, coef = 1, radicand = 2 a] instead of
[exponent = 4, coef = 1, radicand = 4 a^2] ?

- Qian

[Waldek Hebisch]

On Sun, Dec 03, 2023 at 09:24:51PM +0800, Qian Yun wrote:
> Should "froot(4*a^2,4)" return
> [exponent = 2, coef = 1, radicand = 2 a] instead of
> [exponent = 4, coef = 1, radicand = 4 a^2] ?

Probably. AFAICS easy modification to froot would do this:
one needs to separately handle content of polynomials
(numerator and denominator). OTOH there is question
how far should we go. 'zroot' internally used inside
PolynomialRoots needs to factor integers. For large
integers it is very expensive, so rather undesirable.

OTOH we can do essential root simplification by factoring
polynomials. When there is large integer content,
factoring polynomial is much cheaper than factoring
integers. By factoring polynomials we will not
simplify say 'sqrt(2048)' to '32*sqrt(2)', but it
is enough to guarantee independence of roots. And
by factoring polynomials we can simplify 'sqrt(2*sqrt(2) + 3)'
to '1 + sqrt(2)' which 'froot' can not do.

So, 'froot' is fairly incomplete simplifier. Its potential
advantage is that it may be less expensive than factoring.
But then we should avoid really expensive stuff in 'froot'.

--
Waldek Hebisch

[Waldek Hebisch]

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

[Waldek Hebisch]

On Mon, Dec 04, 2023 at 01:52:39AM +0100, Waldek Hebisch wrote:
> 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 fails on integrals like

integrate((z^2+z)^(1/2)/(1+z^2)^2, z)

where the poly is a^4 - (1/64)*a^2 + 1/8192 and
resolvent is 64*u^3 + (1/2)*u^2 - (1/1024)*u which
factors as 64* u*(u^2 + (1/128)*u - 1/65536. Here linear
factor must be discarded and quadratic factor has one
positive root. We should use this positive root, but
currently the routine gives up.

--
Waldek Hebisch

[Qian Yun]

I wonder if the current work on figuring out simplified roots
of degree 4 polynomial is somewhat related to the recent
"radicalSolve" post on sci.math.symbolic?

(aka: will this also improve radicalSolve?)

- Qian

[Waldek Hebisch]

On Tue, Dec 05, 2023 at 09:06:48PM +0800, Qian Yun wrote:
> I wonder if the current work on figuring out simplified roots
> of degree 4 polynomial is somewhat related to the recent
> "radicalSolve" post on sci.math.symbolic?

Well, one of examples came from this post.

>
> (aka: will this also improve radicalSolve?)

Code that I posted only affects integration, more precisely
transformations of logarithms. I suspect that one could also
do some transformation of special functions, but ATM we do
nothing for them.

There are related issues:

- computation of Galois groups, currently not implemented in
FriCAS
- checking for dependence/independence of roots of Trager
double resultant (important when we want to claim nonintegrability)
- limits of 'rootSum'
- 'radicalSolve'
- TranscendentalSolve
- determining sign
- root simplifiction

It would be nice to share code between various uses, however there
are differences, for example 'radicalSolve' may produce complex
valued radical, for definite integration we want real quantities.

Ideally integrator will make sure that all new roots introduced
during integration are independent, for this we probably should
have "renormalization" pass, that is run computations like
normalization, and use result to eliminate dependencies between
new kernels.

--
Waldek Hebisch

[Waldek Hebisch]

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

[Qian Yun]

Looks OK. I only have one question:

"do_rr" returns a list of length 1. Does this mean it can't
handle some parameterized integral?

- Qian

[Waldek Hebisch]

On Sun, Dec 17, 2023 at 08:40:11AM +0800, Qian Yun wrote:
> Looks OK. I only have one question:
>
> "do_rr" returns a list of length 1. Does this mean it can't
> handle some parameterized integral?

Yes, the code can handle some parameterized integrals, but will
fail on majority of parameteric cases. More precisely,
'do_rr' is called only when we get sucessful result
from 'complex_roots' and 'complex_roots' returns sucessfully
only when it can determine actual configuration of roots.

In principle we could call 'do_rr' for each hypotethical
configuration of roots, but already for general degree 4
there is several possible configurations. And wrong configurations
and hard cases are likely to generate complicated expressions.
So, to keep things simple, if there is any doubt, then we fall back
to current code.

Theoretically, we could attach validity conditions to various
configurations, that is produce only configurations that are
valid on a nonempty set of parameters. But again, to keep
things simple I am not doing this. And if there are conditions
for given configuration, then it would be natural to return
list of 'suchThat' constructs. But then the thing would propagate
needing more changes.

--
Waldek Hebisch