a rift in the continuum

[clicl...@freenet.de]

Yesterday I discovered a rift in the continuum of numbers - or perhaps
just another software bug. In order to make sure I need reliable values
of the upper incomplete gamma function (that is of the upper-case
Gamma), namely

Gamma(10^6 + 0.5, 10^6) = ?

and

Gamma(10^6 + 0.5, 1.003*10^6) = ?

to 100 decimal mantissa digits. Independent values (e.g. Maple and MMA)
would be nice for cross-checking.

Thank you,

Martin.

[Nasser M. Abbasi]

Mathematica V9:
In[230]:= N[Gamma[10^6 + 1/2, 10^6], 100]

Out[230]= \
4.13251479962262062103065157181206841728923017264300442213705199608773\
3766056346603731287344905361003*10^5565705


In[234]:= N[Gamma[10^6 + 1/2, 1003/1000 10^6], 100]
Out[234]= \
1.12717938591210252512349905106484944923179280972477237957652285580320\
7306530531823278504921383679816*10^5565703


Maple 17:
> evalf(GAMMA(10^6+1/2,10^6),100);
0.4132514799622620621030651571812068417289230172643004422137051996087733766\
5565706
056346603731287344905361003 10


evalf(GAMMA(10^6+1/2,1003/1000*10^6),100);
1.12717938591210252512349905106484944923179280972477237957652285\
5565703
5803207306530531823278504921383679652 10


--Nasser

[acer]

And it appears that Mathematica 9 got the last three of those digits wrong, for the second example.

And Maple 17's result ending in ...652 appears correct to within 1/2 an ulp (and to a few more places one might get ...652449).

[Miroslaw Kwasniak]

acer <ma...@rogers.com> wrote:
> And it appears that Mathematica 9 got the last three of those digits wrong, for the second example.
>
> And Maple 17's result ending in ...652 appears correct to within 1/2 an ulp (and to a few more places one might get ...652449).

Yes, the same results on wolframalpha ;) and GP/PARI 2.5.1

? default(realprecision,200); printf("%.100g\n%.100g\n",incgam(10^6 + 0.5,10^6),incgam(10^6 + 0.5, 1.003*10^6))
4.132514799622620621030651571812068417289230172643004422137051996087733766056346603731287344905361003 e5565705
1.127179385912102525123499051064849449231792809724772379576522855803207306530531823278504921383679652 e5565703

[acer]

Hmm. This seems a little severe. My (limited) understanding was that Mathematica did some sort of error-tracking with intervals, and that it should be getting all the requested digits correct. I would have expected that this is just the kind of operation for which the approach of manually-raising precision-then-returning-a-rounded-result should not have been necessary.

After all, this is just a single numeric special function call, and both of the input arguments should be getting passed in full (all inbound digits correct).

Maple doesn't promise `Digits` worth of correct digits for numerical evaluation of *compound* expressions, but it is supposed to get within 0.6 ulps on single special function calls like this (where the inbound digits are all passed, retained and not rounded since there are less than Digits' worth of them). And it accomplishes this, when successful, by using guard digits and argument range examination.

So this should be straightforward for Mathematica's error-tracking interval-arithmetic magic, no? But it (well, WolframAlpha) gets the last 3 or so digits wrong even for requests far less precise that in the original post. For example,

Requesting,

N[Gamma[10^6 + 0.5, 1.003*10^6], 20]

gets,

1.1271793859121024758 10^5565703

whereas the first 20 (chopped) correct digits of mantissa are,

1.1271793859121025251

Requesting,

N[Gamma[10^6 + 0.5, 1.003*10^6], 15]

gets,

1.12717938590655×10^5565703

whereas the first 15 (chopped) correct digits or mantissa are,

1.12717938591210


For these examples, as well a the original with 100 digits precision, Maple 17 produces the number of correct digits that matches the specified working precision, which is expected.

What's going on with Mathematica here? Is this wrong usage, or unexpected, or...?

[Nasser M. Abbasi]

