[fricas-devel] interesting error in integration

['Martin R' via FriCAS - computer algebra system]

Nasser found the following interesting bug at
https://trac.sagemath.org/ticket/25905:

(1) -> f := (I*a*tan(d*x + c) + a)^3*tan(d*x + c)

(1)
3 3 4 2 3 3 3 2
I a tan(d x + c) + 3 I a tan(d x + c) + 3 I a tan(d x + c)
+
3
a tan(d x + c)
Type:
Expression(Integer)
(2) -> F := integrate(f, x)

(2)
2 3 1 3 3 3
(9 I - 3)a log(-----------------) + 2 I a tan(d x + c)
2
tan(d x + c) + 1
+
2 3 2 3 3 3 3
9 I a tan(d x + c) + (- 6 I + 18 I)a tan(d x + c) + (6 I - 18 I)a
d x
/
6 d
Type:
Union(Expression(Integer),...)
(3) -> D(eval(F, I =%i), x)-eval(f, I=%i)

(3) 0
Type:
Expression(Complex(Integer))
(4) -> F := integrate(eval(f, I=%i), x)
arithmetic error DIVISION-BY-ZERO signalled
Operation was (/ -2 0).

Any ideas?

Martin

[Waldek Hebisch]

Post by 'Martin R' via FriCAS - computer algebra system
Nasser found the following interesting bug at
(1) -> f := (I*a*tan(d*x + c) + a)^3*tan(d*x + c)
(1)
3 3 4 2 3 3 3 2
I a tan(d x + c) + 3 I a tan(d x + c) + 3 I a tan(d x + c)
+
3
a tan(d x + c)
Expression(Integer)
(2) -> F := integrate(f, x)
(2)
2 3 1 3 3 3
(9 I - 3)a log(-----------------) + 2 I a tan(d x + c)
2
tan(d x + c) + 1
+
2 3 2 3 3 3 3
9 I a tan(d x + c) + (- 6 I + 18 I)a tan(d x + c) + (6 I - 18 I)a
d x
/
6 d
Union(Expression(Integer),...)
(3) -> D(eval(F, I =%i), x)-eval(f, I=%i)
(3) 0
Expression(Complex(Integer))
(4) -> F := integrate(eval(f, I=%i), x)
arithmetic error DIVISION-BY-ZERO signalled
Operation was (/ -2 0).
Any ideas?

This is old confusion between real and complex. We create
Complex(Complex(Integer)) and loose. It would be good to
clean this up. AFAICS proper fix should not require much
code, but it is going to touch several domains and packages.
Basically, old code assumed that functions are real and
allowed complex functions only as a hack. Current
'FunctionSpaceComplexIntegration' integrates _real_
function by passage to complex variables. We need a
function to integrate _complex_ functions.

[Waldek Hebisch]

Post by Waldek Hebisch

Post by 'Martin R' via FriCAS - computer algebra system
Nasser found the following interesting bug at
(1) -> f := (I*a*tan(d*x + c) + a)^3*tan(d*x + c)
(1)
3 3 4 2 3 3 3 2
I a tan(d x + c) + 3 I a tan(d x + c) + 3 I a tan(d x + c)
+
3
a tan(d x + c)
Expression(Integer)
(2) -> F := integrate(f, x)
(2)
2 3 1 3 3 3
(9 I - 3)a log(-----------------) + 2 I a tan(d x + c)
2
tan(d x + c) + 1
+
2 3 2 3 3 3 3
9 I a tan(d x + c) + (- 6 I + 18 I)a tan(d x + c) + (6 I - 18 I)a
d x
/
6 d
Union(Expression(Integer),...)
(3) -> D(eval(F, I =%i), x)-eval(f, I=%i)
(3) 0
Expression(Complex(Integer))
(4) -> F := integrate(eval(f, I=%i), x)
arithmetic error DIVISION-BY-ZERO signalled
Operation was (/ -2 0).
Any ideas?

