From 89664d2921325408b4dedda0a526bab5e1f9068e Mon Sep 17 00:00:00 2001
From: Stephen Heumann <stephenheumann@gmail.com>
Date: Sat, 24 Dec 2022 21:59:52 -0600
Subject: [PATCH] Slightly improve tgamma calculation for x < 8.

Previously, 1-4 low-order bits of the input value were essentially ignored when calculating the numerator, but used to some degree when calculating the denominator. This would lead to the calculated tgamma values decreasing slightly over the range of several consecutive input values (when they should increase). Now, the low-order bits of the input value are effectively just rounded away. This should give slightly more accurate results, and greatly reduces the frequency of cases where consecutive output values go in the wrong direction.
---
 math2.asm | 63 +++++++++++++++++++++++++++++++++++++++----------------
 1 file changed, 45 insertions(+), 18 deletions(-)

diff --git a/math2.asm b/math2.asm
index 05179e3..aec6b73 100644
--- a/math2.asm
+++ b/math2.asm
@@ -3164,30 +3164,57 @@ lp1      lda   one,y
          dey
          bpl   lp1
 
-lp2      ph4   #cutoff                  while x < 9.375 then
-         tdc
-         clc
-         adc   #x
-         pea   0
-         pha
-         FCMPS
-         bvc   stirling
-         
-         tdc                              denom := denom * x
-         clc
-         adc   #x
-         pea   0
-         pha
-         ph4   #denom
-         FMULX
+         ph4   #cutoff
+         ph4   #z
+         FS2X
 
-         ph4   #one                       x := x + 1
+         tdc                            z := 10.375 - x
+         clc
+         adc   #x
+         pea   0
+         pha
+         ph4   #z
+         FSUBX
+         
+         lda   z+8                      if z < 0 (or NAN)
+         jmi   stirling                   just use Stirling series approx.
+         cmp   #$7fff
+         jeq   stirling
+         
+         ph4   #z                       truncate z to integer
+         FTINTX
+         
+         ph4   #z                       x := x + z
          tdc
          clc
          adc   #x
          pea   0
          pha
          FADDX
+         
+         tdc                            y := x
+         clc
+         adc   #x
+         pea   0
+         pha
+         ph4   #y
+         FX2X
+         
+         ph4   #z                       repeat z times :
+         ph4   #z
+         FX2I
+         
+lp2      dec   z
+         bmi   stirling
+
+         ph4   #one                       y := y - 1
+         ph4   #y
+         FSUBX
+         
+         ph4   #y                         denom := denom * y
+         ph4   #denom
+         FMULX
+         
          bra   lp2
          
 * For x >= 9.375, calculate Gamma(x) using a Stirling series approximation
@@ -3296,7 +3323,7 @@ done     FPROCEXIT                      restore env & raise any new exceptions
          sta   x
          creturn 4:x
 
-cutoff   dc    f'9.375'                 cutoff for using Stirling approximation
+cutoff   dc    f'10.375'                cutoff for Stirling approximation (+1)
 
 one_half dc    f'0.5'
 two      dc    i2'2'