(HP71B) Modern GAMMA formula for Forth or ASM?
07-29-2024, 03:46 PM
Post: #1
 floppy Senior Member Posts: 389 Joined: Apr 2021
(HP71B) Modern GAMMA formula for Forth or ASM?
Hello,
Gamma implementation can be different:
- HP41 RPN https://hp41programs.yolasite.com/gamma.php
- to 7 digits.. https://rosettacode.org/wiki/Gamma_function#Forth (we see the use of constants in other assembler programming)
- 2006 year formula https://www.rskey.org/CMS/the-library/?v...icle&id=11
Has anybody something more like anno 2024?
Thanks.

HP71B 4TH/ASM/Multimod, HP41CV/X/Y & Nov64d, PILBOX, HP-IL 821.62A & 64A & 66A, Deb11 64b-PC & PI2 3 4 w/ ILPER, VIDEO80, V41 & EMU71, DM41X
08-03-2024, 05:24 PM (This post was last modified: 08-03-2024 05:34 PM by peacecalc.)
Post: #2
 peacecalc Member Posts: 205 Joined: Dec 2013
RE: (HP71B) Modern GAMMA formula for Forth or ASM?
Hi floppy,

for over 30 years (not a knowledge we have today 2024) I programmed the gamma-function for a computer with a 368 processor and a 387 coprocessor (only for real arguments).

The interesting point is only that I used different expressions for different arguments:

1) for positive x

$x \ge 10: \Gamma(x) \approx \sqrt{\frac{2\pi}{x}}\cdot \exp\left(x*\ln(x) + (S_8 - x)\right)$

This is the expression - Stirling, but the infinite sum is reduced to 8 terms (is possible because x => 10):

$S_8 = \sum_{k=1}^{8} \frac{B_{2k}}{2k(2k-1)x^{2k-1}}$

The B(2k) are the Bernoullinumbers and the sum ist calculated with the Horner-scheme:

$S_8 = (((((((c_8x^{-2} + c_7)x^{-2} + c_6)x^{-2} + c_5)x^{-2} + c_4)x^{-2} + c_3)x^{-2} + c_2)x^{-2} + c_1)x^{-1}$

The c(k) are the precalculated expressions:

$c_k = \frac{B_{2k}}{2k(2k-1)}$

If we want to calculate for smaller positiv arguments let's us say x = 7.32, my former program calculates the value for 10.32 and in a loop it calculates:

$\Gamma(9.32) = \frac{\Gamma(10.32)}{9.32 } \qquad \Gamma(8.32) = \frac{\Gamma(9.32)}{8.32 } \qquad \Gamma(7.32) = \frac{\Gamma(8.32)}{7.32 }$

My test if this works with a good accuracy is:

$\Gamma(0.5) = \sqrt{\pi}$

Unfortunately I didn't document how many digit were correct.

2) And the negative Numbers x < 0 (without the negativ integers) were calculated with the expression:

$\Gamma(x) = \frac{\pi}{\sin(\pi\cdot x)\Gamma(1-x)}$

Let us say we have x = - 3.6 then we calculate:

$\Gamma(-3.6) = \frac{\pi}{\sin(\pi\cdot (-3.6))\Gamma(4.6)}$

This all was only a self educational exercise for learning programming the coprocessor with some aspects. May be someone else remembers or knows a more actual procedure if he reads that.
08-03-2024, 06:41 PM
Post: #3
 floppy Senior Member Posts: 389 Joined: Apr 2021
RE: (HP71B) Modern GAMMA formula for Forth or ASM?
Very interesting. Thanks.
And: using multiprocessing for math calculation (I have already done with python where separate calculation like Result = Series A / Series B where both series are independent) is something I had in my head for HP71 spreading tasks to other HP71 in a loop (just for fun; probably giving the task via HP-IL is too slow but could simulate a parallel processing).

HP71B 4TH/ASM/Multimod, HP41CV/X/Y & Nov64d, PILBOX, HP-IL 821.62A & 64A & 66A, Deb11 64b-PC & PI2 3 4 w/ ILPER, VIDEO80, V41 & EMU71, DM41X
08-04-2024, 11:48 PM (This post was last modified: 08-04-2024 11:58 PM by John Keith.)
Post: #4
 John Keith Senior Member Posts: 1,041 Joined: Dec 2013
RE: (HP71B) Modern GAMMA formula for Forth or ASM?
More information on Stirling's approximation can be found here and here. More details here under "Implementation". Assuming that a reasonable number of terms are sufficient for 12-digit accuracy, it would be faster and easier to use in-line constants for the coefficients in your code.

The HP 49 and 50 have the Gamma function for real and complex arguments. I haven't tried to examine the internal code myself but I would guess that it is rather involved.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)