On 5/9/2013 8:57 AM, acer wrote:

> Requesting,
>
> N[Gamma[10^6 + 0.5, 1.003*10^6], 20]
>
> gets,
>
> 1.1271793859121024758 10^5565703
>
> whereas the first 20 (chopped) correct digits of mantissa are,
>
> 1.1271793859121025251
>

If you write 1.003 in Mathematica, then it will use machine precision.

http://reference.wolfram.com/mathematica/tutorial/NumericalPrecision.html
"
123.4 a machine-precision number
123.45678901234567890 an arbitrary-precision number on some computer systems
123.45678901234567890` a machine-precision number on all computer systems
123.456`200 an arbitrary-precision number with 200 digits of precision
123.456``200 an arbitrary-precision number with 200 digits of accuracy
"

also
http://reference.wolfram.com/mathematica/howto/ControlThePrecisionAndAccuracyOfNumericalResults.html

"Use num`n to directly enter a number you know is correct to
n digits of precision, and num``n to directly enter a number
you know is correct to a digits of accuracy"

So, to get the same numerical value you show up:

In[144]:= N[Gamma[10^6 + 1/2, 1003/1000*10^6], 20]
Out[144]= 1.1271793859121025247*10^5565703
1.1271793859121025251 <---- your value

In[145]:= Precision[%]
Out[145]= 20.

> Requesting,
>
> N[Gamma[10^6 + 0.5, 1.003*10^6], 15]
>
> gets,
>
> 1.12717938590655�10^5565703
>
> whereas the first 15 (chopped) correct digits or mantissa are,
>
> 1.12717938591210
>

In[146]:= N[Gamma[10^6 + 1/2, 1003/1000*10^6], 15]
Out[146]= 1.12717938591210*10^5565703
1.12717938591210 <---- your value
In[147]:= Precision[%]
Out[147]= 15.

>
> For these examples, as well a the original with 100 digits precision,
>Maple 17 produces the number of correct digits that matches the specified
>working precision, which is expected.
>
> What's going on with Mathematica here? Is this wrong usage, or unexpected, or...?
>

--Nasser

[clicl...@freenet.de]

"Nasser M. Abbasi" schrieb:

This settles my original question. The new question is (as already
pointed out by acer) if this happens because of MMA's significance
arithmetic, or in spite of it. :)

The computation becomes still more difficult on the negative axes. Is
PARI/GP also able to handle

