Incomplete gamma function et al. for 35s



#2

Here is a long promised contribution. It is sort of long, and perhaps should eventually be reworked as an article or software library submission.

The gamma function for real numbers is, as we all know, an extension of the factorial function. Gamma(a) is the integral of exp(-t)*t^(a-1) as t goes from 0 to infinity. If we split this interval of integration at some point x, we get two incomplete gamma functions, where the left-sided gamma(a,x) takes the integral from 0 to x and the right sided Gamma(a,x) takes it from x to infinity. Obviously, Gamma(a) = gamma(a,x) + Gamma(a,x).

The incomplete gamma functions may not be of great interest to most here, but certain of their familiar special cases--the error functions and the normal and chi-square probability functions--possibly are.

The following listing computes the incomplete gamma functions gamma(a,x) and Gamma(a,x). The code is adapted from the series and continued fraction routines given in Chapter 6 of Numerical Recipes. I just give the listing here and describe use below:

G001 LBL G
G002 STO X
G003 x<>y
G004 STO A
G005 STO B
G006 1/x
G007 STO C
G008 STO D
G009 STO E
G010 1
G011 STO+ B
G012 RCL X
G013 STO* C
G014 RCL B
G015 STO/ C
G016 RCL C
G017 STO+ D
G018 STO E
G019 RCL D
G020 x#y?
G021 GTO G009
G022 GTO G067

G023 STO X
G024 x<>y
G025 STO A
G026 -
G027 1
G028 =
G029 STO B
G030 1e250
G031 STO E
G032 1/x
G033 +
G034 1/x
G035 STO C
G036 STO D
G037 CLx
G038 STO I
G039 1
G040 STO+ I
G041 STO+ B
G042 STO+ B
G043 RCL A
G044 RCL- I
G045 RCL* I
G046 STO G
G047 RCL* C
G048 RCL+ B
G049 1e-250
G050 +
G051 STO C
G052 RCL G
G053 RCL/ E
G054 RCL+ B
G055 1e-250
G056 +
G057 STO E
G058 RCL C
G059 1/x
G060 x<> C
G061 /
G062 STO* D
G063 1
G064 x#y?
G065 GTO G039
G066 RCL D
G067 RCL X
G068 RCL A
G069 y^x
G070 *
G071 RCL X
G072 +/-
G073 e^x
G074 *
G075 RTN

G076 x^2
G077 0.5
G078 x<>y
G079 XEQ G002
G080 GTO G085

G081 x^2
G082 0.5
G083 x<>y
G084 XEQ G023
G085 PI
G086 SQRT
G087 /
G088 RTN

G089 x^2
G090 2
G091 /
G092 0.5
G093 x<>y
G094 XEQ G002
G095 2
G096 /
G097 PI
G098 SQRT
G099 /
G100 0.5
G101 +
G102 RTN

G103 x^2
G104 2
G105 /
G106 0.5
G107 x<>y
G108 XEQ G023
G109 2
G110 /
G111 PI
G112 SQRT
G113 /
G114RTN

G115 XEQ G122
G116 XEQ G002
G117 GTO G120

G118 XEQ G122
G119 XEQ G023
G120 RCL/ F
G121 RTN

G122 x<>y
G123 2
G124 /
G125 ENTER
G126 ENTER
G127 1
G128 -
G129 !
G130 STO F
G131 Rv
G132 x<>y
G133 2
G134 /
G135 RTN

a ENTER x XEQ G002 gives gamma(a,x). Best used when a < x. For a > x, compute gamma(a,x) = (a-1)! - Gamma(a,x).

a ENTER x XEQ G023 gives Gamma(a,x). Use for a > x usually. For a < x, use Gamma(a,x) = (a-1)! - gamma(a,x).

x XEQ G076 gives erf(x). I use for x < 1.8, and prefer for x > 1.8 erf(x) = 1 - erfc(x).

