in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java [438:539]
private static double evaluateRational(double[] a, int[] b, double x) {
// The choice of algorithm in Boost is based on the compiler
// to suite the available optimisations.
//
// Tests against rational_horner1_13.hpp which uses a first order
// Horner method (no x*x term) show only minor variations in
// error. rational_horner2_13.hpp has the same second order Horner
// method with different code layout of the same sum.
// rational_horner3_13.hpp
if (x <= 1) {
final double x2 = x * x;
double t0 = a[12] * x2 + a[10];
double t1 = a[11] * x2 + a[9];
double t2 = b[12] * x2 + b[10];
double t3 = b[11] * x2 + b[9];
t0 *= x2;
t1 *= x2;
t2 *= x2;
t3 *= x2;
t0 += a[8];
t1 += a[7];
t2 += b[8];
t3 += b[7];
t0 *= x2;
t1 *= x2;
t2 *= x2;
t3 *= x2;
t0 += a[6];
t1 += a[5];
t2 += b[6];
t3 += b[5];
t0 *= x2;
t1 *= x2;
t2 *= x2;
t3 *= x2;
t0 += a[4];
t1 += a[3];
t2 += b[4];
t3 += b[3];
t0 *= x2;
t1 *= x2;
t2 *= x2;
t3 *= x2;
t0 += a[2];
t1 += a[1];
t2 += b[2];
t3 += b[1];
t0 *= x2;
t2 *= x2;
t0 += a[0];
t2 += b[0];
t1 *= x;
t3 *= x;
return (t0 + t1) / (t2 + t3);
}
final double z = 1 / x;
final double z2 = 1 / (x * x);
double t0 = a[0] * z2 + a[2];
double t1 = a[1] * z2 + a[3];
double t2 = b[0] * z2 + b[2];
double t3 = b[1] * z2 + b[3];
t0 *= z2;
t1 *= z2;
t2 *= z2;
t3 *= z2;
t0 += a[4];
t1 += a[5];
t2 += b[4];
t3 += b[5];
t0 *= z2;
t1 *= z2;
t2 *= z2;
t3 *= z2;
t0 += a[6];
t1 += a[7];
t2 += b[6];
t3 += b[7];
t0 *= z2;
t1 *= z2;
t2 *= z2;
t3 *= z2;
t0 += a[8];
t1 += a[9];
t2 += b[8];
t3 += b[9];
t0 *= z2;
t1 *= z2;
t2 *= z2;
t3 *= z2;
t0 += a[10];
t1 += a[11];
t2 += b[10];
t3 += b[11];
t0 *= z2;
t2 *= z2;
t0 += a[12];
t2 += b[12];
t1 *= z;
t3 *= z;
return (t0 + t1) / (t2 + t3);
}