This is old confusion between real and complex. We create
Complex(Complex(Integer)) and loose.

Thinking loudly. In order to integrate trigonometric
functions we need imaginary unit. So, if not present we
need to add it, that is extend the original field. Currently
this is done by changing base ring from R to Complex(R).
OTOH if imaginary unit is present we should use it.

There is extra complication, namely normalization works differently
when imaginary unit is present.

Now, adding imaginary unit when one is already present
is an error: it causes zero divisors which may lead to
wrong answer or crash. Also, algebraic quantities
complicate computations so for efficiency reasons we
would like to avoid adding imaginary unit when it
is not necessary.

One solution could be a wrapper which test for presence
of imaginary unit and produces to new domains: "complex_R"
and "Expression(complex_R)". "complex_R" should
be the same as R when R has imaginary unit, otherwise
"complex_R" should be "Complex(R)". We need convertion
between those domains, which currently use "Complex(Expression(R))",
so we need extra domain to use in itermediate stage.

Main integration routine could take five domains above
and imaginary unit in "complex_R" as parameters and
depending on presence (or absence) of trigonometric
functions convert to "Expression(complex_R)".

[Ralf Hemmecke]

Post by Waldek Hebisch
Thinking loudly. In order to integrate trigonometric
functions we need imaginary unit. So, if not present we
need to add it, that is extend the original field. Currently
this is done by changing base ring from R to Complex(R).
OTOH if imaginary unit is present we should use it.

Hmmm... also thinking loudly...
We cannot and should not have Complex Complex R, but what about a domain

ComplexClosure(R: CommutativeRing): ComplexCategory(R) with ...
== if R has with imaginary: () -> R then R else Complex R

Would this be the wrapper that you want? To me that looks like an
elegant solution.

Of course there can be rings R that don't export "imaginary" and still
have an element x that satisfies x^2+1=0. But I didn't want a condition
that checks whether x^2+1 can be factored over R.

As a first approximation, ComplexClosure would do what is needed, no?
Adding conversion between R and ComplexClosure R shouldn't be difficult.

Ralf

[Riccardo GUIDA]

Post by Waldek Hebisch
This is old confusion between real and complex. We create
Complex(Complex(Integer)) and loose.

Naive question: is it not possible/enough to define Complex(R) to be R if R = Complex (S)?

[Waldek Hebisch]

Post by Riccardo GUIDA

Post by Waldek Hebisch
This is old confusion between real and complex. We create
Complex(Complex(Integer)) and loose.

Naive question: is it not possible/enough to define Complex(R) to be R if R = Complex (S)?

No. There are axioms which link 'real', 'imag' and 'complex' that
essentially force current implementation. Trying to cheat
(ignering axioms) is likely to lead to problems in other places.

[Waldek Hebisch]

Post by Ralf Hemmecke

Post by Waldek Hebisch
Thinking loudly. In order to integrate trigonometric
functions we need imaginary unit. So, if not present we
need to add it, that is extend the original field. Currently
this is done by changing base ring from R to Complex(R).
OTOH if imaginary unit is present we should use it.

Hmmm... also thinking loudly...
We cannot and should not have Complex Complex R, but what about a domain
ComplexClosure(R: CommutativeRing): ComplexCategory(R) with ...
== if R has with imaginary: () -> R then R else Complex R
Would this be the wrapper that you want? To me that looks like an
elegant solution.
Of course there can be rings R that don't export "imaginary" and still
have an element x that satisfies x^2+1=0. But I didn't want a condition
that checks whether x^2+1 can be factored over R.
As a first approximation, ComplexClosure would do what is needed, no?
Adding conversion between R and ComplexClosure R shouldn't be difficult.

Well, ATM I am thinking of doig this, but witout introducing such
global constructor. But this is just small part of the whole problem.

