diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,5 @@ +# Changelog for `numerical-integration` + +## 0.1.0.0 - 2023-08-30 + +First release. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Stéphane Laurent (c) 2023 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Stéphane Laurent nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/LICENSES/LICENSE_Eigen b/LICENSES/LICENSE_Eigen new file mode 100644 --- /dev/null +++ b/LICENSES/LICENSE_Eigen @@ -0,0 +1,26 @@ +/* + Copyright (c) 2011, Intel Corporation. All rights reserved. + + Redistribution and use in source and binary forms, with or without modification, + are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ diff --git a/LICENSES/LICENSE_NumericalIntegration b/LICENSES/LICENSE_NumericalIntegration new file mode 100644 --- /dev/null +++ b/LICENSES/LICENSE_NumericalIntegration @@ -0,0 +1,362 @@ +Mozilla Public License, version 2.0 + +1. Definitions + +1.1. "Contributor" + + means each individual or legal entity that creates, contributes to the + creation of, or owns Covered Software. + +1.2. "Contributor Version" + + means the combination of the Contributions of others (if any) used by a + Contributor and that particular Contributor's Contribution. + +1.3. "Contribution" + + means Covered Software of a particular Contributor. + +1.4. "Covered Software" + + means Source Code Form to which the initial Contributor has attached the + notice in Exhibit A, the Executable Form of such Source Code Form, and + Modifications of such Source Code Form, in each case including portions + thereof. + +1.5. "Incompatible With Secondary Licenses" + means + + a. that the initial Contributor has attached the notice described in + Exhibit B to the Covered Software; or + + b. that the Covered Software was made available under the terms of + version 1.1 or earlier of the License, but not also under the terms of + a Secondary License. + +1.6. "Executable Form" + + means any form of the work other than Source Code Form. + +1.7. "Larger Work" + + means a work that combines Covered Software with other material, in a + separate file or files, that is not Covered Software. + +1.8. "License" + + means this document. + +1.9. "Licensable" + + means having the right to grant, to the maximum extent possible, whether + at the time of the initial grant or subsequently, any and all of the + rights conveyed by this License. + +1.10. "Modifications" + + means any of the following: + + a. any file in Source Code Form that results from an addition to, + deletion from, or modification of the contents of Covered Software; or + + b. any new file in Source Code Form that contains any Covered Software. + +1.11. "Patent Claims" of a Contributor + + means any patent claim(s), including without limitation, method, + process, and apparatus claims, in any patent Licensable by such + Contributor that would be infringed, but for the grant of the License, + by the making, using, selling, offering for sale, having made, import, + or transfer of either its Contributions or its Contributor Version. + +1.12. "Secondary License" + + means either the GNU General Public License, Version 2.0, the GNU Lesser + General Public License, Version 2.1, the GNU Affero General Public + License, Version 3.0, or any later versions of those licenses. + +1.13. "Source Code Form" + + means the form of the work preferred for making modifications. + +1.14. "You" (or "Your") + + means an individual or a legal entity exercising rights under this + License. For legal entities, "You" includes any entity that controls, is + controlled by, or is under common control with You. For purposes of this + definition, "control" means (a) the power, direct or indirect, to cause + the direction or management of such entity, whether by contract or + otherwise, or (b) ownership of more than fifty percent (50%) of the + outstanding shares or beneficial ownership of such entity. + + +2. License Grants and Conditions + +2.1. Grants + + Each Contributor hereby grants You a world-wide, royalty-free, + non-exclusive license: + + a. under intellectual property rights (other than patent or trademark) + Licensable by such Contributor to use, reproduce, make available, + modify, display, perform, distribute, and otherwise exploit its + Contributions, either on an unmodified basis, with Modifications, or + as part of a Larger Work; and + + b. under Patent Claims of such Contributor to make, use, sell, offer for + sale, have made, import, and otherwise transfer either its + Contributions or its Contributor Version. + +2.2. Effective Date + + The licenses granted in Section 2.1 with respect to any Contribution + become effective for each Contribution on the date the Contributor first + distributes such Contribution. + +2.3. Limitations on Grant Scope + + The licenses granted in this Section 2 are the only rights granted under + this License. No additional rights or licenses will be implied from the + distribution or licensing of Covered Software under this License. + Notwithstanding Section 2.1(b) above, no patent license is granted by a + Contributor: + + a. for any code that a Contributor has removed from Covered Software; or + + b. for infringements caused by: (i) Your and any other third party's + modifications of Covered Software, or (ii) the combination of its + Contributions with other software (except as part of its Contributor + Version); or + + c. under Patent Claims infringed by Covered Software in the absence of + its Contributions. + + This License does not grant any rights in the trademarks, service marks, + or logos of any Contributor (except as may be necessary to comply with + the notice requirements in Section 3.4). + +2.4. Subsequent Licenses + + No Contributor makes additional grants as a result of Your choice to + distribute the Covered Software under a subsequent version of this + License (see Section 10.2) or under the terms of a Secondary License (if + permitted under the terms of Section 3.3). + +2.5. Representation + + Each Contributor represents that the Contributor believes its + Contributions are its original creation(s) or it has sufficient rights to + grant the rights to its Contributions conveyed by this License. + +2.6. Fair Use + + This License is not intended to limit any rights You have under + applicable copyright doctrines of fair use, fair dealing, or other + equivalents. + +2.7. Conditions + + Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted in + Section 2.1. + + +3. Responsibilities + +3.1. Distribution of Source Form + + All distribution of Covered Software in Source Code Form, including any + Modifications that You create or to which You contribute, must be under + the terms of this License. You must inform recipients that the Source + Code Form of the Covered Software is governed by the terms of this + License, and how they can obtain a copy of this License. You may not + attempt to alter or restrict the recipients' rights in the Source Code + Form. + +3.2. Distribution of Executable Form + + If You distribute Covered Software in Executable Form then: + + a. such Covered Software must also be made available in Source Code Form, + as described in Section 3.1, and You must inform recipients of the + Executable Form how they can obtain a copy of such Source Code Form by + reasonable means in a timely manner, at a charge no more than the cost + of distribution to the recipient; and + + b. You may distribute such Executable Form under the terms of this + License, or sublicense it under different terms, provided that the + license for the Executable Form does not attempt to limit or alter the + recipients' rights in the Source Code Form under this License. + +3.3. Distribution of a Larger Work + + You may create and distribute a Larger Work under terms of Your choice, + provided that You also comply with the requirements of this License for + the Covered Software. If the Larger Work is a combination of Covered + Software with a work governed by one or more Secondary Licenses, and the + Covered Software is not Incompatible With Secondary Licenses, this + License permits You to additionally distribute such Covered Software + under the terms of such Secondary License(s), so that the recipient of + the Larger Work may, at their option, further distribute the Covered + Software under the terms of either this License or such Secondary + License(s). + +3.4. Notices + + You may not remove or alter the substance of any license notices + (including copyright notices, patent notices, disclaimers of warranty, or + limitations of liability) contained within the Source Code Form of the + Covered Software, except that You may alter any license notices to the + extent required to remedy known factual inaccuracies. + +3.5. Application of Additional Terms + + You may choose to offer, and to charge a fee for, warranty, support, + indemnity or liability obligations to one or more recipients of Covered + Software. However, You may do so only on Your own behalf, and not on + behalf of any Contributor. You must make it absolutely clear that any + such warranty, support, indemnity, or liability obligation is offered by + You alone, and You hereby agree to indemnify every Contributor for any + liability incurred by such Contributor as a result of warranty, support, + indemnity or liability terms You offer. You may include additional + disclaimers of warranty and limitations of liability specific to any + jurisdiction. + +4. Inability to Comply Due to Statute or Regulation + + If it is impossible for You to comply with any of the terms of this License + with respect to some or all of the Covered Software due to statute, + judicial order, or regulation then You must: (a) comply with the terms of + this License to the maximum extent possible; and (b) describe the + limitations and the code they affect. Such description must be placed in a + text file included with all distributions of the Covered Software under + this License. Except to the extent prohibited by statute or regulation, + such description must be sufficiently detailed for a recipient of ordinary + skill to be able to understand it. + +5. Termination + +5.1. The rights granted under this License will terminate automatically if You + fail to comply with any of its terms. However, if You become compliant, + then the rights granted under this License from a particular Contributor + are reinstated (a) provisionally, unless and until such Contributor + explicitly and finally terminates Your grants, and (b) on an ongoing + basis, if such Contributor fails to notify You of the non-compliance by + some reasonable means prior to 60 days after You have come back into + compliance. Moreover, Your grants from a particular Contributor are + reinstated on an ongoing basis if such Contributor notifies You of the + non-compliance by some reasonable means, this is the first time You have + received notice of non-compliance with this License from such + Contributor, and You become compliant prior to 30 days after Your receipt + of the notice. + +5.2. If You initiate litigation against any entity by asserting a patent + infringement claim (excluding declaratory judgment actions, + counter-claims, and cross-claims) alleging that a Contributor Version + directly or indirectly infringes any patent, then the rights granted to + You by any and all Contributors for the Covered Software under Section + 2.1 of this License shall terminate. + +5.3. In the event of termination under Sections 5.1 or 5.2 above, all end user + license agreements (excluding distributors and resellers) which have been + validly granted by You or Your distributors under this License prior to + termination shall survive termination. + +6. Disclaimer of Warranty + + Covered Software is provided under this License on an "as is" basis, + without warranty of any kind, either expressed, implied, or statutory, + including, without limitation, warranties that the Covered Software is free + of defects, merchantable, fit for a particular purpose or non-infringing. + The entire risk as to the quality and performance of the Covered Software + is with You. Should any Covered Software prove defective in any respect, + You (not any Contributor) assume the cost of any necessary servicing, + repair, or correction. This disclaimer of warranty constitutes an essential + part of this License. No use of any Covered Software is authorized under + this License except under this disclaimer. + +7. Limitation of Liability + + Under no circumstances and under no legal theory, whether tort (including + negligence), contract, or otherwise, shall any Contributor, or anyone who + distributes Covered Software as permitted above, be liable to You for any + direct, indirect, special, incidental, or consequential damages of any + character including, without limitation, damages for lost profits, loss of + goodwill, work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses, even if such party shall have been + informed of the possibility of such damages. This limitation of liability + shall not apply to liability for death or personal injury resulting from + such party's negligence to the extent applicable law prohibits such + limitation. Some jurisdictions do not allow the exclusion or limitation of + incidental or consequential damages, so this exclusion and limitation may + not apply to You. + +8. Litigation + + Any litigation relating to this License may be brought only in the courts + of a jurisdiction where the defendant maintains its principal place of + business and such litigation shall be governed by laws of that + jurisdiction, without reference to its conflict-of-law provisions. Nothing + in this Section shall prevent a party's ability to bring cross-claims or + counter-claims. + +9. Miscellaneous + + This License represents the complete agreement concerning the subject + matter hereof. If any provision of this License is held to be + unenforceable, such provision shall be reformed only to the extent + necessary to make it enforceable. Any law or regulation which provides that + the language of a contract shall be construed against the drafter shall not + be used to construe this License against a Contributor. + + +10. Versions of the License + +10.1. New Versions + + Mozilla Foundation is the license steward. Except as provided in Section + 10.3, no one other than the license steward has the right to modify or + publish new versions of this License. Each version will be given a + distinguishing version number. + +10.2. Effect of New Versions + + You may distribute the Covered Software under the terms of the version + of the License under which You originally received the Covered Software, + or under the terms of any subsequent version published by the license + steward. + +10.3. Modified Versions + + If you create software not governed by this License, and you want to + create a new license for such software, you may create and use a + modified version of this License if you rename the license and remove + any references to the name of the license steward (except to note that + such modified license differs from this License). + +10.4. Distributing Source Code Form that is Incompatible With Secondary + Licenses If You choose to distribute Source Code Form that is + Incompatible With Secondary Licenses under the terms of this version of + the License, the notice described in Exhibit B of this License must be + attached. + +Exhibit A - Source Code Form License Notice + + This Source Code Form is subject to the + terms of the Mozilla Public License, v. + 2.0. If a copy of the MPL was not + distributed with this file, You can + obtain one at + http://mozilla.org/MPL/2.0/. + +If it is not possible or desirable to put the notice in a particular file, +then You may include the notice in a location (such as a LICENSE file in a +relevant directory) where a recipient would be likely to look for such a +notice. + +You may add additional accurate notices of copyright ownership. + +Exhibit B - "Incompatible With Secondary Licenses" Notice + + This Source Code Form is "Incompatible + With Secondary Licenses", as defined by + the Mozilla Public License, v. 2.0. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,13 @@ +# numerical-integration + +One-dimensional numerical integration using the 'NumericalIntegration' C++ library. + +___ + +***Example.*** Integrate x² between 0 and 1 with desired relative error 1e-10 and +using 200 subdivisions. + +```haskell +example :: IO (Double, Double, Int) -- value, error estimate, error code +example = integration (\x -> x*x) 0 1 1e-10 200 +``` diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cpp/ComputeGaussKronrodNodesWeights.h b/cpp/ComputeGaussKronrodNodesWeights.h new file mode 100644 --- /dev/null +++ b/cpp/ComputeGaussKronrodNodesWeights.h @@ -0,0 +1,14 @@ +/** +* \file ComputeGaussKronrodNodesWeights.h +* The functions contained in this file calculate the Gauss-Kronrod nodes and weights +* using the Laurie/Gautschi method. +*/ + +#ifndef EIGEN_COMPUTE_GAUSS_KRONROD_NODES_WEIGHTS_H +#define EIGEN_COMPUTE_GAUSS_KRONROD_NODES_WEIGHTS_H + +#include "LaurieGautschi.h" // Stable, slow and requires c++11 tgamma and the latest mpreal.h. +#include "Piessens.h" // Stable for most rules and precisions, faster than LaurieGautschi. +#include "Monegato.h" // Fastest but unstable for rules above 80. + +#endif //EIGEN_COMPUTE_GAUSS_KRONROD_NODES_WEIGHTS_H diff --git a/cpp/Eigen/Eigenvalues b/cpp/Eigen/Eigenvalues new file mode 100644 --- /dev/null +++ b/cpp/Eigen/Eigenvalues @@ -0,0 +1,60 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_EIGENVALUES_MODULE_H +#define EIGEN_EIGENVALUES_MODULE_H + +#include "Core" + +#include "Cholesky" +#include "Jacobi" +#include "Householder" +#include "LU" +#include "Geometry" + +#include "src/Core/util/DisableStupidWarnings.h" + +/** \defgroup Eigenvalues_Module Eigenvalues module + * + * + * + * This module mainly provides various eigenvalue solvers. + * This module also provides some MatrixBase methods, including: + * - MatrixBase::eigenvalues(), + * - MatrixBase::operatorNorm() + * + * \code + * #include + * \endcode + */ + +#include "src/misc/RealSvd2x2.h" +#include "src/Eigenvalues/Tridiagonalization.h" +#include "src/Eigenvalues/RealSchur.h" +#include "src/Eigenvalues/EigenSolver.h" +#include "src/Eigenvalues/SelfAdjointEigenSolver.h" +#include "src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h" +#include "src/Eigenvalues/HessenbergDecomposition.h" +#include "src/Eigenvalues/ComplexSchur.h" +#include "src/Eigenvalues/ComplexEigenSolver.h" +#include "src/Eigenvalues/RealQZ.h" +#include "src/Eigenvalues/GeneralizedEigenSolver.h" +#include "src/Eigenvalues/MatrixBaseEigenvalues.h" +#ifdef EIGEN_USE_LAPACKE +#ifdef EIGEN_USE_MKL +#include "mkl_lapacke.h" +#else +#include "src/misc/lapacke.h" +#endif +#include "src/Eigenvalues/RealSchur_LAPACKE.h" +#include "src/Eigenvalues/ComplexSchur_LAPACKE.h" +#include "src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h" +#endif + +#include "src/Core/util/ReenableStupidWarnings.h" + +#endif // EIGEN_EIGENVALUES_MODULE_H diff --git a/cpp/GaussKronrodNodesWeights.h b/cpp/GaussKronrodNodesWeights.h new file mode 100644 --- /dev/null +++ b/cpp/GaussKronrodNodesWeights.h @@ -0,0 +1,1770 @@ +/** + * \file GaussKronrodNodesWeights.h + * \sa R Piessens, E de Doncker-Kapenger, C Ueberhuber, D Kahaner, QUADPACK, A Subroutine Package + * for Automatic Integration, Springer Verlag, 1983. + */ + +#ifndef EIGEN_QUADRATURE_KRONROD_H +#define EIGEN_QUADRATURE_KRONROD_H + +namespace Eigen +{ + /** + * \brief The abscissae and weights are given for the interval (-1,1). + * Because of symmetry, only the positive abscissae and their + * corresponding weights are given. + * + * \param abscissaeGaussKronrod15 The abscissae of the 15 point kronrod rule. + * \param abscissaeGaussKronrod21 The abscissae of the 21 point kronrod rule. + * \param abscissaeGaussKronrod31 The abscissae of the 31 point kronrod rule. + * \param abscissaeGaussKronrod41 The abscissae of the 41 point kronrod rule. + * \param abscissaeGaussKronrod51 The abscissae of the 51 point kronrod rule. + * \param abscissaeGaussKronrod61 The abscissae of the 61 point kronrod rule. + * \param abscissaeGaussKronrod71 The abscissae of the 71 point kronrod rule. + * \param abscissaeGaussKronrod81 The abscissae of the 81 point kronrod rule. + * \param abscissaeGaussKronrod91 The abscissae of the 91 point kronrod rule. + * \param abscissaeGaussKronrod101 The abscissae of the 101 point kronrod rule. + * \param abscissaeGaussKronrod121 The abscissae of the 121 point kronrod rule. + * \param abscissaeGaussKronrod201 The abscissae of the 201 point kronrod rule. + * + * \param weightsGaussKronrod15 The weights of the 15 point kronrod rule. + * \param weightsGaussKronrod21 The weights of the 21 point kronrod rule. + * \param weightsGaussKronrod31 The weights of the 31 point kronrod rule. + * \param weightsGaussKronrod41 The weights of the 41 point kronrod rule. + * \param weightsGaussKronrod51 The weights of the 51 point kronrod rule. + * \param weightsGaussKronrod61 The weights of the 61 point kronrod rule. + * \param weightsGaussKronrod71 The weights of the 71 point kronrod rule. + * \param weightsGaussKronrod81 The weights of the 81 point kronrod rule. + * \param weightsGaussKronrod91 The weights of the 91 point kronrod rule. + * \param weightsGaussKronrod101 The weights of the 101 point kronrod rule. + * \param weightsGaussKronrod121 The weights of the 121 point kronrod rule. + * \param weightsGaussKronrod201 The weights of the 201 point kronrod rule. + * + * \param abscissaeGauss15 The abscissae of the 7 point gauss rule. + * \param abscissaeGauss21 The abscissae of the 10 point gauss rule. + * \param abscissaeGauss31 The abscissae of the 15 point gauss rule. + * \param abscissaeGauss41 The abscissae of the 20 point gauss rule. + * \param abscissaeGauss51 The abscissae of the 25 point gauss rule. + * \param abscissaeGauss61 The abscissae of the 30 point gauss rule. + * \param abscissaeGauss71 The abscissae of the 35 point gauss rule. + * \param abscissaeGauss81 The abscissae of the 40 point gauss rule. + * \param abscissaeGauss91 The abscissae of the 45 point gauss rule. + * \param abscissaeGauss101 The abscissae of the 50 point gauss rule. + * \param abscissaeGauss121 The abscissae of the 60 point gauss rule. + * \param abscissaeGauss201 The abscissae of the 100 point gauss rule. + * + * \param weightsGauss15 The weights of the 7 point gauss rule. + * \param weightsGauss21 The weights of the 10 point gauss rule. + * \param weightsGauss31 The weights of the 15 point gauss rule. + * \param weightsGauss41 The weights of the 20 point gauss rule. + * \param weightsGauss51 The weights of the 25 point gauss rule. + * \param weightsGauss61 The weights of the 30 point gauss rule. + * \param weightsGauss71 The weights of the 35 point gauss rule. + * \param weightsGauss81 The weights of the 40 point gauss rule. + * \param weightsGauss91 The weights of the 45 point gauss rule. + * \param weightsGauss101 The weights of the 50 point gauss rule. + * \param weightsGauss121 The weights of the 60 point gauss rule. + * \param weightsGauss201 The weights of the 100 point gauss rule. + */ + + template + class QuadratureKronrod + { + public: + static Array abscissaeGaussKronrod15; + static Array weightsGaussKronrod15; + static Array abscissaeGauss15; + static Array weightsGauss15; + + static Array abscissaeGaussKronrod21; + static Array weightsGaussKronrod21; + static Array abscissaeGauss21; + static Array weightsGauss21; + + static Array abscissaeGaussKronrod31; + static Array weightsGaussKronrod31; + static Array abscissaeGauss31; + static Array weightsGauss31; + + static Array abscissaeGaussKronrod41; + static Array weightsGaussKronrod41; + static Array abscissaeGauss41; + static Array weightsGauss41; + + static Array abscissaeGaussKronrod51; + static Array weightsGaussKronrod51; + static Array abscissaeGauss51; + static Array weightsGauss51; + + static Array abscissaeGaussKronrod61; + static Array weightsGaussKronrod61; + static Array abscissaeGauss61; + static Array weightsGauss61; + + static Array abscissaeGaussKronrod71; + static Array weightsGaussKronrod71; + static Array abscissaeGauss71; + static Array weightsGauss71; + + static Array abscissaeGaussKronrod81; + static Array weightsGaussKronrod81; + static Array abscissaeGauss81; + static Array weightsGauss81; + + static Array abscissaeGaussKronrod91; + static Array weightsGaussKronrod91; + static Array abscissaeGauss91; + static Array weightsGauss91; + + static Array abscissaeGaussKronrod101; + static Array weightsGaussKronrod101; + static Array abscissaeGauss101; + static Array weightsGauss101; + + static Array abscissaeGaussKronrod121; + static Array weightsGaussKronrod121; + static Array abscissaeGauss121; + static Array weightsGauss121; + + static Array abscissaeGaussKronrod201; + static Array weightsGaussKronrod201; + static Array abscissaeGauss201; + static Array weightsGauss201; + + static void computeNodesAndWeights() + { + QuadratureKronrod::computeForRule<7>(abscissaeGaussKronrod15, weightsGaussKronrod15, abscissaeGauss15, weightsGauss15); + QuadratureKronrod::computeForRule<10>(abscissaeGaussKronrod21, weightsGaussKronrod21, abscissaeGauss21, weightsGauss21); + QuadratureKronrod::computeForRule<15>(abscissaeGaussKronrod31, weightsGaussKronrod31, abscissaeGauss31, weightsGauss31); + QuadratureKronrod::computeForRule<20>(abscissaeGaussKronrod41, weightsGaussKronrod41, abscissaeGauss41, weightsGauss41); + QuadratureKronrod::computeForRule<25>(abscissaeGaussKronrod51, weightsGaussKronrod51, abscissaeGauss51, weightsGauss51); + QuadratureKronrod::computeForRule<30>(abscissaeGaussKronrod61, weightsGaussKronrod61, abscissaeGauss61, weightsGauss61); + QuadratureKronrod::computeForRule<35>(abscissaeGaussKronrod71, weightsGaussKronrod71, abscissaeGauss71, weightsGauss71); + QuadratureKronrod::computeForRule<40>(abscissaeGaussKronrod81, weightsGaussKronrod81, abscissaeGauss81, weightsGauss81); + QuadratureKronrod::computeForRule<45>(abscissaeGaussKronrod91, weightsGaussKronrod91, abscissaeGauss91, weightsGauss91); + QuadratureKronrod::computeForRule<50>(abscissaeGaussKronrod101, weightsGaussKronrod101, abscissaeGauss101, weightsGauss101); + QuadratureKronrod::computeForRule<60>(abscissaeGaussKronrod121, weightsGaussKronrod121, abscissaeGauss121, weightsGauss121); + QuadratureKronrod::computeForRule<100>(abscissaeGaussKronrod201, weightsGaussKronrod201, abscissaeGauss201, weightsGauss201); + } + + template + static void computeForRule(Array& kronrodAbscissae, Array& kronrodWeights, + Array& gaussAbscissae, Array& gaussWeights) + { + Eigen::Array xGK; + Eigen::Array wGK; + Eigen::Array xG; + Eigen::Array wG; + + Eigen::LaurieGautschi::computeAbscissaeAndWeights(N,xGK,wGK,xG,wG); + //Eigen::Piessens::computeAbscissaeAndWeights(N,xGK,wGK,xG,wG); + //Eigen::Monegato::computeAbscissaeAndWeights(N,xGK,wGK,xG,wG); + + for(size_t i=0; i +Array QuadratureKronrod::abscissaeGaussKronrod15 = + (Array() << + 0.99145537112081263920685469752632851664204433837033, + 0.94910791234275852452618968404785126240077093767062, + 0.86486442335976907278971278864092620121097230707409, + 0.74153118559939443986386477328078840707414764714139, + 0.58608723546769113029414483825872959843678075060436, + 0.40584515137739716690660641207696146334738201409937, + 0.20778495500789846760068940377324491347978440714517, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod15 = + (Array() << + 0.02293532201052922496373200805896959199356081127575, + 0.06309209262997855329070066318920428666507115721155, + 0.10479001032225018383987632254151801744375665421383, + 0.14065325971552591874518959051023792039988975724800, + 0.16900472663926790282658342659855028410624490030294, + 0.19035057806478540991325640242101368282607807545536, + 0.20443294007529889241416199923464908471651760418072, + 0.20948214108472782801299917489171426369776208022370 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss15 = + (Array() << + 0.94910791234275852452618968404785126240077093767062, + 0.74153118559939443986386477328078840707414764714139, + 0.40584515137739716690660641207696146334738201409937, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss15 = + (Array() << + 0.12948496616886969327061143267908201832858740225995, + 0.27970539148927666790146777142377958248692506522660, + 0.38183005050511894495036977548897513387836508353386, + 0.41795918367346938775510204081632653061224489795918 + ).finished(); + +// Nodes and Weights - Rule 21 +template +Array QuadratureKronrod::abscissaeGaussKronrod21 = + (Array() << + 0.99565716302580808073552728068900284792126058721948, + 0.97390652851717172007796401208445205342826994669238, + 0.93015749135570822600120718005950834622516790998194, + 0.86506336668898451073209668842349304852754301496533, + 0.78081772658641689706371757834504237716340752029816, + 0.67940956829902440623432736511487357576929471183481, + 0.56275713466860468333900009927269414084301388194197, + 0.43339539412924719079926594316578416220007183765625, + 0.29439286270146019813112660310386556616268662515696, + 0.14887433898163121088482600112971998461756485942069, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod21 = + (Array() << + 0.01169463886737187427806439606219204839621733248193, + 0.03255816230796472747881897245938976061738893984566, + 0.05475589657435199603138130024458017637372111405833, + 0.07503967481091995276704314091619000939521938200091, + 0.09312545458369760553506546508336634439001882888076, + 0.10938715880229764189921059032580496027181329983435, + 0.12349197626206585107795810983107415951230034952865, + 0.13470921731147332592805400177170683276099191300856, + 0.14277593857706008079709427313871706088597905653191, + 0.14773910490133849137484151597206804552373162548521, + 0.14944555400291690566493646838982120374523631668747 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss21 = + (Array() << + 0.97390652851717172007796401208445205342826994669238, + 0.86506336668898451073209668842349304852754301496533, + 0.67940956829902440623432736511487357576929471183481, + 0.43339539412924719079926594316578416220007183765625, + 0.14887433898163121088482600112971998461756485942069 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss21 = + (Array() << + 0.06667134430868813759356880989333179285786483432016, + 0.14945134915058059314577633965769733240255663966943, + 0.21908636251598204399553493422816319245877187052268, + 0.26926671930999635509122692156946935285975993846088, + 0.29552422471475287017389299465133832942104671702685 + ).finished(); + +// Nodes and Weights - Rule 31 +template +Array QuadratureKronrod::abscissaeGaussKronrod31 = + (Array() << + 0.99800229869339706028517284015227120907340644231556, + 0.98799251802048542848956571858661258114697281712376, + 0.96773907567913913425734797878433722528335733730013, + 0.93727339240070590430775894771020947124399627351530, + 0.89726453234408190088250965645449588283177871149443, + 0.84820658341042721620064832077421685136625617473699, + 0.79041850144246593296764929481794734686214051995698, + 0.72441773136017004741618605461393800963089929458410, + 0.65099674129741697053373589531327469254694822609260, + 0.57097217260853884753722673725391064123838639628275, + 0.48508186364023968069365574023235061286633893089407, + 0.39415134707756336989720737098104546836275277615870, + 0.29918000715316881216678002426638896266160338274382, + 0.20119409399743452230062830339459620781283645446264, + 0.10114206691871749902707423144739233878745105740164, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod31 = + (Array() << + 0.00537747987292334898779205143012764981830804024313, + 0.01500794732931612253837476307580726809463943643739, + 0.02546084732671532018687400101965335939727174504686, + 0.03534636079137584622203794847836004812263067899242, + 0.04458975132476487660822729937327969022325664966792, + 0.05348152469092808726534314723943029677155476094712, + 0.06200956780067064028513923096080293219040000421033, + 0.06985412131872825870952007709914747578604543514067, + 0.07684968075772037889443277748265900672210910116795, + 0.08308050282313302103828924728610378960155418825337, + 0.08856444305621177064727544369377430321226673269066, + 0.09312659817082532122548687274734571856192788132132, + 0.09664272698362367850517990762758933513665656863050, + 0.09917359872179195933239317348460313105956726081671, + 0.10076984552387559504494666261756972191634838013536, + 0.10133000701479154901737479276749254677092627259660 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss31 = + (Array() << + 0.98799251802048542848956571858661258114697281712376, + 0.93727339240070590430775894771020947124399627351530, + 0.84820658341042721620064832077421685136625617473699, + 0.72441773136017004741618605461393800963089929458410, + 0.57097217260853884753722673725391064123838639628275, + 0.39415134707756336989720737098104546836275277615870, + 0.20119409399743452230062830339459620781283645446264, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss31 = + (Array() << + 0.03075324199611726835462839357720441772174814483343, + 0.07036604748810812470926741645066733846670803275433, + 0.10715922046717193501186954668586930341554371575810, + 0.13957067792615431444780479451102832252085027531551, + 0.16626920581699393355320086048120881113090018009841, + 0.18616100001556221102680056186642282450622601227793, + 0.19843148532711157645611832644383932481869255995754, + 0.20257824192556127288062019996751931483866215800948 + ).finished(); + +// Nodes and Weights - Rule 41 +template +Array QuadratureKronrod::abscissaeGaussKronrod41 = + (Array() << + 0.99885903158827766383831557654586300999957020432630, + 0.99312859918509492478612238847132027822264713090166, + 0.98150787745025025919334299472021694456725093981024, + 0.96397192727791379126766613119727722191206032780619, + 0.94082263383175475351998272221244338027429557377965, + 0.91223442825132590586775244120329811304918479742369, + 0.87827681125228197607744299511307846671124526828251, + 0.83911697182221882339452906170152068532962936506564, + 0.79504142883755119835063883327278794295938959911578, + 0.74633190646015079261430507035564159031073067956918, + 0.69323765633475138480549071184593153338642585141021, + 0.63605368072651502545283669622628593674338911679937, + 0.57514044681971031534294603658642513281381264014772, + 0.51086700195082709800436405095525099842549132920243, + 0.44359317523872510319999221349264010784010101082300, + 0.37370608871541956067254817702492723739574632170568, + 0.30162786811491300432055535685859226061539650501373, + 0.22778585114164507808049619536857462474308893768293, + 0.15260546524092267550522024102267752791167622481842, + 0.07652652113349733375464040939883821100479626681350, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod41 = + (Array() << + 0.00307358371852053150121829324603098748803350468825, + 0.00860026985564294219866178795010234725212892276671, + 0.01462616925697125298378796030886835616388105016225, + 0.02038837346126652359801023143275470512283862794019, + 0.02588213360495115883450506709615314299947911804867, + 0.03128730677703279895854311932380073788776928036281, + 0.03660016975820079803055724070721100848745349674750, + 0.04166887332797368626378830593689473804396084315301, + 0.04643482186749767472023188092610751684212707100708, + 0.05094457392372869193270767005034494866483636580926, + 0.05519510534828599474483237241977732919475345622815, + 0.05911140088063957237496722064859421713641936597704, + 0.06265323755478116802587012217425498058581974469890, + 0.06583459713361842211156355696939794314722350634338, + 0.06864867292852161934562341188536780171548970495824, + 0.07105442355344406830579036172321016741291215932221, + 0.07303069033278666749518941765891311276062684523455, + 0.07458287540049918898658141836248752861611649357209, + 0.07570449768455667465954277537661655826336315590041, + 0.07637786767208073670550283503806100180080103676495, + 0.07660071191799965644504990153010174082793250062867 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss41 = + (Array() << + 0.99312859918509492478612238847132027822264713090166, + 0.96397192727791379126766613119727722191206032780619, + 0.91223442825132590586775244120329811304918479742369, + 0.83911697182221882339452906170152068532962936506564, + 0.74633190646015079261430507035564159031073067956918, + 0.63605368072651502545283669622628593674338911679937, + 0.51086700195082709800436405095525099842549132920243, + 0.37370608871541956067254817702492723739574632170568, + 0.22778585114164507808049619536857462474308893768293, + 0.07652652113349733375464040939883821100479626681350 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss41 = + (Array() << + 0.01761400713915211831186196235185281636214310554334, + 0.04060142980038694133103995227493210987909063998995, + 0.06267204833410906356950653518704160635160107657844, + 0.08327674157670474872475814322204620610017782858316, + 0.10193011981724043503675013548034987616669165602339, + 0.11819453196151841731237737771138228700504121954897, + 0.13168863844917662689849449974816313491611051114698, + 0.14209610931838205132929832506716493303451541339202, + 0.14917298647260374678782873700196943669267990408137, + 0.15275338713072585069808433195509759349194864511238 + ).finished(); + +// Nodes and Weights - Rule 51 +template +Array QuadratureKronrod::abscissaeGaussKronrod51 = + (Array() << + 0.99926210499260983419345748654034059370452496042280, + 0.99555696979049809790878494689390161725756264940481, + 0.98803579453407724763733101457740622707248415209161, + 0.97666392145951751149831538647959406774537055531441, + 0.96161498642584251241813003366016724169212642963710, + 0.94297457122897433941401116965847053190520157060899, + 0.92074711528170156174634608454633063157457035996277, + 0.89499199787827536885104200678280495417455484975358, + 0.86584706529327559544899696958834008820284409402824, + 0.83344262876083400142102110869356956946096411382352, + 0.79787379799850005941041090499430656940863230009338, + 0.75925926303735763057728286520436097638752201889833, + 0.71776640681308438818665407977329778059771167555516, + 0.67356636847346836448512063324762217588341672807275, + 0.62681009901031741278812268162451788101954628995069, + 0.57766293024122296772368984161265406739573503929152, + 0.52632528433471918259962377815801017803683252320191, + 0.47300273144571496052218211500919204133181773846163, + 0.41788538219303774885181439459457248709336998140070, + 0.36117230580938783773582173012764066742207834704338, + 0.30308953893110783016747890998033932920041937876655, + 0.24386688372098843204519036279745158640563315632598, + 0.18371893942104889201596988875952841578528447834991, + 0.12286469261071039638735981880803680553220534604978, + 0.06154448300568507888654639236679663128172434803982, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod51 = + (Array() << + 0.00198738389233031592650785188284340988942998042825, + 0.00556193213535671375804023690106552207017692954963, + 0.00947397338617415160720771052365532387164532684837, + 0.01323622919557167481365640584697623807757808499786, + 0.01684781770912829823151666753633631584040265462471, + 0.02043537114588283545656829223593897367875800609767, + 0.02400994560695321622009248916488108139293152820966, + 0.02747531758785173780294845551781107861479601328871, + 0.03079230016738748889110902021522858560087716239329, + 0.03400213027432933783674879522955120322567052825005, + 0.03711627148341554356033062536761987599599780268805, + 0.04008382550403238207483928446707564640141054926659, + 0.04287284502017004947689579243949516110199950419988, + 0.04550291304992178890987058475266039304370776893570, + 0.04798253713883671390639225575691475498359220742327, + 0.05027767908071567196332525943344008444058763060478, + 0.05236288580640747586436671213787271488735155072371, + 0.05425112988854549014454337045987560682607683844126, + 0.05595081122041231730824068638274734682027103511277, + 0.05743711636156783285358269393950647199483285682390, + 0.05868968002239420796197417585678776413979564625483, + 0.05972034032417405997909929193256185383536304547619, + 0.06053945537604586294536026751756542716231236571046, + 0.06112850971705304830585903041629271192267855232196, + 0.06147118987142531666154413196526417758653796287689, + 0.06158081806783293507875982424006455319043693690314 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss51 = + (Array() << + 0.99555696979049809790878494689390161725756264940481, + 0.97666392145951751149831538647959406774537055531441, + 0.94297457122897433941401116965847053190520157060899, + 0.89499199787827536885104200678280495417455484975358, + 0.83344262876083400142102110869356956946096411382352, + 0.75925926303735763057728286520436097638752201889833, + 0.67356636847346836448512063324762217588341672807275, + 0.57766293024122296772368984161265406739573503929152, + 0.47300273144571496052218211500919204133181773846163, + 0.36117230580938783773582173012764066742207834704338, + 0.24386688372098843204519036279745158640563315632598, + 0.12286469261071039638735981880803680553220534604978, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss51 = + (Array() << + 0.01139379850102628794790296411323477360332052629291, + 0.02635498661503213726190181529529914493596328170332, + 0.04093915670130631265562348771164595366084578336410, + 0.05490469597583519192593689154047332416010998555311, + 0.06803833381235691720718718565670796855470949435464, + 0.08014070033500101801323495966911130229022573285368, + 0.09102826198296364981149722070289165338099255895933, + 0.10053594906705064420220689039268582698846609452814, + 0.10851962447426365311609395705011661934007758798672, + 0.11485825914571164833932554586955580864093619166818, + 0.11945576353578477222817812651290104739017670141373, + 0.12224244299031004168895951894585150583505924756306, + 0.12317605372671545120390287307905014243823362751815 + ).finished(); + +// Nodes and Weights - Rule 61 +template +Array QuadratureKronrod::abscissaeGaussKronrod61 = + (Array() << + 0.99948441005049063757132589570581081946887394701851, + 0.99689348407464954027163005091869528334088203811775, + 0.99163099687040459485862836610948572485050033374616, + 0.98366812327974720997003258160566280194031785470971, + 0.97311632250112626837469386842370688488763796428344, + 0.96002186496830751221687102558179766293035921740392, + 0.94437444474855997941583132403743912158564371496498, + 0.92620004742927432587932427708047400408647453682533, + 0.90557330769990779854652255892595831956897536366223, + 0.88256053579205268154311646253022559005668914714648, + 0.85720523354606109895865851065894385682080017062360, + 0.82956576238276839744289811973250191643906869617034, + 0.79972783582183908301366894232268324073569842937778, + 0.76777743210482619491797734097450313169488361723291, + 0.73379006245322680472617113136952764566938172775469, + 0.69785049479331579693229238802664006838235380065395, + 0.66006106412662696137005366814927075303835037480883, + 0.62052618298924286114047755643118929920736469282953, + 0.57934523582636169175602493217254049590705158881215, + 0.53662414814201989926416979331107279416417800693030, + 0.49248046786177857499369306120770879564426564096319, + 0.44703376953808917678060990032285400016240759386142, + 0.40040125483039439253547621154266063361104593297078, + 0.35270472553087811347103720708937386065363100802143, + 0.30407320227362507737267710719925655353115778980946, + 0.25463692616788984643980512981780510788278930330252, + 0.20452511668230989143895767100202470952410426459556, + 0.15386991360858354696379467274325592041855197124434, + 0.10280693796673703014709675131800059247190133296516, + 0.05147184255531769583302521316672257374914145366657, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod61 = + (Array() << + 0.00138901369867700762455159122675969968104884129196, + 0.00389046112709988405126720184451550327851514298489, + 0.00663070391593129217331982636975016813362838821778, + 0.00927327965951776342844114689202436042127002493819, + 0.01182301525349634174223289885325059289626440625061, + 0.01436972950704580481245143244358001019584189989500, + 0.01692088918905327262757228942032209236856670378384, + 0.01941414119394238117340895105012845585142101419143, + 0.02182803582160919229716748573833899340150729605683, + 0.02419116207808060136568637072523202676039137782818, + 0.02650995488233310161060170933507541436651757952275, + 0.02875404876504129284397878535433421114467916054207, + 0.03090725756238776247288425294309227263527045852381, + 0.03298144705748372603181419101685392751059929121386, + 0.03497933802806002413749967073146787509722691279482, + 0.03688236465182122922391106561713596773695516478103, + 0.03867894562472759295034865153228105025092362982155, + 0.04037453895153595911199527975246811421612606212603, + 0.04196981021516424614714754128596975779008865671899, + 0.04345253970135606931683172811707325807460330863170, + 0.04481480013316266319235555161672324375743139279637, + 0.04605923827100698811627173555937358059469287557182, + 0.04718554656929915394526147818109948648288480730063, + 0.04818586175708712914077949229830459260579923610843, + 0.04905543455502977888752816536723817360588740529530, + 0.04979568342707420635781156937994232853920960281370, + 0.05040592140278234684089308565358502890219701825162, + 0.05088179589874960649229747304980469185338491426092, + 0.05122154784925877217065628260494420825114695242525, + 0.05142612853745902593386287921578125982955203486240, + 0.05149472942945156755834043364709930753273688039646 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss61 = + (Array() << + 0.99689348407464954027163005091869528334088203811775, + 0.98366812327974720997003258160566280194031785470971, + 0.96002186496830751221687102558179766293035921740392, + 0.92620004742927432587932427708047400408647453682533, + 0.88256053579205268154311646253022559005668914714648, + 0.82956576238276839744289811973250191643906869617034, + 0.76777743210482619491797734097450313169488361723291, + 0.69785049479331579693229238802664006838235380065395, + 0.62052618298924286114047755643118929920736469282953, + 0.53662414814201989926416979331107279416417800693030, + 0.44703376953808917678060990032285400016240759386142, + 0.35270472553087811347103720708937386065363100802143, + 0.25463692616788984643980512981780510788278930330252, + 0.15386991360858354696379467274325592041855197124434, + 0.05147184255531769583302521316672257374914145366657 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss61 = + (Array() << + 0.00796819249616660561546588347467362245048069658715, + 0.01846646831109095914230213191204726909620653396818, + 0.02878470788332336934971917961129204363958889454629, + 0.03879919256962704959680193644634769203320097676640, + 0.04840267283059405290293814042280751781527180919737, + 0.05749315621761906648172168940205612879712067072176, + 0.06597422988218049512812851511596236123744295365666, + 0.07375597473770520626824385002219073415377052603705, + 0.08075589522942021535469493846052973087589280370844, + 0.08689978720108297980238753071512570257675332874355, + 0.09212252223778612871763270708761876719691323441823, + 0.09636873717464425963946862635180986509640646143016, + 0.09959342058679526706278028210356947652986926366670, + 0.10176238974840550459642895216855404463270628948713, + 0.10285265289355884034128563670541504386837555706493 + ).finished(); + +// Nodes and Weights - Rule 71 +template +Array QuadratureKronrod::abscissaeGaussKronrod71 = + (Array() << + 0.99961929856058782041953501908431046312122344098794, + 0.99770656909960029726016313931209472519327463717101, + 0.99382029303890921225840752201784028500756549314159, + 0.98793576444385149803511708918548598531300829571417, + 0.98013165785134098801985990365260988816688816643335, + 0.97043761603922983321507048258477014252488123931238, + 0.95883869699584307892126219427392217538678025386427, + 0.94534514820782732953872598552997464486281125124787, + 0.93000375305070609922503535864064081717372299765533, + 0.91285426135931761446493706355576423048806101362556, + 0.89391630583904940482411214013268317995163371208619, + 0.87321912502522233152328234914138526664433890227937, + 0.85081354468109158704203280179571293221562523873360, + 0.82674989909222540683405061274855757904689061798231, + 0.80106721312570571470358184335837821672915468266178, + 0.77381025228691255526742300920988748101181007297155, + 0.74503897566640677164430817613858088537216200738750, + 0.71481450155662878326440863122444674696204396852173, + 0.68319041848815657677750427209465627245583162169621, + 0.65022436466589038867579280898455926336976232592837, + 0.61598571048722183053960576573763958060437543282625, + 0.58054534474976450993450200818969007093207020356433, + 0.54396835169625813820629328275772195267411067648783, + 0.50632277324148861502429755583732737658608793863637, + 0.46768618346152964902238333071075817099648473017136, + 0.42813754151781425418762061300147985878471844648853, + 0.38775069602784231260916836932812549811185119580936, + 0.34660155443081394587697983493023823180840876937729, + 0.30477400147105037961997631656077667943504232291863, + 0.26235294120929605797089520045558051348630439864578, + 0.21941825841501800318906027538425963000063532656757, + 0.17605106116598956997430365644505992967542971232976, + 0.13233927061341662561114208274719221355992727169910, + 0.08837134327565926360092943349754880342411743476510, + 0.04423040796047631902497910724648175796711437032383, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod71 = + (Array() << + 0.00102550911074666801006956437298132913634921030706, + 0.00287226001447070188175126346596935559720130883926, + 0.00489809089031614709134489995355345719632799389497, + 0.00685548721878420013489017671059314269258714775681, + 0.00874803476789701226762937017103524588617649633070, + 0.01064412676080364548520974626298521833239102927401, + 0.01255213863161942856328562549320926369826006957553, + 0.01442614862529363429700928603045290106973954647845, + 0.01624977199984979251411432365271640468771985708338, + 0.01804665112955870366183510500927265777645403257829, + 0.01982463073192568430864633416826825982695039151445, + 0.02156072900282074084433582290901980023837200422120, + 0.02324181089546663502483166236397397407862232471778, + 0.02487938986497896142316498112220908754540096768188, + 0.02647872983924452154290308599311750111721575798899, + 0.02802485927048032502869308406000651692913135080711, + 0.02950731294048380503159847508360512949719541602423, + 0.03093298569089254124584577413599566869200123205647, + 0.03230574967486032561801474409387789273907226946487, + 0.03361454962779494244108331458137284946324198895395, + 0.03485077628981657862883631795167583760136994493590, + 0.03601932106443251160286472895539991131861935891999, + 0.03712347803674947651991054179868437562892996972235, + 0.03815455393845179802844670519848738586957355006794, + 0.03910531516466664309273221987791124541581708553489, + 0.03997983486093488779275160637760618088208667334682, + 0.04078134475859293797125394533241037525076355046502, + 0.04150279114110496636106397836894150439113925391415, + 0.04213802289742381614517983543503522074367479896450, + 0.04269093484449388936791278778195848585871259590150, + 0.04316504612011059575016094844755100917095935153618, + 0.04355454541697310298508166546078403157623950944932, + 0.04385415492459730857308029613026789252900261518504, + 0.04406798346693598707683395305061859923249071421063, + 0.04420009752589896947154443219235054970162309288894, + 0.04424566572105622843217879601220347725717507881493 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss71 = + (Array() << + 0.99770656909960029726016313931209472519327463717101, + 0.98793576444385149803511708918548598531300829571417, + 0.97043761603922983321507048258477014252488123931238, + 0.94534514820782732953872598552997464486281125124787, + 0.91285426135931761446493706355576423048806101362556, + 0.87321912502522233152328234914138526664433890227937, + 0.82674989909222540683405061274855757904689061798231, + 0.77381025228691255526742300920988748101181007297155, + 0.71481450155662878326440863122444674696204396852173, + 0.65022436466589038867579280898455926336976232592837, + 0.58054534474976450993450200818969007093207020356433, + 0.50632277324148861502429755583732737658608793863637, + 0.42813754151781425418762061300147985878471844648853, + 0.34660155443081394587697983493023823180840876937729, + 0.26235294120929605797089520045558051348630439864578, + 0.17605106116598956997430365644505992967542971232976, + 0.08837134327565926360092943349754880342411743476510, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss71 = + (Array() << + 0.00588343342044308497575389624011262403075431779823, + 0.01365082834836149226640400292051638397399403784142, + 0.02132297991148358088343798396620507832380964471994, + 0.02882926010889425404871603971448491626152772662959, + 0.03611011586346338053271696964754993830853577896999, + 0.04310842232617021878230645937490820342280658716399, + 0.04976937040135352980519967608499495859449495319457, + 0.05604081621237012857832774716510096444034853490456, + 0.06187367196608018888701413878868865226537122580410, + 0.06722228526908690396430550874814855854709311914960, + 0.07204479477256006466546190978527775614726739302743, + 0.07630345715544205353865853788422622552708008189564, + 0.07996494224232426293266208098504575632883549943009, + 0.08300059372885658837992652821617700348517450381345, + 0.08538665339209912522594398739111756764119098828375, + 0.08710444699718353424332203160554093668505671509385, + 0.08814053043027546297073880759309664142545431727449, + 0.08848679490710429063820738777761566918240915001061 + ).finished(); + +// Nodes and Weights - Rule 81 +template +Array QuadratureKronrod::abscissaeGaussKronrod81 = + (Array() << + 0.99970755925870001652122454213597471358259546611434, + 0.99823770971055920034962270242058649233577038159505, + 0.99525057344607275036560954016723640396135078520144, + 0.99072623869945700645305435222137215496222208135109, + 0.98472283986425006102933341488331281956810416420996, + 0.97725994998377426266337028371290380697866793203799, + 0.96832312685414990090376748802208053929226225295460, + 0.95791681921379165580454099945275928509488349060275, + 0.94607183716250003820183480835310183364449049027233, + 0.93281280827867653336085216684520571643475357528269, + 0.91814954307289887682909109264899900900889565336112, + 0.90209880696887429672825333086849310358448808105767, + 0.88469200870108974596916731568481615247744988567172, + 0.86595950321225950382078180835461996357054655301110, + 0.84592398558731071742075270176063153327242355905011, + 0.82461223083331166319632023066609877390724038424299, + 0.80206056614025212716548248052306700513426840614442, + 0.77830565142651938769497154550649484802069131612689, + 0.75337980343894219817195210861389792731211452767941, + 0.72731825518992710328099645175493054855737867353333, + 0.70016297748732993103067821243443810794745752806774, + 0.67195668461417954837935451496149410997032598138385, + 0.64273952430557996253724399164382458078430653187326, + 0.61255388966798023795261245023069487738012378168315, + 0.58144706582913000652989501202985449166039158983665, + 0.54946712509512820207593130552951797023397510159565, + 0.51666060738638370597736752352686730669850255178962, + 0.48307580168617871290856657424482300459902239553311, + 0.44876451363816376391523113354600061846497447744503, + 0.41377920437160500152487974580371368297409962405292, + 0.37817143547358892456847074204267036554749984414768, + 0.34199409082575847300749248117919431006695362002733, + 0.30530244173524671953945499726037915149778023393627, + 0.26815218500725368114118434480859618342480437323625, + 0.23059852188071949700436103627672748185567094334317, + 0.19269758070137109971551685206514989481409202110521, + 0.15450687937939447709272996582455060888266557312076, + 0.11608407067525520848345128440802411376872853085422, + 0.07748658833128284169115486612617194705087430189926, + 0.03877241750605082193319344402462329467936463438314, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod81 = + (Array() << + 0.00078786332389437149872027155489686731646438112543, + 0.00220748573572677796216880922531121147539447921678, + 0.00376522867934192207419437277688852248174719104712, + 0.00527194271488547391100911400749610401795296325231, + 0.00673181348520739996342079325197287731781780936164, + 0.00819757638675148244956105325366150439252190701325, + 0.00967540148401718791503549171298574769106769313863, + 0.01113132166402750374938616623733449344824035103306, + 0.01255438476851726603177494941601931549647320859678, + 0.01396255986698061404257329276730788872218488725171, + 0.01536132635910245297306719365346221049653776351247, + 0.01673453247500258319616665389692844862618082103791, + 0.01807386840881819058019116416727673799123477315931, + 0.01938764589431774100483071215654747790839959650200, + 0.02067904327352817531538698511576519576908963214656, + 0.02193818733583309346140081267327422472167322967377, + 0.02315893101337702414441542631657291426916138363685, + 0.02434569018227335927008045325065084644250141492898, + 0.02550021760313012760411536558132073210208609486171, + 0.02661573749902468675858406158693142228991791577345, + 0.02768762611106109151543416219022982991880587861249, + 0.02871838684109212328774430242495038364855427947406, + 0.02970892727777659464157767791486382641598503632608, + 0.03065436089141152537823603492997426208913617114202, + 0.03155122361911536248171493595674374505113742718348, + 0.03240098250760594428516927456513374062015435029682, + 0.03320404434125756040053598377665629314697494328677, + 0.03395686283420980625135217300229063782965844230293, + 0.03465693584349753394613505048683965340397811449698, + 0.03530514470862184103888891975079794872037448391518, + 0.03590160278362810442749426164234802328728310552020, + 0.03644382653034092475806451684526800902468673825359, + 0.03693016953404855460457701491128601575795938601292, + 0.03736118002546921808817250104671655660126038378354, + 0.03773680126309354415257237969171128676159927920230, + 0.03805546377885242099072980354540709765269970008461, + 0.03831632400517485967847680297981699802520907959292, + 0.03851974174995072693620904485019618182205811434106, + 0.03866555543914104039741924425967757530593426930461, + 0.03875302937875238614021147436365714157701740560300, + 0.03878210476428280538640259652566010072577182998938 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss81 = + (Array() << + 0.99823770971055920034962270242058649233577038159505, + 0.99072623869945700645305435222137215496222208135109, + 0.97725994998377426266337028371290380697866793203799, + 0.95791681921379165580454099945275928509488349060275, + 0.93281280827867653336085216684520571643475357528269, + 0.90209880696887429672825333086849310358448808105767, + 0.86595950321225950382078180835461996357054655301110, + 0.82461223083331166319632023066609877390724038424299, + 0.77830565142651938769497154550649484802069131612689, + 0.72731825518992710328099645175493054855737867353333, + 0.67195668461417954837935451496149410997032598138385, + 0.61255388966798023795261245023069487738012378168315, + 0.54946712509512820207593130552951797023397510159565, + 0.48307580168617871290856657424482300459902239553311, + 0.41377920437160500152487974580371368297409962405292, + 0.34199409082575847300749248117919431006695362002733, + 0.26815218500725368114118434480859618342480437323625, + 0.19269758070137109971551685206514989481409202110521, + 0.11608407067525520848345128440802411376872853085422, + 0.03877241750605082193319344402462329467936463438314 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss81 = + (Array() << + 0.00452127709853319125847173287818533272783111019971, + 0.01049828453115281361474217106727965237679262131580, + 0.01642105838190788871286348488236392729234229334696, + 0.02224584919416695726150432418420857320703319667935, + 0.02793700698002340109848915750772107730255086205077, + 0.03346019528254784739267818308641084897724178665376, + 0.03878216797447201763997203129044616225345921123203, + 0.04387090818567327199167468604171549581100683717024, + 0.04869580763507223206143416044814638806784302737712, + 0.05322784698393682435499647977226050455532117182201, + 0.05743976909939155136661773091042598560010483585445, + 0.06130624249292893916653799640839859590259376351117, + 0.06480401345660103807455452956675273003269296420849, + 0.06791204581523390382569010823192398598419723837929, + 0.07061164739128677969548363085528683235955910399559, + 0.07288658239580405906051068344251783585755908098580, + 0.07472316905796826420018933626132467319120293442036, + 0.07611036190062624237155807592249482301255955384507, + 0.07703981816424796558830753428381024852443975416394, + 0.07750594797842481126372396295832632696366865278811 + ).finished(); + +// Nodes and Weights - Rule 91 +template +Array QuadratureKronrod::abscissaeGaussKronrod91 = + (Array() << + 0.99976825819812844184050501808521584761778222189509, + 0.99860364518193663815654767690082054809228852721448, + 0.99623648400262042934511873220069067464977892521152, + 0.99264999844720374174861712059773534497500164330364, + 0.98788929755834198035845135712826787436469148879696, + 0.98196871503454056823931847363434148153836524922046, + 0.97487472859984626283011247237626660754074730116341, + 0.96660831039689460473642516089247809711667903060942, + 0.95719162401328142632709071186085980254412868260598, + 0.94664169099562906178472059695383707845046834416714, + 0.93496270170823072849720985525389099772434749839611, + 0.92216393671900038809746736096052731926536287725066, + 0.90826680768346396973605068382241808590437868337987, + 0.89329167175324173846464905149305734636540124041379, + 0.87725158719302201212972841148989357243612789849954, + 0.86016247596066422533907887056713438778341442357702, + 0.84204857691040695822975484057647336894431244265745, + 0.82293422050208633703577526002650202999341525805627, + 0.80283895116683677464879103408690607617541369942643, + 0.78178431259390629131236318809860281766027351683906, + 0.75979816452041089331911717455595240315934812011636, + 0.73690884894549035262373884859489201987960568083762, + 0.71314122427549174186073675931964340779561425372465, + 0.68852168077120052523201982588043823595523300561853, + 0.66308170860223379808352326376326032948162534249687, + 0.63685339445322335927122384590338315069203124346568, + 0.60986605449389570894381986808230584011038649491498, + 0.58215021256935318668096733444417713100008882598323, + 0.55374067215934876332312890149007883417942668039810, + 0.52467282046291606709113410046016240260710827080028, + 0.49497965749810183702993192041027660372919586028089, + 0.46469512391963509857960150230974851083981266723719, + 0.43385684374178294432166447428181289994755173351281, + 0.40250294385854191407797450854834645597029935756940, + 0.37066934235973061895322964085585254593836191642981, + 0.33839265425060216164340410003187298953234144313742, + 0.30571272186623304325858534427814131940857839940069, + 0.27266976975237756060876539161564518189821175427664, + 0.23930185320471224179064130404859896962578658675486, + 0.20564748978326374571978722547154393070323638569943, + 0.17174806594978090965313436541849815218083385269234, + 0.13764520598325302875659004142306555310676555654153, + 0.10337830283214540467276673615057417548074590315259, + 0.06898698016314417249041461410381172969392066693097, + 0.03451344875177766949491534277959965665394030693585, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod91 = + (Array() << + 0.00062429185718321579857663066059662918012169021712, + 0.00174909651748535281285516187173085451060602294185, + 0.00298417191546646641998503281623961554312494201463, + 0.00417984253266509402229000496101711163585669889060, + 0.00533932741808091123159062611497603209976404263396, + 0.00650505471381696245904611633672298845714802006140, + 0.00768283318432149209385511436140106223666562795019, + 0.00884569677730051624106815234208503432536098150154, + 0.00998469672095943272307954145417903769666978042970, + 0.01111490008447190508121416984333781971203603824680, + 0.01224167381140228811054606558412293322915204182020, + 0.01335197006885853998622059924442501930963026478306, + 0.01443876216592827376541528869012420064884302430462, + 0.01550957989105782085098413853805035872321392001603, + 0.01656805497063899757282455227144833435996323736845, + 0.01760592823349179892606175575954299643495333036010, + 0.01861770602562484479718322571024115120208025395210, + 0.01960791451368152118747566973260204511490037731486, + 0.02057913131743379223158230090375461301131538208475, + 0.02152536104976848200630137341864843337072862418308, + 0.02244212099311634372269863847029007990818455855580, + 0.02333248915075203001289733158912301176791164139547, + 0.02419844417232129113332510712740896880597195801619, + 0.02503528547727422745626205933879306013580656266019, + 0.02583923328997513115912223308380120472230411615960, + 0.02661260570454035168063547182183051788811540169958, + 0.02735705869931106756998238581749345123775004557503, + 0.02806873537773934426971004212203570044889225018607, + 0.02874437835191346477732940766106574336017984206002, + 0.02938590583816742107250843684941804435493284536361, + 0.02999482799045285292199349638404121450838986974628, + 0.03056789199447173972249252053845442669853092956157, + 0.03110225031086367084980881406978430961362268833061, + 0.03159963354172971455344564253125185922453333794127, + 0.03206152703486851892940119666539450559009394304134, + 0.03248514058861336432416886222383084559924471982830, + 0.03286795999547159469562119073429122583191588601916, + 0.03321166534480586792126042755694344444575877149100, + 0.03351780561085007806314526562438369890686100054127, + 0.03378396062701886695872739048248957350422138863593, + 0.03400789127506731448126888746030257764555957519386, + 0.03419132258692411201484025624573949944920924709690, + 0.03433593627332217900296412964336557875142219514127, + 0.03443961425268577159747120385314938034943289238764, + 0.03450034175251302542891776315122965078611807559511, + 0.03451995999118589472369263667768431132580391340384 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss91 = + (Array() << + 0.99860364518193663815654767690082054809228852721448, + 0.99264999844720374174861712059773534497500164330364, + 0.98196871503454056823931847363434148153836524922046, + 0.96660831039689460473642516089247809711667903060942, + 0.94664169099562906178472059695383707845046834416714, + 0.92216393671900038809746736096052731926536287725066, + 0.89329167175324173846464905149305734636540124041379, + 0.86016247596066422533907887056713438778341442357702, + 0.82293422050208633703577526002650202999341525805627, + 0.78178431259390629131236318809860281766027351683906, + 0.73690884894549035262373884859489201987960568083762, + 0.68852168077120052523201982588043823595523300561853, + 0.63685339445322335927122384590338315069203124346568, + 0.58215021256935318668096733444417713100008882598323, + 0.52467282046291606709113410046016240260710827080028, + 0.46469512391963509857960150230974851083981266723719, + 0.40250294385854191407797450854834645597029935756940, + 0.33839265425060216164340410003187298953234144313742, + 0.27266976975237756060876539161564518189821175427664, + 0.20564748978326374571978722547154393070323638569943, + 0.13764520598325302875659004142306555310676555654153, + 0.06898698016314417249041461410381172969392066693097, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss91 = + (Array() << + 0.00358266315528355893114302865935139257058563635818, + 0.00832318929621824164573585312223385375615793129295, + 0.01303110499158278432063108246968692602714117158264, + 0.01767753525793759061709254666957708852065127812647, + 0.02223984755057873239395075855216899334487570870130, + 0.02669621396757766480567477879310752778782491056589, + 0.03102537493451546716250793889376806388592115212397, + 0.03520669220160901624769979826157509081008807471271, + 0.03922023672930244756418718534392933821813012459631, + 0.04304688070916497115169111308116694178831923616319, + 0.04666838771837336526776847574165410120498614148785, + 0.05006749923795202979913210247487431739149869293515, + 0.05322801673126895194590404401931040379848221694018, + 0.05613487875978647664392394037486976324148018494616, + 0.05877423271884173857436151763183142253126572546098, + 0.06113350083106652250188637053632557326012423426524, + 0.06320144007381993774996373029066688186630023727582, + 0.06496819575072343085382657035907569508384972548758, + 0.06642534844984252808291471563910373761317454675067, + 0.06756595416360753627091022387364863097947258092395, + 0.06838457737866967453169209933431609707348491305464, + 0.06887731697766132288200284829805578888452413011065, + 0.06904182482923202011079855515940474072332289227500 + ).finished(); + +// Nodes and Weights - Rule 101 +template +Array QuadratureKronrod::abscissaeGaussKronrod101 = + (Array() << + 0.99981190136436471898752764191979778490900714786428, + 0.99886640442007105018545944497421850599624351290408, + 0.99694438701887617830552521284551574074542656874602, + 0.99403196943209071258510820042069472815747797106832, + 0.99016501066968000942392663658159575132877674133477, + 0.98535408404800588230900962563248940401559263094541, + 0.97958729676076942938993428547746899117341045817570, + 0.97286438510669207371334410460625205366917340704995, + 0.96520165106614520049209141719129593627549507978394, + 0.95661095524280794299774564415662209405143412462605, + 0.94709403424493944376335887833533441316436721753676, + 0.93665661894487793378087494727249660215373159809520, + 0.92531354647480188621729501218597896170837009890339, + 0.91307855665579189308973564277165709478418819169778, + 0.89995988756429461789335121594253760342125497503872, + 0.88596797952361304863754098246675363419429031075581, + 0.87111919820410031041116631276794734769410099155684, + 0.85542976942994608461136264393475746765483303948691, + 0.83891258696722439582669453115982171629679423354273, + 0.82158207085933594835625411087393953776074138334349, + 0.80345686805045941104583132344339868468195642520966, + 0.78455583290039926390530519634099120084731627255644, + 0.76489567937235152255213391009405332066417203426695, + 0.74449430222606853826053625268219424287018793132958, + 0.72337276602259251579221941156517307053014928418822, + 0.70155246870682225108954625788365572814971922848470, + 0.67905338940674727147750574623528945542568600842673, + 0.65589646568543936078162486400367981904141052838233, + 0.63210506844506487645226444165317494562655419021859, + 0.60770292718495023918038179639183289360420502067599, + 0.58271281781783653230059341928604868801078591807290, + 0.55715830451465005431552290962580160781589838223775, + 0.53106482452270813956871444365352729284157032257425, + 0.50445814490746420165145913184914119263537867827068, + 0.47736339223304303620238841113348793007606445039659, + 0.44980633497403878914713146777837581731506451346517, + 0.42181415731061309702869639767299001110343390879166, + 0.39341431189756512739422925382381727024613946867266, + 0.36463382886161427721485091875517338779531125640738, + 0.33550024541943735683698825729107169784121859303680, + 0.30604211922918460908352322458238267491364761240808, + 0.27628819377953199032764527852113018571480158713200, + 0.24626694739814402597841233060681906225445876784272, + 0.21600723687604175684728453261710133370575597251013, + 0.18553858172277277270028588015893298081440861045060, + 0.15489058999814590207162862094110950120185022005494, + 0.12409272435916037178300113567216009078143172259769, + 0.09317470156008614085445037763960034788567138392214, + 0.06216656481941616908016236906249366229724267742777, + 0.03109833832718887611232898966594919424729622296000, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod101 = + (Array() << + 0.00050676166803489136800323005679418875182879376563, + 0.00142011023816635705848924160910481934652645245411, + 0.00242310374582057327186801995970143275616529077528, + 0.00339459089289723743085662727982835774647902739813, + 0.00433770360526373578354342360625291962127768748888, + 0.00528690708431070081496631090522651841927406286998, + 0.00624676165289346880527846119294362879366575110822, + 0.00719585591375833177104770901157397246654579768609, + 0.00812753366255938737868037294245956987246699844772, + 0.00905390685313200073482984638318670974026478257534, + 0.00997908051277197271520006896555491352835623156784, + 0.01089303114914542139185373632048379056944269951885, + 0.01179073210054285467626745948807510643296154432320, + 0.01267806183701224440357173722638836784313869359541, + 0.01355761601155914093045815780333977992202586708240, + 0.01442332487228954690873605873699027758803460839370, + 0.01527146381217306258703163983580949905342101469485, + 0.01610536347829330632879586448333389813140260935888, + 0.01692665787573854108211355546392361200532841207090, + 0.01773116846670068736945143328713021401556189039555, + 0.01851604246177747374224081714934718333316987698583, + 0.01928332352180814825301889801162132421117243231650, + 0.02003404583071943486093518979257578870512514064793, + 0.02076512847171916224208023972166701446741771437993, + 0.02147433820180993345171927343976448497995531229818, + 0.02216298492660156610684349832667521094279850151824, + 0.02283171756887400699163339995992033749214277945418, + 0.02347818407891019300458609050334622304513636748516, + 0.02410062587721830597666576014469811068703941778003, + 0.02469990353732573949476020971483749324030400159067, + 0.02527641257483270914647897860737446154499106952106, + 0.02582833680912179520884486377406142381526019962380, + 0.02635430681869013845472576260957721404783898088646, + 0.02685489636290659736893714711520551156209704917336, + 0.02733033189010964383305404350544173245562448595110, + 0.02777922389243983336140685575895975209594524979909, + 0.02820054057947867914629860187746193364776090154299, + 0.02859466838340605773904731791010499100505911124158, + 0.02896171979050103356757141731956685768195060652065, + 0.02930066703511978233684077865154509911023252219986, + 0.02961078494045995503025213539965659277620209307433, + 0.02989233581703979245880464830873980542938987086275, + 0.03014535377457083519407513020913017271701705074418, + 0.03036913320016487318137849566585392442983158643626, + 0.03056323770631003720546053583408408633745212487948, + 0.03072784566941669020593437016581195240625267413048, + 0.03086293412170997376026432332960616247053522064300, + 0.03096809562366843795854210956298413399486076114109, + 0.03104317399694570106876064569864783540375985841535, + 0.03108828778240521910421818490634306613818955150864, + 0.03110336664174957546715464493457708644263854335374 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss101 = + (Array() << + 0.99886640442007105018545944497421850599624351290408, + 0.99403196943209071258510820042069472815747797106832, + 0.98535408404800588230900962563248940401559263094541, + 0.97286438510669207371334410460625205366917340704995, + 0.95661095524280794299774564415662209405143412462605, + 0.93665661894487793378087494727249660215373159809520, + 0.91307855665579189308973564277165709478418819169778, + 0.88596797952361304863754098246675363419429031075581, + 0.85542976942994608461136264393475746765483303948691, + 0.82158207085933594835625411087393953776074138334349, + 0.78455583290039926390530519634099120084731627255644, + 0.74449430222606853826053625268219424287018793132958, + 0.70155246870682225108954625788365572814971922848470, + 0.65589646568543936078162486400367981904141052838233, + 0.60770292718495023918038179639183289360420502067599, + 0.55715830451465005431552290962580160781589838223775, + 0.50445814490746420165145913184914119263537867827068, + 0.44980633497403878914713146777837581731506451346517, + 0.39341431189756512739422925382381727024613946867266, + 0.33550024541943735683698825729107169784121859303680, + 0.27628819377953199032764527852113018571480158713200, + 0.21600723687604175684728453261710133370575597251013, + 0.15489058999814590207162862094110950120185022005494, + 0.09317470156008614085445037763960034788567138392214, + 0.03109833832718887611232898966594919424729622296000 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss101 = + (Array() << + 0.00290862255315514095840072434285548080667299645995, + 0.00675979919574540150277887817798503180187383240647, + 0.01059054838365096926356968149924102233940181908646, + 0.01438082276148557441937890892732434993703178617059, + 0.01811556071348939035125994342235461984466731704973, + 0.02178024317012479298159206906269034122731346235793, + 0.02536067357001239044019487838544272346016125997571, + 0.02884299358053519802990637311323243251784686559354, + 0.03221372822357801664816582732300395344858905883343, + 0.03545983561514615416073461100097579709696000049698, + 0.03856875661258767524477015023638593486477170500052, + 0.04152846309014769742241197896406701780897797548584, + 0.04432750433880327549202228683039419746076129835545, + 0.04695505130394843296563301363498768251406430618605, + 0.04940093844946631492124358075143272869228705096661, + 0.05165570306958113848990529584009527964982544939544, + 0.05371062188899624652345879725566455276802321352992, + 0.05555774480621251762356742561226949759513529998390, + 0.05718992564772838372302931506599316301157537225709, + 0.05860084981322244583512243663084846620976751344403, + 0.05978505870426545750957640531258523079666604207267, + 0.06073797084177021603175001538481100160979927323540, + 0.06145589959031666375640678608391537509726757576401, + 0.06193606742068324338408750978083068857287705669124, + 0.06217661665534726232103310736061343086768246920103 + ).finished(); + +// Nodes and Weights - Rule 121 +template +Array QuadratureKronrod::abscissaeGaussKronrod121 = + (Array() << + 0.99986894039735508603454696805483048308252339364253, + 0.99921012322743602203422958579764926639879897871133, + 0.99787067797281103000055449680255041109802834093833, + 0.99584052511883817387674671337744065277489029497714, + 0.99314406509020969852781424310993230666589953769707, + 0.98978789522222171736727898701609604257739315378694, + 0.98576273403417081406883799910761849326959871083386, + 0.98106720175259818561857679982677003388053857922189, + 0.97571130869063093805424022397979537946475594414538, + 0.96970178876505273372154409891379427254093434486825, + 0.96303809044100194217445769061433368420308650851500, + 0.95572225583999610739723184582969950248009374476661, + 0.94776254337082142416238277974195989092969410130796, + 0.93916627611642324949541901160970509641940342040120, + 0.92993671378354345915802464683245381420494070890654, + 0.92007847617762755285665686251989682878376887601371, + 0.90960021010202694746557389315600630364455686348403, + 0.89851031081004594193778932957262833016445844274707, + 0.88681462876695505446866905561876851436297727385471, + 0.87451992264689831512930809991243580515918229443578, + 0.86163582612712916884010876561015688672720821191562, + 0.84817198478592963249051549499437560129752871650678, + 0.83413634063502318823513056433153895616632797500262, + 0.81953752616214575936851810851972357189949051896959, + 0.80438634842167500661914227180660775048846373673656, + 0.78869373993226405456994479977721573524633292787923, + 0.77246944634325712130468945617575041059983850332073, + 0.75572377530658568686884206660237310592890019956471, + 0.73846873678927334614831784537093050545190770229944, + 0.72071651335573039943602106101352107776477399473324, + 0.70247844071277491341354556676393514423069049151760, + 0.68376632738135543722293023922429715908759460433887, + 0.66459333948891834538628881414033517814480799774455, + 0.64497282848947706781344789642044471716519001934666, + 0.62491753372827970746014412059280947928268675981204, + 0.60444059704851036344420877631120172065042105359106, + 0.58355624965567953080935429206168664732529695996909, + 0.56227890075394453917827258748599872906318823765940, + 0.54062251531844021045684191707599245063719726744878, + 0.51860140005856974741788934848472130132984667628869, + 0.49623073093171768406858727593889830377761134305524, + 0.47352584176170711110816305375279463800786800509984, + 0.45050174562863861734206716711875032408754056360999, + 0.42717374158307838930745285353031192335344213618423, + 0.40355780908279404695333333368401568081423124589053, + 0.37967005657679797715495267052188770405637443316239, + 0.35552636501293175588409210698088763021622142976651, + 0.33114284826844819425235296535055270480012149637631, + 0.30653613242384558698842567017265989136532053439467, + 0.28172293742326169169069486033944159830826306079538, + 0.25671982824253729939215528961771483496009022967204, + 0.23154355137602933801034463134675544074353905193959, + 0.20621121086087678024608930136722267382406551852335, + 0.18073996487342541724087694126185263123177216863905, + 0.15514687531461929153617689855444662362933578843402, + 0.12944913539694500314644416464957577475987352436647, + 0.10366414973805567097692423407501883537630426619863, + 0.07780933394953656941928550708222528757759478289772, + 0.05190205753879922246314126084963723796593693305509, + 0.02595977230124779858917038540034482566031476071728, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod121 = + (Array() << + 0.00035309352008476328873302150506396067820554489853, + 0.00098957111594993012002795471522501001707498606076, + 0.00168881658669538272849089285105177845280274177469, + 0.00236661114697454617075475538520317633315532160502, + 0.00302530527821396938049320055216641290240634905666, + 0.00368912213435713356782149157136472153921088296969, + 0.00436146638924688454590356599869214359782276998026, + 0.00502758855697275303082752008261190854561998819561, + 0.00568300310799120346648467571699154030415431078060, + 0.00633636534216254504799580887837090652934518564671, + 0.00699076245664291472692253153841534937056702016655, + 0.00763935787810120919989669059239195349160759882912, + 0.00827878115796801323966142828447078823241218583717, + 0.00891334185215251454610013327637323498094064454821, + 0.00954507013554903770056092516381352299744866414013, + 0.01016985984860481299714318208572482794015404599887, + 0.01078522284204381631750127741952879057916973310398, + 0.01139368492941710872440593756091481437036855202376, + 0.01199659638647986661732163593813166954215294680801, + 0.01259114041876269638053880248826790662577710230923, + 0.01317540660371342750620081979004402932037003022067, + 0.01375101092307560054089755109416989805354339164989, + 0.01431887215442789937035515803915600785470840146381, + 0.01487690562344235570814655024966329586556913856472, + 0.01542359803641373989397778363148163169047625524956, + 0.01596003744934033355994981103588918423169316792324, + 0.01648685780200353115743507689306805446791494560564, + 0.01700244529575939575670279579554073187580588123302, + 0.01750557806830523681292548656619025120665578599915, + 0.01799701182836498525925640182815936357315182357777, + 0.01847718567886398682948269679412403693476107997848, + 0.01894481804101016688214698965950324848350668749227, + 0.01939891427004672519965813885528064708101401901242, + 0.01984000938527257124526765005244593165383386796636, + 0.02026840565003063810461444077298412851468322177833, + 0.02068307361792965307984840163452720256611289248871, + 0.02108320540912487808345587005640913822849056383719, + 0.02146918421657572995510375908055968070647062620577, + 0.02184121436913142190366087107926959818357547855685, + 0.02219846947126649229671310378395216714738953151697, + 0.02254030223805028212731124333425081265066909380287, + 0.02286698888783978160371241817159477370209744385441, + 0.02317866262411930143640919197308030421619718533173, + 0.02347466875931088521516214970827055070559381432420, + 0.02375450343651776669480062293094101720463883872448, + 0.02401836617317175928431662138798116606010655200640, + 0.02426633771149105650424641365034412183131195279715, + 0.02449791452590019789919446326522050049632493665594, + 0.02471272496709623460174523761205724250739759370788, + 0.02491091263261197247357887006487721264874660274884, + 0.02509251862664474061093109702616181054149913651710, + 0.02525717697351328657763893489291540401679249790887, + 0.02540464114326442690732349325450960796495815128648, + 0.02553501298607099536450151815998273053431323202531, + 0.02564830231130669156727699936154427063030445683730, + 0.02574427177483327835345017378864237099131910116066, + 0.02582279593614917829003082527854324007813305881338, + 0.02588394414600359533197858687343537654901210069360, + 0.02592769956371017513985574241157347171723158109308, + 0.02595394793410606119738320085065067300844960043221, + 0.02596268333963692375153481541246139454955940844745 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss121 = + (Array() << + 0.99921012322743602203422958579764926639879897871133, + 0.99584052511883817387674671337744065277489029497714, + 0.98978789522222171736727898701609604257739315378694, + 0.98106720175259818561857679982677003388053857922189, + 0.96970178876505273372154409891379427254093434486825, + 0.95572225583999610739723184582969950248009374476661, + 0.93916627611642324949541901160970509641940342040120, + 0.92007847617762755285665686251989682878376887601371, + 0.89851031081004594193778932957262833016445844274707, + 0.87451992264689831512930809991243580515918229443578, + 0.84817198478592963249051549499437560129752871650678, + 0.81953752616214575936851810851972357189949051896959, + 0.78869373993226405456994479977721573524633292787923, + 0.75572377530658568686884206660237310592890019956471, + 0.72071651335573039943602106101352107776477399473324, + 0.68376632738135543722293023922429715908759460433887, + 0.64497282848947706781344789642044471716519001934666, + 0.60444059704851036344420877631120172065042105359106, + 0.56227890075394453917827258748599872906318823765940, + 0.51860140005856974741788934848472130132984667628869, + 0.47352584176170711110816305375279463800786800509984, + 0.42717374158307838930745285353031192335344213618423, + 0.37967005657679797715495267052188770405637443316239, + 0.33114284826844819425235296535055270480012149637631, + 0.28172293742326169169069486033944159830826306079538, + 0.23154355137602933801034463134675544074353905193959, + 0.18073996487342541724087694126185263123177216863905, + 0.12944913539694500314644416464957577475987352436647, + 0.07780933394953656941928550708222528757759478289772, + 0.02595977230124779858917038540034482566031476071728 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss121 = + (Array() << + 0.00202681196887375849643171020989232469531166822556, + 0.00471272992695356864089482171407723594808497019372, + 0.00738993116334545553151695602208606088222818020623, + 0.01004755718228798435788576437705725649965174748768, + 0.01267816647681596013149537926951423498833836430263, + 0.01527461859678479930672603809882533517790938411507, + 0.01782990101420772026039626124834856763385921958339, + 0.02033712072945728677503214741710630146201711250534, + 0.02278951694399781986378345819290020781016738421521, + 0.02518047762152124837957096597236127913178083294419, + 0.02750355674992479163522319763862224764840347075416, + 0.02975249150078894524083648467348771033247515812411, + 0.03192121901929632894945889953676047676669729228299, + 0.03400389272494642283491440155525874118727622355286, + 0.03599489805108450306657864628806233859424242002100, + 0.03788886756924344403094079420927603181008428106381, + 0.03968069545238079947012283481171000876280351129417, + 0.04136555123558475561316383680665889688018700852151, + 0.04293889283593564195423122065638280694921177196492, + 0.04439647879578711332778416409137736966388498573401, + 0.04573437971611448664719645529090933300793497425868, + 0.04694898884891220484701315639470156627557490258344, + 0.04803703181997118096366665272873365360250002924692, + 0.04899557545575683538947568685789429835758389599909, + 0.04982203569055018101115923089370331247921559724378, + 0.05051418453250937459823873574165364323995616273742, + 0.05107015606985562740454912073449103536619708311367, + 0.05148845150098093399504439717705431173392382665958, + 0.05176794317491018754380364302882373293234694549920, + 0.05190787763122063973286493836226968018015159556564 + ).finished(); + +// Nodes and Weights - Rule 201 +template +Array QuadratureKronrod::abscissaeGaussKronrod201 = + (Array() << + 0.99995250325234874194558759586872675230201032186296, + 0.99971372677344123367822846934230067671834952730840, + 0.99922816588380125603468689462902570985161537549449, + 0.99849195063959581840016335918634916230485485042057, + 0.99751361312272973925231594907328863008404380861059, + 0.99629513473312514918613173224113103543643128814043, + 0.99483262193692678214229984365392976193623306275107, + 0.99312493703744345965200989284878347073177145886652, + 0.99117498765102584468434111430887815181054026590303, + 0.98898439524299174800441874580773663183933363710695, + 0.98655201560314858543214329463373901430393735434646, + 0.98387754070605701549610015551100816734436701685080, + 0.98096284258061262094794383303858458336270892150820, + 0.97780935848691828855378108842920192863523449426625, + 0.97441692579132830191156952676419699735650628441883, + 0.97078577576370633193089785789750538855055719947821, + 0.96691753657125109640782216846696070235549924264129, + 0.96281365425581552729365932603016638643733150673041, + 0.95847452546444300989234207911968140986589441225372, + 0.95390078292549174284933693089435764464522145101467, + 0.94909405007764308807002661063238399920305087971147, + 0.94405587013625597796277470641521874673972037338210, + 0.93878704353080861159581415143454540129149440268213, + 0.93328853504307954592433366813086250408354607429702, + 0.92756205890488017513109041315467906526813611874056, + 0.92160929814533395266695132848198745912458279773220, + 0.91543138327874017569217880264508697295963148874566, + 0.90902957098252969046712633778914606444327728958459, + 0.90240571020837873902681113376935748186170440444374, + 0.89556164497072698669852102243022776984818176899772, + 0.88849879400991357098633030987871529168085698938796, + 0.88121867938501841557331682542780558244549441102126, + 0.87372330591888454232901721359869700713999812944025, + 0.86601468849716462341073996967624296638031483930559, + 0.85809450656528158256769359279995027785553637603413, + 0.84996452787959128429336259142010465407379077950067, + 0.84162692263159422803428490828685310746728377281023, + 0.83308387988840082354291583384475567990749483030996, + 0.82433731934785511555382066394577351308957189674952, + 0.81538923833917625439398875864925800538255037651370, + 0.80624197522998149898858119947772518581699983645897, + 0.79689789239031447638957288218324598288952685596429, + 0.78735913298215327121915096807989317961453460725592, + 0.77762790964949547562755138683449010653853979993608, + 0.76770672752390789522106835290910988619423163187483, + 0.75759811851970717603566796443840077231310897190004, + 0.74730443440010216722543283418354970398268442664090, + 0.73682808980202070551242771482010100284327844624712, + 0.72617175233016540458785543702154726618162507726741, + 0.71533811757305644645996712270436596408439783859563, + 0.70432973201100171130158275065992755531836865624424, + 0.69314919935580196594864794167543726558700001793073, + 0.68179934329913790750199866457132454892446683289894, + 0.67028301560314101580258701432322661366980568402882, + 0.65860294402143308554822559819231646499591256882672, + 0.64676190851412927983263030445863043501973378424853, + 0.63476288086691287961700686682300223062843045115192, + 0.62260886020370777160419084517231224465381773228982, + 0.61030274223394973052304014680666051705840336769498, + 0.59784747024717872126480654514934063639489919232049, + 0.58524615484722223751330407348246631581348705331813, + 0.57250193262138119131687044352572544896003394967556, + 0.55961785385987287004483883274572350988535188007656, + 0.54659701206509416746799425718174990395624177593753, + 0.53344264632810924391925934260941164781644110851200, + 0.52015801988176305664681574945520853076893769042009, + 0.50674632405992146013910672747973470564589036548717, + 0.49321078920819093356930879344933399099072332535856, + 0.47955477168853305490480433055889014820419073581020, + 0.46578164977335804224921662339575458161165111021221, + 0.45189474208663614330483287246151294430151576448075, + 0.43789740217203151310897804362219596212570176348410, + 0.42379309169388210181506692715279269905002159373794, + 0.40958529167830154252886840005715770149536438916475, + 0.39527743404816671610367673557861642687518028284140, + 0.38087298162462995676336254886958740374970726512371, + 0.36637548876362800231863519558827019975726941799314, + 0.35178852637242172097234382954897056524931809638907, + 0.33711562544519365010380511464074285765053795714571, + 0.32236034390052915172247658239832542740219162302309, + 0.30752631562753657120752788303988100046876429120433, + 0.29261718803847196473755588823549438456153898917258, + 0.27763657675183495616346814290881136285663496860995, + 0.26258812037150347916892933625498214113202269453552, + 0.24747551869798892096614882491742687591540876456389, + 0.23230248184497396964950996320796411069750977150714, + 0.21707269541726899854706722200296707797219389992710, + 0.20178986409573599723604885953039646294369200355905, + 0.18645773958437819065168393518669992652130113816299, + 0.17108008053860327488753237470708980746585972511807, + 0.15566062775725995230028233532247535151338751890639, + 0.14020313723611397320751460468240551661687300626336, + 0.12471139826072812045287156735452075155899823333651, + 0.10918920358006111500342600657938488688489962996916, + 0.09364033428354398098988456682945235155973991630929, + 0.07806858281343663669481737120155257397635002744853, + 0.06247776146923010010357548046587278110387007259524, + 0.04687168242159163161492391293384830953706539908602, + 0.03125415208386678081857435045254965009647374420346, + 0.01562898442154308287221669999742934014775618285556, + 0.00000000000000000000000000000000000000000000000000 + ).finished(); + +template +Array QuadratureKronrod::weightsGaussKronrod201 = + (Array() << + 0.00012796430957024721771296604719777847195024888314, + 0.00035867672428027546451819698699002695000334366552, + 0.00061229953852751686968468056548070834614262353518, + 0.00085841483623935311664167374912297679663629462243, + 0.00109796416132059890255397465155433695320668680377, + 0.00133983780660382341754863019304498885110952126241, + 0.00158539114055099209687612715239439499710635284278, + 0.00182937001345660807993702363601705354484262856545, + 0.00207023170435975873786445531287183441533108495909, + 0.00231123008487149478622574422551013428738043877307, + 0.00255360720035108161107860930247813105889546630505, + 0.00279496754149865554132924610296587377583807221606, + 0.00303416468796449379917327454133990046860593313429, + 0.00327287902049999951685865323103972155152040933937, + 0.00351196946843139471347132584863839017075918804593, + 0.00375002082356796252973315848571211234815922368507, + 0.00398619744931235184550297813723407560532468657842, + 0.00422153240339578242265345682700445192052466529462, + 0.00445663655824494510224150713533978978713271961637, + 0.00469055184548748439481506164708745018600910075929, + 0.00492264171642177778021256522028596859149482883314, + 0.00515360688570317554082069990644640034295769630886, + 0.00538389880957824664615043517317289909958011793608, + 0.00561281343110871597561315793834039501636153255547, + 0.00583984643867452187979812998241983823823136309861, + 0.00606550320973303863921033440183559298376742116472, + 0.00629012776973767504951202123451843743076930245062, + 0.00651317411315535135623979587712921532296014614460, + 0.00673423042383435906546110501890771393375794286618, + 0.00695367647707533778618333444634079607272882365501, + 0.00717178053315253341524750684836420331388733863189, + 0.00738810306985931173687401288774170241141466081296, + 0.00760230005117430829817610049263739864681050864527, + 0.00781466520662034505085711719016354399072460775306, + 0.00802541132258648471308962366547616916200657491446, + 0.00823417501731623690079536275494453055496288986340, + 0.00844066393687280440093304148919420862913744658522, + 0.00864511004867877912656139959329107529390995475094, + 0.00884768427998831020007201725385439489287388704377, + 0.00904808029993671228835071165243690477400271230057, + 0.00924604652391111724687570119177700799209325218608, + 0.00944176898875185552314903500699258800726938238310, + 0.00963538626179462391106485411435148851979409901308, + 0.00982663642373088065571326334283803573470396510601, + 0.01001530101161305053268537450183723383477657205865, + 0.01020153095088495986015919422291444604170764482781, + 0.01038543929706851271926198711586114579021961800284, + 0.01056679983247565098373215943168121363502057278619, + 0.01074542170778693080646692567379385062368419611030, + 0.01092142841500597232724852424011180927323636614697, + 0.01109491257094606883971232155924910988765741030054, + 0.01126567746332956237923986454364445975020263538030, + 0.01143355579967040516909625306002457417492109489982, + 0.01159864925643444521982733110763019057800050447511, + 0.01176103385541580107532667409174481684699820604850, + 0.01192053786750440393445381605750602287206317500254, + 0.01207701451438569994891221026442585729555689105328, + 0.01223054787113266210109559371404225931637605012182, + 0.01238120033058132478094426900828573870330496911285, + 0.01252882177740552594442502987847792106610336206449, + 0.01267328363219971069040254068582595376734875571756, + 0.01281465559402401820734387481377863160543256944343, + 0.01295298874891947508181321347643329979562489424432, + 0.01308815204304083947131249628348670851036112224233, + 0.01322003331459425763799456688719462334189620391816, + 0.01334869039477997219532795119897057181197377482865, + 0.01347416490042760653129384849738331885547322480541, + 0.01359634288062617797837336787649308030170304524353, + 0.01371512721021576932690771893822620264523809273465, + 0.01383056582631592119521912896609004389004651368400, + 0.01394269234050930765052410813271065280007760112975, + 0.01405140838748433908894290183169630886722307263276, + 0.01415663080386553141797847993528261640630526208010, + 0.01425839919822706151742984801432832105691149679315, + 0.01435674034589556885653893978408636354274396518129, + 0.01445157028450311174994356610608770124062851800510, + 0.01454281897629966555537908242478331367093974940807, + 0.01463051894733967823881185958389904589719486669315, + 0.01471469106719224107616070689745858786072765279902, + 0.01479526485272091523273713296861772987622838214603, + 0.01487218274551550851114896175149252660492892913239, + 0.01494547117901571779314143269144247770437604078025, + 0.01501514584992601775929409587145443755032002323588, + 0.01508114903502136704084960300677971411706204873991, + 0.01514343516330660012211995270464269800292539698574, + 0.01520202535514477439125881646888569025211919463334, + 0.01525693069928626681152511683204470629971399732269, + 0.01530810567790395821753937559915116616801682403691, + 0.01535551634411844871530936744265760333339456402763, + 0.01539917910754267008730746329865897219178834017417, + 0.01543910086714820098073247575245749712612887092387, + 0.01547524789208464777691069202608380676296793168412, + 0.01550759760659890369589537474066318545807510752952, + 0.01553616215856523293822391298200717744826237954682, + 0.01556094454341676462938645995315419106586887459681, + 0.01558192251438053678351080364024221646030414904878, + 0.01559908471168710670698089639871094604147763682664, + 0.01561243933403481989468460959619688421126306083370, + 0.01562198563724487346851127872839698584387203112330, + 0.01562771265700169003321136142695117550866696845783, + 0.01562962018460484993210603370275515394376090869371 + ).finished(); + +template +Array QuadratureKronrod::abscissaeGauss201 = + (Array() << + 0.99971372677344123367822846934230067671834952730840, + 0.99849195063959581840016335918634916230485485042057, + 0.99629513473312514918613173224113103543643128814043, + 0.99312493703744345965200989284878347073177145886652, + 0.98898439524299174800441874580773663183933363710695, + 0.98387754070605701549610015551100816734436701685080, + 0.97780935848691828855378108842920192863523449426625, + 0.97078577576370633193089785789750538855055719947821, + 0.96281365425581552729365932603016638643733150673041, + 0.95390078292549174284933693089435764464522145101467, + 0.94405587013625597796277470641521874673972037338210, + 0.93328853504307954592433366813086250408354607429702, + 0.92160929814533395266695132848198745912458279773220, + 0.90902957098252969046712633778914606444327728958459, + 0.89556164497072698669852102243022776984818176899772, + 0.88121867938501841557331682542780558244549441102126, + 0.86601468849716462341073996967624296638031483930559, + 0.84996452787959128429336259142010465407379077950067, + 0.83308387988840082354291583384475567990749483030996, + 0.81538923833917625439398875864925800538255037651370, + 0.79689789239031447638957288218324598288952685596429, + 0.77762790964949547562755138683449010653853979993608, + 0.75759811851970717603566796443840077231310897190004, + 0.73682808980202070551242771482010100284327844624712, + 0.71533811757305644645996712270436596408439783859563, + 0.69314919935580196594864794167543726558700001793073, + 0.67028301560314101580258701432322661366980568402882, + 0.64676190851412927983263030445863043501973378424853, + 0.62260886020370777160419084517231224465381773228982, + 0.59784747024717872126480654514934063639489919232049, + 0.57250193262138119131687044352572544896003394967556, + 0.54659701206509416746799425718174990395624177593753, + 0.52015801988176305664681574945520853076893769042009, + 0.49321078920819093356930879344933399099072332535856, + 0.46578164977335804224921662339575458161165111021221, + 0.43789740217203151310897804362219596212570176348410, + 0.40958529167830154252886840005715770149536438916475, + 0.38087298162462995676336254886958740374970726512371, + 0.35178852637242172097234382954897056524931809638907, + 0.32236034390052915172247658239832542740219162302309, + 0.29261718803847196473755588823549438456153898917258, + 0.26258812037150347916892933625498214113202269453552, + 0.23230248184497396964950996320796411069750977150714, + 0.20178986409573599723604885953039646294369200355905, + 0.17108008053860327488753237470708980746585972511807, + 0.14020313723611397320751460468240551661687300626336, + 0.10918920358006111500342600657938488688489962996916, + 0.07806858281343663669481737120155257397635002744853, + 0.04687168242159163161492391293384830953706539908602, + 0.01562898442154308287221669999742934014775618285556 + ).finished(); + +template +Array QuadratureKronrod::weightsGauss201 = + (Array() << + 0.00073463449050567173040632065833033639067047356248, + 0.00170939265351810523952935837149119524373138549146, + 0.00268392537155348241943959042900112008193111495100, + 0.00365596120132637518234245872752519569920656740515, + 0.00462445006342211935109578908297847665035249529489, + 0.00558842800386551515721194634843921073131869400808, + 0.00654694845084532276415210333149526369993836336648, + 0.00749907325546471157882874401639778316358347894815, + 0.00844387146966897140262083490230100193464445988410, + 0.00938041965369445795141823766081211873078704323867, + 0.01030780257486896958578210172783537797605834384143, + 0.01122511402318597711722157336633358477722641956438, + 0.01213145766297949740774479244874817073696312331126, + 0.01302594789297154228555858375890179013496473584175, + 0.01390771070371877268795414910800463779518081214312, + 0.01477588452744130176887998752035425716938874311460, + 0.01562962107754600272393686595379192555246997980994, + 0.01646808617614521264310498008821078082116766160380, + 0.01729046056832358243934419836674167481162350856517, + 0.01809594072212811666439075142049303134757874495839, + 0.01888373961337490455294116588154323429711127634742, + 0.01965308749443530586538147024544406555526959949125, + 0.02040323264620943276683885165758377060570969930262, + 0.02113344211252764154267230044096968163532972887451, + 0.02184300241624738631395374130439802476534899982325, + 0.02253122025633627270179697093167396234015893534871, + 0.02319742318525412162248885418272728845115448573609, + 0.02384096026596820596256041190228343214470744909262, + 0.02446120270795705271997502334977289064629573239780, + 0.02505754448157958970376422562092326422383855852793, + 0.02562940291020811607564200986215087092697767002027, + 0.02617621923954567634230874175730188501127513119069, + 0.02669745918357096266038466418633635063465575039001, + 0.02719261344657688013649156780217069226698789601200, + 0.02766119822079238829420415587042645529240035866422, + 0.02810275565910117331764833018699455045141809940021, + 0.02851685432239509799093676286445787325984272548397, + 0.02890308960112520313487622813451526531560786805526, + 0.02926108411063827662011902349564095444308419504535, + 0.02959048805991264251175451067883658517280628507137, + 0.02989097959333283091683680666859582765809141426080, + 0.03016226510516914491906868161047923265710232578271, + 0.03040407952645482001650785981882517660560724831012, + 0.03061618658398044849645944326205319285308602378906, + 0.03079837903115259042771390303055976009497083447037, + 0.03095047885049098823406346347074792738298717776694, + 0.03107233742756651658781017024291803484591543634796, + 0.03116383569620990678381832121718665334383636868393, + 0.03122488425484935773237649864809813488180274068218, + 0.03125542345386335694764247438619802878783383672609 + ).finished(); +} +#endif // EIGEN_QUADRATURE_KRONROD_H diff --git a/cpp/Integrator.h b/cpp/Integrator.h new file mode 100644 --- /dev/null +++ b/cpp/Integrator.h @@ -0,0 +1,745 @@ +/** + * \file Integrator.h + * This file contains routines for numerical integration using Gauss-Kronrod quadrature. + * \sa R Piessens, E de Doncker-Kapenger, C Ueberhuber, D Kahaner, QUADPACK, A Subroutine Package + * for Automatic Integration, Springer Verlag, 1983. + */ + +#ifndef EIGEN_INTEGRATOR_H +#define EIGEN_INTEGRATOR_H + +namespace Eigen +{ + /** + * \ingroup Quadrature_Module + * + * \brief This class provides numerical integration functionality. + * + * Memory management and additional information (e.g. number of function calls, estimated error), + * are provided by the class in a way that can support the future porting of additional quadrature + * functions from the QUADPACK library. + * + * \todo Ensure only appropriates types are used for Scalar, e.g. prohibit integers. + */ + template + + class Integrator + { + public: + + /** + * \brief The local Gauss-Kronrod quadrature rule to use. + */ + enum QuadratureRule + { + GaussKronrod15 = 1, /**< Use 7, 15 points. */ + GaussKronrod21 = 2, /**< Use 10, 21 points. */ + GaussKronrod31 = 3, /**< Use 15, 31 points. */ + GaussKronrod41 = 4, /**< Use 20, 41 points. */ + GaussKronrod51 = 5, /**< Use 25, 51 points. */ + GaussKronrod61 = 6, /**< Use 30, 61 points. */ + GaussKronrod71 = 7, /**< Use 35, 71 points. */ + GaussKronrod81 = 8, /**< Use 40, 81 points. */ + GaussKronrod91 = 9, /**< Use 45, 91 points. */ + GaussKronrod101 = 10, /**< Use 50, 101 points. */ + GaussKronrod121 = 11, /**< Use 60, 121 points. */ + GaussKronrod201 = 12 /**< Use 100, 201 points. */ + }; + + /** + * \brief Prepares an Integrator for a call to a quadrature function. + * + * \param[in] maxSubintervals The maximum number of subintervals allowed in the subdivision process + * of quadrature functions. This corresponds to the amount of memory allocated for said + * functions. + */ + Integrator(const int maxSubintervals) + : m_maxSubintervals(maxSubintervals) + { + assert(maxSubintervals >= 1); // \todo use Eigen assert. + + m_errorListIndices.resize(maxSubintervals, 1); + m_lowerList.resize(maxSubintervals, 1); + m_upperList.resize(maxSubintervals, 1); + m_integralList.resize(maxSubintervals, 1); + m_errorList.resize(maxSubintervals, 1); + } + + /** + * \brief This function calculates an approximation I' to a given definite integral I, the + * integral of f from lowerLimit to upperLimit, hopefully satisfying + * abs(I - I') <= max(desiredAbsoluteError, desiredRelativeError * abs(I)). + * + * This function is best suited for integrands without singularities or discontinuities, which + * are too difficult for non-adaptive quadrature, and, in particular, for integrands with + * oscillating behavior of a non-specific type. + * + * \param[in,out] functionType The function type defining the integrand function. + * \param[in] lowerLimit The lower limit of integration. + * \param[in] upperLimit The upper limit of integration. + * \param[in] desiredAbsoluteError The absolute accuracy requested. + * \param[in] desiredRelativeError The relative accuracy requested. + * If desiredAbsoluteError <= 0 and desiredRelativeError < 50 * machinePrecision, + * the routine will end with errorCode = 6. + * \param[in] quadratureRule The local Gauss-Kronrod quadrature rule to use. + * + * \returns The approximation to the integral. + */ + template + Scalar quadratureAdaptive(const FunctionType& functionType, + const Scalar lowerLimit, + const Scalar upperLimit, + const Scalar desiredAbsoluteError = Scalar(0.), + const Scalar desiredRelativeError = Scalar(0.), + const QuadratureRule quadratureRule = Eigen::Integrator::QuadratureRule(1)) + { + using std::abs; + using std::max; + + if ((desiredAbsoluteError <= Scalar(0.) && + desiredRelativeError < NumTraits::epsilon()) || + m_maxSubintervals < 1) + { + m_errorCode = 6; + return Scalar(0.); + } + + m_errorCode = 0; + m_numEvaluations = 0; + m_lowerList[0] = lowerLimit; + m_upperList[0] = upperLimit; + m_integralList[0] = Scalar(0.); + m_errorList[0] = Scalar(0.); + m_errorListIndices[0] = 0; + m_errorListIndices[1] = 1; + + Scalar absDiff = 0.; + Scalar absResult = 0.; + + // First approximation to the integral + Scalar integral = quadratureKronrod(functionType, lowerLimit, upperLimit, m_estimatedError, absDiff, absResult, quadratureRule); + + m_numSubintervals = 1; + m_integralList[0] = integral; + m_errorList[0] = m_estimatedError; + + // Test on accuracy. + Scalar errorBound = max(desiredAbsoluteError, desiredRelativeError * abs(integral)); + + if (m_maxSubintervals == 1) + { + m_errorCode = 1; + } + else if (m_estimatedError <= NumTraits::epsilon() * Scalar(50.) * absDiff + && m_estimatedError > errorBound) + { + m_errorCode = 2; + } + + if (m_errorCode != 0 || + m_estimatedError == Scalar(0.) || + (m_estimatedError <= errorBound && + m_estimatedError != absResult)) + { + if (quadratureRule == GaussKronrod15) + { + m_numEvaluations = m_numEvaluations * 30 + 15; + } + else + { + m_numEvaluations = (quadratureRule * 10 + 1) * (m_numEvaluations * 2 + 1); + } + + return integral; + } + + // The sum of the integrals over the subintervals. + Scalar area = integral; + + Scalar errorSum = m_estimatedError; + + // The maximum interval error. + Scalar errorMax = m_estimatedError; + + // An index into m_errorList at the interval with largest error estimate. + Index maxErrorIndex = 0; + Index nrMax = 0; + + int roundOff1 = 0; + int roundOff2 = 0; + + Scalar error1 = 0.; + Scalar error2 = 0.; + Scalar absDiff1 = 0.; + Scalar absDiff2 = 0.; + + // Main loop for the integration + for (m_numSubintervals = 2; m_numSubintervals <= m_maxSubintervals; ++m_numSubintervals) + { + const Index numSubintervalsIndex = m_numSubintervals - 1; + + // Bisect the subinterval with the largest error estimate. + const Scalar lower1 = m_lowerList[maxErrorIndex]; + const Scalar upper2 = m_upperList[maxErrorIndex]; + + const Scalar upper1 = (lower1 + upper2) * Scalar(.5); + const Scalar lower2 = upper1; + + const Scalar area1 = quadratureKronrod(functionType, lower1, upper1, error1, absResult, absDiff1, quadratureRule); + const Scalar area2 = quadratureKronrod(functionType, lower2, upper2, error2, absResult, absDiff2, quadratureRule); + + // Improve previous approximations to integral and error and test for accuracy. + ++(m_numEvaluations); + + const Scalar area12 = area1 + area2; + const Scalar error12 = error1 + error2; + + errorSum += error12 - errorMax; + area += area12 - m_integralList[maxErrorIndex]; + + if (absDiff1 != error1 && + absDiff2 != error2) + { + if (abs(m_integralList[maxErrorIndex] - area12) <= abs(area12) * Scalar(1.e-5) && + error12 >= errorMax * Scalar(.99)) + { + ++roundOff1; + } + + if (m_numSubintervals > 10 && + error12 > errorMax) + { + ++roundOff2; + } + } + + m_integralList[maxErrorIndex] = area1; + m_integralList[numSubintervalsIndex] = area2; + + errorBound = max(desiredAbsoluteError, desiredRelativeError * abs(area)); + + if (errorSum > errorBound) + { + // Set error flag in the case that the number of subintervals has reached the max allowable. + if (m_numSubintervals == m_maxSubintervals) + { + m_errorCode = 1; + } + // Test for roundoff error and set error flag. + else if (roundOff1 >= 6 || roundOff2 >= 20) + { + m_errorCode = 2; + } + // Set m_error_code in the case of poor integrand behaviour within the integration range. + else if (max(abs(lower1), abs(upper2)) <= + (NumTraits::epsilon() * Scalar(100.) + Scalar(1.)) + * (abs(lower2) + (std::numeric_limits::min)() * Scalar(1.e3))) + { + m_errorCode = 3; + } + } + + // Append the newly-created intervals to the list. + if (error2 > error1) + { + m_lowerList[numSubintervalsIndex] = lower1; + m_lowerList[maxErrorIndex] = lower2; + m_upperList[numSubintervalsIndex] = upper1; + m_integralList[maxErrorIndex] = area2; + m_integralList[numSubintervalsIndex] = area1; + m_errorList[maxErrorIndex] = error2; + m_errorList[numSubintervalsIndex] = error1; + } + else + { + m_lowerList[numSubintervalsIndex] = lower2; + m_upperList[maxErrorIndex] = upper1; + m_upperList[numSubintervalsIndex] = upper2; + m_errorList[maxErrorIndex] = error1; + m_errorList[numSubintervalsIndex] = error2; + } + + // Maintain the descending ordering in the list of error estimates and select the subinterval + // with the largest error estimate, (the next subinterval to be bisected). + quadratureSort(maxErrorIndex, errorMax, nrMax); + + if (m_errorCode != 0 || + errorSum <= errorBound || + m_numSubintervals == m_maxSubintervals) + { + break; + } + } + + integral = Scalar(0.); + + for (Index k = 0; k < m_numSubintervals; ++k) + { + integral += m_integralList[k]; + } + + m_estimatedError = errorSum; + + if (quadratureRule == GaussKronrod15) + { + m_numEvaluations = m_numEvaluations * 30 + 15; + } + else + { + m_numEvaluations = (quadratureRule * 10 + 1) * (m_numEvaluations * 2 + 1); + } + + return integral; + } + + /** + * \brief Returns the estimated absolute error from the last integration. + * + * \returns The value returned will only be valid after calling quadratureAdaptive at least once. + */ + inline Scalar estimatedError() + { + return m_estimatedError; + } + + /** + * \brief Returns the error code. + * + * \returns The value returned will only be valid after calling quadratureAdaptive at least once. + */ + inline Index errorCode() + { + return m_errorCode; + } + + private: + + /** + * \brief This routine maintains the descending ordering in the list of the local error + * estimates resulting from the interval subdivision process. + * + * At each call two error estimates are inserted using the sequential search method, top-down + * for the largest error estimate and bottom-up for the smallest error estimate. + * + * \param[in,out] maxErrorIndex The index to the nrMax-th largest error estimate currently in the list. + * \param[in,out] errorMax The nrMax-th largest error estimate. errorMaxIndex = errorList(maxError). + * \param[in,out] nrMax The integer value such that maxError = errorListIndices(nrMax). + */ + void quadratureSort(Index& maxErrorIndex, + Scalar& errorMax, + Index& nrMax) + { + if (m_numSubintervals <= 2) + { + m_errorListIndices[0] = 0; + m_errorListIndices[1] = 1; + maxErrorIndex = m_errorListIndices[nrMax]; + errorMax = m_errorList[maxErrorIndex]; + return; + } + + // This part of the routine is only executed if, due to a difficult integrand, subdivision has + // increased the error estimate. In the normal case the insert procedure should start after the + // nrMax-th largest error estimate. + Index i = 0; + Index succeed = 0; + const Scalar errorMaximum = m_errorList[maxErrorIndex]; + + if (nrMax != 1) + { + for (i = 1; i < nrMax; ++i) + { + succeed = m_errorListIndices[nrMax - 1]; + + if (errorMaximum <= m_errorList[succeed]) + { + break; + } + + m_errorListIndices[nrMax] = succeed; + --nrMax; + } + } + + // Compute the number of elements in the list to be maintained in descending order. This number + // depends on the number of subdivisions remaining allowed. + Index topBegin = m_numSubintervals - 1; + Index bottomEnd = topBegin - 1; + Index start = nrMax + 1; + + if (m_numSubintervals > m_maxSubintervals / 2 + 2) + { + topBegin = m_maxSubintervals + 3 - m_numSubintervals + 1; + } + + // Insert errorMax by traversing the list top-down, starting comparison from the element + // errorlist(m_errorListIndices(nrMax+1)). + if (start <= bottomEnd) + { + for (i = start; i <= bottomEnd; ++i) + { + succeed = m_errorListIndices[i]; + + if (errorMaximum >= m_errorList[succeed]) + { + break; + } + + m_errorListIndices[i - 1] = succeed; + } + } + + if (start > bottomEnd) + { + m_errorListIndices[bottomEnd] = maxErrorIndex; + m_errorListIndices[topBegin] = m_numSubintervals - 1; + maxErrorIndex = m_errorListIndices[nrMax]; + errorMax = m_errorList[maxErrorIndex]; + return; + } + + // Insert errorMin by traversing the list bottom-up. + m_errorListIndices[i - 1] = maxErrorIndex; + + Index tempIndex = bottomEnd; + + for (Index j = i; j <= bottomEnd; ++j) + { + succeed = m_errorListIndices[tempIndex]; + + if (m_errorList[m_numSubintervals - 1] < m_errorList[succeed]) + { + m_errorListIndices[tempIndex + 1] = m_numSubintervals - 1; + maxErrorIndex = m_errorListIndices[nrMax]; + errorMax = m_errorList[maxErrorIndex]; + return; + } + + m_errorListIndices[tempIndex + 1] = succeed; + --tempIndex; + } + + m_errorListIndices[i] = m_numSubintervals - 1; + maxErrorIndex = m_errorListIndices[nrMax]; + errorMax = m_errorList[maxErrorIndex]; + return; + } + + /** + * \brief This function calculates an approximation I' to a given definite integral I, the + * integral of f from lowerLimit to upperLimit and provides an error estimate. + * + * \param[in] f The variable representing the function f(x) be integrated. + * \param[in] lowerLimit The lower limit of integration. + * \param[in] upperLimit The upper limit of integration. + * \param[in,out] errorEstimate Estimate of the modulus of the absolute error, not to exceed + * abs(I - I'). + * \param[in,out] absIntegral The approximation to the integral of abs(f) from lowerLimit to + * upperLimit. + * \param[in,out] absDiffIntegral The approximation to the integral of + * abs(f - I/(upperLimit - lowerLimit)). + * + * \returns The approximation I' to the integral I. It is computed by applying the 15, 21, 31, + * 41, 51, 61, 71, 81, 91, 101, 121, 201-point kronrod rule obtained by optimal addition + * of abscissae to the 7, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 100-point Gauss rule. + * + * \detail This series of functions represents a priority queue data structure: + * - Apply Gauss-Kronrod on the whole initial interval and estimate the error. + * - Split the interval symmetrically in two and again apply Gauss-Kronrod on the intervals, + * estimate the errors and add a pair (or tuple) "interval, error" to the priority queue. + * - If the total error is smaller than the requested tolerance, or if the maximum number + * of subdivisions is reached, discontinue the process. + * - Otherwise, pop the top element from the queue (highest error), split the interval in two, + * and apply Gauss-Kronrod to the new intervals, add those two new elements to the priority + * queue and repeat from the previous step. + */ + template + Scalar quadratureKronrod( + const FunctionType& functionType, const Scalar lowerLimit, const Scalar upperLimit, + Scalar& estimatedError, Scalar& absIntegral, Scalar& absDiffIntegral, + const QuadratureRule quadratureRule) + { + switch (quadratureRule) + { + case GaussKronrod15: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod15, + QuadratureKronrod::weightsGaussKronrod15, QuadratureKronrod::weightsGauss15, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod21: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod21, + QuadratureKronrod::weightsGaussKronrod21, QuadratureKronrod::weightsGauss21, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod31: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod31, + QuadratureKronrod::weightsGaussKronrod31, QuadratureKronrod::weightsGauss31, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod41: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod41, + QuadratureKronrod::weightsGaussKronrod41, QuadratureKronrod::weightsGauss41, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod51: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod51, + QuadratureKronrod::weightsGaussKronrod51, QuadratureKronrod::weightsGauss51, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod61: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod61, + QuadratureKronrod::weightsGaussKronrod61, QuadratureKronrod::weightsGauss61, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod71: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod71, + QuadratureKronrod::weightsGaussKronrod71, QuadratureKronrod::weightsGauss71, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod81: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod81, + QuadratureKronrod::weightsGaussKronrod81, QuadratureKronrod::weightsGauss81, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod91: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod91, + QuadratureKronrod::weightsGaussKronrod91, QuadratureKronrod::weightsGauss91, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod101: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod101, + QuadratureKronrod::weightsGaussKronrod101, QuadratureKronrod::weightsGauss101, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod121: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod121, + QuadratureKronrod::weightsGaussKronrod121, QuadratureKronrod::weightsGauss121, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + case GaussKronrod201: + return quadratureKronrodHelper( + QuadratureKronrod::abscissaeGaussKronrod201, + QuadratureKronrod::weightsGaussKronrod201, QuadratureKronrod::weightsGauss201, + functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule); + + default: + return Scalar(0.); + } + } + + template + Scalar quadratureKronrodHelper(Array abscissaeGaussKronrod, + Array weightsGaussKronrod, + Array weightsGauss, + const FunctionType& functionType, + const Scalar lowerLimit, + const Scalar upperLimit, + Scalar& estimatedError, + Scalar& absIntegral, + Scalar& absDiffIntegral, + const QuadratureRule quadratureRule) + { + using std::abs; + using std::min; + using std::max; + using std::pow; + + // Half-length of the interval. + const Scalar halfLength = (upperLimit - lowerLimit) * Scalar(.5); + + // Midpoint of the interval. + const Scalar center = (lowerLimit + upperLimit) * Scalar(.5); + const Scalar fCenter = functionType(center); + + Index size1 = weightsGaussKronrod.size() - 1; + Index size2 = weightsGauss.size() - 1; + + Array f1Array; + Array f2Array; + + // The result of the Gauss formula. + Scalar resultGauss = Scalar(0.); + + if (quadratureRule % 2 != 0) + { + resultGauss = weightsGauss[size2] * fCenter; + } + + // The result of the Kronrod formula. + Scalar resultKronrod = weightsGaussKronrod[size1] * fCenter; + + absIntegral = abs(resultKronrod); + + for (Index j = 1; j < weightsGaussKronrod.size() - weightsGauss.size(); ++j) + { + const Index jj = j * 2 - 1; + const Scalar abscissa = halfLength * abscissaeGaussKronrod[jj]; + + const Scalar f1 = functionType(center - abscissa); + const Scalar f2 = functionType(center + abscissa); + + f1Array[jj] = f1; + f2Array[jj] = f2; + + const Scalar funcSum = f1 + f2; + + resultGauss += weightsGauss[j - 1] * funcSum; + resultKronrod += weightsGaussKronrod[jj] * funcSum; + + absIntegral += weightsGaussKronrod[jj] * (abs(f1) + abs(f2)); + } + + for (Index j = 0; j < weightsGauss.size(); ++j) + { + const Index jj = j * 2; + const Scalar abscissa = halfLength * abscissaeGaussKronrod[jj]; + + const Scalar f1 = functionType(center - abscissa); + const Scalar f2 = functionType(center + abscissa); + + f1Array[jj] = f1; + f2Array[jj] = f2; + + const Scalar funcSum = f1 + f2; + + resultKronrod += weightsGaussKronrod[jj] * funcSum; + + absIntegral += weightsGaussKronrod[jj] * (abs(f1) + abs(f2)); + } + + // Approximation to the mean value of f over the interval (lowerLimit, upperLimit), + // i.e. I / (upperLimit - lowerLimit) + Scalar resultMeanKronrod = resultKronrod * Scalar(.5); + + absDiffIntegral = weightsGaussKronrod[size1] * (abs(fCenter - resultMeanKronrod)); + + absDiffIntegral += (((f1Array.head(size1) - resultMeanKronrod).abs() + + (f2Array.head(size1) - resultMeanKronrod).abs()) + * weightsGaussKronrod.head(size1)).sum(); + + Scalar result = resultKronrod * halfLength; + absIntegral *= abs(halfLength); + absDiffIntegral *= abs(halfLength); + estimatedError = abs((resultKronrod - resultGauss) * halfLength); + + if (absDiffIntegral != Scalar(0.) && + estimatedError != Scalar(0.)) + { + estimatedError = absDiffIntegral * min(Scalar(1.), pow((estimatedError * Scalar(200.) / absDiffIntegral), Scalar(1.5))); + } + + if (absIntegral > (std::numeric_limits::min)() / (NumTraits::epsilon() * Scalar(50.))) + { + estimatedError = max(NumTraits::epsilon() * Scalar(50.) * absIntegral, estimatedError); + } + + return result; + } + + + /** + * \brief An Array of dimension m_maxSubintervals for error estimates. + * + * The first k elements are indices to the error estimates over the subintervals, such that + * errorList(errorListIndices(0)), ..., errorList(errorListIndices(k - 1)) forms a decreasing + * sequence, with k = m_numSubintervals if m_numSubintervals <= (m_maxSubintervals/2 + 2), + * otherwise k = m_maxSubintervals + 1 - m_numSubintervals. + */ + Array m_errorListIndices; + + /** + * \brief An Array of dimension m_maxSubintervals for subinterval left endpoints. + * + * The first m_numSubintervals elements are the lower end points of the subintervals in the + * partition of the given integration range (lowerLimit, upperLimit). + */ + Array m_lowerList; + + /** + * \brief An Array of dimension m_maxSubintervals for subinterval upper endpoints. + * + * The first m_numSubintervals elements are the upper end points of the subintervals in the + * partition of the given integration range (lowerLimit, upperLimit). + */ + Array m_upperList; + + /** + * \brief An Array of dimension m_maxSubintervals for integral approximations. + * + * The first m_numSubintervals elements are the integral approximations on the subintervals. + */ + Array m_integralList; + + /** + * \brief An Array of dimension m_maxSubintervals for error estimates. + * + * The first m_numSubintervals elements of which are the moduli of the absolute error estimates + * on the subintervals. + */ + Array m_errorList; + + /** + * \brief Gives an upper bound on the number of subintervals. Must be at least 1. + */ + Index m_maxSubintervals; + + /** + * \brief The number of integrand evaluations. + */ + Index m_numEvaluations; + + /** + * \brief Estimate of the modulus of the absolute error, which should equal or exceed abs(I - I'). + */ + Scalar m_estimatedError; + + /** + * \brief Error messages generated by the routine. + * + * errorCode = 0 Indicates normal and reliable termination of the routine. (It is assumed that + * the requested accuracy has been achieved.) + * errorCode > 0 Any errorCode greater than zero indicates abnormal termination of the routine. + * (The estimates for integral and m_estimatedError are less reliable and the + * requested accuracy has not been achieved.) + * errorCode = 1 The maximum number of subdivisions allowed has been achieved. One can allow more + * subdivisions by increasing the value of m_maxSubintervals. However, if this + * yields no improvement it is advised to analyze the integrand in order to + * determine the integration difficulaties. If the position of a local difficulty + * can be determined, (i.e. singularity or discontinuity within the interval), one + * will probably gain from splitting up the interval at this point and calling the + * integrator on the subranges. If possible, an appropriate special-purpose + * integrator should be used which is designed for handling the type of difficulty + * involved. + * errorCode = 2 The occurrence of roundoff error is detected, preventing the requested tolerance + * from being achieved. + * errorCode = 3 Extremely bad integrand behaviour occurs at points in the integration interval. + * errorCode = 4 Roundoff error on extrapolation + * errorCode = 5 Divergent integral (or very slowly convergent integral) + * errorCode = 6 The input is invalid, because (desiredAbsoluteError <= 0 and + * desiredRealtiveError < 50 * relativeMachineAccuracy, or + * m_maxSubintervals < 1. + * errorCode = 7 Applies to (D)QAWF only - limiting number of cycles has been attained + * + * \todo make relativeMachineAccuracy a member variable. + */ + Index m_errorCode; + + /** + * \brief The number of subintervals actually produced in the subdivision process. + */ + Index m_numSubintervals; + + }; +} +#endif // EIGEN_INTEGRATOR_H diff --git a/cpp/LaurieGautschi.h b/cpp/LaurieGautschi.h new file mode 100644 --- /dev/null +++ b/cpp/LaurieGautschi.h @@ -0,0 +1,481 @@ +#ifndef EIGEN_LAURIEGAUTSCHI_H +#define EIGEN_LAURIEGAUTSCHI_H + +namespace Eigen +{ + + /** + * \ingroup NumericalIntegration_Module + * + * \class LaurieGautschi + * + * \tparam Scalar floating point type + * + * \brief This class computes Kronrod abscissae & weights for arbitrary precision. + * + * Based on work of Dirk Laurie and Walter Gautschi. + * D. P. Laurie (1997). Calculation of Gauss-Kronrod Quadrature Rules. + * Mathematics of Computation, 66(219), 1133-1145 + * Created by Pavel Holoborodko, November 7, 2011. + * Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder, and Matt Beall September 2014 + * + * \TODO Ensure only appropriates types are used for Scalar, e.g. prohibit integers. + */ + template + + class LaurieGautschi + { + public: + + /** + * \brief Recurrence coefficients for monic Jacobi polynomials. + * + * This method generates the first N recurrence + * coefficients for monic Jacobi polynomials with parameters + * alpha and beta. These are orthogonal on [-1,1] relative to the + * weight function w(t)=(1-t)^a(1+t)^b. The N alpha-coefficients + * are stored in \a alphaOut, the n beta-coefficients in \a betaOut. + * http://en.wikipedia.org/wiki/Jacobi_polynomials + * + * Created by Dirk Laurie, 6-22-1998; edited by Walter Gautschi, 4-4-2002. + * Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder and Matt Beall, September 2014 + * + * \param[in] N Number of recurrence coefficients + * \param[in] alpha Alpha parameter of the Jacobi-polynomials + * \param[in] beta Beta parameter of the Jacobi-polynomials + * \param[in,out] alphaOut N alpha-coefficients + * \param[in,out] betaOut N beta-coefficients + */ + static void r_jacobi(const Index N, + const Scalar alpha, + const Scalar beta, + Eigen::Array& alphaOut, + Eigen::Array& betaOut) + { + using std::pow; + using std::tgamma; + + //TODO : make use the eigen assert facilities + assert(alpha > Scalar(-1)); + assert(beta > Scalar(-1)); + assert(alphaOut.rows() == betaOut.rows()); + assert(alphaOut.rows() > 0); + assert(N <= alphaOut.rows()); + + alphaOut(0) = (beta-alpha)/(alpha+beta+Scalar(2.)); + betaOut(0) = pow(Scalar(2.),(alpha+beta+Scalar(1.))) * tgamma(alpha+Scalar(1.)) * + tgamma(beta+Scalar(1.)) / tgamma(alpha+beta+Scalar(2.)); + + for (Index n = 1; n < N; ++n) + { + Scalar nAlphaBeta = Scalar(2.) * n + alpha + beta; + alphaOut(n) = (beta*beta - alpha*alpha) / (nAlphaBeta * (nAlphaBeta+Scalar(2.))); + betaOut(n) = Scalar(4.) * (n+alpha)*(n+beta)*n*(n+alpha+beta) / + (nAlphaBeta*nAlphaBeta*(nAlphaBeta+Scalar(1.)) * (nAlphaBeta-Scalar(1.))); + } + } + + /** + * \brief Recurrence coefficients for monic Jacobi polynomials on [0,1]. + * + * This method generates the first N recurrence + * coefficients for monic Jacobi polynomials on [0,1] with + * parameters alpha and beta. These are orthogonal on [0,1] relative + * to the weight function w(t)=(1-t)^alpha t^beta. The N alpha- + * coefficients are stored in \a alphaOut, the N beta- + * coefficients in \a betaOut. + * http://en.wikipedia.org/wiki/Jacobi_polynomials + * + * Created by Dirk Laurie, 6-22-1998; edited by Walter Gautschi, 4-4-2002. + * Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder and Matt Beall, September 2014 + * + * \param[in] N Number of recurrence coefficients + * \param[in] alpha Alpha parameter of the Jacobi-polynomials + * \param[in] beta Beta parameter of the Jacobi-polynomials + * \param[in,out] alphaOut N alpha-coefficients + * \param[in,out] betaOut N beta-coefficients + */ + static void r_jacobi_01(const Index N, + const Scalar alpha, + const Scalar beta, + Eigen::Array& alphaOut, + Eigen::Array& betaOut) + { + using std::pow; + + //TODO : make use the eigen assert facilities + assert(alpha > Scalar(-1)); + assert(beta > Scalar(-1)); + assert(alphaOut.rows() == betaOut.rows()); + assert(alphaOut.rows() > 0); + assert(N <= alphaOut.rows()); + + r_jacobi(N, alpha, beta, alphaOut, betaOut); + + for (Index n = 0; n < N; ++n) + { + alphaOut(n) = (Scalar(1.)+alphaOut(n)) / Scalar(2.); + } + + betaOut(0) = betaOut(0) / pow(Scalar(2), alpha + beta + Scalar(1.)); + + for (Index n = 1; n < N; ++n) + { + betaOut(n) = betaOut(n) / Scalar(4.); + } + } + + /** + * \brief Jacobi-Kronrod matrix. + * + * This method produces the alpha- and beta-elements in + * the Jacobi-Kronrod matrix of order 2N+1 for the weight + * function (or measure) w. The input data for the weight + * function w are the recurrence coefficients of the associated + * orthogonal polynomials, which are stored in \a alphaIn and \a betaIn . + * At least ceil(3*N/2)+1 coefficients should be provided. + * The 2N+1 alpha- and beta-elements are returned in \a alpha and \a beta + * respectively. + * + * Created by Dirk Laurie, 6-22.1998 + * Edited by Pavel Holoborodko, November 7, 2011 + * Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder and Matt Beall, September 2014 + * + * \param[in] N Number of nodes + * \param[in,out] alphaIn The recurrence coefficients of the associated orthogonal polynomials + * \param[in,out] betaIn The recurrence coefficients of the associated orthogonal polynomials + * \param[in,out] alpha Alpha-elements in the Jacobi-Kronrod matrix of order 2N+1 + * \param[in,out] beta Beta-elements in the Jacobi-Kronrod matrix of order 2N+1 + * + */ + static void r_kronrod(const Index N, + const Eigen::Array& alphaIn, + const Eigen::Array& betaIn, + Eigen::Array& alpha, + Eigen::Array& beta) + { + using std::ceil; + using std::floor; + + //TODO : make use the eigen assert facilities + assert(alphaIn.rows() == betaIn.rows()); + assert(alphaIn.rows() >= ceil(3*N/2) + 1 ); + assert(alphaIn.rows() > 0); + assert(alpha.rows() == 2*N+1); + assert(beta.rows() == 2*N+1); + + alpha = Eigen::Array::Zero(2*N+1); + beta = Eigen::Array::Zero(2*N+1); + + for (Index k = 0; k <= (Index)floor(3*N/2) + 1; ++k) + { + alpha(k) = alphaIn(k); + } + + for (Index k=0; k <= (Index)ceil(3*N/2) + 1; ++k) + { + beta(k) = betaIn(k); + } + + Eigen::Array sigma = Eigen::Array::Zero((Index)floor(N/2) + 2); + Eigen::Array tempVector = Eigen::Array::Zero((Index)floor(N/2) + 2); + + tempVector(1) = beta(N+1); + + for (Index m = 0; m < N-2+1; ++m) + { + Scalar u = 0; + for(Index k = (Index)floor((m+1) / 2); k >= 0;--k) + { + Index l = m-k; + u = u + (alpha(k+N+1)-alpha(l)) * tempVector(k+1) + beta(k+N+1)*sigma(k) - beta(l)*sigma(k+1); + sigma(k+1) = u; + } + + Eigen::Array swap = sigma; + sigma = tempVector; + tempVector = swap; + } + + for (Index j = (Index)floor(N/2); j>=0; --j) + { + sigma(j+1) = sigma(j); + } + + for (Index m = N-1; m < 2*N-3+1; ++m) + { + Index k = m+1-N; + Index j = 0; + Scalar u = 0; + + for (k = m+1-N; k < (Index)floor((m-1) / 2) + 1; ++k) + { + Index l = m-k; + j = N-1-l; + u = u - (alpha(k+N+1)-alpha(l))*tempVector(j+1) - beta(k+N+1)*sigma(j+1) + beta(l)*sigma(j+2); + sigma(j+1) = u; + } + + k = (Index)floor((m+1) / 2); + + if (m % 2 == 0) + { + alpha(k+N+1) = alpha(k) + (sigma(j+1)-beta(k+N+1)*sigma(j+2)) / tempVector(j+2); + } + else + { + beta(k+N+1) = sigma(j+1) / sigma(j+2); + } + + Eigen::Array swap = sigma; + sigma = tempVector; + tempVector = swap; + } + + alpha(2*N) = alpha(N-1)-beta(2*N)*sigma(1) / tempVector(1); + } + + /** + * \brief Gauss-Kronrod quadrature formula. + * + * This method generates the (2N+1)-point Gauss-Kronrod + * quadrature rule for the weight function w encoded by the + * recurrence matrix (alpha,beta) of order [ceil(3*n/2)+1]x2 containing + * in its first and second column respectively the alpha- and + * beta-coefficients in the three-term recurrence relation + * for w. The 2N+1 nodes, in increasing order, are output + * into \a nodes, the corresponding weights into \a weights. + * + * Created by Dirk Laurie, June 22, 1998. + * Edited by Pavel Holoborodko, November 7, 2011: + * Ported to C++/Eigen by Sreekumar Thaithara Balan, Mark Sauder + * and Matt Beall, September 2014 + * + * \param[in] N Number of nodes + * \param[in,out] alpha 2N alpha coefficients (input) + * \param[in,out] beta 2N beta coefficients (input) + * \param[in,out] nodes 2N+1 nodes + * \param[in,out] weights 2N+1 weights corresponding to \a nodes + */ + static void kronrod(const Index N, + const Eigen::Array& alpha, + const Eigen::Array& beta, + Eigen::Array& nodes, + Eigen::Array& weights) + { + using std::abs; + using std::sqrt; + + //TODO : make use the eigen assert facilities + assert(N>0); + assert(alpha.rows() == 2*N); + assert(alpha.rows() == beta.rows()); + assert(nodes.rows() == 2*N+1); + assert(nodes.rows() == weights.rows()); + + Eigen::Array alpha0 = Eigen::Array::Zero(2*N+1); + Eigen::Array beta0 = Eigen::Array::Zero(2*N+1); + + r_kronrod(N, alpha, beta, alpha0, beta0); + + // \TODO : CHECK NEEDED LIKE THE ONE ON LINE 21 IN KRONROD.M + // Do we have an approximately equal function in Eigen? + assert(abs(beta0.sum() - (Scalar) (2*N+1)) > 1e-5); + + Eigen::Matrix J = Eigen::Matrix::Zero(2*N+1, 2*N+1); + + for (Index k = 0; k < 2*N; ++k) + { + J(k,k) = alpha0(k); + J(k,k+1) = sqrt(beta0(k+1)); + J(k+1,k) = J(k,k+1); + } + + J(2*N,2*N) = alpha0(2*N); + + //TODO : Is this assumption of positive definiteness correct? + SelfAdjointEigenSolver< Eigen::Matrix > es(J); + + //TODO : make use the eigen assert facilities + assert(es.info() == Eigen::Success); + + nodes = es.eigenvalues(); + Eigen::Matrix V = es.eigenvectors(); + + weights = beta0(0) * (V.row(0).array() * V.row(0).array()).matrix(); + + } + + /** + * \brief Gauss quadrature rule. + * + * Given a weight function w encoded by (alpha,beta) of the + * first N recurrence coefficients for the associated orthogonal + * polynomials, the first column of (alpha,beta) containing the N alpha- + * coefficients and the second column the N beta-coefficients, + * the method generates the nodes and weights of + * the N-point Gauss quadrature rule for the weight function. + * The Nodes, in increasing order, are stored in \a nodes , + * the N corresponding weights are stored in \a weights . + * + * \param[in] N Number of nodes + * \param[in,out] alpha 2N alpha coefficients (input) + * \param[in,out] beta 2N beta coefficients (input) + * \param[in,out] nodes 2N+1 nodes + * \param[in,out] weights 2N+1 weights corresponding to \a nodes + */ + static void gauss(const Index N, + const Eigen::Array& alpha, + const Eigen::Array& beta, + Eigen::Array& nodes, + Eigen::Array& weights) + { + using std::sqrt; + + //TODO : make use the eigen assert facilities + assert(N > 0); + assert(alpha.rows() == 2*N); + assert(alpha.rows() == beta.rows()); + assert(nodes.rows() == N); + assert(nodes.rows() == weights.rows()); + + Eigen::Matrix J = Eigen::Matrix::Zero(N,N); + + J(0,0) = alpha(0); + + for (Index n=1; n > es(J); + + //TODO : make use the eigen assert facilities + assert(es.info() == Eigen::Success); + + nodes = es.eigenvalues(); + Eigen::Matrix V = es.eigenvectors(); + + weights = beta(0)*(V.row(0).array()*V.row(0).array()).matrix(); + } + + /** + * \brief Arbitrary precision Kronrod abscissae & weights. + * + * This method computes Kronrod points for (-1,1) with any required precision. + * The result is a vector 2N+1 nodes (N points on either side of zero and zero) + * and the corresponding weights. + * + * \param[in] N Number of nodes + * \param[in,out] nodes Returns a vector of 2N+1 nodes + * \param[in,out] w Returns a vector of weights corresponding to \a nodes + */ + static void mpkronrod(const Index N, + Eigen::Array& nodes, + Eigen::Array& weights) + { + //TODO : make use the eigen assert facilities + assert(nodes.rows() == 2*N+1); + assert(weights.rows() == 2*N+1); + assert(N>0); + + Eigen::Array alpha = Eigen::Array::Zero(2*N); + Eigen::Array beta = Eigen::Array::Zero(2*N); + + r_jacobi_01( 2*N, Scalar(0), Scalar(0), alpha, beta); + + kronrod(N,alpha,beta,nodes,weights); + + for (Index i=0; i& nodes, + Eigen::Array& weights) + { + //TODO : make use the eigen assert facilities + assert(nodes.rows() == N); + assert(weights.rows() == N); + assert(N > 0); + + Eigen::Array alpha = Eigen::Array::Zero(2*N); + Eigen::Array beta = Eigen::Array::Zero(2*N); + + r_jacobi_01(2*N, Scalar(0), Scalar(0), alpha, beta); + + gauss(N, alpha, beta, nodes, weights); + + for (Index i = 0;i& abscGaussKronrod, + Eigen::Array& weightGaussKronrod, + Eigen::Array& abscGauss, + Eigen::Array& weightGauss) + { + using std::abs; + + Eigen::Array xGK = Eigen::Array::Zero(2*nNodes+1); + Eigen::Array wGK = Eigen::Array::Zero(2*nNodes+1); + Eigen::Array xG = Eigen::Array::Zero(nNodes); + Eigen::Array wG = Eigen::Array::Zero(nNodes); + + LaurieGautschi::mpkronrod(nNodes,xGK,wGK); + LaurieGautschi::mpgauss(nNodes,xG,wG); + + unsigned int arraySize = nNodes + 1; + + abscGaussKronrod = Eigen::Array::Zero(arraySize); + weightGaussKronrod = Eigen::Array::Zero(arraySize); + abscGauss = Eigen::Array::Zero(arraySize/2); + weightGauss = Eigen::Array::Zero(arraySize/2); + + for (unsigned int i=0; i which is distributed under + * GNU General Public License + */ + template + class Monegato + { + public: + + /** + * \brief Compute the absolute value of a scalar + * \param[in] x inpute scalar value + * \return absolute value of scalar input + */ + static Scalar abs(Scalar x) + { + return (x < Scalar(0)) ? -(x) : x; + } + + /** + * \brief Compute the Legendre polynomials and their error + * \param[in] n degree + * \param[in] x value of the variable + * \param[out] err error to be returned + * \return value of the Legendre polynomial + * + * This function is based on the routine gsl_sf_legendre_Pl_e + * distributed with GSL + */ + static Scalar legendre_err(const int n, const Scalar x, Scalar& err) + { + if (n == 0) + { + err = Scalar(0); + return Scalar(1); + } + else if (n == 1) + { + err = Scalar(0); + return x; + } + + // below Sree modified this to avoid -Wmaybe-uninitialized + Scalar P0 = Scalar(1), P1 = x, P2=x; + Scalar E0 = NumTraits::epsilon(); + Scalar E1 = abs(x) * NumTraits::epsilon(); + for (int k = 1; k < n; ++k) + { + P2 = ((2*k + 1) * x * P1 - k * P0) / (k + 1); + err = ((2*k + 1) * abs(x) * E1 + k * E0) / (2*(k + 1)); + P0 = P1; P1 = P2; + E0 = E1; E1 = err; + } + return P2; + } + + /** + * \brief Three-term recursion identity for the Legendre derivatives + * \param[in] n degree + * \param[in] x value of the variable + * \return value of the Legendre derivative + */ + static Scalar legendre_deriv(int const n, Scalar const x) + { + if (n == 0) + { + return Scalar(0); + } + else if (n == 1) + { + return Scalar(1); + } + + Scalar P0 = Scalar(1); + Scalar P1 = x; + Scalar P2 = P1; + Scalar dP0 = Scalar(0); + Scalar dP1 = Scalar(1); + Scalar dP2 = P1; + + for (int k = 1; k < n; ++k) + { + P2 = ((2*k + 1) * x * P1 - k * P0) / (k + Scalar(1)); + dP2 = (2*k + 1) * P1 + dP0; + P0 = P1; + P1 = P2; + dP0 = dP1; + dP1 = dP2; + } + return dP2; + } + + /** + * \brief Compute derivatives of Chebyshev polynomials + * \param[in] x variable + * \param[in] n_ degree + * \param[in] coefs Chebyshev coefficients + * \return Derivative of Chebyshev polynomials + */ + static Scalar chebyshev_series_deriv(const Scalar x, + const int n_, + const Eigen::Array& coefs) + { + Scalar d1(0), d2(0); + Scalar y2 = 2 * x; // linear term for Clenshaw recursion + + for (int k = n_; k >= 2; --k) + { + Scalar temp = d1; + d1 = y2 * d1 - d2 + k * coefs(k); + d2 = temp; + } + + return y2 * d1 - d2 + coefs(1); + } + + /** + * \breif Evaluation of the Chebyshev polynomial using Clenshaw recursion + * \param[in] x variable + * \param[in] n_ degree + * \param[in] coefs Chebyshev coefficients + * \param[out] err error to be returned + * \return value of Chebyshev polynomial + */ + static Scalar chebyshev_series(const Scalar x, + const int n_, + const Eigen::Array& coefs, + Scalar& err) + { + Scalar d1(0), d2(0); + Scalar absc = abs(coefs(0)); // final term for truncation error + Scalar y2 = 2 * x; // linear term for Clenshaw recursion + + for (int k = n_; k >= 1; --k) + { + Scalar temp = d1; + d1 = y2 * d1 - d2 + coefs(k); + d2 = temp; + absc += abs(coefs(k)); + } + + err = absc * NumTraits::epsilon(); + return x * d1 - d2 + coefs(0)/2.; + } + + /** + * \brief Computes the zeros of the Legendre polynomial + * \param[in] m_ degree + * \param[out] zeros the zeros of the Legendre polynomial + */ + static void legendre_zeros(const int m_, + Eigen::Array& zeros) + { + Eigen::Array temp = Eigen::Array::Zero(m_+1); + zeros(0) = Scalar(-1); + zeros(1) = Scalar(1); + Scalar delta, epsilon; + + for (int k = 1; k <= m_; ++k) + { + // Loop to locate zeros of P_k interlacing z_0,...,z_k + for (int j = 0; j < k; ++j) + { + // Newton's method for P_k : + // initialize solver at midpoint of (z_j, z_{j+1}) + delta = 1; + Scalar x_j = (zeros(j) + zeros(j+1)) / 2.; + Scalar P_k = legendre_err(k, x_j, epsilon); + + while (abs(P_k) > epsilon && + abs(delta) > NumTraits::epsilon()) + { + delta = P_k / legendre_deriv(k, x_j); + x_j -= delta; + P_k = legendre_err(k, x_j, epsilon); + } + temp(j) = x_j; + } + + // Copy roots tmp_0,...,tmp_{k-1} to z_1,...z_k: + zeros(k+1) = zeros(k); + + for (int j = 0; j < k; ++j) + { + zeros(j+1) = temp(j); + } + } + } + + /** + * \brief Computes coefficients of the Chebyshev polynomial. + * \param[in] m_ degree of the Chebyshev polynomial. + * \param[out] coefs coefficients of the Chebyshev polynomial. + */ + static void chebyshev_coefs(const int m_, + Eigen::Array& coefs) + { + int ell = (m_ + 1)/2; + Eigen::Array alpha = Eigen::Array::Zero(ell+1); + Eigen::Array f = Eigen::Array::Zero(ell+1); + + // Care must be exercised in initalizing the constants in the definitions. + // Compilers interpret expressions like "(2*k + 1.0)/(k + 2.0)" as floating + // point precision, before casting to Real. + + f(1) = Scalar(m_+1) / Scalar(2*m_ + 3); + alpha(0) = Scalar(1); // coefficient of T_{m+1} + alpha(1) = -f(1); + + for (int k = 2; k <= ell; ++k) + { + f(k) = f(k-1) * (2*k - 1) * (m_ + k) / (k * (2*m_ + 2*k + 1)); + alpha(k) = -f(k); + + for (int i = 1; i < k; ++i) + { + alpha(k) -= f(i) * alpha(k-i); + } + } + + for (int k = 0; k <= ell; ++k) + { + coefs(m_ + 1 - 2*k) = alpha(k); + + if (m_ >= 2*k) + { + coefs(m_ - 2*k) = Scalar(0); + } + } + } + + /** + * \brief Compute the Gauss/Kronrod abscissae + * \param[in] n_ size of Gauss-Kronrod arrays + * \param[in] Gauss-Legendre degree + * \param[in] coefs Chebyshev coefficients + * \param[out] xgk_ Gauss/Kronrod abscissae + */ + static void gauss_kronrod_abscissae(const int n_, + const int m_, + const Eigen::Array& zeros, + const Eigen::Array& coefs, + Eigen::Array& xgk_) + { + // Now from the function gauss_kronrod_abscissae + Scalar delta, epsilon; + + for (int k = 0; k < n_ / 2; ++k) + { + delta = 1; + + // Newton's method for E_{n+1} : + Scalar x_k = (zeros(m_-k) + zeros(m_+1-k))/Scalar(2); + Scalar E = chebyshev_series(x_k,n_,coefs, epsilon); + + while (abs(E) > epsilon && + abs(delta) > NumTraits::epsilon() ) + { + delta = E / chebyshev_series_deriv(x_k,n_,coefs); + x_k -= delta; + E = chebyshev_series(x_k,n_,coefs, epsilon); + } + + xgk_(2*k) = x_k; + + // Copy adjacent Legendre-zero into the array: + if (2*k+1 < n_) + { + xgk_(2*k+1) = zeros(m_-k); + } + } + } + + /** + * \brief Compute Gauss/Kronrod weights + * \param[in] n_ size of Gauss-Kronrod arrays + * \param[in] m_ Gauss-Legendre degree + * \param[in] xgk_ Gauss/Kronrod abscissae + * \param[out] wg_ Gauss weights + * \param[out] wgk_ Kronrod weights + */ + static void gauss_kronrod_weights(const int& n_, + const int m_, + const Eigen::Array& coefs, + const Eigen::Array& xgk_, + Eigen::Array& wg_, + Eigen::Array& wgk_) + { + Scalar err; + + // Gauss-Legendre weights: + for(int k = 0; k < n_ / 2; ++k) + { + Scalar x = xgk_(2*k + 1); + wg_(k) = (Scalar(-2) / ((m_ + 1) * legendre_deriv(m_, x) * legendre_err(m_+1, x, err))); + } + + // The ratio of leading coefficients of P_n and T_{n+1} is computed + // from the recursive formulae for the respective polynomials. + Scalar F_m = Scalar(2) / Scalar(2*m_ + 1); + + for (int k = 1; k <= m_; ++k) + { + F_m *= (Scalar(2*k) / Scalar(2*k - 1)); + } + + // Gauss-Kronrod weights: + for (int k = 0; k < n_; ++k) + { + Scalar x = xgk_(k); + + if (k % 2 == 0) + { + wgk_(k) = F_m / (legendre_err(m_, x, err) * chebyshev_series_deriv(x,n_,coefs)); + } + else + { + wgk_(k) = (wg_(k/2) + F_m / (legendre_deriv(m_, x) * chebyshev_series(x,n_,coefs, err))); + } + } + + } + + /** + * \brief A method ofr computing Gauss/Kronrod abscissae and weights + * \param[in] m_ Gauss-Legendre degree + * \param[out] xgk_ Gauss/Kronrod abscissae + * \param[out] wgk_ Gauss/Kronrod weights + * \param[out] xk_ Gauss abscissae + * \param[out] wg_ Gauss weights + * + * Note that Gauss abscissae is not calculated here. + */ + static void computeAbscissaeAndWeights(unsigned int m_, + Eigen::Array& xgk_, + Eigen::Array& wgk_, + Eigen::Array& xk_, + Eigen::Array& wg_) + { + const unsigned int n_ = m_ + 1; + + xgk_ = Eigen::Array::Zero(n_); //2*nNodes+1 + wgk_ = Eigen::Array::Zero(n_); //2*nNodes+1 + xk_ = Eigen::Array::Zero(n_/2); //2*nNodes + wg_ = Eigen::Array::Zero(n_/2); //2*nNodes + + // initialise the coefficients to zero + Eigen::Array coefs = Eigen::Array::Zero(n_+1); + Eigen::Array zeros = Eigen::Array::Zero(m_+2); + + legendre_zeros(m_, zeros); + chebyshev_coefs(m_, coefs); + gauss_kronrod_abscissae(n_, m_, zeros,coefs, xgk_); + gauss_kronrod_weights(n_, m_,coefs, xgk_, wg_, wgk_ ); + } + }; + +} // namespace Eigen + +#endif // EIGEN_MONEGATO_H diff --git a/cpp/NumericalIntegration.h b/cpp/NumericalIntegration.h new file mode 100644 --- /dev/null +++ b/cpp/NumericalIntegration.h @@ -0,0 +1,46 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Development work based on work from QUADPACK, Robert Piessens, et al, +// and original work by Dirk Laurie, Walter Gautschi, with support by work +// from John Burkardt. +// +// Multiprecision templating by Pavel Holoborodko and +// code porting, multiprecision templating, and unit tests created by +// Sreekumar Thaithara Balan, Mark Sauder, and Matt Beall 2014. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_NUMERICAL_INTEGRATION_H +#define EIGEN_NUMERICAL_INTEGRATION_H + +namespace Eigen +{ +/** + * \defgroup NumericalIntegration_Module + * \brief This module provides an adaptive quadrature method of numerical integration. + * + * This module provides an adaptive quadrature method of numerical integration of the style + * implemented in the QUADPACK library while offering functionality to calculate nodes/weights + * for Gauss-Kronrod integration, unit tests, and support for multiprecision using mpreal + * precision type. + * + * To use this module, add + * \code + * #include + * \endcode + * at the start of your source file. + */ +} + +//#include +#include +//#include + +#include "ComputeGaussKronrodNodesWeights.h" +#include "GaussKronrodNodesWeights.h" +#include "Integrator.h" + +#endif // EIGEN_NUMERICAL_INTEGRATION_H diff --git a/cpp/Piessens.h b/cpp/Piessens.h new file mode 100644 --- /dev/null +++ b/cpp/Piessens.h @@ -0,0 +1,470 @@ +#ifndef EIGEN_PIESSENS_H +#define EIGEN_PIESSENS_H + +namespace Eigen +{ + + /** + * \ingroup NumericalIntegration_Module + * + * \class Piessens + * + * \brief This class computes Kronrod abscissae & weights for arbitrary precision + * + * \tparam Scalar floating point type + * + * This class is based on the work by R. Piessens, et.al,published in the + * journal Mathematics of Computation, Volume 28, Number 125, January, 1974. + */ + template + class Piessens + { + public: + + /** + * \brief kronrod adds n+1 points to an n-point Gaussian rule. + * + * This function is a C++ implementation of original work by R. Piessens, et.al, + * published in the journal Mathematics of Computation, Volume 28, Number 125, + * January, 1974. Where possible strucutre and variable naming convention + * has been aligned to favor work of the QUADPACK Gaus-Kronrod integration + * routines created by individuals of the same group. + * + * This function calculates the abscissas and weights of the (2n+1)-point + * Gauss-Kronrod quadrature formula which is obtained from the n-point + * Gauss quadrature formula with the optimal addition of (n+1)-points. + * + * The optimally added points are the Kronrod abscissae. The + * abscissas and weights for both the Gauss and Gauss Kronrod rules + * are calculated for integration over the interval (-1, +1). + * + * Because the quadrature formula is symmetric with respect to the origin, + * only the positive abscissas are calculated. Weights corresponding to the + * symetric abscissae are equal. Weights of weightGauss are calculated as well. + * + * Work by Dr. John Burkhardt made note that the code published in Mathematics of + * Computation omitted the definition of the second Chebyshev coefficient (chebCoeff2), + * and Dr. Burkhardt's contributions are reflected here with permission. + * + * The arrays abcsGaussKronrod, weightGaussKronrod and weightGauss contain the + * positive abscissae in decreasing order, and the weights of each abscissa in + * the Gauss-Kronrod and Gauss rules, respectively. + * + * Ported to C++/Eigen and templated for multiprecision by Mark Sauder, + * Sreekumar Thaithara Balan, Matt Beall, and R. Jeff Jenkins - September 2014. + * + * Input Parameters: + * \param[in] n, the order of the Gauss rule. + * + * Return Parameters: + * \param[in,out] abscGaussKronrod[n+1] The Gauss-Kronrod abscissae. + * \param[in,out] weightGaussKronrod[n+1] The weights for the Gauss-Kronrod rule. + * \param[in,out] weightGauss[n+1] The weights for the Gauss rule. + */ + static void kronrod(unsigned int nNodes, + Eigen::Array& abscGaussKronrod, + Eigen::Array& weightGaussKronrod, + Eigen::Array& weightGauss) + { + unsigned int arraySize = nNodes + 1; + abscGaussKronrod = Eigen::Array::Zero(arraySize); + weightGaussKronrod = Eigen::Array::Zero(arraySize); + weightGauss = Eigen::Array::Zero(arraySize / 2); + + Scalar aN(0.0); + Scalar d(2.0); + + for (size_t i = 0; i < nNodes; ++i) + { + aN += Scalar(1.0); + d *= aN / (aN + Scalar(0.5)); + } + + unsigned int m = (nNodes + 1) / 2; + bool even = (nNodes == 2 * m); + + // aK is an index variable to account for calculating only the positive abscissae + Scalar aK = aN; + + // Calculation of the Chebyshev coefficients of the orthogonal polynomial. + Eigen::Array tau(m); + tau(0) = (aN + Scalar(2.0)) / (Scalar(2) * aN + Scalar(3.0)); + + Eigen::Array betaCoeffs(m + 1); + betaCoeffs(m - 1) = tau(0) - Scalar(1.0); + + for (size_t k = 1; k < m; ++k) + { + // This step accounts for both positive and negative abscissae + aK += Scalar(2.0); + + tau(k) = ((aK - Scalar(1.0)) * aK - aN * (aN + Scalar(1.0) )) * (aK + Scalar(2.0) ) * tau(k-1) / + (aK * ((aK + Scalar(3.0) ) * (aK + Scalar(2.0) ) - aN * (aN + Scalar(1.0) ))); + + betaCoeffs(m-k-1) = tau(k); + + for (size_t i = 1; i <= k; ++i) + { + betaCoeffs(m - k - 1) = betaCoeffs(m - k + i - 1) * tau(i - 1) + betaCoeffs(m - k - 1); + } + } + + betaCoeffs(m) = Scalar(1.); + + // Calculation of approximate values for the abscissae as inital values + // for the Newton-Raphson iterative solution. These values are derived + // from Pythagorean identities to the original code to more closely follow + // the mathematics of the 1974 CoM paper. + + // @TODO The usage of constant Pi with fixed precision needs to be changed to the following for multiprecision + //RealScalar pi = NumTraits::Pi(); + + Scalar s1 = sin((M_PI / Scalar(2) ) / (Scalar(2.) * aN + Scalar(1.0) )); + Scalar c1 = cos((M_PI / Scalar(2) ) / (Scalar(2.) * aN + Scalar(1.0) )); + + Scalar s2 = sin((M_PI) / (Scalar(2.) * aN + Scalar(1.0) )); + Scalar c2 = cos((M_PI) / (Scalar(2.) * aN + Scalar(1.0) )); + + // Coefficient for Gauss and Kronrod abscissae and weights + Scalar chebCoeff1 = Scalar(1.0) - Scalar(1.0) / (Scalar(8.0) * aN * aN) + Scalar(1.0) / (Scalar(8.0) * aN * aN * aN); + Scalar chebCoeff2 = Scalar(2.0) / (Scalar(2. * nNodes + 1)); + + for (size_t i = 1; i <= nNodes; ++i) + { + chebCoeff2 = Scalar(4.0) * chebCoeff2 * i / (nNodes + i); + } + + Scalar abscK = chebCoeff1 * c1; + Scalar temp(0.); + + // Calculation of the K-th (Kronrod) abscissa and the corresponding weight. + for (size_t k = 0; k < nNodes; ++k) + { + abscWeightKronrod(nNodes, m, even, chebCoeff2, betaCoeffs, abscK, weightGaussKronrod(k)); + abscGaussKronrod(k) = abscK; + ++k; + + temp = c1; + c1 = temp * c2 - s1 * s2; + s1 = temp * s2 + s1 * c2; + abscK = chebCoeff1 * c1; + + // Calculation of the k+1 (Gauss) abscissa and the corresponding weights. + abscWeightGauss(nNodes, m, even, chebCoeff2, betaCoeffs, abscK, weightGaussKronrod(k), + weightGauss(k/2)); + abscGaussKronrod(k) = abscK; + + temp = c1; + c1 = temp * c2 - s1 * s2; + s1 = temp * s2 + s1 * c2; + abscK = chebCoeff1 * c1; + } + + // Add a Kronrod abscissa at the origin if n is even. + if (even) + { + abscWeightKronrod(nNodes, m, even, chebCoeff2, betaCoeffs, abscK, weightGaussKronrod(nNodes)); + } + + // Set the abscissa value at the origin to zero and exit the function. + abscGaussKronrod(nNodes) = Scalar(0.0); + return; + } + + + /** + * \brief abscWeightKronrod calculates a Kronrod abscissa and weight. + * + * Input Parameters: + * \param[in] betaCoeffs[m+1] The Chebyshev coefficients. + * \param[in] chebCoeff A value needed to compute weights. + * \param[in] even A boolean variable that is TRUE if n is even. + * \param[in] n The order of the Gauss rule. + * \param[in] m The value of ( n + 1 ) / 2. + * + * Input/output: + * \param[in,out] abscGaussKronrod An estimate for the abscissa on input and the computed abscissa on output. + * \param[in,out] weightGaussKronrod The Gauss-Kronrod weight. + */ + static void abscWeightKronrod(unsigned int nNodes, + unsigned int m, + bool even, + Scalar chebCoeff, + Eigen::Array betaCoeffs, + Scalar& abscGaussKronrod, + Scalar& weightGaussKronrod) + { + Scalar ai; + + Scalar b0(0.); + Scalar b1(0.); + Scalar b2(0.); + + Scalar d0(0.); + Scalar d1(0.); + Scalar d2(0.); + + Scalar delta(1.); + Scalar dif(0.); + + Scalar f(0.); + Scalar fd(0.); + + Scalar yy(0.); + + int i = 0; + + size_t iter = 0; + size_t iterationLimit = 50; + + // Iterative process for the computation of a Kronrod abscissa. + while (abs(delta) > machineEpsilon()) + { + ++iter; + + b1 = Scalar(0.0); + b2 = betaCoeffs(m); + + yy = Scalar(4.0) * abscGaussKronrod * abscGaussKronrod - Scalar(2.0); + d1 = Scalar(0.0); + + if (even) + { + ai = Scalar(m + m + 1); + d2 = ai * betaCoeffs(m); + dif = Scalar(2.0); + } + else + { + ai = Scalar(m + 1); + d2 = Scalar(0.0); + dif = Scalar(1.0); + } + + for (size_t k = 0; k < m; ++k) + { + ai -= dif; + i = m - (int)k - 1; + b0 = b1; + b1 = b2; + d0 = d1; + d1 = d2; + b2 = yy * b1 - b0 + betaCoeffs(i); + + if (!even) + { + i += 1; + } + + d2 = yy * d1 - d0 + ai * betaCoeffs(i); + } + + if (even) + { + f = abscGaussKronrod * (b2 - b1); + fd = d2 + d1; + } + else + { + f = Scalar(0.5) * (b2 - b0); + fd = Scalar(4.0) * abscGaussKronrod * d2; + } + + // Newton correction. + delta = f / fd; + abscGaussKronrod -= delta; + + if (abscGaussKronrod == Scalar(0.)) + { + break; + } + + // Identify non-convergence of the iterative solver after 50 iterations + if (iter > iterationLimit) + { + std::cout << "@abscWeightKronrod Newton-Raphson iterative abscissae solver failed."< betaCoeffs, + Scalar& abscGaussKronrod, + Scalar& weightGaussKronrod, + Scalar& weightGauss) + { + Scalar ai(0.); + Scalar delta(1.); + + Scalar p0(0.); + Scalar p1(0.); + Scalar p2(0.); + + Scalar pd0(0.); + Scalar pd1(0.); + Scalar pd2(0.); + + Scalar yy(0.); + + size_t iter = 0; + size_t iterationLimit = 50; + + // Iterative process for the computation of a Gaussian abscissa. + while (abs(delta) > machineEpsilon()) + { + ++iter; + p0 = Scalar(1.); + p1 = abscGaussKronrod; + pd0 = Scalar(0.); + pd1 = Scalar(1.); + + // If nNodes <= 1, initialize p2 and pd2 to avoid problems calculating delta. + if (nNodes <= 1) + { + if (machineEpsilon() < abs(abscGaussKronrod)) + { + p2 = (Scalar(3.0) * (abscGaussKronrod) * (abscGaussKronrod) - Scalar(1.0)) / Scalar(2.0); + pd2 = Scalar(3.0) * (abscGaussKronrod); + } + else + { + p2 = Scalar(3.0) * (abscGaussKronrod); + pd2 = Scalar(3.0); + } + } + + ai = Scalar(0.0); + + for (size_t k = 0; k < nNodes - 1; ++k) + { + ai = ai + Scalar(1.0); + p2 = ((ai + ai + Scalar(1.0) ) * abscGaussKronrod * p1 - ai * p0) / (ai + Scalar(1.0) ); + pd2 = ((ai + ai + Scalar(1.0) ) * (p1 + abscGaussKronrod * pd1) - ai * pd0) / (ai + Scalar(1.0) ); + p0 = p1; + p1 = p2; + pd0 = pd1; + pd1 = pd2; + } + + // Newton correction. + delta = p2 / pd2; + + if (abscGaussKronrod == Scalar(0.0)) + { + abscGaussKronrod -= delta; + break; + } + else + { + abscGaussKronrod -= delta; + } + + // Identify non-convergence of the iterative solver after iteration limit. + if (iter > iterationLimit) + { + std::cout << "@abscWeightGauss Newton-Raphson iterative abscissae solver failed."<::epsilon(); + } + + /** + * \brief A method ofr computing Gauss/Kronrod abscissae and weights + * \param nNodes Gauss-Legendre degree + * \param abscGaussKronrod Gauss/Kronrod abscissae + * \param weightGaussKronrod Gauss/Kronrod weights + * \param abscGauss Gauss abscissae + * \param weightGauss Gauss weights + */ + static void computeAbscissaeAndWeights(unsigned int nNodes, + Eigen::Array &abscGaussKronrod, + Eigen::Array &weightGaussKronrod, + Eigen::Array abscGauss, + Eigen::Array &weightGauss) + { + Piessens::kronrod(nNodes,abscGaussKronrod,weightGaussKronrod,weightGauss); + abscGauss = Eigen::Array::Zero(abscGauss.rows()); + } + }; + +}//namespace Eigen +#endif //EIGEN_PIESSENS_H diff --git a/cpp/integration.cpp b/cpp/integration.cpp new file mode 100644 --- /dev/null +++ b/cpp/integration.cpp @@ -0,0 +1,46 @@ +#define _USE_MATH_DEFINES + +#include +#include +#include +#include +#include +#include +#include + +extern "C" { + +class Integrand { + private: + std::function f; + + public: + Integrand(std::function& f_) : f(f_) {} + double operator()(const double x) const { return f(x); } +}; + +double integration(double f(double), + double lower, + double upper, + double relError, + int subdiv, + double* errorEstimate, + int* errorCode) { + // Define the integrand. + std::function f_ = [&](double x) { return f(x); }; + Integrand integrand(f_); + // Define the integrator. + Eigen::Integrator integrator(subdiv); + // Define a quadrature rule. + Eigen::Integrator::QuadratureRule rule = + Eigen::Integrator::GaussKronrod201; + // Define the desired absolute error. + double absError = 0.0; + // Integrate. + double result = integrator.quadratureAdaptive(integrand, lower, upper, + absError, relError, rule); + *errorEstimate = integrator.estimatedError(); + *errorCode = integrator.errorCode(); + return result; +} +} diff --git a/numerical-integration.cabal b/numerical-integration.cabal new file mode 100644 --- /dev/null +++ b/numerical-integration.cabal @@ -0,0 +1,50 @@ +cabal-version: 2.2 +name: numerical-integration +version: 0.1.0.0 +synopsis: Numerical integration. +description: One-dimensional numerical integration using the 'NumericalIntegration' C++ library. +homepage: https://github.com/stla/numerical-integration#readme +license: BSD-3-Clause +license-file: LICENSE +author: Stéphane Laurent +maintainer: laurent_step@outlook.fr +copyright: 2023 Stéphane Laurent +category: Numerical +build-type: Simple +extra-source-files: README.md + CHANGELOG.md + LICENSES/LICENSE_Eigen + LICENSES/LICENSE_NumericalIntegration + +library + hs-source-dirs: src + exposed-modules: Numerical.Integration + build-depends: base >= 4.7 && < 5 + other-extensions: ForeignFunctionInterface + include-dirs: cpp + C-sources: cpp/integration.cpp + install-includes: cpp/ComputeGaussKronrodNodesWeights.h + , cpp/Eigen/Eigenvalues + , cpp/GaussKronrodNodesWeights.h + , cpp/Integrator.h + , cpp/LaurieGautschi.h + , cpp/Monegato.h + , cpp/NumericalIntegration.h + , cpp/Piessens.h + extra-libraries: stdc++ + default-language: Haskell2010 + ghc-options: -Wall + -Wcompat + -Widentities + -Wincomplete-record-updates + -Wincomplete-uni-patterns + -Wmissing-export-lists + -Wmissing-home-modules + -Wpartial-fields + -Wredundant-constraints + -optcxx-std=c++11 + cxx-options: -fPIC -std=c++11 + +source-repository head + type: git + location: https://github.com/stla/numerical-integration diff --git a/src/Numerical/Integration.hs b/src/Numerical/Integration.hs new file mode 100644 --- /dev/null +++ b/src/Numerical/Integration.hs @@ -0,0 +1,35 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +module Numerical.Integration + (integration) + where +import Foreign.Marshal.Alloc (free, mallocBytes) +import Foreign.Ptr (FunPtr, Ptr, freeHaskellFunPtr) +import Foreign.Storable (peek, sizeOf) + +foreign import ccall safe "wrapper" funPtr + :: (Double -> Double) -> IO(FunPtr (Double -> Double)) + +foreign import ccall safe "integration" c_integration + :: FunPtr (Double -> Double) -> Double -> Double -> Double -> Int + -> Ptr Double -> Ptr Int -> IO Double + +-- | Numerical integration. +integration :: (Double -> Double) -- ^ integrand + -> Double -- ^ lower bound + -> Double -- ^ upper bound + -> Double -- ^ desired relative error + -> Int -- ^ number of subdivisions + -> IO (Double, Double, Int) -- ^ value, error estimate, error code +integration f lower upper relError subdiv = do + errorEstimatePtr <- mallocBytes (sizeOf (0 :: Double)) + errorCodePtr <- mallocBytes (sizeOf (0 :: Int)) + fPtr <- funPtr f + result <- + c_integration fPtr lower upper relError subdiv errorEstimatePtr errorCodePtr + errorEstimate <- peek errorEstimatePtr + errorCode <- peek errorCodePtr + let out = (result, errorEstimate, errorCode) + free errorEstimatePtr + free errorCodePtr + freeHaskellFunPtr fPtr + return out