x XEQ G081 gives erfc(x). I use it for x > 1.8, and prefer erfc(x) = 1 - erf(x) for smaller values.

z XEQ G089 gives the lower tailed normal probability associated with a standard normal variate z--i.e., the probability that a normal variate takes on a value less z. I use it for z < 2.3. For larger values, I prefer to compute 1 - the upper tailed probability.

z XEQ G103 gives the upper tailed probability associated with a standard normal variate z--i.e., the probability that variate takes on a value exceeding z. I use it for z > 2.3, and for smaller values compute 1 - lower-tailed probability.

nu ENTER x XEQ G115 gives the lower tailed probability associated with a chi-squared variate x at nu degrees of freedom. Works best for nu < x. For nu > x, compute 1 - upper-tailed probability.

nu ENTER x XEQ G118 gives the upper tailed probability associated with a chi-squared variate x at nu degrees of freedom. Use for nu > x. For nu < x compute 1 - lower-tailed probability.

Here are some examples:

5 ENTER 4 XEQ G002 gives gamma(5,4) as 8.90791255566. The last digit should be an 8.

5 ENTER 6 XEQ G023 gives Gamma(5,6) as 6.84135600761. The last digit should be a 0.

1 XEQ G076 gives erf(1) = 8.42700792945e-1. The last two digits should be 50.

2 XEQ G081 gives erfc(2) = 4.67773498102e-3. The last digit should be a 5.

1.5 XEQ G089 gives 9.33192798730e-1 as the probability that a standard normal variate will be less than 1.5. The last digit should be a 1.

2.5 XEQ G103 gives 6.20966532587e-3 as the probability that a standard normal variate will exceed 2.5. The last two digits should be 77.

4 ENTER 3 XEQ G115 gives 4.42174599629e-1 as the probability that a chisquare variate at 4 df is less than 3. All digits are correct here.

2 ENTER 7 XEQ G118 gives 3.01973834223e-2 as the probability that a chisquare variate at 2 df exceeds 7. All digits are correct here.

A few caveats are in order. First of all, I don't do sign checking, so make sure that input for all functions is positive and meaningful. I am assuming that anyone with an interest in these functions is going to be aware of the usual reflection and complementary formulae. Second, the user needs to know about the basics of the computation to judge whether one launches a routine that enters a series computation or a continued fraction computation. This is very much a matter of judgement and personal taste--I give only rough guidelines. Finally, keep in mind that due to rounding and digit loss due to subtraction in some cases it is uncommon to expect all twelve digits to be correct. Indeed, I find that in general the twelfth digit is usually off a bit, and occasionally even the eleventh. In general, I think at least 10 digits are good the vast majority of the time, at least within 1 or 2 ULP.

Sadly, these routines, while pretty fast, and certainly much faster than equivalents I have for the 41CX and 42S, are slower than on the 33S. Regrettably, they don't fit the 33S paradigm as neatly, since they would gobble up oodles of those precious labels with the various entry points and GTOs etc.

Enjoy this work, and let me know what you think.

Les

Edited: 26 Aug 2007, 4:00 p.m.


#3

Les,

Many Thanks! or as we say in Germany "Ein herzliches Dankeschön!" for bringing back online the work you spent so much time on.

As another saying goes: "There are no Monuments for Critics!" - people like you set their own monuments with their contributions for and to an ever interested community of still learning readers!

Respectfully,

Peter A. Gebhardt

#4

Quote:
Sadly, these routines, while pretty fast.

How fast, do you have any timings?

Using the native INTEG/EQN support I get the following answers from your examples (FIX ALL):

Quote:
5 ENTER 4 XEQ G002 gives gamma(5,4) as 8.90791255566. The last digit should be an 8.

0
4
EQN T^(Z-1)/EXP(T)
INTEGRATE
INTEGRATE FN dT
Z? 5
=8.90791255568 (38 sec) (FIX 9: 8.907912556 (19 sec), FIX 4: 8.9079 (5 sec))
Quote:
5 ENTER 6 XEQ G023 gives Gamma(5,6) as 6.84135600761. The last digit should be a 0.

