if (increment <= PG_SUM_RELATIVE_ERROR * sum) {
break;
}
}
if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
}
ret = sum * FastMath.sqrt(2 * FastMath.PI) / z;
// K_1(z)
// Sum is -inf to inf, but k term is always (k + 1/2) ^ 2, so really have
// twice the sum from k = 0 to inf (k = -1 is same as 0, -2 same as 1, ...)
final double twoZ2 = 2 * z2;
sum = 0;
kTerm = 0;
double kTerm2 = 0;
for (k = 0; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
kTerm = k + 0.5;
kTerm2 = kTerm * kTerm;
increment = (PI_SQUARED * kTerm2 - z2) * FastMath.exp(-PI_SQUARED * kTerm2 / twoZ2);
sum += increment;
if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum)) {
break;
}
}
if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
}
final double sqrtHalfPi = FastMath.sqrt(PI / 2);
// Instead of doubling sum, divide by 3 instead of 6
ret += sum * sqrtHalfPi / (3 * z4 * sqrtN);
// K_2(z)
// Same drill as K_1, but with two doubly infinite sums, all k terms are even powers.
final double z4Term = 2 * z4;
final double z6Term = 6 * z6;
z2Term = 5 * z2;
final double pi4 = PI_SQUARED * PI_SQUARED;
sum = 0;
kTerm = 0;
kTerm2 = 0;
for (k = 0; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
kTerm = k + 0.5;
kTerm2 = kTerm * kTerm;
increment = (z6Term + z4Term + PI_SQUARED * (z4Term - z2Term) * kTerm2 +
pi4 * (1 - twoZ2) * kTerm2 * kTerm2) * FastMath.exp(-PI_SQUARED * kTerm2 / twoZ2);
sum += increment;
if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum)) {
break;
}
}
if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
}
double sum2 = 0;
kTerm2 = 0;
for (k = 1; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
kTerm2 = k * k;
increment = PI_SQUARED * kTerm2 * FastMath.exp(-PI_SQUARED * kTerm2 / twoZ2);
sum2 += increment;
if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum2)) {
break;
}
}
if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
}
// Again, adjust coefficients instead of doubling sum, sum2
ret += (sqrtHalfPi / n) * (sum / (36 * z2 * z2 * z2 * z) - sum2 / (18 * z2 * z));
// K_3(z) One more time with feeling - two doubly infinite sums, all k powers even.
// Multiply coefficient denominators by 2, so omit doubling sums.
final double pi6 = pi4 * PI_SQUARED;
sum = 0;
double kTerm4 = 0;
double kTerm6 = 0;
for (k = 0; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
kTerm = k + 0.5;
kTerm2 = kTerm * kTerm;
kTerm4 = kTerm2 * kTerm2;
kTerm6 = kTerm4 * kTerm2;
increment = (pi6 * kTerm6 * (5 - 30 * z2) + pi4 * kTerm4 * (-60 * z2 + 212 * z4) +
PI_SQUARED * kTerm2 * (135 * z4 - 96 * z6) - 30 * z6 - 90 * z8) *
FastMath.exp(-PI_SQUARED * kTerm2 / twoZ2);
sum += increment;
if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum)) {
break;
}
}
if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
}
sum2 = 0;
for (k = 1; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
kTerm2 = k * k;
kTerm4 = kTerm2 * kTerm2;
increment = (-pi4 * kTerm4 + 3 * PI_SQUARED * kTerm2 * z2) *
FastMath.exp(-PI_SQUARED * kTerm2 / twoZ2);
sum2 += increment;
if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum2)) {
break;
}
}
if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
}
return ret + (sqrtHalfPi / (sqrtN * n)) * (sum / (3240 * z6 * z4) +
+ sum2 / (108 * z6));
}