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 */
017 package org.apache.commons.math3.analysis.integration.gauss;
018
019 import java.math.MathContext;
020 import java.math.BigDecimal;
021 import org.apache.commons.math3.util.Pair;
022 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
023 import org.apache.commons.math3.exception.util.LocalizedFormats;
024
025 /**
026 * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
027 * In this implementation, the lower and upper bounds of the natural interval
028 * of integration are -1 and 1, respectively.
029 * The Legendre polynomials are evaluated using the recurrence relation
030 * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"
031 * Abramowitz and Stegun, 1964</a>.
032 *
033 * @since 3.1
034 * @version $Id: LegendreHighPrecisionRuleFactory.java 1374187 2012-08-17 09:50:54Z erans $
035 */
036 public class LegendreHighPrecisionRuleFactory extends BaseRuleFactory<BigDecimal> {
037 /** Settings for enhanced precision computations. */
038 private final MathContext mContext;
039 /** The number {@code 2}. */
040 private final BigDecimal two;
041 /** The number {@code -1}. */
042 private final BigDecimal minusOne;
043 /** The number {@code 0.5}. */
044 private final BigDecimal oneHalf;
045
046 /**
047 * Default precision is {@link MathContext#DECIMAL128 DECIMAL128}.
048 */
049 public LegendreHighPrecisionRuleFactory() {
050 this(MathContext.DECIMAL128);
051 }
052
053 /**
054 * @param mContext Precision setting for computing the quadrature rules.
055 */
056 public LegendreHighPrecisionRuleFactory(MathContext mContext) {
057 this.mContext = mContext;
058 two = new BigDecimal("2", mContext);
059 minusOne = new BigDecimal("-1", mContext);
060 oneHalf = new BigDecimal("0.5", mContext);
061 }
062
063 /**
064 * {@inheritDoc}
065 *
066 * @throws NotStrictlyPositiveException if {@code numberOfPoints < 1}.
067 */
068 @Override
069 protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints) {
070 if (numberOfPoints <= 0) {
071 throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS,
072 numberOfPoints);
073 }
074
075 if (numberOfPoints == 1) {
076 // Break recursion.
077 return new Pair<BigDecimal[], BigDecimal[]>(new BigDecimal[] { BigDecimal.ZERO },
078 new BigDecimal[] { two });
079 }
080
081 // Get previous rule.
082 // If it has not been computed yet it will trigger a recursive call
083 // to this method.
084 final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
085
086 // Compute next rule.
087 final BigDecimal[] points = new BigDecimal[numberOfPoints];
088 final BigDecimal[] weights = new BigDecimal[numberOfPoints];
089
090 // Find i-th root of P[n+1] by bracketing.
091 final int iMax = numberOfPoints / 2;
092 for (int i = 0; i < iMax; i++) {
093 // Lower-bound of the interval.
094 BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1];
095 // Upper-bound of the interval.
096 BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i];
097 // P[j-1](a)
098 BigDecimal pma = BigDecimal.ONE;
099 // P[j](a)
100 BigDecimal pa = a;
101 // P[j-1](b)
102 BigDecimal pmb = BigDecimal.ONE;
103 // P[j](b)
104 BigDecimal pb = b;
105 for (int j = 1; j < numberOfPoints; j++) {
106 final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
107 final BigDecimal b_j = new BigDecimal(j, mContext);
108 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
109
110 // Compute P[j+1](a)
111 // ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1);
112
113 BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext);
114 tmp1 = pa.multiply(tmp1, mContext);
115 BigDecimal tmp2 = pma.multiply(b_j, mContext);
116 // P[j+1](a)
117 BigDecimal ppa = tmp1.subtract(tmp2, mContext);
118 ppa = ppa.divide(b_j_p_1, mContext);
119
120 // Compute P[j+1](b)
121 // ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1);
122
123 tmp1 = b.multiply(b_two_j_p_1, mContext);
124 tmp1 = pb.multiply(tmp1, mContext);
125 tmp2 = pmb.multiply(b_j, mContext);
126 // P[j+1](b)
127 BigDecimal ppb = tmp1.subtract(tmp2, mContext);
128 ppb = ppb.divide(b_j_p_1, mContext);
129
130 pma = pa;
131 pa = ppa;
132 pmb = pb;
133 pb = ppb;
134 }
135 // Now pa = P[n+1](a), and pma = P[n](a). Same holds for b.
136 // Middle of the interval.
137 BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext);
138 // P[j-1](c)
139 BigDecimal pmc = BigDecimal.ONE;
140 // P[j](c)
141 BigDecimal pc = c;
142 boolean done = false;
143 while (!done) {
144 BigDecimal tmp1 = b.subtract(a, mContext);
145 BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext);
146 done = tmp1.compareTo(tmp2) <= 0;
147 pmc = BigDecimal.ONE;
148 pc = c;
149 for (int j = 1; j < numberOfPoints; j++) {
150 final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
151 final BigDecimal b_j = new BigDecimal(j, mContext);
152 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
153
154 // Compute P[j+1](c)
155 tmp1 = c.multiply(b_two_j_p_1, mContext);
156 tmp1 = pc.multiply(tmp1, mContext);
157 tmp2 = pmc.multiply(b_j, mContext);
158 // P[j+1](c)
159 BigDecimal ppc = tmp1.subtract(tmp2, mContext);
160 ppc = ppc.divide(b_j_p_1, mContext);
161
162 pmc = pc;
163 pc = ppc;
164 }
165 // Now pc = P[n+1](c) and pmc = P[n](c).
166 if (!done) {
167 if (pa.signum() * pc.signum() <= 0) {
168 b = c;
169 pmb = pmc;
170 pb = pc;
171 } else {
172 a = c;
173 pma = pmc;
174 pa = pc;
175 }
176 c = a.add(b, mContext).multiply(oneHalf, mContext);
177 }
178 }
179 final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
180 BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext);
181 tmp1 = tmp1.multiply(nP);
182 tmp1 = tmp1.pow(2, mContext);
183 BigDecimal tmp2 = c.pow(2, mContext);
184 tmp2 = BigDecimal.ONE.subtract(tmp2, mContext);
185 tmp2 = tmp2.multiply(two, mContext);
186 tmp2 = tmp2.divide(tmp1, mContext);
187
188 points[i] = c;
189 weights[i] = tmp2;
190
191 final int idx = numberOfPoints - i - 1;
192 points[idx] = c.negate(mContext);
193 weights[idx] = tmp2;
194 }
195 // If "numberOfPoints" is odd, 0 is a root.
196 // Note: as written, the test for oddness will work for negative
197 // integers too (although it is not necessary here), preventing
198 // a FindBugs warning.
199 if (numberOfPoints % 2 != 0) {
200 BigDecimal pmc = BigDecimal.ONE;
201 for (int j = 1; j < numberOfPoints; j += 2) {
202 final BigDecimal b_j = new BigDecimal(j, mContext);
203 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
204
205 // pmc = -j * pmc / (j + 1);
206 pmc = pmc.multiply(b_j, mContext);
207 pmc = pmc.divide(b_j_p_1, mContext);
208 pmc = pmc.negate(mContext);
209 }
210
211 // 2 / pow(numberOfPoints * pmc, 2);
212 final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
213 BigDecimal tmp1 = pmc.multiply(nP, mContext);
214 tmp1 = tmp1.pow(2, mContext);
215 BigDecimal tmp2 = two.divide(tmp1, mContext);
216
217 points[iMax] = BigDecimal.ZERO;
218 weights[iMax] = tmp2;
219 }
220
221 return new Pair<BigDecimal[], BigDecimal[]>(points, weights);
222 }
223 }