    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    }