Gamma(- 10^6 + 0.5, - 10^6 + #i, 0)

to 100 decimal digits by incrementing the precision and then truncating?

Martin.

[acer]

Le jeudi 9 mai 2013 15:25:49 UTC-4, Nasser M. Abbasi a écrit :
[deleted]

> So, to get the same numerical value you show up:
>
>
>
> In[144]:= N[Gamma[10^6 + 1/2, 1003/1000*10^6], 20]
>
> Out[144]= 1.1271793859121025247*10^5565703
>
> 1.1271793859121025251 <---- your value
>
>
>
> In[145]:= Precision[%]
>
> Out[145]= 20.
>
>
>
> > Requesting,
>
> >
>
> > N[Gamma[10^6 + 0.5, 1.003*10^6], 15]
>
> >
>
> > gets,
>
> >
>

> > 1.12717938590655�10^5565703

>
> >
>
> > whereas the first 15 (chopped) correct digits or mantissa are,
>
> >
>
> > 1.12717938591210
>
> >
>
>
>
> In[146]:= N[Gamma[10^6 + 1/2, 1003/1000*10^6], 15]
>
> Out[146]= 1.12717938591210*10^5565703
>
> 1.12717938591210 <---- your value
>
> In[147]:= Precision[%]
>
> Out[147]= 15.
>

Could you please show (using that forward-tick ` notation, or otherwise) how to get Mathematica 9 or current WolframAlpha to produce 100 correct decimal digits of Gamma[10^6 + 1/2, (1003/1000)*10^6] without having to temporarily bump working precision to some arbitrary value higher than 100?

I don't see where you have shown that explicitly, though you seem to be implying that it's straightforwardly possible.

For example, what simple modification of syntax to,
N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 20]
would needed in order to obtain 20 (and only 20) correct decimal digits of a floating-point approximation?

And what simple modification of syntax to,
N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 100]
would be needed in order to obtain 100 (and only 100) correct decimal digits of a floating-point approximation?

In either case, an acceptable syntax should of course *not* require any temporary raise to working-precision by a problem-specific amount. What's the cleanest and preferably problem independent syntax?

Manually doubling of the working precision, followed by rounding down to the requested accuracy, is not ok since (if that were really necessary as technique) then there would exist some other example which would require trebling the working precision, etc. And given an arbitrary problem you wouldn't know how high to raise it, a priori. That's the whole point about Mathematica's significance arithmetic and error tracking -- that's its automatic.

If there is a straightforward syntax to obtain this clear objective, than hopefully someone will share it here.

[Nasser M. Abbasi]

On 5/9/2013 3:40 PM, acer wrote:

>
> Could you please show (using that forward-tick ` notation, or otherwise) how
>to get Mathematica 9 or current WolframAlpha to produce 100 correct decimal
>digits of Gamma[10^6 + 1/2, (1003/1000)*10^6] without having to temporarily
>bump working precision to some arbitrary value higher than 100?
>

Because integers have infinite accuracy, N[computation_wit_integers,n]
will always(*) manage to get the required n precision:

http://reference.wolfram.com/mathematica/tutorial/ArbitraryPrecisionNumbers.html

"If you start with an expression that contains only integers
and other exact numeric quantities, then N[expr, n] will
in almost all cases succeed in giving you a result to n digits
of precision. You should realize, however, that to do this
Mathematica sometimes has to perform internal intermediate
calculations to much higher precision. The global variable
$MaxExtraPrecision specifies how many additional digits
should be allowed in such intermediate calculations. "

But If you want to use floating points instead of integers, then
one needs to tell M the digits of precision of each explicitly.

Gamma[10^6 + 0.5`21, 1.003`21*10^6]
1.12717938591210247577*10^5565703
1.1271793859121025251 <---- your value

Precision[%]
20.699

> I don't see where you have shown that explicitly, though you seem to be implying
>that it's straightforwardly possible.
>
> For example, what simple modification of syntax to,
> N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 20]
> would needed in order to obtain 20 (and only 20) correct decimal
> digits of a floating-point approximation?
>

If you do not want to use exact values for input, but want to
use floating points in the input, then I do not know other than
using `nn. In the above for example, using `21 gave
atleast the 20 digits precision in the output.

From above link:

" Mathematica tries to give you results which have the
highest possible precision, given the precision of the
input you provided. "

"The precision of the output from a function can depend in
a complicated way on the precision of the input."

"If you give input only to a few digits of precision,
Mathematica cannot give you such high-precision output."

> And what simple modification of syntax to,
> N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 100]
> would be needed in order to obtain 100 (and only 100) correct
>decimal digits of a floating-point approximation?
>

---------------------------
Gamma[10^6 + 0.5`101, 1.003`101*10^6] (*<------ Precision[%] ==> 100.699*)
1.12717938591210252512349905106484944923179280972477237957652285580320\
73065305318232785049213836798160 * 10^5565703

Maple
1.12717938591210252512349905106484944923179280972477237957652285580320\
7306530531823278504921383679652 * 10^5565703
-------------------------

