Digit loss mystery solved Message #26 Posted by Les Wright on 27 Aug 2009, 2:51 a.m., in response to message #25 by Les Wright
Steve, I ported your routine to Mathematica, set it up so that the argument was handled not as an exact number but an arbitrary precision decimal number to the desired number of digits. I worked with 100 just for the sake of speed.
Guess what? If we call "sum" the sum of all the positive terms and "sumneg" the sum of all the negative terms, these numbers happen to be almost identical, and with 100 digits of arbitrary precision the first 16 to 20 or so digits often look identical depending on the argument.
For illustration, here is my Mathematica code. Even if you don't speak Mathematica, it is pretty clear:
gamma_steve[zin_, n_] := Module[
{k, den, sum, sumneg, term, z},
z = N[Rationalize[zin], n];
a = n;
den = 1;
sum = Sqrt[2 Pi];
sumneg = 0;
For[k = 1, k < n, k++,
term = Sqrt[(a - k)^(2 k - 1)] Exp[a - k]/(den (z + k));
If[OddQ[k], sum += term, sumneg += term];
den *= k
];
{sum, sumneg, N[(sum - sumneg) (a + z)^(z + 1/2)/Exp[z + a], n]}
]
Here is output of sum, sumneg, and the final result for Pi and 100 digits. BTW, your code actually computes Gamma[z+1], hence the shift:
In[102]:= gamma_steve[Pi - 1, 100]
Out[102]=
{1.3383358053937979747418612286965372395777129410743133382110955900822128710027816852366593774862467173*10^54,
1.3383358053937953954626266996265458635993940570656707504694486183871951053643816586975959442344028621*10^54,
2.2880377953400324179595889090602339228896881533562224411993807454704710066085042825007}
Mathematica has a built in function that estimates the precision of the final result:
In[103]:= Precision[%]
Out[103]= 86.0209
The greater the desired precision, the closer together sum and sumneg are, and the more leading digits get lost to subtraction. Hence the need for about 20% more guard digits.
I am thinking Spouge ain't all that good for this application unless one is on a fast machine and can pay the price for so many more guard digits.
FYI, for those who know Mathematica, if I change the definition of z to z=Rationalize[zin], and pass it to the routine as an exact number, Mathematica uses it as such and carries out the computations exactly before converting it to a 100 digit decimal number at the end, maintaining full precision. The N[] command converts it to an approximate 100 digit number, and a cost is paid for that.
Well, at least we know know we aren't going crazy!
Hope you find this helpful. If you have access to Mathematica please email and I will send the notebooks.
Les
|