12C+ Gamma



#10

Gamma(x+1)   (0 <= x <= 39.4)

01 STO 0 24 RCL 2
02 . 25 y^x
03 6 26 /
04 y^x 27 STO+ 3
05 1 28 1
06 9 29 STO+ 2
07 * 30 GTO 12
08 STO 1 31 RCL 3
09 0 32 2
10 STO 2 33 SQRTx
11 STO 3 34 LASTx
12 RCL 2 35 /
13 RCL 1 36 /
14 x<=y 37 2
15 GTO 31 38 ENTER
16 x<>y 39 LN
17 2 40 /
18 * 41 RCL 0
19 1 42 1
20 + 43 +
21 RCL 0 44 y^x
22 y^x 45 /
23 2 46 GTO 00

12.12 R/S => 648,976,950.1 after 3.5 seconds (HP-12C+)

On the 15C,

12.12 x! => 648,976,950.9

This is based on the following approximation I've come up with:


Gamma(x+1) ~ SUM(k=0,inf,(2*k+1)^x/2^k)/(sqrt2/2*(2/ln(2))^(x+1))

Obviously, this is not the proper way to compute the Gamma function. It's just an example of what can be done on the fast 12C+. Just for an idea, it would take about 7 minutes to compute Gamma(40) on a normal 12C using this method.

Gerson.

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

The approximation is being used reversely, that is, it is meant to approximate
this integer sequence, not the Gamma function:

an ~ n!/(2*sqrt2)*(2/ln(2))^(n+1) 

I found it when thinking about John Smitherman's math joke turned into a problem. The numbers 2, 6, 34, 294... (which showed as the problem was getting progressively more complex) ought to mean something... I couldn't find it at OEIS but I found 1, 3, 17, 147..., the halved sequence. Eventually I found out the original sequence is the row sums of Sierpinski-Pascal-MacMahon triangle, whatever this might be.
Next time someone starts a joke I'll try to just laugh, as everybody else :-)

I just intended to put a computationally intensive programming problem to the 12C+ to see how it performed. For computing Gamma function on the 12C+ or the normal 12C and 12c Platinum there isn't anything better than Egan Ford's program:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv016.cgi?read=100275

See message #6.

Edited: 26 June 2009, 10:30 a.m.


#11

Just for fun I thought I'd try it on the new 35s to see how it did (considering it's slighly(?) old CPU).



12.12 XEQ is about 30.5s with an answer of 648,976,950.899 (FIX = ALL).

40 XEQ is about 1 min.

Bart


#12

13 and 23 seconds, respectively, on the 33s which used to be my fastest RPN calculator. I wonder how nice it would be if they released a new 15C+. If it were as fast as the new 12C+, 60.5 CHS x! would take about 200 ms to display -1.527756e-97 instead of current 13 seconds.

Gerson.


#13

A fast 15C would be a dream come true.



I tried a HP-20S implementation:

01 LBL A	23 2		45 1/x
02 STO 0 24 + 46 *
03 y^x 25 1 47 2
04 . 26 = 48 =
05 6 27 y^x 49 y^x
06 = 28 RCL 0 50 (
07 * 29 = 51 RCL 0
08 1 30 2 52 +
09 9 31 y^x 53 1
10 = 32 RCL 2 54 )
11 STO 1 33 = 55 =
12 0 34 1/x 56 *
13 STO 2 35 * 57 2
14 STO 3 36 LAST 58 SQRT
15 LBL B 37 = 59 /
16 RCL 2 38 STO+ 3 60 2
17 INPUT 39 1 61 =
18 RCL 1 40 STO+ 2 62 1/x
19 x<=y? 41 GTO B 63 *
20 GTO C 42 LBL C 64 RCL 3
21 RCL 2 43 2 65 =
22 * 44 LN 66 RTN


Compared to the original post, one can see the advantage of RPN here. Although this is probably not the most optimum.



12.12 XEQ A, 15.5s; 648,976,950.9

40 XEQ A, 30s



Considering the program is 50% longer, that makes the 20S 3x faster than the 35s. Possibly the advantage of simple data types?

Bart


edited to try and make listing more readable

Edited: 30 June 2009, 11:09 a.m.

#14

HP-71B with Math Pack

a=time @ GAMMA(13.12); @ time-a

results in

648976950.903 .23

thus 0.23 seconds ....


The good old HP-71B ...


BR Raymund


Edited: 29 June 2009, 6:37 p.m.


#15

Quote:
648976950.903 .23
thus 0.23 seconds ....

Haven't you missed a "1" before ".23"?
1.22 seconds on mine...

...and 2.04 seconds for GAMMA(-68.5).

Gerson.


#16

a=time @ GAMMA(13.12); @ time-a
execution time between 0.22 up to 0.26


a=time @ GAMMA(-68.5); @ time-a
result = -1.527757e-97 execution time 1.01 (10 runs)





FYI VER$ => HP71:1BBBB HPIL:1A MATH:1A FTH:1A EDT:A KBD:B



I've tested it on three 71-B's with different memory sizes (16K up to 144K Ram), all about the same.


Maybe you have a Plug-In or Lex File causing a lot of Interrupts or Polls?

BR

Ray


#17

Mine must be on steroids:

.09 seconds for GAMMA(13.12).

It appears to be a very inaccurate measurement because several tries resulted in vastly disparate timing results, a few outcomes even with a negative sign.

#18

I was doing a=time @ GAMMA(13.12) @ time-a (without semicollon).

Now I get 0.21 to 0.22 and 1.02 to 1.05 for GAMMA(-68.5).

1BBBB version also with a 32K RAM module and card reader. They have no effect on timing though.

Thanks!

Gerson.


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP15c continued fraction for Ln(Gamma) Tom Grydeland 0 350 09-30-2013, 05:48 AM
Last Post: Tom Grydeland
  HP Prime emulator: Gamma function Stephan Matthys 28 2,231 08-21-2013, 04:52 PM
Last Post: Namir
  What is the Gamma approximation you use? Namir 21 1,791 08-05-2013, 07:14 AM
Last Post: Namir
  SandMath routine of the week: Inverse Gamma Function Ángel Martin 39 3,198 03-24-2013, 08:19 AM
Last Post: peacecalc
  [wp34s] Incomplete Gamma on the wp34s Les Wright 18 1,669 12-06-2011, 11:07 AM
Last Post: Namir
  HP 12C, 12C Platinum & 15C iOS App Walter Lam 2 479 06-02-2011, 01:25 PM
Last Post: Andrés C. Rodríguez (Argentina)
  gamma approximation in 48/49/50? Eric Smith 19 1,564 09-25-2009, 08:06 PM
Last Post: hugh steers
  Complex Gamma Function... revisited Ángel Martin 13 1,052 09-22-2009, 08:51 AM
Last Post: Ángel Martin
  Gamma (and Factorial) for 12C/12CP Les Wright 9 901 02-02-2009, 06:18 PM
Last Post: Egan Ford
  Gamma and Factorial for 10C/12C/12CP (thread hijack) Les Wright 1 389 01-28-2009, 08:19 PM
Last Post: Egan Ford

Forum Jump: