Gamma Function Using Spouge's Method

[Namir]

Gamma Function using Spouge's method.

Memory Map

R00 = x and = x-1
R01 = a
R02 = CHS
R03 = Sum
R04 = I
R05 = sqrt(2*pi)
R06 = integer part of I

Implementation

Code:

1 LBL "GAMMA"
2 LBL A
3 "X?"
4 PROMPT
5 1
6 -
7 STO 00
8 12.5
9 STO 01        # a = 12.5
10 1
11 STO 02        # CHS = 1
12 STO 03        # Sum = 1
13 2
14 PI
15 *
16 SQRT
17 STO 05        # sqrt(2*pi)
18 1.012
19 STO 04        # set up loop control variable
20 LBL 00        # start the loop
21 RCL 02
22 RCL 05
23 /
24 RCL 04
25 INT
26 STO 06
27 1
28 -
29 FACT
30 /
31 RCL 01
32 RCL 06
33 -
34 RCL 06
35 0.5
36 -
37 Y^X
38 *
39 RCL 01
40 RCL 07
41 -
42 EXP
43 *
44 RCL 06
45 RCL 00
46 +
47 /
48 STO+ 03            # Sum = Sum + C/(X+I)
49 RCL 02
50 CHS
51 STO 02        # CHS = -CHS
52 ISG 04        # end of loop
53 GTO 00
54 RCL 00
55 RCL 01
56 +
57 STO 06
58 RCL 00
59 0.5
60 +
61 Y^X
62 RCL 06
63 EXP
64 /
65 RCL 05
66 *
67 RCL 03
68 *
69 "GAMA="
70 ARCL X
71 PROMPT
72 GTO A

[lcwright1964]

This is great work, Namir. The Spouge approximation isn't quite as good as the related Lanczos approximation for for the same number of terms in the series, but it is much easier to program, especially on a calculator as the coefficients are simpler.

Les

P.S. I have used JMB's GAM+ program in the "old" program listings as my Gamma of choice for years. Will be interesting to see if this performs as well. Jean-Michel used the first few terms of standard Stirling series, cleverly factored Horner-style to save steps, IIRC. He also stuck to synthetic registers.

[lcwright1964]

Namir,

I looked at the original Spouge paper from 1994 (attached). There is heavy theoretical math in the proofs, and I don't follow it all, but the max relative error is give by a fairly simple formula in Theorem 1.3.1 at the top of page 934. Basically, the larger a is, the smaller the relative error, and the number of terms N = Ceiling(a) - 1. You use a =12.5, but a = 13 would do better (in theory) than your choice of 12.5, and you still only do the 12 loops.

To compare, Spouge's relative error for a = 12.5 is about 1.2e-11, whereas it is about 4.7e-12 for a =13.

That said it may not make a lot of difference in HP41 or the HP67. I understand that once you start wracking up the terms the promised theoretical precision can get eclipsed by rounding error if one doesn't have a lot of guard digits on hand. I know this is a consideration in arbitrary precision environments (e.g., I read someplace that if you want 45 digits accurate you need to work with 70), but maybe the impact is less with the 13 internal digits of our vintage machines.

I don't know if you made your choice of a based on the original paper or another source, but I thought this was interesting.

Les

[Dieter]

(08-12-2015 04:03 AM)lcwright1964 Wrote:  To compare, Spouge's relative error for a = 12.5 is about 1.2e-11, whereas it is about 4.7e-12 for a =13.

I assume that's the calculated upper bound of the error.
But things are different on real calculators. ;-)

(08-12-2015 04:03 AM)lcwright1964 Wrote:  That said it may not make a lot of difference in HP41 or the HP67. I understand that once you start wracking up the terms the promised theoretical precision can get eclipsed by rounding error if one doesn't have a lot of guard digits on hand.

That's exactly the crucial point. The theoretical error that would be possible if (!) all calculations were done with infinite (or at least sufficient) precision is much lower than what we can expect from a real calculator with, here, 10 digit working precision. The three additional guard digits you mentioned are only used internally, so they do not matter here. On the '41 or '67, all we got are ten digits.

Now let's take a look at a few different values of the variable "a" and see what we get:

Mathematically, a value as low as a=9 is sufficient for a 10-digit result over the whole (positive) 41/67 working range, i.e. for input between 0 and approx. 71 (overflow limit). Starting at x=0 the error increases with x, reaches a maximum and then decreases again as x gets larger. The largest relative error is 5,33 E–11 near x=32,8. Which is fine for 10 digits (maybe +1 ULP).

However, on a real '41 or '67 the last two or even three digits are off. Increasing the value of "a" does not neccessarily help here. For instance, a=13 would be good for a largest relative error of 2,72 E–15 (near x=77,9). But in real life the error gets even larger (!) than for a=9.

A value of a=15 theoretically yields an error less than 4 E–17, i.e. full 15–16 digit IEEE double precision. I tried this in Excel VBA as well as on a WP34s in SP mode and got about 10 valid digits near the error maximum.

So indeed the working precision has to be greater than the desired accuracy of the result. I assume that 13 digits would suffice for a 10 digit result. But this cannot be accomplished in a simple user program. On the '41, a 13-digit MCode implementation would be possible.

As long as no additional guard digits are available, a careful choice of "a" is crucial. I would suggest a=9 for the 10-digit machines. Larger values do not seem to provide better accuracy. About 8 digits is all we get.

Addendum: This also means that a value as high as a=12,5 does not make much sense here. Try Namir's program with a=7 (replace 12.5 with 7 and 1.012 with 1.006). My tests returned results at least as accurate as with a=12,5 (mostly even better) while the program runs about twice as fast. BTW there is an error in the original listing: line 40 should read RCL 06.

Dieter

[lcwright1964]

(08-12-2015 05:39 PM)Dieter Wrote:  So indeed the working precision has to be greater than the desired accuracy of the result. I assume that 13 digits would suffice for a 10 digit result. But this cannot be accomplished in a simple user program. On the '41, a 13-digit MCode implementation would be possible.

Indeed. I believe this is exactly what we have in Angel Martin's SANDMATH.

Some years back I was quite engaged with approximating gamma. I learned that Lanczos and Spouge were attractive to some because they were convergent and as such one didn't get into the shift and divide thing for small arguments. Despite this, I think that programmers to this day go with the good old-fashioned Stirling's series, since in practice Spouge and Lanczos don't offer much more in the way of efficiency or accuracy.

That said, I really can see that calculator programmers are drawn to the Spouge formula, as the coefficients are computed simply and on the fly. Nonetheless, JM Baillairds excellent GAM+ program in the old HP41 program archive here cleverly rearranges a truncated Stirling approximation in a way that maximizes efficiency and storage of coefficients.

Les

[Dieter]

(08-12-2015 05:39 PM)Dieter Wrote:  Addendum: This also means that a value as high as a=12,5 does not make much sense here. Try Namir's program with a=7 (replace 12.5 with 7 and 1.012 with 1.006). My tests returned results at least as accurate as with a=12,5 (mostly even better) while the program runs about twice as fast.

FTR, here is my implementation for the '41, based on a "hard-wired" a=7.
Like the factorial function on various other HPs, it actually calculates Γ(x+1).

Code:

01 LBL "GAM+1"
02 STO 00
03 6
04 +
05 STO 04
06 5,5
07 STO 02
08 INT
09 FACT
10 STO 03
11 SIGN
12 STO 01
13 CLx
14 LBL 01
15 RCL 01
16 RCL 02
17 y^x
18 RCL 01
19 e^x
20 *
21 RCL 03
22 /
23 RCL 04
24 /
25 +
26 RCL 02
27 INT
28 X=0?
29 GTO 02
30 CHS
31 ST/ 03
32 SIGN
33 ST- 01
34 ST+ 02
35 ST+ 04
36 RDN
37 GTO 01
38 LBL 02
39 RDN
40 ABS
41 Pi
42 ST+ X
43 SQRT
44 +
45 RCL 00
46 7
47 +
48 STO 03
49 RCL 00
50 ,5
51 +
52 y^x
53 *
54 RCL 03
55 e^x
56 /
57 END

