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.
This uses an approximation based on the Stirling series for large enough x (for which it is highly accurate). For smaller x, identities are used to express gamma(x) in terms of gamma(x+1) or gamma(1-x), ultimately letting the Stirling series approximation be used.
This implementation is based on the approximations given in the following paper:
W. J. Cody, Rational Chebyshev Approximations for the Error Function, Mathematics of Computation, Vol. 23, No. 107 (Jul., 1969), pp. 631-637.
Per the paper, the approximations have maximal relative error of 6e-19 or lower (although I have not verified what the figure is for this actual implementation).
See also Cody's FORTRAN implementation based on the same approach:
https://netlib.org/specfun/erf
This is similar to the approach recommended in Apple Numerics Manual Ch. 9, except that there is an added case for large values that would otherwise cause an overflow or spuriously report underflow.
The previous simple one could be wrong in several low-order digits due to the inaccuracy in the representation of the exponent (1/3). This version effectively breaks the number up into the form a*8^b, computes the cube root of 8^b exactly (i.e. 2^b), and uses the slightly inaccurate exponentiation only for a.
This uses the obvious calculation, except with scaling to avoid unnecessary overflow/underflow.
There is a discussion of hypot implementations in C. Borges, An Improved Algorithm for hypot(a,b) (https://arxiv.org/pdf/1904.09481.pdf). This implementation is similar to the "Naive (Unfused)" version discussed in that paper. As the paper notes, it is possible to get better accuracy by adding a correction term, but the "naive" version is already reasonably good, so we skip the correction in the interest of code size and speed.
The lrint functions could give the wrong result for negative numbers in upward/downward rounding modes. Casts to comp could also have different rounding behavior.
This is basically the implementation recommended in Apple Numerics Manual Ch. 9, except that there is an added case for large values that would otherwise cause an overflow.
These SANE operations can sometimes return incorrect values for certain negative integers such as -2147483648 and -53021371269120 (numbers with at least 16 low-order zero bits in their two's-complement representation). To work around this, we now avoid calling FX2C or FX2L on negative numbers, generally by saving and restoring the sign separately.
These workarounds are used in several of the new <math.h> rounding functions, and also for code that converts floating-point values to comp or long long. There are some places in SysFloat that should be patched similarly, so we may still hit this problem in certain situations until that is done.
This should work, and mostly does. However, it is affected by a bug in FX2L (and FX2C) which can sometimes give the wrong results for certain negative integers (such as -2147483648). I believe this can occur when at least the lower 16 bits if the integer (in two's-complement representation) are zeros.
This is a bit more complex than other rounding functions, because it rounds to nearest but always away from zero in halfway cases, which is not a rounding direction directly supported by SANE.
This differs from scalbn in that the exponent has type long. When scaling an extended value, exponents slightly outside the range of int can actually be used meaningfully. We address this by doing multiple SCALBX calls (at most 2) in a loop.
Most of these actually just jump to the existing functions, since they really use extended precision anyway. The exception is the modf functions, which need a separate implementation for each type because they store a value through a pointer to that type.
The functions implemented so far are largely the ones that map (nearly) directly to SANE calls.
Note that C99 specifies separate float/double/long double versions of each of these functions, but under ORCA/C they generally use the same code.
These macros differ from the normal comparison operators in that they will not signal invalid due to the operands being unordered. However, this implementation will still signal invalid for SNaNs. That is clearly OK according to the wording in draft C23. C17 and earlier do not mention that possibility, but they do not really specify the behavior of SNaNs in general.
Stephen Heumann