# isn't that funny? -Infinity is +Infinity (with double precision)

Discussion in 'Symbolic Algebra' started by Peter Pein, Mar 15, 2008.

1. ### Peter PeinGuest

Dear group,

the following has been enterd into Maple 10. Is there a fix in later versions?
sqrt(x) GAMMA(x + 1/2)
f := x -> ----------------------
GAMMA(x + 1)
exp(-1) exp(1) exp(-1) exp(1)
------------------ - 1/8 --------------------
exp(-1/2) exp(1/2) exp(-1/2) exp(1/2) x

exp(-1) exp(1) 1
+ 1/128 --------------------- + O(----)
2 3
exp(-1/2) exp(1/2) x x
1 1 1
1 - --- + ------ + O(----)
8 x 2 3
128 x x

the following integral should evaluate to a finite value (Mathematica gives
0.00979687)
infinity
/ 1/2
| x GAMMA(x + 1/2) 1
i0 := | ------------------- - 1 + --- dx
| GAMMA(x + 1) 8 x
/
1
Float(infinity)
/ 1/2
| x GAMMA(x + 0.500000000000000000000000000000)
| ------------------------------------------------
| GAMMA(x + 1.)
/
1.

0.125000000000000000000000000000
- 1. + -------------------------------- dx
x
Classic Worksheet Interface, Maple 10.00, Windows, May 13 2005 B\
uild ID 190196

And why do I get 0.500... but not 1.000... in the evalf(i0,3*Digits)?

I've got difficulties to understand the error messages which Maxima 5.14 gives me:
(%i1) f(x):=sqrt(x)*gamma(x+1/2)/gamma(x+1)-1+1/8/x;
(%o1) f(x):=(sqrt(x)*gamma(x+1/2))/gamma(x+1)-1+(1/8)/x
(%i2) integrate(f(x),x,1,inf);
`quotient' by `zero'
-- an error. To debug this try debugmode(true);

where can there be a division by zero?

Maxima encountered a Lisp error:
Error in MACSYMA-TOP-LEVEL [or a callee]: 4000000.0 is not of type SEQUENCE.
Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.

Well, I know by myself, that 4000000.0 is not of type SEQUENCE. But how does
it come to this error?

And without specifications for accuracy and maximum steps I get the even more
cryptic:

***MESSAGE FROM ROUTINE DQAGI IN LIBRARY SLATEC.
***INFORMATIVE MESSAGE, PROG CONTINUES, TRACEBACK REQUESTED
* ABNORMAL RETURN
* ERROR NUMBER = 2
*
***END OF MESSAGE

Maxima encountered a Lisp error:

Error in CONDITIONS::CLCS-UNIVERSAL-ERROR-HANDLER [or a callee]: Bind stack
overflow.
Broken at PCL: RINT-STD-INSTANCE.
1 (Continue) Maxima top-level
dbl:MAXIMA>>

What is going on here?

Peter Pein, Mar 15, 2008

2. Dear Peter,

Here come three Maple 11.02 defects for your convergent integral
Int(x^(1/2)*GAMMA(x+.5000000000)/GAMMA(x+1.)-1.+.1250000000/x,
x = 1. .. Float(infinity))
Int(x^(1/2)*GAMMA(x+.50000000000000000000)/GAMMA(x+1.)-1.+
.12500000000000000000/x,x = 1. .. Float(infinity))

#
# After a couple of minutes...
#
Float(undefined) # HOWLER!

-------------------------------------------------------------

Contrast with the instant perfect Mathematica 6.0.2 results

NIntegrate[Sqrt[x] Gamma[x + 1/2]/Gamma[x + 1] - 1 + 1/(8 x),
{x, 1, Infinity}]

0.00979687

NIntegrate[Sqrt[x] Gamma[x + 1/2]/Gamma[x + 1] - 1 + 1/(8 x),
{x, 1, Infinity}, WorkingPrecision -> 30]

0.00979686520191666594405859200437

-------------------------------------------------------------

PP> What is going on here?

You did not upgrade to Mathematica 6.0.2 yet Cheers,

3. ### Peter PeinGuest

And I do not want to while I'm reading about problems in the HelpBrowser and
other problems.

But at least the newer Maple does not imply that some finite value is -inf and
calculated with double precision si +inf. It simply fails, which is OK to me,
because no program of this complexity can be perfect.

Cheerio!
Peter

Peter Pein, Mar 15, 2008
4. ### Axel VogtGuest

This is when Maple finds, that the results has relative errors larger
then requested (would be fine to have a mode, where it talks about it
with suggestions ...)