We need convertions between Expression(R) and Expression(Complex(R))
which while not very hard is more complicated than convertion
between R and Complex(R). In particular imaginary unit may appear
in final result so convertion from Expression(Complex(R)) to
Expression(R) must properly handle imaginary unit (replace %i by
sqrt(-1) if R is real). Extra factor is that convertions are
supposed to be isomorphisms of differential fields so we need
to remember enough information to be able to compute the inverse.
More precisely, we are embedding smaller field in a bigger one
and getting representation in smaller field requires some
effort.

ATM all involved convertions are split between several steps.
There is some possibility to rationalize things, for example
we first convert trigonometric functions to tangents and
later to complex eponentials. It would be more natural
to do this in one step. But limit code needs real functions
(that is tangents) and currently normalization is shared
between integrator and limit code.

We could probably get simpler code by adding imaginary unit
before anything else. OTOH unneeded imaginary unit is an
extra algebraic extention which may increase cost of other
operations. In the past it could lead to other problems,
but now I hope that the only trouble will be loss of
efficiency. But current code makes some effort to introduce
imaginary unit only when needed and it would be nice
to preserve this property.

[Ralf Hemmecke]

Post by Waldek Hebisch
in final result so convertion from Expression(Complex(R)) to
Expression(R) must properly handle imaginary unit (replace %i by
sqrt(-1) if R is real).

Oh... maybe I'm wrong, but AFAIU, the Expression domain is somewhat
tailored towards a nice representation for integration. It is, in fact,
something like

K(t1, t2, ..., tn) (*)

where the ti may have relations with the tj for j<i. Unfortunately, our
Expression domain, does not explicitly show this relations. Waldek, you
know certainly better than me that the integration code builds this
tower of fields somehow, but in fact it's always in Expression(R) or
maybe Expression(Complex(R)). Am I wrong?

It's probably too much of rewriting to actually get rid of Expression(R)
in the integration code and let it work with a tower of fields that is
explicitly constructed for the expression that is to be integrated.
What I mean is something like this

given integrand -->
build tower (*) and map integrand to that field -->
apply Risch's algorithm -->
transform the result back to Expression(X)

I.e. during Risch's algorith there would be no involvement of
Expression(X), but only some rational function fields. Isn't that all
what Risch's algorithm cares about?

Ralf

[Kurt Pagani]

Post by Waldek Hebisch

Post by Waldek Hebisch
Thinking loudly. In order to integrate trigonometric
functions we need imaginary unit. So, if not present we

Well, ATM I am thinking of doig this, but witout introducing such
global constructor.
But this is just small part of the whole problem.

Indeed! It would be nice to have full conformity with the (group) of roots of
unity (makes sense for any unital ring).

(Un)fortunately - depending on the pint of view - we have with X==>EXPR
INT,XC==>EXPR COMPLEX INT


(1) -> I:=rootsOf(x::X^2+1)

(1) [%x0, - %x0]
Type: List(Expression(Integer))

(2) -> J:=rootsOf(x::XC^2+1)

(2) [%x8, - %x8]
Type: List(Expression(Complex(Integer)))


but neither test(J.?=%i::XC) -> true. It's even unclear to me how one could
"save" XC so that exp(2*%pi*%i/N) is "known" as a principal N-th root of unity?

https://en.wikipedia.org/wiki/Principal_root_of_unity
In an integral domain, every primitive n-th root of unity is also a principal
n-th root of unity. In any ring,...


It works at least in the context of primitive roots, e.g. N=5

(3) -> rootsOf(x::XC^5+1)

(3)
2 3 3 2
[%x0, %x0 %x1, %x0 %x1 , %x0 %x1 , - %x0 %x1 - %x0 %x1 - %x0 %x1 - %x0]
Type: List(Expression(Complex(Integer)))

but simply adding an "imaginary" might lead to more confusion.

[Waldek Hebisch]

Post by Ralf Hemmecke

Post by Waldek Hebisch
in final result so convertion from Expression(Complex(R)) to
Expression(R) must properly handle imaginary unit (replace %i by
sqrt(-1) if R is real).

