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 org.apache.commons.math3.exception.DimensionMismatchException;
020import org.apache.commons.math3.util.Pair;
021
022/**
023 * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
024 * In this implementation, the lower and upper bounds of the natural interval
025 * of integration are -1 and 1, respectively.
026 * The Legendre polynomials are evaluated using the recurrence relation
027 * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">
028 * Abramowitz and Stegun, 1964</a>.
029 *
030 * @since 3.1
031 */
032public class LegendreRuleFactory extends BaseRuleFactory<Double> {
033    /** {@inheritDoc} */
034    @Override
035    protected Pair<Double[], Double[]> computeRule(int numberOfPoints)
036        throws DimensionMismatchException {
037
038        if (numberOfPoints == 1) {
039            // Break recursion.
040            return new Pair<Double[], Double[]>(new Double[] { 0d },
041                                                new Double[] { 2d });
042        }
043
044        // Get previous rule.
045        // If it has not been computed yet it will trigger a recursive call
046        // to this method.
047        final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
048
049        // Compute next rule.
050        final Double[] points = new Double[numberOfPoints];
051        final Double[] weights = new Double[numberOfPoints];
052
053        // Find i-th root of P[n+1] by bracketing.
054        final int iMax = numberOfPoints / 2;
055        for (int i = 0; i < iMax; i++) {
056            // Lower-bound of the interval.
057            double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue();
058            // Upper-bound of the interval.
059            double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue();
060            // P[j-1](a)
061            double pma = 1;
062            // P[j](a)
063            double pa = a;
064            // P[j-1](b)
065            double pmb = 1;
066            // P[j](b)
067            double pb = b;
068            for (int j = 1; j < numberOfPoints; j++) {
069                final int two_j_p_1 = 2 * j + 1;
070                final int j_p_1 = j + 1;
071                // P[j+1](a)
072                final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1;
073                // P[j+1](b)
074                final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1;
075                pma = pa;
076                pa = ppa;
077                pmb = pb;
078                pb = ppb;
079            }
080            // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b).
081            // Middle of the interval.
082            double c = 0.5 * (a + b);
083            // P[j-1](c)
084            double pmc = 1;
085            // P[j](c)
086            double pc = c;
087            boolean done = false;
088            while (!done) {
089                done = b - a <= Math.ulp(c);
090                pmc = 1;
091                pc = c;
092                for (int j = 1; j < numberOfPoints; j++) {
093                    // P[j+1](c)
094                    final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1);
095                    pmc = pc;
096                    pc = ppc;
097                }
098                // Now pc = P[n+1](c) and pmc = P[n](c).
099                if (!done) {
100                    if (pa * pc <= 0) {
101                        b = c;
102                        pmb = pmc;
103                        pb = pc;
104                    } else {
105                        a = c;
106                        pma = pmc;
107                        pa = pc;
108                    }
109                    c = 0.5 * (a + b);
110                }
111            }
112            final double d = numberOfPoints * (pmc - c * pc);
113            final double w = 2 * (1 - c * c) / (d * d);
114
115            points[i] = c;
116            weights[i] = w;
117
118            final int idx = numberOfPoints - i - 1;
119            points[idx] = -c;
120            weights[idx] = w;
121        }
122        // If "numberOfPoints" is odd, 0 is a root.
123        // Note: as written, the test for oddness will work for negative
124        // integers too (although it is not necessary here), preventing
125        // a FindBugs warning.
126        if (numberOfPoints % 2 != 0) {
127            double pmc = 1;
128            for (int j = 1; j < numberOfPoints; j += 2) {
129                pmc = -j * pmc / (j + 1);
130            }
131            final double d = numberOfPoints * pmc;
132            final double w = 2 / (d * d);
133
134            points[iMax] = 0d;
135            weights[iMax] = w;
136        }
137
138        return new Pair<Double[], Double[]>(points, weights);
139    }
140}