Doubts about correctness of `integrate(floor(x)^2,x)`

[Georgi Guninski]

I have theoretical reasons to doubt the correctness
of integrals involving `floor`.

The smallest testcases:

sage: integrate(floor(x)^2,x)
// Giac share root-directory:/usr/share/giac/
// Giac share root-directory:/usr/share/giac/
Added 0 synonyms
x*floor(x)^2

sage: integrate(2**floor(x),x)
2^floor(x)*x

Would someone check with another CAS or prove/disprove by hand?

[Jonathan Thornburg]

Maple reports the same result for your first testcase:

tux% maple
|\^/| Maple 2022 (X86 64 LINUX)
._|\| |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2022
\ MAPLE / All rights reserved. Maple is a trademark of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
> int(floor(x)^2, x);
2
floor(x) x

>

But, I think Maple and Sage/Giac are both wrong: consider the *definite*
integral (latex notation) $I = \int_0^{3/2} \lfloor x \rfloor^2 \, dx$:

Then the Sage/Giac indefinite integral implies that
$I = \left. x \lfloor x \rfloor^2 \right|^{3/2}_0
= (3/2) \lfloor (3/2) \rfloor^2 - 0 \lfloor 0 \rfloor^2
= 3/2$

But it seems to me that the correct result should be
$I = \int_0^{3/2} \lfloor x \rfloor^2 \, dx$
= \int_0^1 \lfloor x \rfloor^2 \, dx
+ \int_1^{3/2} \lfloor x \rfloor^2 \, dx
= \int_0^1 0 \, dx + \int_1^{3/2} 1 \, dx
= 0 + 1/2
= 1/2$

--
-- "Jonathan Thornburg [remove -color to reply]" <dr.j.th...@gmail-pink.com>
on the west coast of Canada
"Now back when I worked in banking, if someone went to Barclays,
pretended to be me, borrowed UKP10,000 and legged it, that was
`impersonation', and it was the bank's money that had been stolen,
not my identity. How did things change?" -- Ross Anderson

[Georgi Guninski]

On Sat, Jan 21, 2023 at 8:43 AM Jonathan Thornburg <jthor...@gmail.com> wrote:
>
>
> Maple reports the same result for your first testcase:
>

> But, I think Maple and Sage/Giac are both wrong: consider the *definite*
> integral (latex notation) $I = \int_0^{3/2} \lfloor x \rfloor^2 \, dx$:
>

Lol, a cross-CAS exploit.
I doubt sage and maple share source code, could
they share a buggy paper?

Here is a geometric disprove of the closed form.
Plot floor(x)^2 from 1 to 3. The area of the definite
integral is 1^2+2^2=5, which agrees with
integrate(floor(x)^2,x,1,3)
but disagrees with the indefinite integral.

[Oscar Benjamin]

The fundamental theorem of calculus requires that the integrand be
continuous over the interval of integration but in this case the
integrand is discontinuous at 1 so you cannot apply the theorem over
[0, 3/2]. The question is: what do you want to give as the indefinite
integral of a discontinuous function f given that its antiderivative
does not exist everywhere?

The expression F = floor(x)^2 * x is a valid antiderivative for f =
floor(x)^2 for all values of x such that f is continuous. At integer
values of x the function f is discontinuous and no (two-sided)
antiderivative exists.

It can still be possible to define a function F so that the integral
of f over [a, b] is equal to F(b) - F(a) but then is it worth doing
that rather than just returning a locally valid antiderivative that is
most likely a simpler expression. In this case a function F could be
something along the lines of (I haven't checked this carefully):

sum(n^2, (n, 1, floor(x)-1)) + floor(x)^2*(x - floor(x))

In context you should consider whether this is actually a better
expression or not though. Also keep in mind that it won't be possible
to do this for all discontinuous functions (e.g. 1/x). Any CAS that
promises to return indefinite integrals or antiderivatives of
potentially discontinuous functions will at least in some cases need
to return expressions that are only locally valid antiderivatives
rather than always returning a function F that can be used to compute
definite integrals blindly without checking continuity of the
integrand.

Given that you can't compute a suitable F in general, is it better to
try some of the time or just to state clearly in the documentation
what the limitations are and return the simplest locally valid
antiderivative expression?

--
Oscar

[parisse]

The antiderivative returned by giac (and by maple) for floor(x)^2 is only piecewise continuous and this is expected. But both CAS implement additional code to check for non continuous antiderivative (in simple situations for giac), and they correctly evaluate integrate(floor(x)^2,x,0,3/2) to 1/2. For definite integration without parameters, giac also checks that the exact answer matches a numeric value using approximate algorithms.

[Georgi Guninski]

[Georgi Guninski]

Thanks. sage appears to be slow on `integrate(floor(x)^2,x,0,n0)`,
is Maple fast?

Session:
sage: n0=10^6
sage: time integrate(floor(x)^2,x,0,n0)
CPU times: user 6.71 s, sys: 6.32 ms, total: 6.71 s
Wall time: 6.87 s
333332833333500000
sage: time integrate(x^3+x,x,0,n0)
CPU times: user 48.6 ms, sys: 2 µs, total: 48.6 ms
Wall time: 69.8 ms
250000000000500000000000

[parisse]

Maple seems (much) slower than giac on this example

giac:

0>> int(floor(x)^2,x=0..10000);
333283335000
// Time 0.2

maple:

int(floor(x)^2,x=0..10000);                                                  
memory used=142.1MB, alloc=150.1MB, time=1.12
memory used=230.5MB, alloc=182.1MB, time=1.81
memory used=267.3MB, alloc=182.1MB, time=2.21
memory used=383.1MB, alloc=215.8MB, time=3.02
memory used=455.8MB, alloc=215.8MB, time=3.59
memory used=519.3MB, alloc=215.8MB, time=4.08
memory used=577.8MB, alloc=215.8MB, time=4.50
memory used=633.2MB, alloc=215.8MB, time=4.91
memory used=685.8MB, alloc=215.8MB, time=5.29
memory used=736.9MB, alloc=215.8MB, time=5.66
memory used=786.5MB, alloc=215.8MB, time=6.07
memory used=834.8MB, alloc=215.8MB, time=6.44
memory used=882.3MB, alloc=215.8MB, time=6.77
memory used=928.9MB, alloc=215.8MB, time=7.12
memory used=974.6MB, alloc=215.8MB, time=7.51
memory used=1019.7MB, alloc=215.8MB, time=7.82
                                 333283335000

[Emmanuel Charpentier]

BTW :

```
sage: a, b = var("a, b")
sage: f(x) = floor(x)^2
sage: f(x).integrate(x, a, b)
// Giac share root-directory:/usr/local/sage-9/local/share/giac/
// Giac share root-directory:/usr/local/sage-9/local/share/giac/
Added 0 synonyms
No checks were made for singular points of antiderivative floor(sageVARa)^2*sageVARx for definite integration in [sageVARa,sageVARb]
-a*floor(a)^2 + b*floor(a)^2

