in commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java [63:211]
protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints) {
if (numberOfPoints == 1) {
// Break recursion.
return new Pair<>(new BigDecimal[] { BigDecimal.ZERO },
new BigDecimal[] { two });
}
// Get previous rule.
// If it has not been computed yet it will trigger a recursive call
// to this method.
final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
// Compute next rule.
final BigDecimal[] points = new BigDecimal[numberOfPoints];
final BigDecimal[] weights = new BigDecimal[numberOfPoints];
// Find i-th root of P[n+1] by bracketing.
final int iMax = numberOfPoints / 2;
for (int i = 0; i < iMax; i++) {
// Lower-bound of the interval.
BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1];
// Upper-bound of the interval.
BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i];
// P[j-1](a)
BigDecimal pma = BigDecimal.ONE;
// P[j](a)
BigDecimal pa = a;
// P[j-1](b)
BigDecimal pmb = BigDecimal.ONE;
// P[j](b)
BigDecimal pb = b;
for (int j = 1; j < numberOfPoints; j++) {
final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
final BigDecimal b_j = new BigDecimal(j, mContext);
final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
// Compute P[j+1](a)
// ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1);
BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext);
tmp1 = pa.multiply(tmp1, mContext);
BigDecimal tmp2 = pma.multiply(b_j, mContext);
// P[j+1](a)
BigDecimal ppa = tmp1.subtract(tmp2, mContext);
ppa = ppa.divide(b_j_p_1, mContext);
// Compute P[j+1](b)
// ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1);
tmp1 = b.multiply(b_two_j_p_1, mContext);
tmp1 = pb.multiply(tmp1, mContext);
tmp2 = pmb.multiply(b_j, mContext);
// P[j+1](b)
BigDecimal ppb = tmp1.subtract(tmp2, mContext);
ppb = ppb.divide(b_j_p_1, mContext);
pma = pa;
pa = ppa;
pmb = pb;
pb = ppb;
}
// Now pa = P[n+1](a), and pma = P[n](a). Same holds for b.
// Middle of the interval.
BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext);
// P[j-1](c)
BigDecimal pmc = BigDecimal.ONE;
// P[j](c)
BigDecimal pc = c;
boolean done = false;
while (!done) {
BigDecimal tmp1 = b.subtract(a, mContext);
BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext);
done = tmp1.compareTo(tmp2) <= 0;
pmc = BigDecimal.ONE;
pc = c;
for (int j = 1; j < numberOfPoints; j++) {
final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
final BigDecimal b_j = new BigDecimal(j, mContext);
final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
// Compute P[j+1](c)
tmp1 = c.multiply(b_two_j_p_1, mContext);
tmp1 = pc.multiply(tmp1, mContext);
tmp2 = pmc.multiply(b_j, mContext);
// P[j+1](c)
BigDecimal ppc = tmp1.subtract(tmp2, mContext);
ppc = ppc.divide(b_j_p_1, mContext);
pmc = pc;
pc = ppc;
}
// Now pc = P[n+1](c) and pmc = P[n](c).
if (!done) {
if (pa.signum() * pc.signum() <= 0) {
b = c;
pmb = pmc;
pb = pc;
} else {
a = c;
pma = pmc;
pa = pc;
}
c = a.add(b, mContext).multiply(oneHalf, mContext);
}
}
final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext);
tmp1 = tmp1.multiply(nP);
tmp1 = tmp1.pow(2, mContext);
BigDecimal tmp2 = c.pow(2, mContext);
tmp2 = BigDecimal.ONE.subtract(tmp2, mContext);
tmp2 = tmp2.multiply(two, mContext);
tmp2 = tmp2.divide(tmp1, mContext);
points[i] = c;
weights[i] = tmp2;
final int idx = numberOfPoints - i - 1;
points[idx] = c.negate(mContext);
weights[idx] = tmp2;
}
// If "numberOfPoints" is odd, 0 is a root.
// Note: as written, the test for oddness will work for negative
// integers too (although it is not necessary here), preventing
// a FindBugs warning.
if ((numberOfPoints & 1) != 0) {
BigDecimal pmc = BigDecimal.ONE;
for (int j = 1; j < numberOfPoints; j += 2) {
final BigDecimal b_j = new BigDecimal(j, mContext);
final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
// pmc = -j * pmc / (j + 1);
pmc = pmc.multiply(b_j, mContext);
pmc = pmc.divide(b_j_p_1, mContext);
pmc = pmc.negate(mContext);
}
// 2 / pow(numberOfPoints * pmc, 2);
final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
BigDecimal tmp1 = pmc.multiply(nP, mContext);
tmp1 = tmp1.pow(2, mContext);
BigDecimal tmp2 = two.divide(tmp1, mContext);
points[iMax] = BigDecimal.ZERO;
weights[iMax] = tmp2;
}
return new Pair<>(points, weights);
}