[
https://issues.apache.org/jira/browse/NUMBERS-179?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17459272#comment-17459272
]
Alex Herbert commented on NUMBERS-179:
--------------------------------------
The Boost threshold to switch to an alternative computation is set using:
{code:java}
double agh = a + Lanczos.g() - 0.5;
double d = ((z - a) - Lanczos.g() + 0.5) / agh;
if (Math.abs(d * d * a) <= 100 && a > 150) {
// ...
{code}
The threshold of a > 150 is too high.
However the current code in Numbers has an addition to the Boost code. It
detects when the regularised gamma prefix (which is the gamma P derivative
multiplied by z) can be computed directly using:
{code:java}
if (a <= 171) {
// Check for overflow or underflow ...
// if OK
return Math.pow(z, a) * Math.exp(-z) / gamma(a);
} {code}
This was added as it increased accuracy when the gamma function can be computed
without a combined formula using the Lanczos approximation. Since this has been
added it is rare that the check for a ~ z is made for small a. Most of the
small a computations use the direct formula. The alternative formula is used
when the direct formula would overflow. This occurs when:
{noformat}
pow(z, a) > MAX_VALUE
z = exp(ln(MAX_VALUE) / a)
{noformat}
Using ln(MAX_VALUE) as 709 (as is done in the current code) the smallest value
of z that will require detection of a ~ z is shown in the following table. Also
listed is 2a which would be a conservative upper limit for any z where a ~ z,
the value d for min z and the condition value abs(d*d*a):
||a||min z||2a||d||abs(d*d*a)||
|120|368.1|240|1.93E+00|4.48E+02|
|121|350.6|242|1.77E+00|3.79E+02|
|122|334.1|244|1.62E+00|3.20E+02|
|123|318.7|246|1.48E+00|2.69E+02|
|124|304.2|248|1.35E+00|2.26E+02|
|125|290.6|250|1.23E+00|1.88E+02|
|126|277.8|252|1.11E+00|1.56E+02|
|127|265.8|254|1.01E+00|1.28E+02|
|128|254.4|256|9.06E-01|1.05E+02|
|129|243.7|258|8.12E-01|8.50E+01|
|130|233.7|260|7.24E-01|6.82E+01|
|131|224.1|262|6.42E-01|5.39E+01|
|132|215.1|264|5.64E-01|4.20E+01|
|133|206.6|266|4.91E-01|3.21E+01|
|134|198.6|268|4.23E-01|2.40E+01|
|135|190.9|270|3.59E-01|1.74E+01|
|136|183.7|272|2.98E-01|1.21E+01|
|137|176.8|274|2.41E-01|7.94E+00|
|138|170.3|276|1.87E-01|4.81E+00|
|139|164.1|278|1.36E-01|2.56E+00|
|140|158.3|280|8.76E-02|1.07E+00|
|141|152.7|282|4.20E-02|2.49E-01|
|142|147.4|284|-1.04E-03|1.53E-04|
|143|142.3|286|-4.18E-02|2.50E-01|
|144|137.5|288|-8.04E-02|9.32E-01|
|145|132.9|290|-1.17E-01|1.99E+00|
|146|128.5|292|-1.52E-01|3.36E+00|
|147|124.4|294|-1.85E-01|5.01E+00|
|148|120.4|296|-2.16E-01|6.90E+00|
|149|116.6|298|-2.46E-01|8.99E+00|
|150|112.9|300|-2.74E-01|1.13E+01|
|151|109.4|302|-3.01E-01|1.37E+01|
|152|106.1|304|-3.26E-01|1.62E+01|
|153|102.9|306|-3.51E-01|1.88E+01|
|154|99.9|308|-3.74E-01|2.15E+01|
|155|96.9|310|-3.96E-01|2.43E+01|
|156|94.1|312|-4.17E-01|2.71E+01|
|157|91.5|314|-4.37E-01|3.00E+01|
|158|88.9|316|-4.56E-01|3.29E+01|
|159|86.4|318|-4.75E-01|3.58E+01|
The values of interest are those where min z < 2a. This is the range for a
where a ~ z is possible. This is a >= 128.
Note that the smallest value of d is when a=142. This is the shape parameter
used in the test in Commons Statistics that discovered an error in this
threshold (see STATISTICS-39). Thus it may be assumed that the test data was
created to specifically target this condition.
The previous threshold was a > 150. This is well into the range where a ~ z is
possible. The current level of 100 is not in the table. It is not possible for
a to be similar in magnitude to z using the overflow condition for the regular
computation.
There is an underflow condition too:
{noformat}
pow(z, a) < MIN_NORMAL
z = exp(ln(MIN_NORMAL) / a)
{noformat}
Here are the maximum values of z before underflow:
||a||max z||
|1|3.31E-308|
|10|1.79E-31|
|100|8.42E-04|
|150|8.92E-03|
Note: Any 'a' less than 1 is handled using a different computation. Here the
values of z are small compared to a. It is not possible for a ~ z. So the
condition must be robust to avoid choosing the wrong computation when z is
small.
If z is near zero the condition creates a value of d close to but greater than
-1:
{code:java}
double agh = a + Lanczos.g() - 0.5;
double d = ((z-a) - Lanczos.g() + 0.5) / agh;
// d => (-a - g + 0.5) / (a + g - 0.5) -> -1
if (Math.abs(d * d * a) <= 100 && a > threshold) {
// ...
{code}
To eliminate very small z from triggering the computation the threshold for a
should be large enough that {{(d*d*a) <= 100}} is false when d is close to -1.
So the threshold should be higher than 100.
Looking at the previous table the value a=128 creates min z around twice the
value of a. It is the largest a where {{abs(d*d*a) <= 100}} is false. It would
eliminate the condition being true when z is tiny and the direct computation
would underflow in the power function.
This analysis suggests changing the Boost threshold for large a from 150 to
128. It does not alter the second threshold using {{{}d*d*a <= 100{}}}.
> Optimise the choice of function for the gamma P derivative
> ----------------------------------------------------------
>
> Key: NUMBERS-179
> URL: https://issues.apache.org/jira/browse/NUMBERS-179
> Project: Commons Numbers
> Issue Type: Improvement
> Components: gamma
> Affects Versions: 1.1
> Reporter: Alex Herbert
> Priority: Minor
>
> The incomplete gamma function is evaluated using the Lanczos approximation.
> This is evaluated using two alternate expressions when the terms a and z are
> large:
> See [Incomplete Gamma Eq. 15 and Eq.
> 16|https://www.boost.org/doc/libs/1_78_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html]
> The defaults used in the Boost C++ implementation required updating to
> increase the precision of the function when a=142 for [STATISTICS-39]. The
> threshold for large a was lowered from 150 to 100. The switch between the two
> functions requires optimisation to ensure the threshold is applicable to
> other parameterisations, e.g. a from 50 to 150.
--
This message was sent by Atlassian Jira
(v8.20.1#820001)
