#2207 Incorrect integral of x^(.5)*%e^(-x-x^2)

[Anonymous]

This is based on Sage bug report at http://trac.sagemath.org/sage_trac/ticket/11493.

(%i16) integrate(x^(1/2)*%e^(-x-x^2),x,0,inf);
inf
/ 2
[ - x - x
(%o16) I sqrt(x) %e dx
]
/
0

(%i18) integrate(x^(.5)*%e^(-x-x^2),x,0,inf);
<snip lots of these>
rat: replaced 0.5 by 1/2 = 0.5
sqrt(%pi)
(%o18) ---------
2

The second one is quite wrong, as quad_qags shows:

(%i28) quad_qags(x^(1/2)*%e^(-x-x^2),x,0,1000);
(%o28) [.3201570903601473, 1.5628860428051894e-9, 609, 0]
(%i29) quad_qags(x^(1/2)*%e^(-x-x^2),x,0,10000);
***MESSAGE FROM ROUTINE DQAGS IN LIBRARY SLATEC.
***INFORMATIVE MESSAGE, PROG CONTINUES, TRACEBACK REQUESTED
* ABNORMAL RETURN
* ERROR NUMBER = 4
*
***END OF MESSAGE

(%o29) [.3201753870288189, .03110853789056672, 7707, 4]

Separately, what's up with that error message?

[Robert Dodier]

See also https://sourceforge.net/p/maxima/bugs/2499/ for another example of a float causing trouble in an integral. There are probably other open bug reports on this topic.

[Raymond Toy]

Some tracing shows that Maxima eventually calls rischint with the given integrand. It returns sqrt(%pi)*erf(x)/2. This is clearly wrong. Also risch(sqrt(x)*exp(-x-x^2),x) returns a noun form. So that value of 0.5 confuses the Risch integrator somehow.

[Raymond Toy]

Oh, also:

(%i26) specint(sqrt(x)*exp(-x-x^2),x);
                  1/8
(%o26) (sqrt(2) %e    sqrt(%pi)
                   1           1  1               3  1        1
  sqrt(%pi) (gamma(─) bessel_i(─, ─) + bessel_i(- ─, ─) gamma(─))
                   4           4  8               4  8        4
 (───────────────────────────────────────────────────────────────
                            7/4       1
                           2    gamma(─)
                                      4
                    3           3  1               1  1        3
   sqrt(%pi) (gamma(─) bessel_i(─, ─) + bessel_i(- ─, ─) gamma(─))
                    4           4  8               4  8        4      3/4

 - ───────────────────────────────────────────────────────────────))/8
                             7/4       3
                            2    gamma(─)
                                       4
(%i27) float(%);
(%o27)                        0.3201570903601469
(%i28) quad_qagi(sqrt(x)*exp(-x-x^2),x,0,inf);
(%o28)        [0.3201570903601456, 8.176410659643807e-10, 285, 0]

And specint works with x^(0.5) too.