/*********************************************************************** Test a higher-precision Pochhammer function against an fp_t version, or against a recursive fp_t version with claimed lower average error, as given by Johansson:2023:APC, page 6, in elefunt.bib. Usage: cc [-DMAXTEST=nnn] [-DMAXN=nnn] [-DMAXX=d.dd] [-DTEST1 | -DTEST2] \ [-DUSE_FLOAT] poch.c && ./a.out > poch.dat TEST1 selects the standard version, while TEST2 selects the divide-and-conquer version. Testing is over the range n in [2, MAXN], and x uniformly sampled on [0, MAXX]. The output file contains triples n x ulps, and is suitable for input to gnuplot, with coomands like this: % gnuplot set yrange [-10:10] set term x11 1; plot "poch.dat.1" using ($1):($3) set term x11 2; plot "poch.dat.2" using ($1):($3) set term x11 3; plot "poch.dat.1" using ($2):($3) set term x11 4; plot "poch.dat.2" using ($2):($3) set term x11 5; plot "poch.dat.3" using ($1):($3) set term x11 6; plot "poch.dat.4" using ($1):($3) set term x11 7; plot "poch.dat.3" using ($2):($3) set term x11 8; plot "poch.dat.4" using ($2):($3) set term x11 9; plot "poch.dat.5" using ($1):($3) set term x11 10; plot "poch.dat.6" using ($1):($3) set term x11 11; plot "poch.dat.5" using ($2):($3) set term x11 12; plot "poch.dat.6" using ($2):($3) set term x11 13; plot "poch.dat.7" using ($1):($3) set term x11 14; plot "poch.dat.8" using ($1):($3) set term x11 15; plot "poch.dat.7" using ($2):($3) set term x11 16; plot "poch.dat.8" using ($2):($3) Plots for x over [0, 1/16], [0,1], [0,8] and [0,1024] show that the accuracy of the divide-and-conquer version is INFERIOR to the standard method, despite a statement to the contrary in the cited paper. Errors increase as MAXX decreases. The reason for that is that, for |x| >= 1, forming x + m in the Pochhammer product is an EXACT operation, but for |x| < 1, each x + m value introduces a further rounding error into the product, and that error increases with m. The data files were made like this with binary80 on x86_64; plots for binary128 on s/390x are similar: % cc -O3 -DMAXN=1000 -DTEST2 poch.c && time ./a.out > poch.dat.2 % cc -O3 -DMAXX=8.0 -DMAXN=1000 -DTEST1 poch.c && time ./a.out > poch.dat.3 % cc -O3 -DMAXX=8.0 -DMAXN=1000 -DTEST2 poch.c && time ./a.out > poch.dat.4 % cc -O3 -DMAXX=1.0 -DMAXN=1000 -DTEST1 poch.c && time ./a.out > poch.dat.5 % cc -O3 -DMAXX=1.0 -DMAXN=1000 -DTEST2 poch.c && time ./a.out > poch.dat.6 % cc -O3 -DMAXX=0.0625 -DMAXN=1000 -DTEST1 poch.c && time ./a.out > poch.dat.7 % cc -O3 -DMAXX=0.0625 -DMAXN=1000 -DTEST2 poch.c && time ./a.out > poch.dat.8 [06-Jul-2026] ***********************************************************************/ #include #include #include #include #include #if !defined(TEST1) && !defined(TEST2) #define TEST1 #endif #if !defined(MAXTEST) #define MAXTEST 250 #endif /* !defined(MAXTEST) */ #if !defined(MAXN) #define MAXN 251 #endif /* !defined(MAXN) */ #if !defined(MAXX) #define MAXX 1024.0 #endif /* !defined(MAXX) */ #if defined(USE_FLOAT) typedef float fp_t; typedef long double hp_t; const fp_t EPS = FLT_EPSILON; static fp_t ONE = 1.0F; static fp_t ZERO = 0.0F; #else typedef double fp_t; typedef long double hp_t; const fp_t EPS = DBL_EPSILON; static fp_t ONE = 1.0; static fp_t ZERO = 0.0; #endif #define iseven(n) ((n & 1) == 0) #define qnan nan static fp_t poch_std(fp_t x, int n) { /* ** Pochhammer symbol, with standard mathematical notation (x)_n, ** defined by: ** ** (x)_n = (x + 0) * (x + 1) * (x + 2) * ... * (x + n - 1) ** ** Special cases: (+/-0)_n = +/-0, (x)_0 = 1, and (x)_1 = x ** ** Identities: (1)_n = n! and (2)_n = (n + 1)! */ fp_t result; if (isnan(x)) /* generate run-time NaN, preserve payload */ { errno = ERANGE; (void)nan(""); result = x; } else if (n < 0) { errno = ERANGE; result = qnan("0xfac50000"); } else if (x == ZERO) result = x; /* NB: preserve sign of zero! */ else if (n == 0) result = ONE; else if (n == 1) result = x; else if (isinf(x)) result = (x > (fp_t)0.0) ? x : (iseven(n) ? -x : x); else { int k; fp_t prod; prod = ONE; for (k = 0; k < n; ++k) prod *= (x + (fp_t)k); result = prod; } return (result); } static hp_t poch_hp(fp_t x, int n) { /* ** Pochhammer symbol, with standard mathematical notation (x)_n, ** defined by: ** ** (x)_n = (x + 0) * (x + 1) * (x + 2) * ... * (x + n - 1) ** ** Special cases: (+/-0)_n = +/-0, (x)_0 = 1, and (x)_1 = x ** ** Identities: (1)_n = n! and (2)_n = (n + 1)! */ hp_t result; if (isnan(x)) /* generate run-time NaN, preserve payload */ { errno = ERANGE; (void)nan(""); result = x; } else if (n < 0) { errno = ERANGE; result = qnan("0xfac50000"); } else if (x == ZERO) result = x; /* NB: preserve sign of zero! */ else if (n == 0) result = ONE; else if (n == 1) result = x; else if (isinf(x)) result = (x > (fp_t)0.0) ? x : (iseven(n) ? -x : x); else { int k; hp_t prod, xx; prod = (hp_t)ONE; xx = (hp_t)x; for (k = 0; k < n; ++k) prod *= (xx + (hp_t)k); result = prod; } return (result); } #if defined(TEST2) static const int __NPOCH_DAQ__ = 10; static fp_t poch_daq(fp_t x, int n) { fp_t result; if (n < __NPOCH_DAQ__) result = poch_std(x, n); else { int m; m = n / 2; result = poch_daq(x, m) * poch_daq(x + m, n - m); } return (result); } #endif /* defined(TEST2) */ int main(void) { fp_t a, ulp, ulpmax, ulpsum, x; hp_t e; int k, n, nsum; srand48(0xcafefeed); nsum = 0; ulpsum = ZERO; ulpmax = ZERO; (void)fprintf(stderr, "hp_t is binary%d\n", (LDBL_MANT_DIG == 64) ? 80 : (int)sizeof(hp_t) * 8); #if defined(TEST1) (void)fprintf(stderr, "Test of poch_std() in binary%d\n", (int)sizeof(fp_t) * 8); for (k = 0; k < MAXTEST; ++k) { x = drand48() * MAXX; for (n = 2; n < MAXN; ++n) { e = poch_hp(x, n); a = poch_std(x, n); if (isinf((fp_t)e)) break; ulp = (fp_t)(((hp_t)a - e) / e) / EPS; if (isinf(ulp)) continue; ++nsum; ulpsum += fabs(ulp); if (fabs(ulpmax) < fabs(ulp)) ulpmax = ulp; (void)printf("%7d\t%15.7f\t%15.7f\n", n, x, ulp); } } (void)fprintf(stderr, "Average ulps for poch_std() = %.7f\n", (double)ulpsum / (double)nsum); (void)fprintf(stderr, "Maximum ulps for poch_std() = %.7f\n", (double)ulpmax); (void)fprintf(stderr, "Sample count for poch_std() = %7d\n", nsum); #elif defined(TEST2) (void)fprintf(stderr, "Test of poch_daq() in binary%d\n", (int)sizeof(fp_t) * 8); for (k = 0; k < MAXTEST; ++k) { x = drand48() * MAXX; for (n = 2; n < MAXN; ++n) { e = poch_hp(x, n); a = poch_daq(x, n); if (isinf((fp_t)e)) break; ulp = (fp_t)(((hp_t)a - e) / e) / EPS; if (isinf(ulp)) continue; ++nsum; ulpsum += fabs(ulp); if (fabs(ulpmax) < fabs(ulp)) ulpmax = ulp; (void)printf("%7d\t%15.7f\t%15.7f\n", n, x, ulp); } } (void)fprintf(stderr, "Average ulps for poch_daq() = %.7f\n", (double)ulpsum / (double)nsum); (void)fprintf(stderr, "Maximum ulps for poch_daq() = %.7f\n", (double)ulpmax); (void)fprintf(stderr, "Sample count for poch_daq() = %7d\n", nsum); #endif /* defined(TEST1) */ return (EXIT_SUCCESS); }