```

Even accepting `x*floor(x)^2` as an antiderivative of `floor(x)`, this *definite* integral is wrong, *wrong*, **wrong**. One could expect :

```

sage: F(x) = f(x).integrate(x) ; F
x |--> x*floor(x)^2
sage: F(b) - F(a)
-a*floor(a)^2 + b*floor(b)^2

```

Something is amiss in Giac's definite integration. Is thois already known ?

[Emmanuel Charpentier]

BTW :

```

sage: a, b = var("a, b")
sage: f(x) = floor(x)^2

sage: F(x) = f(x).integrate(x) ; F
x |--> x*floor(x)^2

sage: def G(x): return numerical_integral(f, 0, x)[0]
sage: plot([F, G], (0, 3))
Launched png viewer for Graphics object consisting of 2 graphics primitives
```

Giac's "antiderivative" implicitly adds `heaviside(x-u)^2` terms... Again, this is wrong...

[Oscar Benjamin]

On Fri, 3 Feb 2023 at 09:31, Emmanuel Charpentier
<emanuel.c...@gmail.com> wrote:
>
> BTW :
>
> ```
> sage: a, b = var("a, b")
> sage: f(x) = floor(x)^2
> sage: f(x).integrate(x, a, b)
> // Giac share root-directory:/usr/local/sage-9/local/share/giac/
> // Giac share root-directory:/usr/local/sage-9/local/share/giac/
> Added 0 synonyms
> No checks were made for singular points of antiderivative floor(sageVARa)^2*sageVARx for definite integration in [sageVARa,sageVARb]
> -a*floor(a)^2 + b*floor(a)^2
> ```
>
> Even accepting `x*floor(x)^2` as an antiderivative of `floor(x)`, this *definite* integral is wrong, *wrong*, **wrong**. One could expect :
>
> ```
> sage: F(x) = f(x).integrate(x) ; F
> x |--> x*floor(x)^2
> sage: F(b) - F(a)
> -a*floor(a)^2 + b*floor(b)^2
> ```

I guess somewhere it is assumed that indefinite integration can be
used in FTOC without checking for discontinuities. The result is then
valid as a definite integral only if both end points lie within an
interval in which the function is continuous. In this particular case
that might not seem like a very useful guarantee but there are
probably other cases where it means that integrate can return
something useful.

Simple cases like floor(x) can be patched up but it's much harder to
solve this in general. The question then is really what guarantee you
want the integrate function to uphold when there are symbolic
parameters because most of the time the choice will be between
returning something like the above or just not returning anything.

A case that is unambiguously wrong is this definite integral that has
no symbolic parameters:

>>> integrate(atan2(sin(x), cos(x)), x)
1/2*arctan(sin(x)/cos(x))^2
>>> integrate(atan2(sin(x), cos(x)), (x,0,pi))
0

Here the integrand actually is continuous within the interior of the
region of integration and really does have an antiderivative but the
antiderivative expression returned by integrate is not differentiable
at pi/2 even though the integrand is continuous there. In this case
the result from SymPy is better:

>>> integrate(atan2(sin(x), cos(x)), x))
atan2(sin(x), cos(x))**2/2
>>> integrate(atan2(sin(x), cos(x)), (x, 0, pi))
pi**2/2

That antiderivative expression is continuous and differentiable
throughout the *interior* of (-pi, pi) just like integrand itself. You
still have to be careful about the end points though so if the
integral goes to pi then it needs to be known whether the
antiderivative expression is valid when approaching pi from the left
or the right so basically a limit should be computed. An example where
SymPy gets this kind of thing wrong is
(https://github.com/sympy/sympy/issues/24004):

>>> integrate(sin(x)*atan2(sin(x), cos(x)), (x, -pi, pi))
0
>>> integrate(sin(x)*atan2(sin(x), cos(x)), (x, -pi+0.001, pi)).n()
6.28318373671672
>>> integrate(sin(x)*atan2(sin(x), cos(x)), x)
sin(x) - cos(x)*atan2(sin(x), cos(x))

Here the integrand is nonnegative and continuous although not
differentiable at the end points. The antiderivative expression is
discontinuous though. This happens because even though the integrand
is continuous our symbolic expression for it contains discontinuous
functions. The issue is not with the fundamental theorem of calculus
but rather specifically with the limitations of *symbolic*
differentiation and integration with expressions involving non-smooth
functions.

Sage also gets this example wrong but in a different way:

>>> integrate(sin(x)*atan2(sin(x), cos(x)), (x, -pi, pi))
0
>>> integrate(sin(x)*atan2(sin(x), cos(x)), x)
-2*arctan(-2*sin(x)/((sin(x)^2/(cos(x) + 1)^2 - 1)*(cos(x) + 1)))
/(sin(x)^2/(cos(x) + 1)^2 + 1)
+ 2*sin(x)/((sin(x)^2/(cos(x) + 1)^2 + 1)*(cos(x) + 1))
+ 2*arctan(sin(x)/(cos(x) + 1))

That expression is definitely wrong because its derivative is negative
for half of every 2*pi interval.

It could be possible to patch up SymPy's antiderivative by adding
something like 2*pi*floor((x-pi)/(2*pi)) which would give a valid
antiderivative in a formal sense but then that means using more
discontinuous functions in antiderivative expressions. I am not sure
that it actually is better to try to make indefinite integration
always return complete antiderivatives where discontinuous functions
are involved. The same applies in reverse when thinking about
differentiation: adding in delta functions is probably not helpful.
Rather I think that it just needs to be understood and documented that
there are limits to
symbolic differentiation and (indefinite) integration where
discontinuous functions are involved.

The Maple docs explicitly discuss this giving floor(x) as an example.
It says there that:

> Note that the indefinite integral in Maple is defined up to a piecewise constant
> Hence, the results returned by int may be discontinuous at some points.
> For symbolic definite integration, two options control how discontinuities are handled.
> By default, int checks for discontinuities, and computes the integral as a sum of
> independent definite integrals, each of which involves an integrand which has
> no discontinuities in the interior of the interval of integration

https://www.maplesoft.com/support/help/maple/view.aspx?path=int%2Fdetails

--
Oscar

[Georgi Guninski]

One of the reasons I asked this is to get correct closed form
for stuff like sum/int 2^2^floor(x).
Judging by the discussions, this won't work.

[Oscar Benjamin]

The integration part of that is easy and sage can do it for a single
integer interval:

>>> a = var('a')
>>> assume(a, 'integer')
>>> integrate(2^2^floor(x), (x, a, a+1))
2^(2^a)

Now you just need the summation of that over the range of integer
intervals you are interested in but Sage's sum can't do anything with
that and I'm not sure if anything else can either.

--
Oscar

[Georgi Guninski]

> >>> a = var('a')
> >>> assume(a, 'integer')
> >>> integrate(2^2^floor(x), (x, a, a+1))
> 2^(2^a)
>
> Now you just need the summation of that over the range of integer
> intervals you are interested in but Sage's sum can't do anything with
> that and I'm not sure if anything else can either.
>

Thanks. I expected the definite integral to give solution.