#3235 ECL lisp arithmetic error in definite integration with large limits

[Dima Pasechnik]

With Maxima 5.38.1 compiled with ECL 16.1.2 on x86_64 Linux, one gets a lisp arithmetic error in an easy integral that does not lead to problems with the same ECL, but Maxima 5.35.1, or with Maxima 5.38.1 compiled with SBCL.

Maxima 5.38.1 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) integrate((x^2)*exp(x) / (1 + exp(x))^2,x,-1000,1000);

Maxima encountered a Lisp error:

 #<a ARITHMETIC-ERROR>

Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.

[Dima Pasechnik]

this popped up on https://trac.sagemath.org/ticket/18920

[Robert Dodier]

I see the result should be

(-(2000*%e^1000*log(%e^-1000*(%e^1000+1)))/(%e^1000+1))
-(2000*log(%e^-1000*(%e^1000+1)))/(%e^1000+1)
 -((2000*%e^1000+2000)*log(%e^1000+1)
 +(2*%e^1000+2)*li[2](-%e^1000)-1000000*%e^1000)
 /(%e^1000+1)+(2*%e^1000*li[2](-%e^-1000))/(%e^1000+1)
+(2*li[2](-%e^-1000))/(%e^1000+1)-1000000/(%e^1000+1)

My first guess is that this has something to do with factoring sums of exponentials which has caused trouble in various ways. To try to narrow the problem, can you enter trace(?factor); (the question mark means it's a Lisp function) and then try the same input again. I believe we'll be able to see whether the error occurs in a call to FACTOR. In any event you can append the output to this bug report and we'll go from there.

[Dima Pasechnik]

Here is the trace:

[dima@hilbert logs]$ ../sage --maxima
;;; Loading #P"/home/dima/Sage/sage-dev/local/lib/ecl/sb-bsd-sockets.fas"
;;; Loading #P"/home/dima/Sage/sage-dev/local/lib/ecl/sockets.fas"
;;; Loading #P"/home/dima/Sage/sage-dev/local/lib/ecl/defsystem.fas"
;;; Loading #P"/home/dima/Sage/sage-dev/local/lib/ecl/cmp.fas"
Maxima 5.38.1 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) trace(?factor);
(%o1)                              [factor]
(%i2) integrate((x^2)*exp(x) / (1 + exp(x))^2,x,-1000,1000);
1 Enter factor [1.0e-8 - epsilon]
                 100000000 epsilon - 1
1 Exit  factor - ---------------------
                       100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [1.0e-8 - g125]
                 100000000 g125 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g125 + 1.0e-8]
               100000000 g125 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g147]
                 100000000 g147 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g147 + 1.0e-8]
               100000000 g147 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g174]
                 100000000 g174 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g174 + 1.0e-8]
               100000000 g174 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g201]
                 100000000 g201 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g201 + 1.0e-8]
               100000000 g201 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [x + 1000]
1 Exit  factor x + 1000
1 Enter factor [1000 - x]
1 Exit  factor - (x - 1000)
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g257]
                 100000000 g257 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g257 + 1.0e-8]
               100000000 g257 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g284]
                 100000000 g284 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g284 + 1.0e-8]
               100000000 g284 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g311]
                 100000000 g311 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g311 + 1.0e-8]
               100000000 g311 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g338]
                 100000000 g338 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g338 + 1.0e-8]
               100000000 g338 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [x + 1000]
1 Exit  factor x + 1000
1 Enter factor [1000 - x]
1 Exit  factor - (x - 1000)
                   2   x
                  x  %e
1 Enter factor [----------]
                   x     2
                (%e  + 1)
                  2   x
                 x  %e
1 Exit  factor ----------
                  x     2
               (%e  + 1)
1 Enter factor [g390 + 1]
1 Exit  factor g390 + 1
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g482]
                 100000000 g482 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g482 + 1.0e-8]
               100000000 g482 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g509]
                 100000000 g509 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g509 + 1.0e-8]
               100000000 g509 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g536]
                 100000000 g536 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g536 + 1.0e-8]
               100000000 g536 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g563]
                 100000000 g563 - 1
1 Exit  factor - ------------------
                     100000000
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [g563 + 1.0e-8]
               100000000 g563 + 1
1 Exit  factor ------------------
                   100000000
1 Enter factor [x + 1000]
1 Exit  factor x + 1000
1 Enter factor [1000 - x]
1 Exit  factor - (x - 1000)
1 Enter factor [1.0e-8 - lim-epsilon]
                 100000000 lim-epsilon - 1
1 Exit  factor - -------------------------
                         100000000
1 Enter factor [prin-inf - 1.0e+8]
1 Exit  factor prin-inf - 100000000
1 Enter factor [1.0e-8 - g609]
                 100000000 g609 - 1
1 Exit  factor - ------------------
                     100000000
Maxima encountered a Lisp error:

 #<a ARITHMETIC-ERROR>

Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
1 Enter factor [x + 1000]
1 Exit  factor x + 1000
1 Enter factor [1000 - x]
1 Exit  factor - (x - 1000)
(%i3) 

[Robert Dodier]