Oh... maybe I'm wrong, but AFAIU, the Expression domain is somewhat
tailored towards a nice representation for integration. It is, in fact,
something like
K(t1, t2, ..., tn) (*)
where the ti may have relations with the tj for j<i. Unfortunately, our
Expression domain, does not explicitly show this relations. Waldek, you
know certainly better than me that the integration code builds this
tower of fields somehow, but in fact it's always in Expression(R) or
maybe Expression(Complex(R)). Am I wrong?

Well, we need to look at kernels, so in some places we require domain
of category FunctionSpace(R). ATM HyperDoc says me that there 3 domains
of that category: Expression, JetBundleExpression, JetBundleXExpression.
I do not know if jet bundle stuff would work for integration, but
formally it seem to have needed properties.

Post by Ralf Hemmecke
It's probably too much of rewriting to actually get rid of Expression(R)
in the integration code and let it work with a tower of fields that is
explicitly constructed for the expression that is to be integrated.
What I mean is something like this
given integrand -->
build tower (*) and map integrand to that field -->
apply Risch's algorithm -->
transform the result back to Expression(X)
I.e. during Risch's algorith there would be no involvement of
Expression(X), but only some rational function fields. Isn't that all
what Risch's algorithm cares about?

I looked shortly at this and I not sure if this is good idea. Parts
of integration code are quite general, for example 'TranscendentalIntegration'
need just a field and few special operations. Algebraic integration
builds its own types which are in some sense independent of Expression.
So actually relatively small part of integration code depends on
special properties of Expression. OTOH in part depending on Expression
(dealing with kernels) rigid structure of type tower would be
inconvenient. In particular, we handle some cases by re-arranging
towers. At several places we perform various ad-hoc extentions.
And we need to see structure of the whole field.

Maybe two examples will explain this better. As I wrote, algebraic
integration uses its own types. The effect is that to use those
routines from other parts of integrator one needs several lines
of setup code (to create types and perform needed convertions).
This makes reuse less convenient and potentially inefficient.
PolynomialFactorizationExplicit is potentially much worse than
old multivariate factorizer, because choice of good evaluation
points is really a global problem, one should choose best of
few evaluation points (with number depending on size of polynomial).
But in recursive structure like our implementation of
PolynomialFactorizationExplicit no sigle layer has all needed
information. Worse, the "best" type for intermediate results
can be determined only after call requesting given result.

[Ralf Hemmecke]

Post by Waldek Hebisch

Post by Ralf Hemmecke
Oh... maybe I'm wrong, but AFAIU, the Expression domain is
somewhat tailored towards a nice representation for integration. It
is, in fact, something like
K(t1, t2, ..., tn) (*)
where the ti may have relations with the tj for j<i.

Well, we need to look at kernels, so in some places we require
domain of category FunctionSpace(R).

Wait a minute... maybe you did not understand my transformation business

Post by Waldek Hebisch

Post by Ralf Hemmecke
given integrand --> build tower (*) and map integrand to that field
--> apply Risch's algorithm --> transform the result back to
Expression(X)

I looked shortly at this and I not sure if this is good idea.

I certainly cannot compete with your insight into the whole Risch
business, but from what I understand from integration or summation
algorithms, the problem is always first to analyse the integrand/summand
an build a tower of respective field extensions. In FriCAS that should
be particularly be easy to do. Thus, in the whole integration there are
no Expression(X) elements anymore. I don't know why it makes sense to
keep a Expression(X)-element in the middle of the integration algorithm
although the only thing that is needed are the properties of that
element. So whether it is sin or cos or acrtan or whatever function, is
unimportant. They should become indeterminates (plus properties) during
the translation phase "Expression(X) --> tower of fields".

