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    }