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