search

rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

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

Closed rtoy closed 4 months ago

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:15:49 Created by dimpase on 2016-10-29 00:03:54 Original: https://sourceforge.net/p/maxima/bugs/3235

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.
rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:15:51 Created by dimpase on 2016-10-29 00:06:56 Original: https://sourceforge.net/p/maxima/bugs/3235/#3fcd

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

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:15:53 Created by robert_dodier on 2016-10-29 18:18:39 Original: https://sourceforge.net/p/maxima/bugs/3235/#a0a4

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.

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:15:55 Created by dimpase on 2016-10-29 18:53:50 Original: https://sourceforge.net/p/maxima/bugs/3235/#2fbe

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)
rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:15:57 Created by robert_dodier on 2016-10-30 00:43:06 Original: https://sourceforge.net/p/maxima/bugs/3235/#bbc4

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>
rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:15:59 Created by robert_dodier on 2016-10-30 00:44:04 Original: https://sourceforge.net/p/maxima/bugs/3235/#83c2

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:16:01 Created by robert_dodier on 2016-10-30 01:57:24 Original: https://sourceforge.net/p/maxima/bugs/3235/#06d2

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.

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:16:03 Created by robert_dodier on 2016-10-30 05:08:35 Original: https://sourceforge.net/p/maxima/bugs/3235/#6c8b

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

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:16:06 Created by robert_dodier on 2016-11-03 14:37:39 Original: https://sourceforge.net/p/maxima/bugs/3235/#72dd

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:16:08 Created by robert_dodier on 2016-11-03 14:37:40 Original: https://sourceforge.net/p/maxima/bugs/3235/#427f

Fixed by commit e916dd0. Closing this report.

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-02 08:16:10 Created by dgildea on 2016-11-07 10:04:57 Original: https://sourceforge.net/p/maxima/bugs/3235/#9c3d

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