![[ICO]](/icons/blank.gif){kind=icon}
![[PARENTDIR]](/icons/back.gif){kind=icon}
![[TXT]](/icons/text.gif){kind=icon}
![[ ]](/icons/unknown.gif){kind=icon}
![[DIR]](/icons/folder.gif){kind=icon}
| Name | Last modified | Size |
|---|---|---|---|
| Parent Directory | - | |
| README.html | 2026-10-30 13:24 | 80K |
| README.sok | 2026-08-17 10:19 | 2.6K |
| gnuplot/ | 2026-08-11 12:14 | - |
| images/ | 2026-08-11 12:26 | - |
| src/ | 2026-08-17 10:13 | - |
| zip/ | 2026-08-17 10:18 | - |
Original version:
Sat Jul 11 14:20:20 2026
Last updates:
…
Fri Oct 30 13:24:22 2026
Current time in Salt Lake City, Utah, USA:
Saturday, 14-Nov-2026 11:20:52 MST
The Pochhammer symbol, also called the rising factorial, is an n-term product defined for complex and real x, and integer n ≥ 0, by the notation
(x)n = x * (x + 1) * (x + 2) * … * (x + n - 1),
with these special relations:
(−0)n = −0,
(+0)n = +0,
(1)n = n!,
(2)n = (n + 1)!,
(x)0 = 1,
(x)1 = x,
(x)n = Γ(x + n) / Γ(x),
(x)n = (x)k (x + k)n - k,
binomial(m, n) = (m - n)! / n! = (m - n + 1)n / n!.
Notice that the relation to the Gamma function allows us to recur upward from the small-argument region where sums or polynomials might be used to compute Γ(x), or downward from the large argument region Γ(x + n), where a faster converging asymptotic expansion might be usable.
Once the special cases have been dealt with, the Pochhammer symbol value is apparently simple to compute, such as with C code like this:
#include <errno.h>
#define iseven(n) ((n & 1) == 0)
static const fp_t ONE = (fp_t)1.;
static const fp_t ZERO = (fp_t)0.;
fp_t /* fp_t = float, double, long double, … */
poch(fp_t x, int n)
{
fp_t result;
if (isnan(x))
{ /* generate run-time NaN, and preserve payload */
errno = ERANGE;
(void)qnan("");
result = x;
}
else if (n < 0)
{
errno = ERANGE;
result = qnan("");
}
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 > ZERO) ? x : (iseven(n) ? -x : x);
else
{
int k;
result = x;
for (k = 1; k < n; ++k)
result *= (x + (fp_t)k);
}
return (result);
}
The initial six checks are essential to guarantee correct behavior with respect to the mathematical properties of the Pochhammer symbol. In the rest of this document, we investigate alternatives in the final else block with the loop that computes the n−1 additions and multiplications to form the Pochhammer product.
A major goal in post-2000 numerical analysis has been to produce provably correctly rounded floating-point results, differing by no more than a half-unit in the last place (a half-ulp) from the infinitely precise result. That can, for example, guarantee computational monotonicity when the mathematical function is known to be monotonic in a given region.
That also matches the default behavior of IEEE 754 floating-point arithmetic which, in the default rounding mode of round to nearest, with ties to even, produces correctly rounded results for add, subtract, multiply, divide, square root, reciprocal square root, and fused multiply-add (FMA) operations.
With guaranteed correct rounding, it is possible to build up more complicated operations, such as complex arithmetic, multiple sums or products, dot products, and matrix multiplications, that are as accurate as possible without having to resort to higher intermediate working precision.
If rounding to within a half-ulp is not possible, then one should strive for faithful rounding, differing by no more than one ulp, such as might be obtained by truncation of an exact result.
From the definition of the Pochhammer symbol, we can immediately see a significant problem: the n-term computation requires n−1 additions, and n−1 multiplications, so in IEEE 754 arithmetic, the worst case error could reach almost n ulps in the default rounding mode, and 2n ulps in other modes. For example, in the computation of the Gamma function in the IEEE 754 binary128 format, overflow is reached just above x = 1755. In the decimal128 format, that limit is about x = 2124. Thus, we might need to compute with n values of a few thousand, implying a loss of about 4 decimal digits from a result that we would like to be correct to about 34 decimal digits.
One obvious improvement is to compute the product in the next higher available precision, but that is not possible if we are already at the highest precision, or have deficient C compiler implementations, or programming languages like C#, Go, Java, and Rust, that fail to provide anything beyond the IEEE 754 binary32 and binary64 formats, even when the host computer may offer a longer type such as binary80 (11 more significand bits) in hardware, or binary128 (60 more significand bits) in hardware or software. Programmers should be extremely annoyed by a compiler or a programming language that denies them access to hardware that has already been paid for, and is frequently needed to improve computational accuracy.
Because poch(x,n) is a bivariate function, plotting its errors would require a 3D point plot. Instead, we choose one coordinate as the horizontal axis of a 2D point plot, and project all values of the other coordinate onto that axis. The resulting error plots for computation of the Pochhammer symbol in working precision look like this:
The errors grow roughly linearly with n, but have similar variation for any x. However, the error growth with increasing n is not as severe as our discussion above predicted: our measurements show an average error of 1.17 ulps, with a maximum error of −8.23 ulps.
In those, and following, plots, three horizontal lines in blue mark the −1/2, 0, and +1/2 ulps boundaries. A correctly rounded result would always lie in that region. Here, we needed a longer vertical scale to contain the measured errors, but later plots may reduce the error range.
The plot titles are suffixed with −n and −x to indicate the horizontal axis name, and are prefixed with the name (here, poch) of the basename of the C source file that generated the plot data. All of the tests, or subsets thereof, are run by a single Unix shell script that compiles, links, and runs the test programs, captures the output data, creates parametrized input files for the gnuplot utility, and then runs that tool to produce plots in bitmap PNG (Portable Network Graphics) and typeset EPS (Encapsulated PostScript) formats. The former allow fast Web browser display, and the latter are suitable for inclusion in typeset documents that preserve image quality under magnification. The final section of this document supplies links to all of the needed files, to allow others to reproduce this work, and experiment with code and data variations.
The C programs for the various implementations of the Pochhammer symbol computation allow compile-time specification of the maximum values of n and x, and the number of data points to be generated. The n values run from 0 to the specified maximum, while the x values are sampled via two different distributions from a POSIX standard pseudo-random number generator that guarantees identical sequences across systems. The default (n,x) pairs produce Pochhammer symbol values from zero up to the function overflow limit to provide thorough test coverage for how the functions behave in practical applications. Experiments show that increasing the point counts beyond the defaults shown in this document do not noticeably change the plot appearances, but just make them slower to display.
Testing with uniformly distributed x values is insufficient, and deceptive, as William Cody and William Waite pointed out in their influential 1980 book Software Manual for the Elementary Functions (ISBN 0-13-822064-6). The reason is that floating-point numbers are not themselves uniformly distributed: they are so when they share the same exponent, but not when viewed over their entire representable range. Cody and Waite therefore used a logarithmic distribution derived from a uniform one, with code like these three functions:
fp_t
Rand(void)
{
return (drand48()); /* uniformly distributed */
}
fp_t
randl(fp_t x)
{
return (EXP(x * Rand()));
}
fp_t
randlab(fp_t a, fp_t b)
{
return (a * randl(LOG(b / a))); /* log distributed on [a, b] */
}
To avoid overflow in forming b/a (avoidable if we instead wrote slower code LOG(b) − LOG(a)), it is necessary to split the representable range into two parts, which our test programs named poch*randlab*.c do with code like this:
for (k = 0; k < MAXTEST; ++k)
{
x = (k % 2) ? randlab(MINNORMAL, ONE) : randlab(ONE, MAXX);
⋮
}
Ideally, we should also use x values from randlab(MINSUBNORMAL, MINNORMAL), but we ignore that range because its numbers have reduced precision.
Here is a view of the errors found from logarithmically distributed x values, where the tests report an average error of 4.16 ulps and a maximum error of 48.35 ulps:
The vertical scale has been significantly expanded to match the new error measurements, and the plot with the x axis shows that the largest errors happen with small x values: that is expected because forming x + k loses one or more bits in floating-point arithmetic when |x| < 1.
In the companion plot with the n axis, errors drop for large-n and tiny x values, as confirmed by interactive manipulation of a 3D point plot.
We revisit the accuracy issue in a later section on worst-case arguments for the Pochhammer symbol computation.
When x is large, many of the Pochhammer terms (x + k) lie in the same binade where the addition has zero rounding error. However, it is expected that each multiplication in the product contributes a rounding error. It may therefore improve accuracy by computing the product with error compensation, where a double-length product x * y is represented as a high part xy = fl(x * y), plus a small correction given by err = fma(x, y, −xy). Here, fl() is the common designation of a floating-point result in working precision.
In a later section, we investigate inclusion of error compensation for both addition and multiplication.
The Pochhammer loop code with error compensation of the product looks like this:
{
fp_t err, prod;
int k;
err = ZERO;
prod = x;
for (k = 1; k < n; ++k)
{
fp_t e, f, p;
f = x + (fp_t)k; /* fl(new factor) */
p = prod * f; /* fl(new product) */
e = FMA(prod, f, -p); /* error in new product */
prod = p; /* current product */
err = FMA(err, f, e); /* current error */
}
result = prod + err;
}
A test of that idea finds that errors lie in [−1/2, +1/2], with an average error of about 0.18 ulps. Here are the error plots:
A test with logarithmically distributed x values has an average error of 3.11 ulps, and a maximum error of 45.67 ulps:
Good-to-the-last bit computation of the Pochhammer symbol may be possible with the help of a multiple-precision arithmetic package. We investigate several alternatives in this section.
With the MPFR library, -lmpfr, the loop block can be written like this:
#include <mpfr.h>
⋮
{
int k;
mpfr_set_default_prec(2 * LDBL_MANT_DIG);
mpfr_set_ld(x_plus_k, (hp_t)x, MPFR_RNDN);
mpfr_set_ld(u, (hp_t)x, MPFR_RNDN);
for (k = 1; k < n; ++k)
{
mpfr_add_si(t, x_plus_k, (long int)k, MPFR_RNDN);
mpfr_mul(u, u, t, MPFR_RNDN);
}
result = (fp_t)mpfr_get_ld(u, MPFR_RNDN);
}
That sample code uses a working precision of 106, 128, or 226 bits, which is almost certainly larger than needed. French researchers reported in late 2025 that the worst cases for rounding of all single-argument binary64 functions in the Standard C math library, -lm, have now been found. The hardest to round correctly is sin(0x1.6ac5b262ca1ffp+849), which has 68 consecutive identical bits after the round bit. The Pochhammer symbol is far simpler, and requires no transcendental functions, so eliminating a 2 n-ulps error for n < 2125 should need only about 12 additional bits of precision, plus a few more for hard-to-round cases: 20 extra bits beyond working precision is likely to suffice.
Point plots of the errors in Pochhammer values computed with the MPFR routines show visually uniform errors in ulps in [−1/2, +1/2], with an average error of about 0.18 ulps. They look like this with sampling from a uniform distribution on [0, MAXX]:
A test with logarithmically distributed x values has an average error of 0.18 ulps, and a maximum error of −0.499 ulps:
The GNU Multiple Precision library, -lgmp, can be used for the Pochhammer loop, with code like this:
#include <gmp.h>
⋮
{
const int gmp_prec = 2 * LDBL_MANT_DIG;
int k;
mpf_t t, u, xx;
mpf_set_default_prec(gmp_prec);
mpf_init(t);
mpf_init_set_d(u, x);
mpf_init_set_d(xx, x);
for (k = 1; k < n; ++k)
{
mpf_add_ui(t, xx, (unsigned long int)k);
mpf_mul(u, u, t);
}
/* NB: GMP provides no float or long double types! */
result = mpf_get_d(u);
mpf_clear(t);
mpf_clear(u);
mpf_clear(xx);
}
The GMP library is more limited than the MPFR library, in that the only native floating-point type supported is the binary64 format supplied by the C type double. Error plots show a distinct difference from the MPFR version: errors lie in the interval [−1, 0], indicating truncating arithmetic. Increasing the precision to even several hundred bits does not change the average error from its default value of about 0.36 ulps, and the worst-case error remains near −1 ulp. with user-specified rounding direction.
The documentation of mpf_get_d() notes that it returns a truncated result, and there does not appear to be an alternate GMP function for conversion to a floating-point value with user-specified rounding direction. Presumably, by delving into GMP internals, this could readily be done by the user, but the MPFR alternative seems a better choice. Indeed, the GMP Info node documentation in Chapter 7, Floating-point Functions, reports
New projects should consider using the GMP extension library MPFR (http://mpfr.org) instead. MPFR provides well-defined precision and accurate rounding, and thereby naturally extends IEEE P754.
In a binary32 version that uses a cast of float on the return value of mpf_get_d(), the GMP interface produces errors in ulps in [−1/2, +1/2], with an average error of about 0.18 ulps.
A test with logarithmically distributed x values has an average error of 0.36 ulps, and a maximum error of −1.00 ulps:
Recent versions of FLINT (Fast Library for Number Theory) have absorbed several earlier libraries for multiple-precision arithmetic, among them, Fredrik Johansson's Arb library for interval arithmetic. The -lflint library is huge, with almost 14_400 function in mid-2026. We only require six of its functions to implement the Pochhammer loop, with code that looks like this:
#include <flint/arf.h>
⋮
{
arf_t t, u, xx;
const int arf_prec = 2 * LDBL_MANT_DIG;
int k;
arf_init(t);
arf_init(u);
arf_init(xx);
arf_set_d(u, x);
arf_set_d(xx, x);
for (k = 1; k < n; ++k)
{
arf_add_ui(t, xx, (unsigned long int)k,
arf_prec, ARF_RND_NEAR);
arf_mul(u, u, t, arf_prec, ARF_RND_NEAR);
}
/* NB: FLINT provides no float or long double types! */
result = arf_get_d(u, ARF_RND_NEAR);
arf_clear(t);
arf_clear(u);
arf_clear(xx);
}
Error plots are visually similar to those for the MPFR version: errors in ulps are in [−0.499, +0.499], with an average error of about 0.18 ulps.
As with the GMP library, the only native floating-point type supported by FLINT for conversions to and from multiple-precision results is the C type double, for the IEEE 754 binary64 format. In this author's view, that is an extremely serious deficiency of those libraries. Of course, suitable get and set conversion functions could be created in user code, but they would then necessarily have to use library implementation details that might change in a future library version.
With uniform sampling of x, errors in ulps are in [−0.499, +0.499], with an average error of about 0.18 ulps:
With logarithmic sampling of x, errors in ulps are in [−0.499, +0.499], with an average error of about 0.18 ulps:
One of the motivations for the development of this author's fixed-point library, -lfx, is that it provides a sufficiently wide format that addition, subtraction, and multiplication of two values in any of the Standard C floating-point types are always exact, unlike floating-point arithmetic, where all of the basic operations may be subject to rounding. Rounding of values of type fx_t to one of the standard floating-point types is then a precise operation, and several different rounding modes, including the traditional four modes of IEEE 754 arithmetic, are available. The default mode is, of course, round to nearest, with ties to even, just as in modern floating-point arithmetic.
That library can be used to implement the Pochhammer loop code like this:
#include <fx.h>
⋮
{
int k;
fx_t t, u, xx;
fxone(&t);
fp64tofx(&u,x);
fp64tofx(&xx,x);
for (k = 1; k < n; ++k)
{
fxaddint(&t, &xx, (int64_t)k);
fxmul(&u, &u, &t);
}
result = fxtofp64(&u);
}
Error plots are visually similar to those for the MPFR version: with uniform sampling of x, errors in ulps are in [−1/2, +1/2], with an average error of about 0.18 ulps.
With logarithmic sampling, errors in ulps are in [−0.501, +0.501], with an average error of about 0.18 ulps.
Other multiple precision arithmetic packages usually define a type that is a one-element array containing a complicated data structure. Because the C language family passes an array argument as a pointer to its first element, no ampersand operator is then needed to take the address of the argument. However, C does not support functions returning arrays, so functions in those libraries are usually declared with type void.
By contrast, the FX library defines the fx_t type as a scalar data structure, and calls to its functions therefore require explicit pointer arguments. However, the functions uniformly return their first (result) argument, so FX function calls can be nested. That is often convenient for coding. For example, our four-line loop body could be reduced to the single line fxmul(&u, &u, fxaddint(&t, &xx, (int64_t)k)), preferably prefixed with (void) to indicate that the outer return value is discarded.
If general multiple-precision arithmetic libraries prove too slow for accurate computation of the Pochhammer symbol, there is an old technique that appears to have been first suggested independently in 1965 by William Kahan (doi: 10.1145/363707.363723) and Ole Møller (doi: 10.1007/BF01975722), and further improved in 1971 by Theodorus Jozef (Dirk) Dekker (1927—2021) (doi: 10.1007/BF01397083). Their idea was to represent numbers as sums of nonoverlapping high and low parts, and then to carry out arithmetic operations with those pairs, thereby simulating arithmetic of roughly twice the working precision. Dekker showed how to implement those operations in reasonably portable high-level programming language code.
The phrase doubled-double has often been used in the numerical literature for this technique, but that is a poor designation, because it has nothing to do with double precision data. This author's book, The Mathematical-Function Computation Handbook (referenced later), devotes two long chapters to what he prefers to call pair-precision arithmetic. The book describes numerous pitfalls with that arithmetic, which are mainly due to limited exponent range, and aberrant behavior near underflow and overflow limits.
The deficient floating-point designs in computers up to the late 1970s made Dekker's techniques tricky to implement, and usually, some machine dependence had to be introduced into his code.
Fortunately, the introduction of IEEE 754 arithmetic, developed from 1977 to 1980, and first implemented in hardware in the Intel 8087 coprocessor chip that reached the market in the summer of 1980, has repaired those blemishes, and become the worldwide standard for computer arithmetic for more than four decades. Consequently, it is now possible to implement pair-precision arithmetic with perfect portability, and without having to delve into the machine representation of floating-point numbers using assembly code instructions, or bit operations on data overlays.
Pair-precision arithmetic has never been standardized in programming languages, and is likely to be unfamiliar to most programmers, so we present it here in several subsections, before applying it to the Pochhammer symbol computation.
The MathCW library, -lmcw, that accompanies this author's book is broadly targeted at multiple historical and current floating-point designs, including those for binary, decimal, and hexadecimal arithmetic, and for use from multiple programming languages. Consequently, it had to recognize that the only common interlanguage argument-passing convention for data structures is the array, so it uses a data type for pairs defined as a two-element array.
Prior to the first ISO C Standard of 1989, one could pass struct arguments to functions, but they could not be returned as values; instead, functions had to return the result via an argument that is a pointer to a struct.
Elimination of that restriction in the 1989 Standard means that we can now write code for pair-precision arithmetic more compactly, and more conveniently, and without an extra argument for the function result.
Here is a four-line preamble of what we need:
typedef double fp_t;
typedef struct { fp_t hi; fp_t lo; } pair_t;
#define FABS(x) fabs(x)
#define FMA(x,y,z) fma((x), (y), (z))
That block can be readily adjusted to support other floating-point types, and the rest of our code is then completely general for every such type.
Just as with complex arithmetic, it is convenient to have a function to construct a pair from its parts. Here is code to do that:
pair_t
makepair(fp_t hi, fp_t lo)
{
pair_t result;
result.hi = hi;
result.lo = lo;
return (result);
}
It is the user's responsibility to ensure that |lo| < |hi| * xxx_EPSILON.
A product of two n-bit numbers requires at most 2n bits, as long as the result does not overflow, or underflow, or descend into the region of IEEE 754 subnormal numbers. That italicized clause is a critical condition for accurate, and correct, operation of pair-precision arithmetic!
Computing such a product was challenging on older designs, but the introduction of the fused multiply-add instruction in the IBM POWER CPU in early 1990, and standardized in later versions of the IEEE 754 Standard, makes it easy, with code like this:
pair_t
twomul(fp_t x, fp_t y)
{
fp_t xy;
xy = x * y;
return (makepair(xy, FMA(x, y, -xy)));
}
In a hardware design for higher-precision arithmetic, there would normally be a wider exponent range, and a single significand with somewhat more than the twice the number of digits in normal precision. That can ensure exact representation of products, and holds, except at the highest precision, for all IEEE 754 binary and decimal formats. The binary formats have 24, 53, 64, and 113 significand bits, and the decimal formats have 7, 16, and 34 significand decimal digits. A future 256-bit format proposed and supported in this author's book will have 70 decimal digits, following a uniform 2n + 2 digit significand growth relation.
Unlike a hardware design, in pair-precision arithmetic the exponents of the high and low parts are not linked: in the extreme, such a number could have a high part of MAXNORMAL, and a low part of MINSUBNORMAL, with a large number of unknown, and omitted, digits in between. That is another anomaly of pair-precision arithmetic that can lead to surprises, and can even derail computations.
Here is the code for splitting a number exactly into high and low parts, an operation that has often been called the 2sum in recent research publications:
pair_t
twosum(fp_t x, fp_t y)
{
volatile fp_t err, sum, tmp;
if (FABS(x) >= FABS(y))
{
sum = x + y;
tmp = sum - x;
err = y - tmp;
}
else
{
sum = x + y;
tmp = sum - y;
err = x - tmp;
}
return (makepair(sum, err));
}
The arithmetic in twosum() must be done in the IEEE 754 default rounding mode. We could guarantee that with outer calls to fegetround() and fesetround(), but have not done that here: a robust library version would have to do so.
A subtle, and critical, point in this code is the volatile qualifier: without it, a compiler code optimizer might erroneously reduce the err value to zero.
There is a longer variant of twosum that avoids a conditional branch and two function calls, but requires six numeric statements instead of three. It was published in the cited 1965 article by Ole Møller, and further clarified in Donald Knuth's three editions of The Art of Computer Programming: Volume 2: Seminumerical Algorithms (1969, 1981, 1998). It is also described in the Handbook of Floating-Point Arithmetic (referenced later). It can replace our first version, because both produce identical results. The longer version is also needed for any base greater than 3, in particular, for decimal arithmetic. Here is its code:
pair_t
twosum(fp_t x, fp_t y)
{
fp_t err, sum;
volatile fp_t t1, t2, t3, t4;
sum = x + y;
t1 = sum - y;
t2 = sum - t1;
t3 = x - t1;
t4 = y - t2;
err = t3 + t4;
return (makepair(sum, err));
}
The four temporary variables require the volatile qualifier, but the err and sum variables do not. As in the previous version, IEEE 754 default rounding is required.
In both versions of twosum(), the arguments should remain away from the overflow, underflow, and subnormal regions, to avoid generating NaN values, or deceptively lowered precision of sum + err.
When we implement the operations of the following subsections, we need one more critical function to provide accurate summation of an arbitrary number of scalar values. Its code looks like this:
pair_t
vsum(fp_t v[/* n */], int n)
{
fp_t err, sum;
int k;
pair_t the_sum;
err = ZERO;
sum = ZERO;
for (k = 0; k < n; ++k)
{
the_sum = twosum(sum, v[k]);
sum = the_sum.hi;
err += the_sum.lo;
}
return (makepair(sum, err));
}
Notice that this code cannot deal properly with an exact sum that spans many digits: it can only produce the closest pair to that value. That problem can only be solved by higher precision computation, such as with the four libraries described earlier.
As an example, consider binary64 arithmetic with x = 1 and e = 0x1p-54. The square (x + e)2 is exactly 1 + 0x2p-54 + 0x1p-108, or +0x1.0000_0000_0000_0800_0000_0000_001p+0, which rounds to +0x1.0000_0000_0000_1p+0. However, if x + e were first evaluated as a binary64 value, it would reduce to 1. Its square would be the same, and incorrect by one ulp.
The factors in the Pochhammer symbol product are each represented by a pair, and we need to multiply and accumulate them componentwise, with this code:
pair_t
pairmul(pair_t x, pair_t y)
{
pair_t the_prod;
fp_t v[8];
the_prod = twomul(x.hi, y.hi);
v[0] = the_prod.hi;
v[1] = the_prod.lo;
the_prod = twomul(x.hi, y.lo);
v[2] = the_prod.hi;
v[3] = the_prod.lo;
the_prod = twomul(x.lo, y.hi);
v[4] = the_prod.hi;
v[5] = the_prod.lo;
the_prod = twomul(x.lo, y.lo);
v[6] = the_prod.hi;
v[7] = the_prod.lo;
return (vsum(v, 8));
}
The last helper function that we need sums two pairs, and looks like this:
pair_t
pairsum(pair_t x, pair_t y)
{
fp_t v[4];
v[0] = x.hi;
v[1] = y.hi;
v[2] = x.lo;
v[3] = y.lo;
return (vsum(v, 4));
}
We have now described all of the support functions that we need, and they could be moved into a common library, with variant name suffixes for each floating-point format. At this point, they are all black boxes, and the Pochhammer loop code looks similar to that for the four multiple-precision variants given earlier, but is shorter:
{
prod = makepair(x, ZERO);
for (k = 1; k < n; ++k)
prod = pairmul(prod, twosum(x, (fp_t)k));
/* NB: intermediate overflow produces NaN parts,
NOT Infinity! */
if (isnan(prod.hi) || isnan(prod.lo))
prod = makepair(infty(), ZERO); /* NB: this assumes x >= 0 */
result = prod.hi;
}
The last statement discards the tiny low part, but in some applications, it might be useful to return the Pochhammer symbol value, prod, as a pair that supplies roughly twice working precision.
Tests of a Pochhammer function that uses pair-precision arithmetic show that its errors are in the range [−0.97, +0.97] ulps, with an average error of about 0.24 ulps. The error plots look like this:
A test with logarithmically distributed x values has an average error of 0.25 ulps, and a maximum error of −1.28 ulps:
Small modifications of four of our test programs created versions that compute the Pochhammer symbol in decimal floating-point arithmetic. Even though our first version of twosum() is not strictly valid for decimal arithmetic, the Pochhammer test results are essentially identical for both versions of that function, so we show only one set of results here.
With uniform sampling of x, errors in ulps are in [−1.27, +1.27], with an average error of about 0.13 ulps:
With logarithmic sampling of x, errors in ulps are in [−1.27, +1.27], with an average error of about 0.13 ulps, similar to the results for uniform sampling:
Compared to the results of the preceding section for binary arithmetic, more of the decimal results are correctly rounded: 88.6% versus 98.2%.
The prolific numerical analyst Fredrik Johansson recently suggested in an article, Arbitrary-precision computation of the gamma function, that our second-from-last Pochhammer relation, (x)n = (x)k (x + k)n − k, given at the beginning of this document might be usefully exploited in a divide-and-conquer method. For large n, split that number in half, compute the Pochhammer value for each half by the same function, and multiply the pairs. Near some small value of n (he suggested 10), stop subdividing, and use the normal Pochhammer loop.
Because each term in the Pochhammer sequence is unique, that approach still requires n−1 additions and multiplications, so unlike many applications of divide-and-conquer methods that reduce the total work to a fraction (log2n) / n of the original method, there is no time savings for the Pochhammer computation. However, the computed sequences are shorter, so perhaps error growth can be reduced.
Each of our test programs implements this alterative, selected by a compile-time definition of the symbol TEST2. The code for one of them looks like this:
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_mpfr_std(x, n);
else
{
int m;
m = n / 2;
result = poch_daq(x, m) * poch_daq(x + m, n - m);
}
return (result);
}
We expect about log2n executions of the else block, contributing about that many rounding errors (unless they are done at somewhat higher intermediate precision). Here is a graphical comparison of the simple Pochhammer version against a multiple-precision version using MPFR, both using Johansson's idea:
As we suggested, the uncorrected extra multiplications increase the total error, and Johansson's divide-and-conquer idea is not to be recommended. The only apparent area where it still might be useful is multiple precision computation, where a somewhat smaller intermediate precision might allow a slight speedup over the standard Pochhammer loop.
For computing any function, f(x), it is helpful to know in advance the worst-case argument x. There are at least two of these arguments. One of them is the x value for which the relative error in f(x) has the largest magnitude. The other is a value for which high-precision computation of f(x) gives a result that is almost exactly halfway between two representable floating-point numbers, and thus, minor perturbations in its calculation could cause incorrect rounding, producing a result that is at least one ulp in error. The maximum relative error gives us an idea of how much extra precision is needed, but the minimum nonzero error identifies the smallest computational precision needed to obtain a correctly rounded result. The worst cases depend on the working precision, and the difficulty of finding them grows exponentially with increasing precision: there are about 2**(53 − 24) = 536_870_912 as many numbers to test when we move from the binary32 to the binary64 format, 10**12 as many in the binary80 format, and almost 10**27 as many in the binary128 format.
That problem has long been known as the Table Maker's Dilemma, a term introduced in a 1974 conference talk by William Kahan. For most floating-point functions in programming language libraries, finding the worst case essentially requires exhaustive evaluation of f(x) at higher precision, for all x in working precision. Book treatments of the problem are in Chapter 12 in the first edition of the Handbook of Floating-Point Arithmetic, the long Section 10.5 in the second edition, and Chapter 12 of the third edition of Elementary Functions: Algorithms and Implementation (DOI 10.1007/978-1-4899-7983-4).
Once the worst cases are known, we can bound the accuracy needed to compute f(x) with guaranteed correct rounding, just as we have with the basic IEEE 754 operations.
The research reported earlier found the worst-case arguments for the sine, cosine, and tangent in binary64 arithmetic. They have, respectively, runs of 68, 66, and 62 consecutive 0 or 1 bits after the rounding bit. Thus, for the sine, the result that is hardest to round correctly must be computed to a precision at least 53 + 68 = 121 bits, which is more than even the 113 bits of the binary128 format can provide. The report required an extremely complex analysis to drastically reduce the number of x values to be tested, and even then, the search required parallel computation and about 6000 CPU-core-years at a large national computer facility.
For the binary32 and decimal32 formats, there are only about 2**32 = 4_294_967_296 possible arguments, and a brute force search can try them all. For the Pochhammer symbol, we can make a further reduction in the argument count by noting that when x is sufficiently small that fl(1 + x) == 1, then poch(x,n) reduces to x * (n − 1)!, and that can be computed accurately with a single multiplication against a precomputed table of factorials stored as sums of exact high, and accurate low, parts.
Similarly, when x is sufficiently large that fl(x + n) == x, then poch(x,n) reduces to x**n, which can be evaluated with about log2(n) multiplications. The power function, x**y or pow(x,y), is the most difficult of all of the elementary functions to compute accurately, essentially because its error magnification is proportional to the power y. Thus, math libraries may have rather complicated implementations of that function to achieve results with small errors.
For decimal32 arithmetic, those constraints mean that for each n value, we only have to examine all significands from 1.000_000 to 9.999_999, with power-of-ten exponents from −7 to +8. That is just 16 decades, instead of the 192 decades that are representable, a 12× speedup. In each decade, there are about 9_000_000 significands to test. The total work for all needed (x,n) pairs is about 288_000_000 floating-point operations, plus some overhead for recording the worst cases for each decade and n value.
Test programs poch-dec-worst*.c and poch-worst*.c find the worst cases for n in [2,1000], where poch(x,n) near the upper limit on n overflows.
For the binary32 format, the complete test on a 2025-vintage server with a 3.15GHz AMD EPYC 9455 CPU took only about 2 minutes, and found the maximum relative error at x = 0x1.f066_96p+4 ~= 31.025 and n = 23, for which the computed Pochhammer result is 0x1.897b_8cp+123 ~= 1.634e37. An accurate binary64 result is 0x1.897b_6410_e025_8p+123, and a binary128 result is +0x1.897b_6410_e025_8ed8_82b5_1b83_addbp+123. The binary32 result is wrong in its 7 low-order bits. The worst-case error is 12.99 ulps, which suggests that we need to compute a binary32 Pochhammer result with at least 7 + 4 = 11 extra significand bits. That is easily done by computing in the binary64 format, and then rounding.
The minimum nonzero error in binary32 is found for x = 0x1.ff40_08p+10 ~= 2_045 and n = 7, for which poch(x,n) = 0x1.ffff_c8p+76. The binary64 result is 0x1.ffff_c800_0000_1p+76, with a run of 29 0-bits after the last binary32 bit. The binary128 result is +0x1.ffff_c800_0000_0003_0000_0cp+76, which has a run of 39 0-bits after the last binary32 bit. A 100-decimal-digit computation in Maple produces a result of +0x1.ffff_c800_0000_0003_0000_0bff_ffff_ffff_2000_008p+76, which rounds to our binary128 result, and has the same bit run. In this case, however, the long 0-bit runs do not affect the rounding, and the computed binary32 Pochhammer result is correct.
Tests with the software implementation of decimal32 arithmetic were notably slower: 11 hours. The worst case for maximum relative error is for poch(9.987_065e+01, 46) which has an error of −22.40 ulps, with a decimal32 result of 8.209_397e+95, a decimal64 result of 8.209_580_839_234_025e+95, and a decimal128 result of 8.209_580_839_234_030_463_532_429_805_542_56e+95. We need at least 3 or 4 more decimal digits to get correct rounding; decimal64 computation is more than ample.
The minimum nonzero error in decimal32 is found for poch(9.999_940e+05, 11), which returns 9.999_890e+65. The decimal64 result is 9.999_890_000_000_003e+65, and the decimal128 result is 9.999_890_000_000_003_299_993_729_969_970_074e+65. For both cases, that is a run of 8 decimal digits, and decimal64 computation suffices.
In summary, we have now presented eight independent implementations of the Pochhammer symbol computation in IEEE 754 binary floating-point arithmetic, a simple one with large rounding errors, one with faithful rounding, one with errors below 1.30 ulps, one with divide-and-conquer, and four with results that are almost always correctly rounded, lying within a half ulp of an exact infinitely precise result. However, we have certainly not proved that experimental finding mathematically, or with a computer algebra proof system.
The test code and plots, and worst-case searches for the 32-bit formats, suggest that a reasonable approach for many applications of the Pochhammer symbol is to use higher internal precision. Otherwise, use pair-precision computation of the sums and product; error compensation of the product alone is not sufficient. If correct rounding with high probability is required, then any one of three multiple-precision packages (FLINT, FX, and MPFR) are recommended.
A production version of poch(x,n) should also include two additional tests for the cases discussed in the previous section where the Pochhammer symbol reduces computationally to x * (n − 1)! or pow(x,n), provided that the implementation of the power function is sufficiently accurate, which historically has often not been the case.
The C implementations are available here, including a test script to run all them and produce gnuplot data files that can be tweaked to plot other data ranges, and create output in other graphics formats.
The gnuplot files are available here, and the graphics output files are here. For convenience, there is also an InfoZip file with the complete set of files for this document.
For more information, see the MathWorld article Pochhammer Symbol, and the encyclopedia article Falling and rising factorials. The former discusses confusion from alternate historical notations (ours follows the well-known NBS Handbook of Mathematical Functions (1964) and the NIST Handbook of Mathematical Functions (2010), ISBN 0-521-19225-0), and also treats the falling factorial, which just labels the Pochhammer symbol product from its other end, with a possible sign change.
For background on floating-point arithmetic, see the Handbook of Floating-Point Arithmetic (DOI 10.1007/978-3-319-76526-6), and this author's The Mathematical-Function Computation Handbook (DOI 10.1007/978-3-319-64110-2).
The bibliographies on floating-point arithmetic and elementary and special functions together provide more than 11_000 literature references, including ones for more than 600 books. There are bibliographies of about 30 journals on numerical computation in the TeX User Group archive, and there are bibliographies of about 70 important computer scientists and numerical analysts in the BibNet Project archive.
Those bibliographies supply valuable resources for the history of the progress in numerical computation, much of which began largely in the 1930s and 1940s, and led to developments of expensive mainframe computers, and with later advances in semiconductor manufacturing, and immense production volumes, billions, and possibly even trillions, of processors in desktop and embedded devices for communication, personal, industrial, military, and transportation applications. Humanity has been changed enormously, but there are older roots that predate electronic computers: the bibliographies record publications back to 1624 in function computation, to 1703 in binary arithmetic, and to 1404 in cryptography.