Some debugging shows the problem is in 'limit'.

(%i2) trace (limit) $

(%i3) integrate((x^2)*exp(x) / (1 + exp(x))^2,x,-1000,1000);

1 Enter limit [2000]
1 Exit  limit 2000
1 Enter limit [-1000]
1 Exit  limit -1000
1 Enter limit [1000]
1 Exit  limit 1000
1 Enter limit [2000]
1 Exit  limit 2000
1 Enter limit [2000]
1 Exit  limit 2000
1 Enter limit [-1000]
1 Exit  limit -1000
1 Enter limit [1000]
1 Exit  limit 1000
1 Enter limit [2000]
1 Exit  limit 2000
1 Enter limit [2000]
1 Exit  limit 2000
1 Enter limit [-1000]
1 Exit  limit -1000
1 Enter limit [1000]
1 Exit  limit 1000
1 Enter limit [2000]
1 Exit  limit 2000
1 Enter limit 
[(-((2*x*%e^x+2*x)*log(%e^x+1))/(%e^x+1))
  -((2*%e^x+2)*li[2](-%e^x))/(%e^x+1)+(x^2*%e^x)/(%e^x+1),x,-1000,plus]
Maxima encountered a Lisp error:

 #<a ARITHMETIC-ERROR>

The simplest related example I can find:

Maxima branch_5_38_base_223_gcf9cbb2 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
(%i1) limit((-((2*x*%e^x+2*x))),x,-1000,plus);

                           - 1000         1000
(%o1)                    %e       (2000 %e     + 2000)
(%i2) limit(log(%e^x+1),x,-1000,plus);
                                - 1000    1000
(%o2)                     log(%e       (%e     + 1))
(%i3) limit(exp(x)*log(exp(x)+1),x,-1000,plus);
                        - 1000       - 1000    1000
(%o3)                 %e       log(%e       (%e     + 1))
(%i4) limit(x*exp(x)*log(exp(x)+1),x,-1000,plus);

Maxima encountered a Lisp error:

 #<a ARITHMETIC-ERROR>

[Robert Dodier]

The problem appears to stem from this:

$ ecl
ECL (Embeddable Common-Lisp) 16.1.2 (git:UNKNOWN)
> (> #.ext::double-float-negative-infinity (/ 1 100000000))

Condition of type: ARITHMETIC-ERROR
#<a ARITHMETIC-ERROR>

I'm not sure where the float infinity came from, but anyway if I work around this error, I get something that seems to be correct. Stay tuned.

I'll report the ECL bug too.

[Robert Dodier]

Bug reported to ECL issue tracker: https://gitlab.com/embeddable-common-lisp/ecl/issues/299

[Robert Dodier]

Fixed by commit e916dd0. Closing this report.

[Dan Gildea]

After this commit I get the following errors with both gcl and ecl.

The error seems to be caused by exp(1000) => floating point inf,
and seems to happen only after executing earlier files in the test suite.

Running tests in rtest_limit: 
********************** Problem 212 (line 801) ***************
Input:
block([actual, expected], actual : limit(x exp(x) log(exp(x) + 1), x, - 1000, 

                             - 1000       - 1000    1000
plus), expected : (- 1000) %e       log(%e       (%e     + 1)), 
if ev(equal(actual, expected), logexpand = 'super) then true
 else [actual, expected])


Result:
bfloat: attempted conversion of floating-point infinity.
error-catch

This differed from the expected result:
true

********************** Problem 213 (line 813) ***************
Input:
                                                2
                                               x  exp(x)
block([actual, expected], actual : integrate(-------------, x, - 1000, 1000), 
                                                         2
                                             (1 + exp(x))
                    1000       - 1000    1000

           - 2000 %e     log(%e       (%e     + 1))
expected : ----------------------------------------
                            1000
                          %e     + 1
              - 1000    1000
   2000 log(%e       (%e     + 1))               1000               1000
 - ------------------------------- + (- ((2000 %e     + 2000) log(%e     + 1)
               1000
             %e     + 1
        1000              1000                  1000      1000
 + (2 %e     + 2) li (- %e    ) + (- 1000000) %e    ))/(%e     + 1)
                    2
       1000         - 1000              - 1000
   2 %e     li (- %e      )   2 li (- %e      )
              2                   2               - 1000000
 + ------------------------ + ----------------- + ----------, 
            1000                   1000             1000
          %e     + 1             %e     + 1       %e     + 1
if ev(equal(actual, expected), logexpand = 'super) then true
 else [actual, expected])


Result:
bfloat: attempted conversion of floating-point infinity.
error-catch

This differed from the expected result:
true

211/213 tests passed

[Robert Dodier]

The problem appears to be that if there is a bigfloat in the assume database for some reason, Maxima tries to compare negative float infinity to the bigfloat, and then tries to promote the float to a bigfloat, and gets the error. I'll open a new bug report about that.

[Robert Dodier]

I've opened #3244 to work on this problem.

[Dima Pasechnik]

The patch provided on https://gitlab.com/embeddable-common-lisp/ecl/issues/299 does fix the original issue, so perhaps a workaround is just not necessary (with their master branch, at least).