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