Since it doesn't seem to get much better than 8 valid digits in the final result, a=7 is sufficient here. Evaluated exactly, the largest relative error is approx. 2,3 E–9. Actually the error may be slightly higher. But it should be less than in the original a=12,5 version.

Some examples:

Code:

 3   XEQ"GAM+1"  =>  5,999999994
 5   XEQ"GAM+1"  =>  120,0000000
 7   XEQ"GAM+1"  =>  5040,000009
8,35 XEQ"GAM+1"  =>  85899,08333      (exact: ...363)
17,3 XEQ"GAM+1"  =>  8,416039696 E+14 (exact: ...738)
50,5 XEQ"GAM+1"  =>  2,166683759 E+65 (exact: ...771)

The current version throws an OUT OF RANGE error for x > 54,9. Maybe this can be improved somehow.

Edit: the following modification may be a solution. However, the resulting accuracy is slightly less than in the first version.

Code:

...
45 RCL 00
46 7
47 +
48 STO 03
49 LN
50 RCL 00
51 ,5
52 +
53 *
54 RCL 03
55 -
56 e^x
57 *
58 END

This way x=69,9575744 returns 9,999997034 E+99 (exact: ...7560).

Dieter

[lcwright1964]

Thanks, Dieter, this is excellent--and very compact too.

Now that I am revisiting the gamma/factorial for calculators issue, it is coming back to me why Spouge and Lanczos, though theoretically attractive because they are convergent and accurate for small real arguments, have seemed less attractive for our calculator work. It is that the just don't do better than the old fashioned asymptotic formulae, even if in theory they should. In Spouge especially, on computes coefficients on the fly and it is no more accurate.

I don't think I ever managed to improve on J-M Baillard's GAM+ implementation in the old HP41 program archives. His 70-step program has some key features--it uses synthetic registers when storage is needed, does most everything on the stack, cleverly rearranges the Stirling series to simplify the constants embedded in the code, and uses the calculator's FACT function if the argument is an integer. It does the Pochhammer shift for arguments less than 5, which seems to suffice in a 10-digit environment to give at least eight good digits and usually nine for most anything.

I also rediscovered Peter Luschny's superb page:

http://www.luschny.de/math/factorial/app...Cases.html