For your case best is to call with specific integration methods like

Int(g(x),x=1..infinity, method = _NCrule);
Int(g(x),x=1..infinity, method = _Dexp);

and use evalf on it. The Newton Cotes method usually works in even
very bad cases, the double exponential needs certain behavior in
infinity (if I remember correctly).

I always work with 14 Digits first, not 10 and get

.97968652014535e-2

This works in M10 and M11.

Axel Vogt, Mar 15, 2008
5. AV> .97968652014535e-2

Could you please show the exact input in
Maple 11/10 for which you got this number?

6. ### Axel VogtGuest

restart: interface(version); Digits:=14:

Classic Worksheet Interface, Maple 11.02, Windows, Nov 10 2007 Build ID 330022

f:= x -> sqrt(x)*GAMMA(x+1/2)/GAMMA(x+1):
g:= x -> f(x)-1+1/8/x:

Int(g(x),x=1..infinity, method = _NCrule): evalf(%);

0.0097968652014535

Int(g(x),x=1..infinity, method = _Dexp): evalf(%);

0.0097968652014535

Axel Vogt, Mar 15, 2008
7. ......................................................

Mathematica 6.0.2

NIntegrate[Sqrt[x] Gamma[x + 1/2]/Gamma[x + 1] - 1 + 1/(8 x),
{x, 1, Infinity}, WorkingPrecision -> 75]

0.00979686520191666594405859200435067912845184025280088677877823557479\
752620948

......................................................

Maple 11.02

restart; f:= x -> sqrt(x)*GAMMA(x+1/2)/GAMMA(x+1):g:= x -> f(x)-1+1/8/
x:
for j from 25 to 30 do
print(j, evalf(Int(g(x),x=1..infinity, method = _Dexp,j))); od;

25, Int(x^(1/2)*GAMMA(x+.5000000000000000000000000)/GAMMA(x+1.)-1.+.
1250000000000000000000000/x,x = 1. .. Float(infinity))

26, Int(x^(1/2)*GAMMA(x+.50000000000000000000000000)/GAMMA(x+1.)-1.+.
12500000000000000000000000/x,x = 1. .. Float(infinity))

27, Int(x^(1/2)*GAMMA(x+.500000000000000000000000000)/GAMMA(x+1.)-1.+.
125000000000000000000000000/x,x = 1. .. Float(infinity))

28, Int(x^(1/2)*GAMMA(x+.5000000000000000000000000000)/GAMMA(x+1.)-1.+.
1250000000000000000000000000/x,x = 1. .. Float(infinity))

29, Int(x^(1/2)*GAMMA(x+.50000000000000000000000000000)/GAMMA(x
+1.)-1.+.12500000000000000000000000000/x,x = 1. .. Float(infinity))

30, Int(x^(1/2)*GAMMA(x+.500000000000000000000000000000)/GAMMA(x
+1.)-1.+.125000000000000000000000000000/x,x = 1. .. Float(infinity))

restart; f:= x -> sqrt(x)*GAMMA(x+1/2)/GAMMA(x+1):g:= x -> f(x)-1+1/8/
x:
for j from 25 to 30 do
print(j, evalf(Int(g(x),x=1..infinity, method = _NCrule,j))); od;

25, Int(x^(1/2)*GAMMA(x+.5000000000000000000000000)/GAMMA(x+1.)-1.+.
1250000000000000000000000/x,x = 1. .. Float(infinity))

26, Int(x^(1/2)*GAMMA(x+.50000000000000000000000000)/GAMMA(x+1.)-1.+.
12500000000000000000000000/x,x = 1. .. Float(infinity))

27, Int(x^(1/2)*GAMMA(x+.500000000000000000000000000)/GAMMA(x+1.)-1.+.
125000000000000000000000000/x,x = 1. .. Float(infinity))

28, Int(x^(1/2)*GAMMA(x+.5000000000000000000000000000)/GAMMA(x+1.)-1.+.
1250000000000000000000000000/x,x = 1. .. Float(infinity))

29, Int(x^(1/2)*GAMMA(x+.50000000000000000000000000000)/GAMMA(x
+1.)-1.+.12500000000000000000000000000/x,x = 1. .. Float(infinity))

30, Int(x^(1/2)*GAMMA(x+.500000000000000000000000000000)/GAMMA(x
+1.)-1.+.125000000000000000000000000000/x,x = 1. .. Float(infinity))

8. ### Axel VogtGuest

....

As quite often you are a bit noisy. Instead try to read what I have
written before and you will find it out yourself.

Axel Vogt, Mar 15, 2008