Gamma[10^6 + 0.5`102, 10^6] (*Precision[%]==>100.627*)
4.13251479962262062103065157181206841728923017264300442213705199608773\
37660563466037312873449053610031*10^5565705

Maple 17:
evalf(GAMMA(10^6+1/2,10^6),100);
0.413251479962262062103065157181206841728923017264300442213705199608773\
3766056346603731287344905361003*10^5565706

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

> In either case, an acceptable syntax should of course *not* require
>any temporary raise to working-precision by a problem-specific amount.
>What's the cleanest and preferably problem independent syntax?
>
> Manually doubling of the working precision, followed by rounding
>down to the requested accuracy, is not ok since (if that were really
>necessary as technique) then there would exist some other example
>which would require trebling the working precision, etc.
>And given an arbitrary problem you wouldn't know how
>high to raise it, a priori. That's the whole point about
>Mathematica's significance arithmetic and error tracking -- that's its automatic.
>
> If there is a straightforward syntax to obtain this clear objective,
>than hopefully someone will share it here.
>

--Nasser

[acer]

Le jeudi 9 mai 2013 22:20:23 UTC-4, Nasser M. Abbasi a écrit :
> On 5/9/2013 3:40 PM, acer wrote:
>
> > Could you please show (using that forward-tick ` notation, or otherwise) how
>
> >to get Mathematica 9 or current WolframAlpha to produce 100 correct decimal
>
> >digits of Gamma[10^6 + 1/2, (1003/1000)*10^6] without having to temporarily
>
> >bump working precision to some arbitrary value higher than 100?
>
> Because integers have infinite accuracy, N[computation_wit_integers,n]
>
> will always(*) manage to get the required n precision:
>
>
> http://reference.wolfram.com/mathematica/tutorial/ArbitraryPrecisionNumbers.html
>
>
> "If you start with an expression that contains only integers
>
> and other exact numeric quantities, then N[expr, n] will
>
> in almost all cases succeed in giving you a result to n digits
>
> of precision. You should realize, however, that to do this
>
> Mathematica sometimes has to perform internal intermediate
>
> calculations to much higher precision. The global variable
>
> $MaxExtraPrecision specifies how many additional digits
>
> should be allowed in such intermediate calculations. "

I do indeed want to use exact rational input, and not floating-point input, so let's focus on the above statements. If the claims are true then,

N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 20]

ought to produce 20 correct decimal digits of answer unless we are in the purportedly rare case that $MaxExtraPrecision needs to be bumped up manually, in a way that is example specific. (And similarly for the request for 100 digits of result, but let's focus on the shorter example.)

The output given by WolframAlpha today for,

N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 20]

is

1.1271793859121024758×10^5565703

And the last four digits of that mantissa are wrong.

Mathematica makes certain claims about obtaining the highest accuracy possible. But the example being considered is not some super difficult compound expression with nested function calls and approximate float input, or what have you. Rather, this is a single call to a well-known special function and the two inputs are supplied as exact rationals (which, incidentally, each require less than ten digits to approximate as radix-10 floats).

The idea that $MaxExtraPrecision need be set to some magic value that depends on the specific input values is quite contrary to notions of Mathematica's obtaining all requested digits correct automatically.

Why doesn't Mathematica know enough about Gamma[] to be able to set enough guard digits (ie. extra working precision, via $MaxExtraPrecision or whatever) to obtain the 20 correct requested digits of result? Ok, range reduction techniques can get difficult. But shouldn't Mathematica know enough about common special functions such as Gamma[] to be able to gauge the magnitude of the arguments and internally add enough guard digits? Note that current Maple gets both these 20-digit and 100-digit examples' computations fully correct, without any need to manually adjust settings.

This looks like a bug in the arbitrary precision Gamma[] implementation to me, judging by what's on that Mathematica help-page you've cited.

But maybe things are worse still. The page you cited has an example which claims that Mathematica cannot obtain all 30 correct digits, using default settings, of the floating-point approximation of N[Sin[10^100],30]. Getting 30 correct digits, automatically and with default settings, is obtained in current Maple with the comamnd,

> evalf[30](sin(10^100));

-0.372376123661276688262086695553

And the online help-page (cited above) shows that Mathematica can obtain the above result if one first sets $MaxExtraPrecision to 200. But why 200? How does one know whether that should be adequate? Why might not it need to be set to 250, or 30000, or...? This is quite unsatisfactory. Range reduction techniques for Sin[] ought to allow the internal working precision adjustment to be fully automatic for such exact rational input.