I think I was once drawn to this pseudocode (edited by me to include Luschny's suggested simplification in the formula) which implements the Stieltjes continued fraction:

Stieltjes3Factorial(x)
y = x+1; p = 1;
while y < 7 do p = p*y; y = y+1; od;
r = exp(y*(log(y)-1)+1/(2*(6*y+1/(5*(y+1/(4*y))))));
if x < 7 then r = r/p fi;
return r*sqrt(2*Pi / y) end;

However, despite it's compactness and elegance and guarantee of at least 8 digits, I don't think I was ever able to better my favourite Baillard program.

Fun times!

Les

[Ángel Martin]

(08-12-2015 09:11 PM)Dieter Wrote:  This way x=69,9575744 returns 9,999997034 E+99 (exact: ...7560).

Yes I used 13-digit routines for the MCODE version in the SandMath, which uses the Lanczos formula - and indeed the result for x=70.9575744 is exactly:

9.9999975 99,

and using VMANT to see the mantissa:
9.999997560

So it appears to hold up nicely in that region.

[Dieter]

(08-19-2015 08:56 PM)Ángel Martin Wrote:  Yes I used 13-digit routines for the MCODE version in the SandMath, which uses the Lanczos formula - and indeed the result for x=70.9575744 is exactly:

I wouldn't have expected anything less. ;-)

The accuracy of the Lanczos method, like many others, increases with the argument. The larger x, the more valid digits you get, c.f. Peter Luschny's Gamma/factorial pages. So large values are the easy part.

Now the trick is to set a minimum x with sufficient accuracy for which the Lanczos (or another) approximation is applied, and evaluate Gamma for smaller arguments by the usual shift/divide trick.

BTW, Sandmath's GAMMA nicely overflows at x=70,95757446 while x=70,95757445 returns 9,999999688 E+99 (exact: ...687).

Dieter

[Dieter]

(08-19-2015 05:32 PM)lcwright1964 Wrote:  I think I was once drawn to this pseudocode (edited by me to include Luschny's suggested simplification in the formula) which implements the Stieltjes continued fraction:

Stieltjes3Factorial(x)
y = x+1; p = 1;
while y < 7 do p = p*y; y = y+1; od;
r = exp(y*(log(y)-1)+1/(2*(6*y+1/(5*(y+1/(4*y))))));
if x < 7 then r = r/p fi;
return r*sqrt(2*Pi / y) end;

As far as I can see the last x<7 test can be omitted. I tried this in VBA for Excel:

Code:

Function Gamma(x)
y = x
p = 1
Do While y < 7
   p = p * y
   y = y + 1
Loop
r = Exp(y * (Log(y) - 1) + 1 / (2 * (6 * y + 1 / (5 * (y + 1 / (4 * y))))))
Gamma = 2.506628274631 * r / p / Sqr(y)
End Function

The constant in the last line is sqrt(2*pi). It turned out that the results seem to have not just eight but 10 valid digits, at least when evaluated with 15-digit precision. Only in some close cases I noticed that the 10th digit may be rounded down instead of up (e.g. for x=1,7 or 7,8). I'll try this one on my '41. ;-)

Thank you very much for reminding me of Peter Luschny's excellent Gamma/factorial pages.

Dieter

[lcwright1964]

(08-19-2015 08:56 PM)Ángel Martin Wrote:  Yes I used 13-digit routines for the MCODE version in the SandMath, which uses the Lanczos formula

May I ask which coefficients? Lanczos offers some in his original paper, and I believe that rendering is quoted in Numerical Recipes verbatim. But there is also the unpublished work of Paul Godfrey, discussed by Viktor Toth on his old http://www.rskey.org website, that reports a matrix algebra method of generating the Lanczos formula coefficients. I think I observed years ago that Godfrey's reported errors in a Matlab implementation were somewhat overly optimistic, as he only seemed to test his formulae with integral and half-integral arguments.

Finally, I think anyone serious about the Lanczos approximation really needs to look at the PhD dissertation of Glen Pugh. Most of the theory and proofs are beyond me, BUT Professor Pugh did some extensive arbitrary precision computer algebra experiments and concluded that the selection of parameters as presented in the original Lanczos paper isn't the best way to go about things, and it is possible to get much less relative error than that originally reported with the same number of terms. For programmers with a mastery of elementary algebra, that part of his work is completely accessible. The Pugh thesis is freely available in a few places--simply Google something like "Pugh Lanczos Dissertation pdf" and you will find it. Dr. Pugh was also receptive to inquiries--he was very friendly with me when I wrote.

Hope this is of interest.

Les

EDIT: Top Google hit for my search:

http://web.viu.ca/pughg/phdThesis/phdThesis.pdf

Enjoy!

[Thomas Klemm]

(08-20-2015 07:49 PM)lcwright1964 Wrote:  Enjoy!

Thanks for sharing.

Kind regards
Thomas

[Ángel Martin]

(08-20-2015 07:49 PM)lcwright1964 Wrote:  

(08-19-2015 08:56 PM)Ángel Martin Wrote:  Yes I used 13-digit routines for the MCODE version in the SandMath, which uses the Lanczos formula

May I ask which coefficients? Lanczos offers some in his original paper, and I believe that rendering is quoted in Numerical Recipes verbatim. But there is also the unpublished work of Paul Godfrey, discussed by Viktor Toth on his old http://www.rskey.org website, that reports a matrix algebra method of generating the Lanczos formula coefficients. I think I observed years ago that Godfrey's reported errors in a Matlab implementation were somewhat overly optimistic, as he only seemed to test his formulae with integral and half-integral arguments.

Indeed it was the formula with the coefficients published on Viktor's page. You can also see it in the SandMath manual - the most recent version is available at TOS, search for SandMath - and an older version is available here:
http://systemyde.com/pdf/SANDMATH-44_3x3_Manual.pdf

Incidentally I also implemented the same approach on the 41Z module using complex math routines - not 13-digit accurate unfortunately but still holding very nicely given the task to perform. It is an ALL-MCODE version as well.

Thanks for sharing the link to Dr. Pugh's work, very enjoyable.

[lcwright1964]

(08-21-2015 09:48 AM)Ángel Martin Wrote:  Indeed it was the formula with the coefficients published on Viktor's page.

Do you mean the seven q's and the formula rearranged for calculators?

If that is the case, you might want to know that the selection of the free parameter (Viktor calls it g, but it is a elsewhere) is from the original Lanczos paper and is quite arbitrary. One very practical result of Pugh's work is that there is an optimal value of this free parameter for each number of terms one takes in the series (n or n + 1, depending how you index). Indeed, Pugh published the list as an appendix to his thesis. The optimal value for n = 5 is about 5.58 and is about 6.78 for n = 6. He reports a max relative error of 1.2e-10 for n = 5 and 2.7e-12 for n = 6, for 9.9 and 11.6 edd (exact decimal digits) respectively. Keep in mind that these errors cover complex arguments as well, and the performance for eligible real values tends to be somewhat better.

I can swear that I did up a Maple worksheet of this ages ago, that accepted these optimal parameters as input, solved Godfrey's matrix equations for the associated coefficients, doing it all to very generous arbitrary precision, and transformed them into the form recommended by Victor for calculator programming. But can I find the darn thing? Alas, no--that's a couple of hard drives ago.

If I do redo this I can compute and send you "improved" coefficients, in the event you ever tweak the code. It might not make much difference unless you can shorten your series by a term, but I sure would like to ponder the question.

I want to continue this discussion elsewhere in the forum as I have some broader newbie questions about the newest SandMath, Lib#4, and programming Nov-64 with this stuff (if possible).

Les

[Ángel Martin]

(08-22-2015 02:09 AM)lcwright1964 Wrote:  

(08-21-2015 09:48 AM)Ángel Martin Wrote:  Indeed it was the formula with the coefficients published on Viktor's page.

Do you mean the seven q's and the formula rearranged for calculators?

If that is the case, you might want to know that the selection of the free parameter (Viktor calls it g, but it is a elsewhere) is from the original Lanczos paper and is quite arbitrary.

I can swear that I did up a Maple worksheet of this ages ago, that accepted these optimal parameters as input, solved Godfrey's matrix equations for the associated coefficients, doing it all to very generous arbitrary precision, and transformed them into the form recommended by Victor for calculator programming. But can I find the darn thing? Alas, no--that's a couple of hard drives ago.

If I do redo this I can compute and send you "improved" coefficients, in the event you ever tweak the code. It might not make much difference unless you can shorten your series by a term, but I sure would like to ponder the question.

I want to continue this discussion elsewhere in the forum as I have some broader newbie questions about the newest SandMath, Lib#4, and programming Nov-64 with this stuff (if possible).

Les

Well, yes I tweak the code all the time (in fact just recently since Dieter found a bug in FLOOR) - but first of all, the reference I use doesn't have any free parameter that I can see. This is the link, and the coefficients are already shown in the more efficient way, I think:

http://www.rskey.org/CMS/index.php/the-library/11

I'm sure I'm missing something so will appreciate the enlightenment.

SandMath questions are always welcome, as a matter of fact I'm goinf to burn my own NoVRAM today with the updated version. It's been a while I did this, I mostly work on my CL units of course - much more comfortable and super fast.

\AM

[lcwright1964]

(08-22-2015 08:35 AM)Ángel Martin Wrote:  Well, yes I tweak the code all the time (in fact just recently since Dieter found a bug in FLOOR) - but first of all, the reference I use doesn't have any free parameter that I can see. This is the link, and the coefficients are already shown in the more efficient way, I think:

http://www.rskey.org/CMS/index.php/the-library/11

I'm sure I'm missing something so will appreciate the enlightenment.

The free, or arbitrary, parameter used in the Numerical Recipes Lanczos formula and rearranged by Viktor is embedded in the 5.5. It is little gamma in the original paper, g on Viktor's page, and r in the Pugh thesis. In the formula you are using it happens to be 5. Read further down on that very page, to the discuss of the Godfrey matrix decomposition, where Viktor identifies this g as "an arbitrary real parameter."

The key to minimizing relative error to choose this parameter--and yes it is chosen, and the coefficients then computed in the complex way according to it and the number of desired terms--just so. Both Lanczos and Godfrey, who inspired Viktor's calculations, were more or less arbitrary in the selection. Pugh did empirical tests and theoretical analysis to find the best value for a desired number of terms. Indeed, in subsequent years he extended this work up to n = 100 (the thesis only goes to n = 60) and I just found that page that he sent me in 2006. In it he gives the parameters to more digits and has improved the error estimates by using a more rigourous calculation. And indeed, it seems that in the very situation you use (n = 6, or 7 coefficients with the index starting at 0), the optimal value of g is around 6.78, not 5.

That's what I mean by "free parameter." Sorry for the confusion.

As for SandMath use, I am curious about two things:

1. What is Library 4 and why do I need it?
2. Given the switching going on in SandMath (as in the Advantage Pac), how do I correctly configure ClonixConfig to burn it into my Nov-64?

Thanks again,


Les

[lcwright1964]

Well I did find that Maple worksheet. Looks like I basically ported Viktor's C++ code for Godfrey's method and took advantage of Maple's easy to use linear algebra capabilities.

Unfortunately I did NOT expand things out and collect terms as Viktor did (and you used) so my result looks a bit different.

My preliminary conclusion is that going with Pugh's recommendation for the arbitrary parameter one can keep the relative error the same (at least 10 digits good for positive real arguments as long as one computes in 13-digit precision), with one fewer terms in the series.

Here is my result with the coefficients presented to 13 significant digits:

(0.9446965718570e-2 + 1.367185253712/z - 1.839535625335/(z+1) + .6838535775175/(z+2) - 0.6542905498252e-1/(z+3) + 0.6572071053721e-3/(z+4)) * ((z+5.081000210223)/exp(1))^(z-1/2)

Input has been adjust to compute Gamma(z), not Gamma(z+1) or z!.

The scaling by e in the second factor before exponentiation should prevent overflow when other formulas tend to for larger arguments, but I have not rigourously examined that.

In this case the arbitrary parameter recommended by Pugh is 5.58100021022. In the original paper Lanczos does not present an example with n = 5 (so 6 coefficients) for comparison.

In the version you use now you embed 7 12-digit coefficients and one 2-digit parameter--5.5 (from 5 + .5). Here you embed 6 13-digit coefficients and one 13-digit parameter. And there is one less series term. So this should be faster and smaller, at least a little.

I should work on doing what Viktor does so as to produce a polynomial divided by a product, etc. In the meantime, this is what I have.

I work in Mathematica more these days so I should see what happens there.

Les

[Ángel Martin]

(08-22-2015 08:58 AM)lcwright1964 Wrote:  As for SandMath use, I am curious about two things:

1. What is Library 4 and why do I need it?
2. Given the switching going on in SandMath (as in the Advantage Pac), how do I correctly configure ClonixConfig to burn it into my Nov-64?

Thanks for the additional explanation on the parameter's role in the formulas, I think I can follow now.

The Library#4 is a ROM that contains a library of subroutines used by several other modules of mine. This saved lots of space and made the code/operation much more consistent across all modules.

It sits in page#4, where typically there's nothing plugged - exceptions being the diagnostics ROM and a disabled printer on the HP-IL module.

The modules that use it are as follows: (I think I'm not missing anything but it's from memory)

- OSX/3 and OSX/4 - one page w/ 4 banks in it = 4 ROMs
-CLUTILS and POWERCL - one page w/ 4 banks in it = 4 ROMs
-ALPHA44
-RAMPAGE 4L
-TOOLBOX 4
-41Z Complex Numbers - 2 pages = 2 ROMS
-SANDMATH - 2 pages of 4 banks each = 8 ROMS
-SANDMatrix - 2 pages of 2 banks each = 4 ROMs
-TotalREKALL
-HP-16C Emulator (about to be released) - one page w/ 4 banks = 4 ROMs

ClonixConfig has a very flexible U/I that allows you to burn bank-switched modules like the OS/X, SandMath & SandMatrix. You simply need to select the right combination of banks/pages more or less in contiguous sequence.

The attached pic shows an example of the SandMath + SandMatrix burned on a Clonix-D with JOINED blocks. A Novo64 would use the same setup, once you disable the HEPAX RAM checkbox.

This doesn't get you there yet - as it's still missing the Library#4. So you'd probably need another Clonix41 with the Lib#4, and (highly recommended) the OS/X4. You can throw in there the TotalRekall to fill it up to 6 ROMS.

Hope this clarifies, let me know if you have more questions.

[Ángel Martin]

forgot to mention that all those modules (Lib#4 included) are available at TOS so you can do your dry runs on V41 to get a pretty good idea of how things work and what's in there.

ps. see two suggestions for Clonix41 configurations as per previous post

[lcwright1964]

Thanks for all this. I have to admit I have obviously been out of touch. The last time I paid attention to SandMath I think it consisted of just lower and upper halves, and that's about it--basically, what is still included with Diego's most recent ClonixConfig software. I took a peak at the most recent offerings on TOS, including the simpler SandMath-IV module that doesn't use Library #4, and modfile.exe tells me that there is a lot more going on with these things in the three or four years since.

And thanks for the configuration suggestion. I would probably just burn the SandMath 4x4 constituents with Library #4, and leave SandMatrix and the others be for now.

Of course, you have tempted me to order another of the Clonix family from Diego

Getting back to the original topic, I have used that very same Maple worksheet and Pugh's recommendations to give this n = 4 (i.e. 5 coefficient) version for a 10-digit FOCAL environment:

(0.3264892413e-1 + 1.188688981/z - 1.068410309/(z+1) + .1887119902/(z+2) - 0.2745021709e-2/(z+3)) * ((z+ 3.840881909)/exp(1))^(z-1/2)

In an ideal situation (i.e., a couple of guard digits kicking around) this is supposed to give between 9.5 and 11 edd in the z = 0 to 70 range, but of course we don't have guard digits so 8 to 10 is about the best I think it would do. Given this I think JMB's implementation of the Stirling series in GAM+, with its simpler coefficients, will still do the same or better. That said, I would like to see how this fares in FOCAL code for the HP41. It would even work for the 67/97 where one could store the program on one side of the card (it would be very short) and the six 10-digit constants as register data on the other side . Without Pochhammer shifting for small arguments this should be faster than a Stirling approach, and of course the Spouge approach uses more terms and calculates more.

If I actually turn this into program steps for the 41 or the 67/97 I will post.

Les

P.S. I notice that the coefficients (q0 to q6) from Viktor's page are only presented to 12 digits, not 13. Is that what you are using, or did he or you recompute 13-digit versions for use in your MCODE? That does seem possible given that the original Lanczos coefficients, p0 to p6, are given just above them to up to 16 digits.