View Javadoc
1   /*
2    * Licensed to the Hipparchus project under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.hipparchus.analysis.integration.gauss;
18  
19  import java.util.Arrays;
20  import java.util.SortedMap;
21  import java.util.TreeMap;
22  
23  import org.hipparchus.analysis.UnivariateFunction;
24  import org.hipparchus.exception.LocalizedCoreFormats;
25  import org.hipparchus.exception.MathIllegalArgumentException;
26  import org.hipparchus.util.FastMath;
27  import org.hipparchus.util.Incrementor;
28  import org.hipparchus.util.Pair;
29  
30  /**
31   * Base class for rules that determines the integration nodes and their
32   * weights.
33   * Subclasses must implement the {@link #computeRule(int) computeRule} method.
34   *
35   * @since 2.0
36   */
37  public abstract class AbstractRuleFactory implements RuleFactory {
38  
39      /** List of points and weights, indexed by the order of the rule. */
40      private final SortedMap<Integer, Pair<double[], double[]>> pointsAndWeights = new TreeMap<>();
41  
42      /** Empty constructor.
43       * <p>
44       * This constructor is not strictly necessary, but it prevents spurious
45       * javadoc warnings with JDK 18 and later.
46       * </p>
47       * @since 3.0
48       */
49      public AbstractRuleFactory() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
50          // nothing to do
51      }
52  
53      /** {@inheritDoc} */
54      @Override
55      public Pair<double[], double[]> getRule(int numberOfPoints)
56          throws MathIllegalArgumentException {
57  
58          if (numberOfPoints <= 0) {
59              throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
60                                                     numberOfPoints);
61          }
62          if (numberOfPoints > 1000) {
63              throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
64                                                     numberOfPoints, 1000);
65          }
66  
67          Pair<double[], double[]> rule;
68          synchronized (pointsAndWeights) {
69              // Try to obtain the rule from the cache.
70              rule = pointsAndWeights.get(numberOfPoints);
71  
72              if (rule == null) {
73                  // Rule not computed yet.
74  
75                  // Compute the rule.
76                  rule = computeRule(numberOfPoints);
77  
78                  // Cache it.
79                  pointsAndWeights.put(numberOfPoints, rule);
80              }
81          }
82  
83          // Return a copy.
84          return new Pair<>(rule.getFirst().clone(), rule.getSecond().clone());
85  
86      }
87  
88      /**
89       * Computes the rule for the given order.
90       *
91       * @param numberOfPoints Order of the rule to be computed.
92       * @return the computed rule.
93       * @throws MathIllegalArgumentException if the elements of the pair do not
94       * have the same length.
95       */
96      protected abstract Pair<double[], double[]> computeRule(int numberOfPoints)
97          throws MathIllegalArgumentException;
98  
99      /** Computes roots of the associated orthogonal polynomials.
100      * <p>
101      * The roots are found using the <a href="https://en.wikipedia.org/wiki/Aberth_method">Aberth method</a>.
102      * The guess points for initializing search for degree n are fixed for degrees 1 and 2 and are
103      * selected from n-1 roots of rule n-1 (the two extreme roots are used, plus the n-1 intermediate
104      * points between all roots).
105      * </p>
106      * @param n number of roots to search for
107      * @param ratioEvaluator function evaluating the ratio Pₙ(x)/Pₙ'(x)
108      * @return sorted array of roots
109      */
110     protected double[] findRoots(final int n, final UnivariateFunction ratioEvaluator) {
111 
112         final double[] roots  = new double[n];
113 
114         // set up initial guess
115         if (n == 1) {
116             // arbitrary guess
117             roots[0] = 0;
118         } else if (n == 2) {
119             // arbitrary guess
120             roots[0] = -1;
121             roots[1] = +1;
122         } else {
123 
124             // get roots from previous rule.
125             // If it has not been computed yet it will trigger a recursive call
126             final double[] previousPoints = getRule(n - 1).getFirst();
127 
128             // first guess at previous first root
129             roots[0] = previousPoints[0];
130 
131             // intermediate guesses between previous roots
132             for (int i = 1; i < n - 1; ++i) {
133                 roots[i] = (previousPoints[i - 1] + previousPoints[i]) * 0.5;
134             }
135 
136             // last guess at previous last root
137             roots[n - 1] = previousPoints[n - 2];
138 
139         }
140 
141         // use Aberth method to find all roots simultaneously
142         final double[]    ratio       = new double[n];
143         final Incrementor incrementor = new Incrementor(1000);
144         double            tol;
145         double            maxOffset;
146         do {
147 
148             // safety check that triggers an exception if too much iterations are made
149             incrementor.increment();
150 
151             // find the ratio P(xᵢ)/P'(xᵢ) for all current roots approximations
152             for (int i = 0; i < n; ++i) {
153                 ratio[i] = ratioEvaluator.value(roots[i]);
154             }
155 
156             // move roots approximations all at once, using Aberth method
157             maxOffset = 0;
158             for (int i = 0; i < n; ++i) {
159                 double sum = 0;
160                 for (int j = 0; j < n; ++j) {
161                     if (j != i) {
162                         sum += 1 / (roots[i] - roots[j]);
163                     }
164                 }
165                 final double offset = ratio[i] / (1 - ratio[i] * sum);
166                 maxOffset = FastMath.max(maxOffset, FastMath.abs(offset));
167                 roots[i] -= offset;
168             }
169 
170             // we set tolerance to 1 ulp of the largest root
171             tol = 0;
172             for (final double r : roots) {
173                 tol = FastMath.max(tol, FastMath.ulp(r));
174             }
175 
176         } while (maxOffset > tol);
177 
178         // sort the roots
179         Arrays.sort(roots);
180 
181         return roots;
182 
183     }
184 
185     /** Enforce symmetry of roots.
186      * @param roots roots to process in place
187      */
188     protected void enforceSymmetry(final double[] roots) {
189 
190         final int n = roots.length;
191 
192         // enforce symmetry
193         for (int i = 0; i < n / 2; ++i) {
194             final int idx = n - i - 1;
195             final double c = (roots[i] - roots[idx]) * 0.5;
196             roots[i]   = +c;
197             roots[idx] = -c;
198         }
199 
200         // If n is odd, 0 is a root.
201         // Note: as written, the test for oddness will work for negative
202         // integers too (although it is not necessary here), preventing
203         // a FindBugs warning.
204         if (n % 2 != 0) {
205             roots[n / 2] = 0;
206         }
207 
208     }
209 
210 }