Post by Waldek Hebisch
Parts of integration code are quite general, for example
'TranscendentalIntegration' need just a field and few special
operations. Algebraic integration builds its own types which are in
some sense independent of Expression. So actually relatively small
part of integration code depends on special properties of Expression.
OTOH in part depending on Expression (dealing with kernels) rigid
structure of type tower would be inconvenient. In particular, we
handle some cases by re-arranging towers. At several places we
perform various ad-hoc extentions. And we need to see structure of
the whole field.

As I said, you know certainly better how the FriCAS integrator works,
but if I had to implement Risch, I would probably avoid the type
Expression. Maybe I can only see the problems with this
no-Expression(X)-approach, if I actually do such an implementation.

Ralf

[Waldek Hebisch]

Post by Ralf Hemmecke

Post by Waldek Hebisch

Post by Ralf Hemmecke
Oh... maybe I'm wrong, but AFAIU, the Expression domain is
somewhat tailored towards a nice representation for integration. It
is, in fact, something like
K(t1, t2, ..., tn) (*)
where the ti may have relations with the tj for j<i.

Well, we need to look at kernels, so in some places we require
domain of category FunctionSpace(R).

Wait a minute... maybe you did not understand my transformation business

Post by Waldek Hebisch

Post by Ralf Hemmecke
given integrand --> build tower (*) and map integrand to that field
--> apply Risch's algorithm --> transform the result back to
Expression(X)

I looked shortly at this and I not sure if this is good idea.

I certainly cannot compete with your insight into the whole Risch
business, but from what I understand from integration or summation
algorithms, the problem is always first to analyse the integrand/summand
an build a tower of respective field extensions. In FriCAS that should
be particularly be easy to do.

Well, tell me what type hipotethical 'to_tower' should have?
At this point I gave up.

Post by Ralf Hemmecke
Thus, in the whole integration there are
no Expression(X) elements anymore. I don't know why it makes sense to
keep a Expression(X)-element in the middle of the integration algorithm
although the only thing that is needed are the properties of that
element. So whether it is sin or cos or acrtan or whatever function, is
unimportant. They should become indeterminates (plus properties) during
the translation phase "Expression(X) --> tower of fields".

Post by Waldek Hebisch
Parts of integration code are quite general, for example
'TranscendentalIntegration' need just a field and few special
operations. Algebraic integration builds its own types which are in
some sense independent of Expression. So actually relatively small
part of integration code depends on special properties of Expression.
OTOH in part depending on Expression (dealing with kernels) rigid
structure of type tower would be inconvenient. In particular, we
handle some cases by re-arranging towers. At several places we
perform various ad-hoc extentions. And we need to see structure of
the whole field.

As I said, you know certainly better how the FriCAS integrator works,
but if I had to implement Risch, I would probably avoid the type
Expression. Maybe I can only see the problems with this
no-Expression(X)-approach, if I actually do such an implementation.

Let me put my point of view differently. It should not matter
what types we have as long as the types have needed properties.
And this is reasonably close to truth for current code. It
just happens that Expression has needed properties and there
are no other suitable type. Now the question is what would we
gain and at what cost by using different type? Some costs
are clear: we need to implement new type. Even if this is
just delegating real work to existing code we still need to
write something. More subtle cost is loss of flexibility:
presumably "tower type" would intoduce some contraints and
exclude other types. There is cost of convertions and
instantiation for new types. In particular, very
specific types require substantial effort.

Now, what would we gain?
Maybe better type checking -- but ATM I do not know how
to write tighter type constrants without going into
truble. One would like types which well correspond to
theory. Which means that we need essentially arbitrary
differential fields. Such fields can be represented
as subfields of Expression and we do this in current
code. In other words current code corresponds quite
well with theory, I see little room for improvemnent
here. So essentially I see no gains.

BTW: It is instructive to compare ExponentialExpansion
with MrvLimitPackage. ExponentialExpansion due to
type constraints can do only very special case,
while MrvLimitPackage do not have such limitation.
In principle one could introduce something like
GeneralExpansion and try to rewrite MrvLimitPackage
in terms of this. But we can do what is needed
using existing series types so I see no gain
from such change.