By the way, WolframAlpha falls down on this example, at the current date. For input,

N[Sin[10^100], 30]

it doesn't even show the error message like that web-page, of, "N::meprec: Internal precision limit $MaxExtraPrecision=50. reached while evaluating...." or similar. Instead, today, it just displays "0." as result.

[clicl...@freenet.de]

"Nasser M. Abbasi" schrieb:

>
> Because integers have infinite accuracy, N[computation_wit_integers,n]
> will always(*) manage to get the required n precision:
>
> http://reference.wolfram.com/mathematica/tutorial/ArbitraryPrecisionNumbers.html
>
> "If you start with an expression that contains only integers
> and other exact numeric quantities, then N[expr, n] will
> in almost all cases succeed in giving you a result to n digits
> of precision. You should realize, however, that to do this
> Mathematica sometimes has to perform internal intermediate
> calculations to much higher precision. The global variable
> $MaxExtraPrecision specifies how many additional digits
> should be allowed in such intermediate calculations. "
>

As you have shown in your first post, my two incomplete gamma cases do
not belong to MMA's "almost all cases" where the result is claimed to
have n digits precision. By the way, it is not necessary in these two
cases to perform intermediate calculations to higher precision, provided
the algorithm is well chosen. Thus a Derive calculation with just 100
digits gives

Gamma(10^6 + 0.5, 10^6) =

4.13251479962262062103065157181206841728923017264300442213705199~
6087733766056346603731287344905361004*10^(5.565705*10^6)

Gamma(10^6 + 0.5, 1.003*10^6, 0) =
1.12717938591210252512349905106484944923179280972477237957652285~
5803207306530531823278504921383679651*10^(5.565703*10^6)

Some erosion of the final digit(s) has then to be expected, of course.
The rift in the continuum, however, was caused by a bug in the Derive
6.10 numerics (a new one to me, but I know of another already).

Martin.

[Nasser M. Abbasi]

On 5/9/2013 11:01 PM, acer wrote:

> I do indeed want to use exact rational input, and not floating-point input,
>so let's focus on the above statements. If the claims are true then,
>
> N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 20]
>
> ought to produce 20 correct decimal digits of answer
>unless we are in the purportedly rare case that $MaxExtraPrecision needs
>to be bumped up manually, in a way that is example specific. (And similarly
>for the request for 100 digits of result, but let's focus on the shorter example.)
>
> The output given by WolframAlpha today for,
>
> N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 20]
>
> is
>
> 1.1271793859121024758�10^5565703
>
> And the last four digits of that mantissa are wrong.
>

I do not know how WAlpha works. But I just did a quick test. The above
seems to be a function some other effect, or itmight be using
different version on the web than the one that can be accessed
from inside M via the WolframAlpha[] command or some other formatting
issue.

http://reference.wolfram.com/mathematica/ref/WolframAlpha.html

Because when I tried the same command but via `=` inside Mathematica,
or when using WolframAlpha[] command from inside M, but
using the "Result" option, I get now the same result as Mathematica,
which is different from what you show above.

Here is a screen shot

http://12000.org/tmp/050913/gamma.jpg

-----------------------------
WolframAlpha["N[Gamma[10^6+1/2,(1003/1000)*10^6],20]","Result"]
1.127179385912102524742433496936420685852660370169`20.*^5565703

N[Gamma[10^6 + 1/2, (1003/1000)*10^6], 20]

1.1271793859121025247*10^5565703

Maple
evalf(GAMMA(10^6+1/2,1003/1000*10^6),20);
1.1271793859121025251 * 10^5565703

1.1271793859121024758 * 10^5565703 ----> WAlpha from the web
------------------------------

I do not know why they are different. May be someone knows, (i.e why
the web result shows different from M) But the first result above
is the same as direct M command as you can see and compares
to Maple's as well.

You are right. needed to increase the $MaxExtraPrecision:

In[214]:= $MaxExtraPrecision = 100;

N[Sin[10^100], 30]

Out[215]= -0.372376123661276688262086695553

--Nasser

[Nasser M. Abbasi]

On 5/10/2013 12:13 AM, Nasser M. Abbasi wrote:

>
> You are right. needed to increase the $MaxExtraPrecision:
>
> In[214]:= $MaxExtraPrecision = 100;
> N[Sin[10^100], 30]
> Out[215]= -0.372376123661276688262086695553
>

Make that 101 not 100 for $MaxExtraPrecision. With 100
I was still getting

N::meprec: Internal precision limit $MaxExtraPrecision = 100.` reached while evaluating

