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}