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