But I did not see the message on the console window.

In[242]:= $MaxExtraPrecision = 101;

N[Sin[10^100], 30]

Out[243]= -0.372376123661276688262086695553

--Nasser

[Dave Linder]

On Thursday, May 9, 2013 10:20:23 PM UTC-4, Nasser M. Abbasi wrote:
> On 5/9/2013 3:40 PM, acer wrote:
>
>
>
> >
>
> > Could you please show (using that forward-tick ` notation, or otherwise) how
>
> >to get Mathematica 9 or current WolframAlpha to produce 100 correct decimal
>
> >digits of Gamma[10^6 + 1/2, (1003/1000)*10^6] without having to temporarily
>
> >bump working precision to some arbitrary value higher than 100?
>
> >
>
>
>
> Because integers have infinite accuracy, N[computation_wit_integers,n]
>
> will always(*) manage to get the required n precision:
>
>
>
> http://reference.wolfram.com/mathematica/tutorial/ArbitraryPrecisionNumbers.html
>
>
>
> "If you start with an expression that contains only integers
>
> and other exact numeric quantities, then N[expr, n] will
>
> in almost all cases succeed in giving you a result to n digits
>
> of precision. You should realize, however, that to do this
>
> Mathematica sometimes has to perform internal intermediate
>
> calculations to much higher precision. The global variable
>
> $MaxExtraPrecision specifies how many additional digits
>
> should be allowed in such intermediate calculations. "

We can try it using Mathematica 9.0 instead of WolframAlpha.

In this next instance no warning about $MaxExtraPrecision is shown. The last four digits are incorrect.

Mathematica 9.0 for Linux x86 (64-bit)
Copyright 1988-2012 Wolfram Research, Inc.

In[1]:= N[Gamma[2000001/2, 1003000], 20]

5565703
Out[1]= 1.1271793859121024758 10

In[2]:= Quit


In a fresh session we can try increasing $MaxExtraPrecision anyway. The same result is returned, still with the last four digits incorrect.

Mathematica 9.0 for Linux x86 (64-bit)
Copyright 1988-2012 Wolfram Research, Inc.

In[1]:= $MaxExtraPrecision = 10000;

In[2]:= N[Gamma[2000001/2, 1003000], 20]

5565703
Out[2]= 1.1271793859121024758 10

In[3]:= Precision[%]

Out[3]= 20.


Using $MaxExtraPrecision = 10000000 doesn't help; the same result as above is returned.

By requesting 40 digits in the result the first 20 digits become correctly calculated and can be obtained by further processing.

In[4]:= N[Gamma[2000001/2, 1003000], 40];

In[5]:= N[%, 20]

5565703
Out[5]= 1.1271793859121025251 10

In[6]:= Quit

But it's not ideal if we can't immediately distinguish that it's necessary to take such an approach, and if we cannot tell how many more digits to request in order to obtain 20 correct significant figures.

This looks like an issue in the arbitrary precision floating-point calculation of `Gamma[..]` for exact rational arguments, because not only are several trailing digits incorrect but also no warning "N::meprec:..." is given and raising $MaxExtraPrecision doesn't appear to help either. These aspects are quite different than for the Sin[10^100] example.

[Richard Fateman]

If you think there is a bug in Mathematica, it should be reported
to Wolfram.

If you are surprised that Mathematica got fewer
correct digits than it claims, you could read
the documentation, and then you might not be so surprised.

Maybe disappointed that you didn't find a bug, just a feature.