Incomplete gamma function et al. for 35s « Next Oldest | Next Newest »

 ▼ Les Wright Unregistered Posts: 1,368 Threads: 212 Joined: Dec 2006 08-26-2007, 04:00 PM 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. ▼ Peter A. Gebhardt Unregistered Posts: 174 Threads: 20 Joined: Sep 2006 08-26-2007, 05:25 PM 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 Egan Ford Unregistered Posts: 1,619 Threads: 147 Joined: May 2006 08-26-2007, 06:28 PM 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. ▼ Les Wright Unregistered Posts: 1,368 Threads: 212 Joined: Dec 2006 08-27-2007, 04:00 AM 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 ▼ Egan Ford Unregistered Posts: 1,619 Threads: 147 Joined: May 2006 08-27-2007, 05:14 AM 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. ▼ Gerson W. Barbosa Unregistered Posts: 2,761 Threads: 100 Joined: Jul 2005 08-27-2007, 07:31 AM 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? ▼ Arne Halvorsen (Norway) Unregistered Posts: 151 Threads: 18 Joined: Aug 2007 08-27-2007, 07:42 AM 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. ▼ Gerson W. Barbosa Unregistered Posts: 2,761 Threads: 100 Joined: Jul 2005 08-27-2007, 08:49 PM Hei, Arne! Even older cell phones offered optional data cables. Well, at least the HP-35s has Continuous Memory :-) Regards, Gerson. Karl Schneider Unregistered Posts: 1,792 Threads: 62 Joined: Jan 2005 08-27-2007, 10:21 PM 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 2,117 12-10-2013, 06:49 PM Last Post: Han HP15c continued fraction for Ln(Gamma) Tom Grydeland 0 1,109 09-30-2013, 05:48 AM Last Post: Tom Grydeland HP Prime emulator: Gamma function Stephan Matthys 28 7,948 08-21-2013, 04:52 PM Last Post: Namir What is the Gamma approximation you use? Namir 21 5,297 08-05-2013, 07:14 AM Last Post: Namir SandMath routine of the week: Inverse Gamma Function Ángel Martin 39 9,465 03-24-2013, 08:19 AM Last Post: peacecalc [WP34s et al.] Solving the TVM equation for the interest rate Dieter 24 5,983 12-01-2012, 05:53 AM Last Post: Paul Dale [WP34S] Improving the Incomplete Gammas Les Wright 19 4,908 05-16-2012, 05:15 PM Last Post: Marcus von Cube, Germany [WP34S] Eagerly Anticipating Non-Regularized Incomplete Gammas Les Wright 6 2,106 05-13-2012, 06:40 PM Last Post: Paul Dale Dave Ramsey et al with early HP 41C conversions to CL Geoff Quickfall 4 1,819 01-01-2012, 12:51 AM Last Post: David Ramsey [wp34s] Incomplete Gamma on the wp34s Les Wright 18 5,117 12-06-2011, 11:07 AM Last Post: Namir

Forum Jump: