    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 java.math.MathContext;
    020    import java.math.BigDecimal;
    021    import org.apache.commons.math3.util.Pair;
    022    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    023    import org.apache.commons.math3.exception.util.LocalizedFormats;
    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     * @version $Id: LegendreHighPrecisionRuleFactory.java 1374187 2012-08-17 09:50:54Z erans $
    035     */
    036    public class LegendreHighPrecisionRuleFactory extends BaseRuleFactory<BigDecimal> {
    037        /** Settings for enhanced precision computations. */
    038        private final MathContext mContext;
    039        /** The number {@code 2}. */
    040        private final BigDecimal two;
    041        /** The number {@code -1}. */
    042        private final BigDecimal minusOne;
    043        /** The number {@code 0.5}. */
    044        private final BigDecimal oneHalf;
    045
    046        /**
    047         * Default precision is {@link MathContext#DECIMAL128 DECIMAL128}.
    048         */
    049        public LegendreHighPrecisionRuleFactory() {
    050            this(MathContext.DECIMAL128);
    051        }
    052
    053        /**
    054         * @param mContext Precision setting for computing the quadrature rules.
    055         */
    056        public LegendreHighPrecisionRuleFactory(MathContext mContext) {
    057            this.mContext = mContext;
    058            two = new BigDecimal("2", mContext);
    059            minusOne = new BigDecimal("-1", mContext);
    060            oneHalf = new BigDecimal("0.5", mContext);
    061        }
    062
    063        /**
    064         * {@inheritDoc}
    065         *
    066         * @throws NotStrictlyPositiveException if {@code numberOfPoints < 1}.
    067         */
    068        @Override
    069        protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints) {
    070            if (numberOfPoints <= 0) {
    071                throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS,
    072                                                       numberOfPoints);
    073            }
    074
    075            if (numberOfPoints == 1) {
    076                // Break recursion.
    077                return new Pair<BigDecimal[], BigDecimal[]>(new BigDecimal[] { BigDecimal.ZERO },
    078                                                            new BigDecimal[] { two });
    079            }
    080
    081            // Get previous rule.
    082            // If it has not been computed yet it will trigger a recursive call
    083            // to this method.
    084            final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
    085
    086            // Compute next rule.
    087            final BigDecimal[] points = new BigDecimal[numberOfPoints];
    088            final BigDecimal[] weights = new BigDecimal[numberOfPoints];
    089
    090            // Find i-th root of P[n+1] by bracketing.
    091            final int iMax = numberOfPoints / 2;
    092            for (int i = 0; i < iMax; i++) {
    093                // Lower-bound of the interval.
    094                BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1];
    095                // Upper-bound of the interval.
    096                BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i];
    097                // P[j-1](a)
    098                BigDecimal pma = BigDecimal.ONE;
    099                // P[j](a)
    100                BigDecimal pa = a;
    101                // P[j-1](b)
    102                BigDecimal pmb = BigDecimal.ONE;
    103                // P[j](b)
    104                BigDecimal pb = b;
    105                for (int j = 1; j < numberOfPoints; j++) {
    106                    final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
    107                    final BigDecimal b_j = new BigDecimal(j, mContext);
    108                    final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
    109
    110                    // Compute P[j+1](a)
    111                    // ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1);
    112
    113                    BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext);
    114                    tmp1 = pa.multiply(tmp1, mContext);
    115                    BigDecimal tmp2 = pma.multiply(b_j, mContext);
    116                    // P[j+1](a)
    117                    BigDecimal ppa = tmp1.subtract(tmp2, mContext);
    118                    ppa = ppa.divide(b_j_p_1, mContext);
    119
    120                    // Compute P[j+1](b)
    121                    // ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1);
    122
    123                    tmp1 = b.multiply(b_two_j_p_1, mContext);
    124                    tmp1 = pb.multiply(tmp1, mContext);
    125                    tmp2 = pmb.multiply(b_j, mContext);
    126                    // P[j+1](b)
    127                    BigDecimal ppb = tmp1.subtract(tmp2, mContext);
    128                    ppb = ppb.divide(b_j_p_1, mContext);
    129
    130                    pma = pa;
    131                    pa = ppa;
    132                    pmb = pb;
    133                    pb = ppb;
    134                }
    135                // Now pa = P[n+1](a), and pma = P[n](a). Same holds for b.
    136                // Middle of the interval.
    137                BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext);
    138                // P[j-1](c)
    139                BigDecimal pmc = BigDecimal.ONE;
    140                // P[j](c)
    141                BigDecimal pc = c;
    142                boolean done = false;
    143                while (!done) {
    144                    BigDecimal tmp1 = b.subtract(a, mContext);
    145                    BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext);
    146                    done = tmp1.compareTo(tmp2) <= 0;
    147                    pmc = BigDecimal.ONE;
    148                    pc = c;
    149                    for (int j = 1; j < numberOfPoints; j++) {
    150                        final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
    151                        final BigDecimal b_j = new BigDecimal(j, mContext);
    152                        final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
    153
    154                        // Compute P[j+1](c)
    155                        tmp1 = c.multiply(b_two_j_p_1, mContext);
    156                        tmp1 = pc.multiply(tmp1, mContext);
    157                        tmp2 = pmc.multiply(b_j, mContext);
    158                        // P[j+1](c)
    159                        BigDecimal ppc = tmp1.subtract(tmp2, mContext);
    160                        ppc = ppc.divide(b_j_p_1, mContext);
    161
    162                        pmc = pc;
    163                        pc = ppc;
    164                    }
    165                    // Now pc = P[n+1](c) and pmc = P[n](c).
    166                    if (!done) {
    167                        if (pa.signum() * pc.signum() <= 0) {
    168                            b = c;
    169                            pmb = pmc;
    170                            pb = pc;
    171                        } else {
    172                            a = c;
    173                            pma = pmc;
    174                            pa = pc;
    175                        }
    176                        c = a.add(b, mContext).multiply(oneHalf, mContext);
    177                    }
    178                }
    179                final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
    180                BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext);
    181                tmp1 = tmp1.multiply(nP);
    182                tmp1 = tmp1.pow(2, mContext);
    183                BigDecimal tmp2 = c.pow(2, mContext);
    184                tmp2 = BigDecimal.ONE.subtract(tmp2, mContext);
    185                tmp2 = tmp2.multiply(two, mContext);
    186                tmp2 = tmp2.divide(tmp1, mContext);
    187
    188                points[i] = c;
    189                weights[i] = tmp2;
    190
    191                final int idx = numberOfPoints - i - 1;
    192                points[idx] = c.negate(mContext);
    193                weights[idx] = tmp2;
    194            }
    195            // If "numberOfPoints" is odd, 0 is a root.
    196            // Note: as written, the test for oddness will work for negative
    197            // integers too (although it is not necessary here), preventing
    198            // a FindBugs warning.
    199            if (numberOfPoints % 2 != 0) {
    200                BigDecimal pmc = BigDecimal.ONE;
    201                for (int j = 1; j < numberOfPoints; j += 2) {
    202                    final BigDecimal b_j = new BigDecimal(j, mContext);
    203                    final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
    204
    205                    // pmc = -j * pmc / (j + 1);
    206                    pmc = pmc.multiply(b_j, mContext);
    207                    pmc = pmc.divide(b_j_p_1, mContext);
    208                    pmc = pmc.negate(mContext);
    209                }
    210
    211                // 2 / pow(numberOfPoints * pmc, 2);
    212                final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
    213                BigDecimal tmp1 = pmc.multiply(nP, mContext);
    214                tmp1 = tmp1.pow(2, mContext);
    215                BigDecimal tmp2 = two.divide(tmp1, mContext);
    216
    217                points[iMax] = BigDecimal.ZERO;
    218                weights[iMax] = tmp2;
    219            }
    220
    221            return new Pair<BigDecimal[], BigDecimal[]>(points, weights);
    222        }
    223    }