001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math4.legacy.analysis.integration.gauss; 018 019import java.math.BigDecimal; 020import java.math.MathContext; 021 022import org.apache.commons.math4.legacy.core.Pair; 023 024/** 025 * Factory that creates Gauss-type quadrature rule using Legendre polynomials. 026 * In this implementation, the lower and upper bounds of the natural interval 027 * of integration are -1 and 1, respectively. 028 * The Legendre polynomials are evaluated using the recurrence relation 029 * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"> 030 * Abramowitz and Stegun, 1964</a>. 031 * 032 * @since 3.1 033 */ 034public class LegendreHighPrecisionRuleFactory extends BaseRuleFactory<BigDecimal> { 035 /** Settings for enhanced precision computations. */ 036 private final MathContext mContext; 037 /** The number {@code 2}. */ 038 private final BigDecimal two; 039 /** The number {@code -1}. */ 040 private final BigDecimal minusOne; 041 /** The number {@code 0.5}. */ 042 private final BigDecimal oneHalf; 043 044 /** 045 * Default precision is {@link MathContext#DECIMAL128 DECIMAL128}. 046 */ 047 public LegendreHighPrecisionRuleFactory() { 048 this(MathContext.DECIMAL128); 049 } 050 051 /** 052 * @param mContext Precision setting for computing the quadrature rules. 053 */ 054 public LegendreHighPrecisionRuleFactory(MathContext mContext) { 055 this.mContext = mContext; 056 two = new BigDecimal("2", mContext); 057 minusOne = new BigDecimal("-1", mContext); 058 oneHalf = new BigDecimal("0.5", mContext); 059 } 060 061 /** {@inheritDoc} */ 062 @Override 063 protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints) { 064 if (numberOfPoints == 1) { 065 // Break recursion. 066 return new Pair<>(new BigDecimal[] { BigDecimal.ZERO }, 067 new BigDecimal[] { two }); 068 } 069 070 // Get previous rule. 071 // If it has not been computed yet it will trigger a recursive call 072 // to this method. 073 final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst(); 074 075 // Compute next rule. 076 final BigDecimal[] points = new BigDecimal[numberOfPoints]; 077 final BigDecimal[] weights = new BigDecimal[numberOfPoints]; 078 079 // Find i-th root of P[n+1] by bracketing. 080 final int iMax = numberOfPoints / 2; 081 for (int i = 0; i < iMax; i++) { 082 // Lower-bound of the interval. 083 BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1]; 084 // Upper-bound of the interval. 085 BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i]; 086 // P[j-1](a) 087 BigDecimal pma = BigDecimal.ONE; 088 // P[j](a) 089 BigDecimal pa = a; 090 // P[j-1](b) 091 BigDecimal pmb = BigDecimal.ONE; 092 // P[j](b) 093 BigDecimal pb = b; 094 for (int j = 1; j < numberOfPoints; j++) { 095 final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext); 096 final BigDecimal b_j = new BigDecimal(j, mContext); 097 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); 098 099 // Compute P[j+1](a) 100 // ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1); 101 102 BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext); 103 tmp1 = pa.multiply(tmp1, mContext); 104 BigDecimal tmp2 = pma.multiply(b_j, mContext); 105 // P[j+1](a) 106 BigDecimal ppa = tmp1.subtract(tmp2, mContext); 107 ppa = ppa.divide(b_j_p_1, mContext); 108 109 // Compute P[j+1](b) 110 // ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1); 111 112 tmp1 = b.multiply(b_two_j_p_1, mContext); 113 tmp1 = pb.multiply(tmp1, mContext); 114 tmp2 = pmb.multiply(b_j, mContext); 115 // P[j+1](b) 116 BigDecimal ppb = tmp1.subtract(tmp2, mContext); 117 ppb = ppb.divide(b_j_p_1, mContext); 118 119 pma = pa; 120 pa = ppa; 121 pmb = pb; 122 pb = ppb; 123 } 124 // Now pa = P[n+1](a), and pma = P[n](a). Same holds for b. 125 // Middle of the interval. 126 BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext); 127 // P[j-1](c) 128 BigDecimal pmc = BigDecimal.ONE; 129 // P[j](c) 130 BigDecimal pc = c; 131 boolean done = false; 132 while (!done) { 133 BigDecimal tmp1 = b.subtract(a, mContext); 134 BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext); 135 done = tmp1.compareTo(tmp2) <= 0; 136 pmc = BigDecimal.ONE; 137 pc = c; 138 for (int j = 1; j < numberOfPoints; j++) { 139 final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext); 140 final BigDecimal b_j = new BigDecimal(j, mContext); 141 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); 142 143 // Compute P[j+1](c) 144 tmp1 = c.multiply(b_two_j_p_1, mContext); 145 tmp1 = pc.multiply(tmp1, mContext); 146 tmp2 = pmc.multiply(b_j, mContext); 147 // P[j+1](c) 148 BigDecimal ppc = tmp1.subtract(tmp2, mContext); 149 ppc = ppc.divide(b_j_p_1, mContext); 150 151 pmc = pc; 152 pc = ppc; 153 } 154 // Now pc = P[n+1](c) and pmc = P[n](c). 155 if (!done) { 156 if (pa.signum() * pc.signum() <= 0) { 157 b = c; 158 pmb = pmc; 159 pb = pc; 160 } else { 161 a = c; 162 pma = pmc; 163 pa = pc; 164 } 165 c = a.add(b, mContext).multiply(oneHalf, mContext); 166 } 167 } 168 final BigDecimal nP = new BigDecimal(numberOfPoints, mContext); 169 BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext); 170 tmp1 = tmp1.multiply(nP); 171 tmp1 = tmp1.pow(2, mContext); 172 BigDecimal tmp2 = c.pow(2, mContext); 173 tmp2 = BigDecimal.ONE.subtract(tmp2, mContext); 174 tmp2 = tmp2.multiply(two, mContext); 175 tmp2 = tmp2.divide(tmp1, mContext); 176 177 points[i] = c; 178 weights[i] = tmp2; 179 180 final int idx = numberOfPoints - i - 1; 181 points[idx] = c.negate(mContext); 182 weights[idx] = tmp2; 183 } 184 // If "numberOfPoints" is odd, 0 is a root. 185 // Note: as written, the test for oddness will work for negative 186 // integers too (although it is not necessary here), preventing 187 // a FindBugs warning. 188 if ((numberOfPoints & 1) != 0) { 189 BigDecimal pmc = BigDecimal.ONE; 190 for (int j = 1; j < numberOfPoints; j += 2) { 191 final BigDecimal b_j = new BigDecimal(j, mContext); 192 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); 193 194 // pmc = -j * pmc / (j + 1); 195 pmc = pmc.multiply(b_j, mContext); 196 pmc = pmc.divide(b_j_p_1, mContext); 197 pmc = pmc.negate(mContext); 198 } 199 200 // 2 / pow(numberOfPoints * pmc, 2); 201 final BigDecimal nP = new BigDecimal(numberOfPoints, mContext); 202 BigDecimal tmp1 = pmc.multiply(nP, mContext); 203 tmp1 = tmp1.pow(2, mContext); 204 BigDecimal tmp2 = two.divide(tmp1, mContext); 205 206 points[iMax] = BigDecimal.ZERO; 207 weights[iMax] = tmp2; 208 } 209 210 return new Pair<>(points, weights); 211 } 212}