WP34s: Gamma function for x very close to zero

[Dieter]

I just found what I think is a glitch in the WP34s Gamma function: arguments very close to zero can cause the calculator to hang and/or throw an error.

Code:

 1E-10     Γ     OK
 1E-20     Γ     OK
 1E-37     Γ     OK
 1E-38     Γ     +∞ Error
 1E-39     Γ     +∞ Error
 1E-40     Γ     calculator hangs and may eventually reset
 below     Γ     calculator hangs and may eventually reset
-1E-30     Γ     OK
-1E-38     Γ     OK
-1E-39     Γ     Domain Error
 below     Γ     Domain Error

I assume this is because somewhere in the code (using 39 digit precision) 1+some_small_value rounds to 1.

But there is a simple fix: If |x|<1E–20 simply evaluate \(\Gamma(x) \approx \frac{1}{x}-\gamma_{EM}\).
The absolute value of the relative error here is < 1E–40, so the result is as least as accurate as the internal 39 digit routine.

Dieter

[Paul Dale]

Yeah, there was a 1-x in there causing problems.
I've fixed the real versions of gamma and lngamma.
I've not touched the complex flavours at this stage.


Pauli

[Dieter]

(10-05-2015 02:25 AM)Paul Dale Wrote:  Yeah, there was a 1-x in there causing problems.
I've fixed the real versions of gamma and lngamma.
I've not touched the complex flavours at this stage.

Thank you very much. I see you set the threshold to 1E–24, both in Gamma and ln Gamma.

But there seems to be another problem with the regularized Gamma functions Γp and Γq for large arguments. Try e.g. 1E12 ENTER Γp or Γq. This should yield roughly 0,5 while it actually returns –22,6673 (!) resp. 1 – this = 23,6673.

I discovered this problem while posting in the Prime forum (cf. "Inverse Poisson" thread) where the errors in Γq caused problems in my Poisson quantile program.

Dieter

[Paul Dale]

(10-05-2015 07:32 PM)Dieter Wrote:  I see you set the threshold to 1E–24, both in Gamma and ln Gamma.

Reusing a constant I already had rather than introducing a new one.

Quote:But there seems to be another problem with the regularized Gamma functions Γp and Γq for large arguments. Try e.g. 1E12 ENTER Γp or Γq. This should yield roughly 0,5 while it actually returns –22,6673 (!) resp. 1 – this = 23,6673.

I'm about ready to give on on these functions and the incomplete beta


- Pauli

[brouhaha]

(10-06-2015 03:07 AM)Paul Dale Wrote:  I'm about ready to give on on these functions and the incomplete beta

And yet everyone bitched about Talking Barbie saying "Math is hard."

[Dieter]

(10-06-2015 03:07 AM)Paul Dale Wrote:  I'm about ready to give on on these functions and the incomplete beta

I am sorry to say this, but it's getting even worse:

Code:

-1,5  Γ     2,3632718...
-1,5  lnΓ   calculator hangs

-2,5  Γ    -0,9453087...
-2,5  lnΓ   should throw an error, but calculator crashes,
            requires hardware reset

Strange: Γ seems to work, while lnΓ crashes. The lnΓ(–2,5) case (should throw an error) caused the display of my hardware 34s to slowly fade away until it got inresponsive. Pressing the reset button in the battery compartment was required to bring it back to life (→ "Restored"). The batteries were not new, but the low battery symbol did not show up. The emulator (wp34sgui.exe) simply hangs and needs an Alt+F4 to close and restart.

Dieter

[emece67]

(10-07-2015 06:34 PM)Dieter Wrote:  

Code:

-1,5  Γ     2,3632718...
-1,5  lnΓ   calculator hangs

-2,5  Γ    -0,9453087...
-2,5  lnΓ   should throw an error, but calculator crashes,
            requires hardware reset

Strange: Γ seems to work, while lnΓ crashes. The lnΓ(–2,5) case (should throw an error) caused the display of my hardware 34s to slowly fade away until it got inresponsive. Pressing the reset button in the battery compartment was required to bring it back to life (→ "Restored"). The batteries were not new, but the low battery symbol did not show up. The emulator (wp34sgui.exe) simply hangs and needs an Alt+F4 to close and restart.

Dieter

A few weeks ago I modified (my local copy of) the source code for the Γ and LnΓ in order to save a few bytes combining both in a single function. An unexpected side effect of this change is that now it does not show this issue (Γp and Γq still fail for large arguments, though).

Perhaps Paul might like to take a look to this code. If so I'll PM it.

Regards.

[Paul Dale]

Weird, I can reproduce on my 34S (v3225) but not on the console emulator which is producing the correct results for these two cases.

No time to track down the issue at present...


Pauli

[Paul Dale]

(10-07-2015 11:21 PM)emece67 Wrote:  Perhaps Paul might like to take a look to this code. If so I'll PM it.

Can't hurt. Don't know when I'll get around to looking at it and integrating it.
Gamma, log gamma and the complex version could do with all being unified.


Pauli

[fhub]

(10-07-2015 06:34 PM)Dieter Wrote:  The emulator (wp34sgui.exe) simply hangs and needs an Alt+F4 to close and restart.

Not with the latest emulator version SVN 3813, here everything works as expected:

Code:


-1.5  Γ     2.3632718...
-1.5  lnΓ   0.8600470...

-2.5  Γ    -0.9453087...
-2.5  lnΓ   Domain Error

Franz

[Dieter]

(10-08-2015 11:19 AM)Paul Dale Wrote:  Weird, I can reproduce on my 34S (v3225) but not on the console emulator which is producing the correct results for these two cases.

(10-08-2015 01:06 PM)fhub Wrote:  Not with the latest emulator version SVN 3813, here everything works as expected:

FTR: My findings refer to version 3.3 3742 with printer support, both for the emulator and the physical device.

Dieter

[Marcus von Cube]

Pauli dis some work on the unified code. It did pay off! We have gained at least one library page (even two in the standard build without xtal or ir). The total savings is some 350 bytes. Changes like this allow for more GUI enhancements. I love it.

As usual, some testing is required!

[Dieter]

(10-11-2015 01:32 PM)Marcus von Cube Wrote:  Pauli dis some work on the unified code. It did pay off! We have gained at least one library page (even two in the standard build without xtal or ir). The total savings is some 350 bytes.

Great!

(10-11-2015 01:32 PM)Marcus von Cube Wrote:  Changes like this allow for more GUI enhancements. I love it.

There are other functions that require some updates, for instance the Poisson and Binomial CDF, PDF and Inverses. I have some (experimental) code for the quantile functions that seems to work fine (provided the Gamma functions and the PDFs work correctly), but this will require more memory than the current versions. But it will also return correct fractional results. ;-)

I would also appreciate if the one or other key function could be switched from XROM to C-code, for instance the Normal quantile. In the days of the original HP35 a design goal was to return results for all functions within a second. This goal is not met by the Q^–1 function, at least not for x or 1–x < 0,04. Which again is caused by the way the Gamma function is coded... #-)

Dieter