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.FastMath;
020import org.apache.commons.math4.util.Pair;
021
022/**
023 * Factory that creates a
024 * <a href="http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature">
025 * Gauss-type quadrature rule using Hermite polynomials</a>
026 * of the first kind.
027 * Such a quadrature rule allows the calculation of improper integrals
028 * of a function
029 * <p>
030 *  \(f(x) e^{-x^2}\)
031 * </p><p>
032 * Recurrence relation and weights computation follow
033 * <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">
034 * Abramowitz and Stegun, 1964</a>.
035 * </p><p>
036 * The coefficients of the standard Hermite polynomials grow very rapidly.
037 * In order to avoid overflows, each Hermite polynomial is normalized with
038 * respect to the underlying scalar product.
039 * The initial interval for the application of the bisection method is
040 * based on the roots of the previous Hermite polynomial (interlacing).
041 * Upper and lower bounds of these roots are provided by </p>
042 * <blockquote>
043 *  I. Krasikov,
044 *  <em>Nonnegative quadratic forms and bounds on orthogonal polynomials</em>,
045 *  Journal of Approximation theory <b>111</b>, 31-49
046 * </blockquote>
047 *
048 * @since 3.3
049 */
050public class HermiteRuleFactory extends BaseRuleFactory<Double> {
051    /** &pi;<sup>1/2</sup> */
052    private static final double SQRT_PI = 1.77245385090551602729;
053    /** &pi;<sup>-1/4</sup> */
054    private static final double H0 = 7.5112554446494248286e-1;
055    /** &pi;<sup>-1/4</sup> &radic;2 */
056    private static final double H1 = 1.0622519320271969145;
057
058    /** {@inheritDoc} */
059    @Override
060    protected Pair<Double[], Double[]> computeRule(int numberOfPoints) {
061        if (numberOfPoints == 1) {
062            // Break recursion.
063            return new Pair<>(new Double[] { 0d },
064                                                new Double[] { SQRT_PI });
065        }
066
067        // Get previous rule.
068        // If it has not been computed yet it will trigger a recursive call
069        // to this method.
070        final int lastNumPoints = numberOfPoints - 1;
071        final Double[] previousPoints = getRuleInternal(lastNumPoints).getFirst();
072
073        // Compute next rule.
074        final Double[] points = new Double[numberOfPoints];
075        final Double[] weights = new Double[numberOfPoints];
076
077        final double sqrtTwoTimesLastNumPoints = FastMath.sqrt(2 * lastNumPoints);
078        final double sqrtTwoTimesNumPoints = FastMath.sqrt(2 * numberOfPoints);
079
080        // Find i-th root of H[n+1] by bracketing.
081        final int iMax = numberOfPoints / 2;
082        for (int i = 0; i < iMax; i++) {
083            // Lower-bound of the interval.
084            double a = (i == 0) ? -sqrtTwoTimesLastNumPoints : previousPoints[i - 1].doubleValue();
085            // Upper-bound of the interval.
086            double b = (iMax == 1) ? -0.5 : previousPoints[i].doubleValue();
087
088            // H[j-1](a)
089            double hma = H0;
090            // H[j](a)
091            double ha = H1 * a;
092            // H[j-1](b)
093            double hmb = H0;
094            // H[j](b)
095            double hb = H1 * b;
096            for (int j = 1; j < numberOfPoints; j++) {
097                // Compute H[j+1](a) and H[j+1](b)
098                final double jp1 = j + 1;
099                final double s = FastMath.sqrt(2 / jp1);
100                final double sm = FastMath.sqrt(j / jp1);
101                final double hpa = s * a * ha - sm * hma;
102                final double hpb = s * b * hb - sm * hmb;
103                hma = ha;
104                ha = hpa;
105                hmb = hb;
106                hb = hpb;
107            }
108
109            // Now ha = H[n+1](a), and hma = H[n](a) (same holds for b).
110            // Middle of the interval.
111            double c = 0.5 * (a + b);
112            // P[j-1](c)
113            double hmc = H0;
114            // P[j](c)
115            double hc = H1 * c;
116            boolean done = false;
117            while (!done) {
118                done = b - a <= Math.ulp(c);
119                hmc = H0;
120                hc = H1 * c;
121                for (int j = 1; j < numberOfPoints; j++) {
122                    // Compute H[j+1](c)
123                    final double jp1 = j + 1;
124                    final double s = FastMath.sqrt(2 / jp1);
125                    final double sm = FastMath.sqrt(j / jp1);
126                    final double hpc = s * c * hc - sm * hmc;
127                    hmc = hc;
128                    hc = hpc;
129                }
130                // Now h = H[n+1](c) and hm = H[n](c).
131                if (!done) {
132                    if (ha * hc < 0) {
133                        b = c;
134                        hmb = hmc;
135                        hb = hc;
136                    } else {
137                        a = c;
138                        hma = hmc;
139                        ha = hc;
140                    }
141                    c = 0.5 * (a + b);
142                }
143            }
144            final double d = sqrtTwoTimesNumPoints * hmc;
145            final double w = 2 / (d * d);
146
147            points[i] = c;
148            weights[i] = w;
149
150            final int idx = lastNumPoints - i;
151            points[idx] = -c;
152            weights[idx] = w;
153        }
154
155        // If "numberOfPoints" is odd, 0 is a root.
156        // Note: as written, the test for oddness will work for negative
157        // integers too (although it is not necessary here), preventing
158        // a FindBugs warning.
159        if (numberOfPoints % 2 != 0) {
160            double hm = H0;
161            for (int j = 1; j < numberOfPoints; j += 2) {
162                final double jp1 = j + 1;
163                hm = -FastMath.sqrt(j / jp1) * hm;
164            }
165            final double d = sqrtTwoTimesNumPoints * hm;
166            final double w = 2 / (d * d);
167
168            points[iMax] = 0d;
169            weights[iMax] = w;
170        }
171
172        return new Pair<>(points, weights);
173    }
174}