0
6
EQN T^(Z-1)/EXP(T)
INTEGRATE
INTEGRATE FN dT
Z? 5
=17.1586439924 (38 sec)
RCL Z
1
-
!
X<>Y
-
=6.8413560076

Edited: 26 Aug 2007, 6:33 p.m.


#5

A couple, three seconds tops in most cases, as far as I can tell.

Computing these functions by their series or CF expansions is typically much faster than integrating under the curve.

Les


#6

Quote:
Computing these functions by their series or CF expansions is typically much faster than integrating under the curve.

True. Its all about cost-benefit. Given how infrequent I would use the incomplete gamma function on my 35s, 38 seconds is not too bad. If I used it frequently I would invest the time to key in your wonderful routine. Its a pity that the 35s has no I/O, imagine all the code we could all share.

#7

Quote:
Its a pity that the 35s has no I/O, imagine all the code we could all share.

Yet there were some cheap SHARP organizers that used to come with a simple serial interface. The cable was optional, but it was easy to build: one 9-pin RS-232 plug in one end and one stereo plug in the other. I wonder adding a mini-USB port to the 35s would not cost so much. Also, I don't think this could be a reason for banning the calculator during tests.
Who is willing to fill up the 32KB RAM with programs and data just to have to enter them again in case of a crash?


#8

They propably thought that this only would sell to hardcore oldtimers that aint even using cellphones :-) Worked on me... (got a phone but its a ancient nokia)

Edited: 27 Aug 2007, 7:43 a.m.


#9

Hei, Arne!

Even older cell phones offered optional data cables. Well, at least the HP-35s has Continuous Memory :-)

Regards,

Gerson.

#10

Quote:
(edited for grammar and spelling)

They probably thought that this would sell only to hardcore old-timers who ain't even using cell phones :-)

Worked on me... (got a phone, but it's an ancient Nokia)


Several years ago, I got an "ancient" Nokia (model with documentation copyrighted 2002) pay-as-you-go cell phone for infrequent use. The youngsters at a recent family gathering remarked how it brought back memories...

Recently, network upgrades soon to be implemented have compelled its replacement, so now I must learn a new one and transfer the stored numbers.

This wasn't a problem with our HP calc's...

-- KS

Edited: 27 Aug 2007, 10:23 p.m.


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP50g: Writing a function that returns a function Chris de Castro 2 448 12-10-2013, 06:49 PM
Last Post: Han
  HP15c continued fraction for Ln(Gamma) Tom Grydeland 0 214 09-30-2013, 05:48 AM
Last Post: Tom Grydeland
  HP Prime emulator: Gamma function Stephan Matthys 28 1,431 08-21-2013, 04:52 PM
Last Post: Namir
  What is the Gamma approximation you use? Namir 21 1,180 08-05-2013, 07:14 AM
Last Post: Namir
  SandMath routine of the week: Inverse Gamma Function Ángel Martin 39 2,159 03-24-2013, 08:19 AM
Last Post: peacecalc
  [WP34s et al.] Solving the TVM equation for the interest rate Dieter 24 1,399 12-01-2012, 05:53 AM
Last Post: Paul Dale
  [WP34S] Improving the Incomplete Gammas Les Wright 19 1,091 05-16-2012, 05:15 PM
Last Post: Marcus von Cube, Germany
  [WP34S] Eagerly Anticipating Non-Regularized Incomplete Gammas Les Wright 6 460 05-13-2012, 06:40 PM
Last Post: Paul Dale
  Dave Ramsey et al with early HP 41C conversions to CL Geoff Quickfall 4 355 01-01-2012, 12:51 AM
Last Post: David Ramsey
  [wp34s] Incomplete Gamma on the wp34s Les Wright 18 1,080 12-06-2011, 11:07 AM
Last Post: Namir

Forum Jump: