numerical-integration (empty) → 0.1.0.0

raw patch · 17 files changed

+4535/−0 lines, 17 filesdep +basesetup-changed

Dependencies added: base

Files

• CHANGELOG.md +5/−0
• LICENSE +30/−0
• LICENSES/LICENSE_Eigen +26/−0
• LICENSES/LICENSE_NumericalIntegration +362/−0
• README.md +13/−0
• Setup.hs +2/−0
• cpp/ComputeGaussKronrodNodesWeights.h +14/−0
• cpp/Eigen/Eigenvalues +60/−0
• cpp/GaussKronrodNodesWeights.h +1770/−0
• cpp/Integrator.h +745/−0
• cpp/LaurieGautschi.h +481/−0
• cpp/Monegato.h +380/−0
• cpp/NumericalIntegration.h +46/−0
• cpp/Piessens.h +470/−0
• cpp/integration.cpp +46/−0
• numerical-integration.cabal +50/−0
• src/Numerical/Integration.hs +35/−0

@@ -0,0 +1,5 @@+# Changelog for `numerical-integration`++## 0.1.0.0 - 2023-08-30++First release.

@@ -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.

@@ -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.+*/

@@ -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.

@@ -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+```

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -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

@@ -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 <Eigen/Eigenvalues>+  * \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

@@ -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 <typename Scalar>+    class QuadratureKronrod+    {+    public:+        static Array<Scalar, 8, 1> abscissaeGaussKronrod15;+        static Array<Scalar, 8, 1> weightsGaussKronrod15;+        static Array<Scalar, 4, 1> abscissaeGauss15;+        static Array<Scalar, 4, 1> weightsGauss15;++        static Array<Scalar, 11, 1> abscissaeGaussKronrod21;+        static Array<Scalar, 11, 1> weightsGaussKronrod21;+        static Array<Scalar, 5, 1> abscissaeGauss21;+        static Array<Scalar, 5, 1> weightsGauss21;++        static Array<Scalar, 16, 1> abscissaeGaussKronrod31;+        static Array<Scalar, 16, 1> weightsGaussKronrod31;+        static Array<Scalar, 8, 1> abscissaeGauss31;+        static Array<Scalar, 8, 1> weightsGauss31;++        static Array<Scalar, 21, 1> abscissaeGaussKronrod41;+        static Array<Scalar, 21, 1> weightsGaussKronrod41;+        static Array<Scalar, 10, 1> abscissaeGauss41;+        static Array<Scalar, 10, 1> weightsGauss41;++        static Array<Scalar, 26, 1> abscissaeGaussKronrod51;+        static Array<Scalar, 26, 1> weightsGaussKronrod51;+        static Array<Scalar, 13, 1> abscissaeGauss51;+        static Array<Scalar, 13, 1> weightsGauss51;++        static Array<Scalar, 31, 1> abscissaeGaussKronrod61;+        static Array<Scalar, 31, 1> weightsGaussKronrod61;+        static Array<Scalar, 15, 1> abscissaeGauss61;+        static Array<Scalar, 15, 1> weightsGauss61;++        static Array<Scalar, 36, 1> abscissaeGaussKronrod71;+        static Array<Scalar, 36, 1> weightsGaussKronrod71;+        static Array<Scalar, 18, 1> abscissaeGauss71;+        static Array<Scalar, 18, 1> weightsGauss71;++        static Array<Scalar, 41, 1> abscissaeGaussKronrod81;+        static Array<Scalar, 41, 1> weightsGaussKronrod81;+        static Array<Scalar, 20, 1> abscissaeGauss81;+        static Array<Scalar, 20, 1> weightsGauss81;++        static Array<Scalar, 46, 1> abscissaeGaussKronrod91;+        static Array<Scalar, 46, 1> weightsGaussKronrod91;+        static Array<Scalar, 23, 1> abscissaeGauss91;+        static Array<Scalar, 23, 1> weightsGauss91;++        static Array<Scalar, 51, 1> abscissaeGaussKronrod101;+        static Array<Scalar, 51, 1> weightsGaussKronrod101;+        static Array<Scalar, 25, 1> abscissaeGauss101;+        static Array<Scalar, 25, 1> weightsGauss101;++        static Array<Scalar, 61, 1> abscissaeGaussKronrod121;+        static Array<Scalar, 61, 1> weightsGaussKronrod121;+        static Array<Scalar, 30, 1> abscissaeGauss121;+        static Array<Scalar, 30, 1> weightsGauss121;++        static Array<Scalar, 101, 1> abscissaeGaussKronrod201;+        static Array<Scalar, 101, 1> weightsGaussKronrod201;+        static Array<Scalar, 50, 1> abscissaeGauss201;+        static Array<Scalar, 50, 1> 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 <Index N>+        static void computeForRule(Array<Scalar, N+1, 1>& kronrodAbscissae, Array<Scalar, N+1, 1>& kronrodWeights,+                                   Array<Scalar, (N+1)/2, 1>& gaussAbscissae, Array<Scalar, (N+1)/2, 1>& gaussWeights)+        {+            Eigen::Array<Scalar, Eigen::Dynamic, 1> xGK;+            Eigen::Array<Scalar, Eigen::Dynamic, 1> wGK;+            Eigen::Array<Scalar, Eigen::Dynamic, 1> xG;+            Eigen::Array<Scalar, Eigen::Dynamic, 1> wG;++            Eigen::LaurieGautschi<Scalar>::computeAbscissaeAndWeights(N,xGK,wGK,xG,wG);+            //Eigen::Piessens<Scalar>::computeAbscissaeAndWeights(N,xGK,wGK,xG,wG);+            //Eigen::Monegato<Scalar>::computeAbscissaeAndWeights(N,xGK,wGK,xG,wG);++            for(size_t i=0; i<N+1; ++i)+            {+                kronrodAbscissae(i) = xGK(i);+                kronrodWeights(i) =  wGK(i);+            }++            for(size_t i=0; i<(N+1)/2; ++i)+            {+                gaussAbscissae(i) = xG(i);+                gaussWeights(i) = wG(i);+            }+        }+    };+++// Nodes and Weights - Rule 15+template <typename Scalar>+Array<Scalar, 8, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod15 =+    (Array<Scalar, 8, 1>() <<+        0.99145537112081263920685469752632851664204433837033,+        0.94910791234275852452618968404785126240077093767062,+        0.86486442335976907278971278864092620121097230707409,+        0.74153118559939443986386477328078840707414764714139,+        0.58608723546769113029414483825872959843678075060436,+        0.40584515137739716690660641207696146334738201409937,+        0.20778495500789846760068940377324491347978440714517,+        0.00000000000000000000000000000000000000000000000000+    ).finished();++template <typename Scalar>+Array<Scalar, 8, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod15 =+    (Array<Scalar, 8, 1>() <<+        0.02293532201052922496373200805896959199356081127575,+        0.06309209262997855329070066318920428666507115721155,+        0.10479001032225018383987632254151801744375665421383,+        0.14065325971552591874518959051023792039988975724800,+        0.16900472663926790282658342659855028410624490030294,+        0.19035057806478540991325640242101368282607807545536,+        0.20443294007529889241416199923464908471651760418072,+        0.20948214108472782801299917489171426369776208022370+    ).finished();++template <typename Scalar>+Array<Scalar, 4, 1> QuadratureKronrod<Scalar>::abscissaeGauss15 =+    (Array<Scalar, 4, 1>() <<+        0.94910791234275852452618968404785126240077093767062,+        0.74153118559939443986386477328078840707414764714139,+        0.40584515137739716690660641207696146334738201409937,+        0.00000000000000000000000000000000000000000000000000+    ).finished();++template <typename Scalar>+Array<Scalar, 4, 1> QuadratureKronrod<Scalar>::weightsGauss15 =+    (Array<Scalar, 4, 1>() <<+        0.12948496616886969327061143267908201832858740225995,+        0.27970539148927666790146777142377958248692506522660,+        0.38183005050511894495036977548897513387836508353386,+        0.41795918367346938775510204081632653061224489795918+    ).finished();++// Nodes and Weights - Rule 21+template <typename Scalar>+Array<Scalar, 11, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod21 =+    (Array<Scalar, 11, 1>() <<+        0.99565716302580808073552728068900284792126058721948,+        0.97390652851717172007796401208445205342826994669238,+        0.93015749135570822600120718005950834622516790998194,+        0.86506336668898451073209668842349304852754301496533,+        0.78081772658641689706371757834504237716340752029816,+        0.67940956829902440623432736511487357576929471183481,+        0.56275713466860468333900009927269414084301388194197,+        0.43339539412924719079926594316578416220007183765625,+        0.29439286270146019813112660310386556616268662515696,+        0.14887433898163121088482600112971998461756485942069,+        0.00000000000000000000000000000000000000000000000000+    ).finished();++template <typename Scalar>+Array<Scalar, 11, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod21 =+    (Array<Scalar, 11, 1>() <<+        0.01169463886737187427806439606219204839621733248193,+        0.03255816230796472747881897245938976061738893984566,+        0.05475589657435199603138130024458017637372111405833,+        0.07503967481091995276704314091619000939521938200091,+        0.09312545458369760553506546508336634439001882888076,+        0.10938715880229764189921059032580496027181329983435,+        0.12349197626206585107795810983107415951230034952865,+        0.13470921731147332592805400177170683276099191300856,+        0.14277593857706008079709427313871706088597905653191,+        0.14773910490133849137484151597206804552373162548521,+        0.14944555400291690566493646838982120374523631668747+    ).finished();++template <typename Scalar>+Array<Scalar, 5, 1> QuadratureKronrod<Scalar>::abscissaeGauss21 =+    (Array<Scalar, 5, 1>() <<+        0.97390652851717172007796401208445205342826994669238,+        0.86506336668898451073209668842349304852754301496533,+        0.67940956829902440623432736511487357576929471183481,+        0.43339539412924719079926594316578416220007183765625,+        0.14887433898163121088482600112971998461756485942069+    ).finished();++template <typename Scalar>+Array<Scalar, 5, 1> QuadratureKronrod<Scalar>::weightsGauss21 =+    (Array<Scalar, 5, 1>() <<+        0.06667134430868813759356880989333179285786483432016,+        0.14945134915058059314577633965769733240255663966943,+        0.21908636251598204399553493422816319245877187052268,+        0.26926671930999635509122692156946935285975993846088,+        0.29552422471475287017389299465133832942104671702685+    ).finished();++// Nodes and Weights - Rule 31+template <typename Scalar>+Array<Scalar, 16, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod31 =+    (Array<Scalar, 16, 1>() <<+        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 <typename Scalar>+Array<Scalar, 16, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod31 =+    (Array<Scalar, 16, 1>() <<+        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 <typename Scalar>+Array<Scalar, 8, 1> QuadratureKronrod<Scalar>::abscissaeGauss31 =+    (Array<Scalar, 8, 1>() <<+        0.98799251802048542848956571858661258114697281712376,+        0.93727339240070590430775894771020947124399627351530,+        0.84820658341042721620064832077421685136625617473699,+        0.72441773136017004741618605461393800963089929458410,+        0.57097217260853884753722673725391064123838639628275,+        0.39415134707756336989720737098104546836275277615870,+        0.20119409399743452230062830339459620781283645446264,+        0.00000000000000000000000000000000000000000000000000+    ).finished();++template <typename Scalar>+Array<Scalar, 8, 1> QuadratureKronrod<Scalar>::weightsGauss31 =+    (Array<Scalar, 8, 1>() <<+        0.03075324199611726835462839357720441772174814483343,+        0.07036604748810812470926741645066733846670803275433,+        0.10715922046717193501186954668586930341554371575810,+        0.13957067792615431444780479451102832252085027531551,+        0.16626920581699393355320086048120881113090018009841,+        0.18616100001556221102680056186642282450622601227793,+        0.19843148532711157645611832644383932481869255995754,+        0.20257824192556127288062019996751931483866215800948+    ).finished();++// Nodes and Weights - Rule 41+template <typename Scalar>+Array<Scalar, 21, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod41 =+    (Array<Scalar, 21, 1>() <<+        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 <typename Scalar>+Array<Scalar, 21, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod41 =+    (Array<Scalar, 21, 1>() <<+        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 <typename Scalar>+Array<Scalar, 10, 1> QuadratureKronrod<Scalar>::abscissaeGauss41 =+    (Array<Scalar, 10, 1>() <<+        0.99312859918509492478612238847132027822264713090166,+        0.96397192727791379126766613119727722191206032780619,+        0.91223442825132590586775244120329811304918479742369,+        0.83911697182221882339452906170152068532962936506564,+        0.74633190646015079261430507035564159031073067956918,+        0.63605368072651502545283669622628593674338911679937,+        0.51086700195082709800436405095525099842549132920243,+        0.37370608871541956067254817702492723739574632170568,+        0.22778585114164507808049619536857462474308893768293,+        0.07652652113349733375464040939883821100479626681350+    ).finished();++template <typename Scalar>+Array<Scalar, 10, 1> QuadratureKronrod<Scalar>::weightsGauss41 =+    (Array<Scalar, 10, 1>() <<+        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 <typename Scalar>+Array<Scalar, 26, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod51 =+    (Array<Scalar, 26, 1>() <<+        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 <typename Scalar>+Array<Scalar, 26, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod51 =+    (Array<Scalar, 26, 1>() <<+        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 <typename Scalar>+Array<Scalar, 13, 1> QuadratureKronrod<Scalar>::abscissaeGauss51 =+    (Array<Scalar, 13, 1>() <<+        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 <typename Scalar>+Array<Scalar, 13, 1> QuadratureKronrod<Scalar>::weightsGauss51 =+    (Array<Scalar, 13, 1>() <<+        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 <typename Scalar>+Array<Scalar, 31, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod61 =+    (Array<Scalar, 31, 1>() <<+        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 <typename Scalar>+Array<Scalar, 31, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod61 =+    (Array<Scalar, 31, 1>() <<+        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 <typename Scalar>+Array<Scalar, 15, 1> QuadratureKronrod<Scalar>::abscissaeGauss61 =+    (Array<Scalar, 15, 1>() <<+        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 <typename Scalar>+Array<Scalar, 15, 1> QuadratureKronrod<Scalar>::weightsGauss61 =+    (Array<Scalar, 15, 1>() <<+        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 <typename Scalar>+Array<Scalar, 36, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod71 =+    (Array<Scalar, 36, 1>() <<+        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 <typename Scalar>+Array<Scalar, 36, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod71 =+    (Array<Scalar, 36, 1>() <<+        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 <typename Scalar>+Array<Scalar, 18, 1> QuadratureKronrod<Scalar>::abscissaeGauss71 =+    (Array<Scalar, 18, 1>() <<+        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 <typename Scalar>+Array<Scalar, 18, 1> QuadratureKronrod<Scalar>::weightsGauss71 =+    (Array<Scalar, 18, 1>() <<+        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 <typename Scalar>+Array<Scalar, 41, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod81 =+    (Array<Scalar, 41, 1>() <<+        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 <typename Scalar>+Array<Scalar, 41, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod81 =+    (Array<Scalar, 41, 1>() <<+        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 <typename Scalar>+Array<Scalar, 20, 1> QuadratureKronrod<Scalar>::abscissaeGauss81 =+    (Array<Scalar, 20, 1>() <<+        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 <typename Scalar>+Array<Scalar, 20, 1> QuadratureKronrod<Scalar>::weightsGauss81 =+    (Array<Scalar, 20, 1>() <<+        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 <typename Scalar>+Array<Scalar, 46, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod91 =+    (Array<Scalar, 46, 1>() <<+        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 <typename Scalar>+Array<Scalar, 46, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod91 =+    (Array<Scalar, 46, 1>() <<+        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 <typename Scalar>+Array<Scalar, 23, 1> QuadratureKronrod<Scalar>::abscissaeGauss91 =+    (Array<Scalar, 23, 1>() <<+        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 <typename Scalar>+Array<Scalar, 23, 1> QuadratureKronrod<Scalar>::weightsGauss91 =+    (Array<Scalar, 23, 1>() <<+        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 <typename Scalar>+Array<Scalar, 51, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod101 =+    (Array<Scalar, 51, 1>() <<+        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 <typename Scalar>+Array<Scalar, 51, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod101 =+    (Array<Scalar, 51, 1>() <<+        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 <typename Scalar>+Array<Scalar, 25, 1> QuadratureKronrod<Scalar>::abscissaeGauss101 =+    (Array<Scalar, 25, 1>() <<+        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 <typename Scalar>+Array<Scalar, 25, 1> QuadratureKronrod<Scalar>::weightsGauss101 =+    (Array<Scalar, 25, 1>() <<+        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 <typename Scalar>+Array<Scalar, 61, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod121 =+    (Array<Scalar, 61, 1>() <<+        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 <typename Scalar>+Array<Scalar, 61, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod121 =+    (Array<Scalar, 61, 1>() <<+        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 <typename Scalar>+Array<Scalar, 30, 1> QuadratureKronrod<Scalar>::abscissaeGauss121 =+    (Array<Scalar, 30, 1>() <<+        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 <typename Scalar>+Array<Scalar, 30, 1> QuadratureKronrod<Scalar>::weightsGauss121 =+    (Array<Scalar, 30, 1>() <<+        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 <typename Scalar>+Array<Scalar, 101, 1> QuadratureKronrod<Scalar>::abscissaeGaussKronrod201 =+    (Array<Scalar, 101, 1>() <<+        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 <typename Scalar>+Array<Scalar, 101, 1> QuadratureKronrod<Scalar>::weightsGaussKronrod201 =+    (Array<Scalar, 101, 1>() <<+        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 <typename Scalar>+Array<Scalar, 50, 1> QuadratureKronrod<Scalar>::abscissaeGauss201 =+    (Array<Scalar, 50, 1>() <<+        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 <typename Scalar>+Array<Scalar, 50, 1> QuadratureKronrod<Scalar>::weightsGauss201 =+    (Array<Scalar, 50, 1>() <<+        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

@@ -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 <typename Scalar>++    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 <typename FunctionType>+        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<Scalar>::QuadratureRule(1))+        {+            using std::abs;+            using std::max;++            if ((desiredAbsoluteError <= Scalar(0.) &&+                desiredRelativeError < NumTraits<Scalar>::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<Scalar>::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<Scalar>::epsilon() * Scalar(100.) + Scalar(1.))+                             * (abs(lower2) + (std::numeric_limits<Scalar>::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 <typename FunctionType>+        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<Scalar>::abscissaeGaussKronrod15,+                QuadratureKronrod<Scalar>::weightsGaussKronrod15, QuadratureKronrod<Scalar>::weightsGauss15,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod21:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod21,+                QuadratureKronrod<Scalar>::weightsGaussKronrod21, QuadratureKronrod<Scalar>::weightsGauss21,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod31:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod31,+                QuadratureKronrod<Scalar>::weightsGaussKronrod31, QuadratureKronrod<Scalar>::weightsGauss31,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod41:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod41,+                QuadratureKronrod<Scalar>::weightsGaussKronrod41, QuadratureKronrod<Scalar>::weightsGauss41,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod51:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod51,+                QuadratureKronrod<Scalar>::weightsGaussKronrod51, QuadratureKronrod<Scalar>::weightsGauss51,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod61:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod61,+                QuadratureKronrod<Scalar>::weightsGaussKronrod61, QuadratureKronrod<Scalar>::weightsGauss61,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod71:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod71,+                QuadratureKronrod<Scalar>::weightsGaussKronrod71, QuadratureKronrod<Scalar>::weightsGauss71,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod81:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod81,+                QuadratureKronrod<Scalar>::weightsGaussKronrod81, QuadratureKronrod<Scalar>::weightsGauss81,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod91:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod91,+                QuadratureKronrod<Scalar>::weightsGaussKronrod91, QuadratureKronrod<Scalar>::weightsGauss91,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod101:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod101,+                QuadratureKronrod<Scalar>::weightsGaussKronrod101, QuadratureKronrod<Scalar>::weightsGauss101,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod121:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod121,+                QuadratureKronrod<Scalar>::weightsGaussKronrod121, QuadratureKronrod<Scalar>::weightsGauss121,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            case GaussKronrod201:+              return quadratureKronrodHelper(+                QuadratureKronrod<Scalar>::abscissaeGaussKronrod201,+                QuadratureKronrod<Scalar>::weightsGaussKronrod201, QuadratureKronrod<Scalar>::weightsGauss201,+                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);++            default:+              return Scalar(0.);+            }+        }++        template <typename FunctionType, int numKronrodRows, int numGaussRows, int alignment>+        Scalar quadratureKronrodHelper(Array<Scalar, numKronrodRows, 1, alignment, numKronrodRows, 1> abscissaeGaussKronrod,+                                       Array<Scalar, numKronrodRows, 1, alignment, numKronrodRows, 1> weightsGaussKronrod,+                                       Array<Scalar, numGaussRows, 1, alignment, numGaussRows, 1> 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<Scalar, numKronrodRows - 1, 1> f1Array;+            Array<Scalar, numKronrodRows - 1, 1> 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<Scalar>::min)() / (NumTraits<Scalar>::epsilon() * Scalar(50.)))+            {+                estimatedError = max(NumTraits<Scalar>::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<Index, Dynamic, 1> 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<Scalar, Dynamic, 1> 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<Scalar, Dynamic, 1> 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<Scalar, Dynamic, 1> 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<Scalar, Dynamic, 1> 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

@@ -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 <typename Scalar>++    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<Scalar, Eigen::Dynamic, 1>& alphaOut,+                             Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<Scalar, Eigen::Dynamic, 1>& alphaOut,+                                Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<Scalar, Eigen::Dynamic, 1>& alphaIn,+                              const Eigen::Array<Scalar, Eigen::Dynamic, 1>& betaIn,+                              Eigen::Array<Scalar, Eigen::Dynamic, 1>& alpha,+                              Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<Scalar, Eigen::Dynamic, 1>::Zero(2*N+1);+            beta = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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<Scalar, Eigen::Dynamic, 1> sigma = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero((Index)floor(N/2) + 2);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> tempVector = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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<Scalar, Eigen::Dynamic, 1> 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<Scalar, Eigen::Dynamic, 1> 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<Scalar, Eigen::Dynamic, 1>& alpha,+                            const Eigen::Array<Scalar, Eigen::Dynamic, 1>& beta,+                            Eigen::Array<Scalar, Eigen::Dynamic, 1>& nodes,+                            Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<Scalar, Eigen::Dynamic, 1> alpha0 = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N+1);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> beta0 = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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<Scalar,Eigen::Dynamic, Eigen::Dynamic> J = Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic>::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<Scalar,Eigen::Dynamic, Eigen::Dynamic> > es(J);++            //TODO : make use the eigen assert facilities+            assert(es.info() == Eigen::Success);++            nodes = es.eigenvalues();+            Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> 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<Scalar, Eigen::Dynamic, 1>& alpha,+                          const Eigen::Array<Scalar, Eigen::Dynamic, 1>& beta,+                          Eigen::Array<Scalar, Eigen::Dynamic, 1>& nodes,+                          Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<Scalar,Eigen::Dynamic, Eigen::Dynamic> J = Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic>::Zero(N,N);++            J(0,0) = alpha(0);++            for (Index n=1; n<N; ++n)+            {+                J(n,n) = alpha(n);+                J(n,n-1) = sqrt(beta(n));+                J(n-1,n) = J(n,n-1);+            }++            //TODO : Is this assumption of positive definiteness correct?+            SelfAdjointEigenSolver< Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> > es(J);++            //TODO : make use the eigen assert facilities+            assert(es.info() == Eigen::Success);++            nodes = es.eigenvalues();+            Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> 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<Scalar, Eigen::Dynamic, 1>& nodes,+                              Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<Scalar, Eigen::Dynamic, 1> alpha = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> beta = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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.rows(); ++i)+            {+                nodes(i) = Scalar(2.)*nodes(i) - Scalar(1.);+                weights(i) = Scalar(2.)*weights(i);+            }+        }++        /**+         * \brief Arbitrary precision Gauss abscissae & weights.+         *+         * This method computes Kronrod points for (-1,1) with any required precision.+         * The result is a vector of N nodes the corresponding weights.+         *+         * \param[in] N Number of nodes+         * \param[in,out] nodes Returns a vector of 2N+1 nodes+         * \param[in,out] weights Returns a vector of weights corresponding to \a nodes+         */+        static void mpgauss(const Index N,+                            Eigen::Array<Scalar, Eigen::Dynamic, 1>& nodes,+                            Eigen::Array<Scalar, Eigen::Dynamic, 1>& weights)+        {+            //TODO : make use the eigen assert facilities+            assert(nodes.rows() == N);+            assert(weights.rows() == N);+            assert(N > 0);++            Eigen::Array<Scalar, Eigen::Dynamic, 1> alpha = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*N);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> beta = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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<nodes.rows();++i)+            {+                nodes(i) = Scalar(2.)*nodes(i) - Scalar(1.);+                weights(i) = Scalar(2.)*weights(i);+            }+        }++        /**+         * \brief A method ofr computing Gauss/Kronrod abscissae and weights+         * \param[in] nNodes Gauss-Legendre degree+         * \param[out] abscGaussKronrod Gauss/Kronrod abscissae+         * \param[out] weightGaussKronrod Gauss/Kronrod weights+         * \param[out] abscGauss Gauss abscissae+         * \param[out] weightGauss Gauss weights+         */+        static void computeAbscissaeAndWeights(unsigned int nNodes,+            Eigen::Array<Scalar, Eigen::Dynamic, 1>& abscGaussKronrod,+            Eigen::Array<Scalar, Eigen::Dynamic, 1>& weightGaussKronrod,+            Eigen::Array<Scalar, Eigen::Dynamic, 1>& abscGauss,+            Eigen::Array<Scalar, Eigen::Dynamic, 1>& weightGauss)+        {+            using std::abs;++            Eigen::Array<Scalar, Eigen::Dynamic, 1> xGK = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*nNodes+1);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> wGK = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(2*nNodes+1);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> xG  = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(nNodes);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> wG  = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(nNodes);++            LaurieGautschi::mpkronrod(nNodes,xGK,wGK);+            LaurieGautschi::mpgauss(nNodes,xG,wG);++            unsigned int arraySize = nNodes + 1;++            abscGaussKronrod = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize);+            weightGaussKronrod = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize);+            abscGauss = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize/2);+            weightGauss = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize/2);++            for (unsigned int i=0; i<arraySize; ++i)+            {+                abscGaussKronrod(i) = abs(xGK(i));+                weightGaussKronrod(i) = wGK(i);+            }++            abscGaussKronrod(arraySize-1) = Scalar(0);++            for (unsigned int i=0; i<arraySize/2; ++i)+            {+                abscGauss(i) = abs(xG(i));+                weightGauss(i) = wG(i);+            }+        }+    };++}//namespace Eigen++#endif //EIGEN_LAURIEGAUTSCHI_H

@@ -0,0 +1,380 @@+#ifndef EIGEN_MONEGATO_H+#define EIGEN_MONEGATO_H++namespace Eigen+{++    /**+     * \ingroup NumericalIntegration_Module+     *+     * \class Monegato+     *+     * \brief A class for computing the Gauss/Kronrod weights.+     *+     * \tparam Scalar floating point type+     *+     * This class is based on "Some remarks on the construction of+     * extended Gaussian quadrature rules", Giovanni Monegato,+     * Math. Comp., Vol. 32 (1978) pp. 247-252. http://www.jstor.org/stable/2006272 .+     *+     * The code is based on quadpackcpp (http://quadpackpp.sourceforge.net/)+     * library written by Jerry Gagelman <jerry@os-scientific.org> which is distributed under+     * GNU General Public License+     */+    template <typename Scalar>+    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<Scalar>::epsilon();+            Scalar E1 = abs(x) * NumTraits<Scalar>::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<Scalar,+                                             Eigen::Dynamic, 1>& 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<Scalar, Eigen::Dynamic, 1>& 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<Scalar>::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<Scalar, Eigen::Dynamic, 1>& zeros)+        {+            Eigen::Array<Scalar, Eigen::Dynamic, 1> temp = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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<Scalar>::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<Scalar, Eigen::Dynamic, 1>& coefs)+        {+            int ell = (m_ + 1)/2;+            Eigen::Array<Scalar, Eigen::Dynamic, 1> alpha = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(ell+1);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> f = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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<Scalar, Eigen::Dynamic, 1>& zeros,+                                            const Eigen::Array<Scalar, Eigen::Dynamic, 1>& coefs,+                                            Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<Scalar>::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<Scalar, Eigen::Dynamic, 1>& coefs,+                                          const Eigen::Array<Scalar, Eigen::Dynamic, 1>& xgk_,+                                          Eigen::Array<Scalar, Eigen::Dynamic, 1>& wg_,+                                          Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<Scalar, Eigen::Dynamic, 1>& xgk_,+                                               Eigen::Array<Scalar, Eigen::Dynamic, 1>& wgk_,+                                               Eigen::Array<Scalar, Eigen::Dynamic, 1>& xk_,+                                               Eigen::Array<Scalar, Eigen::Dynamic, 1>& wg_)+        {+            const unsigned int n_ = m_ + 1;++            xgk_ = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_);   //2*nNodes+1+            wgk_ = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_);   //2*nNodes+1+            xk_ = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_/2);  //2*nNodes+            wg_ = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_/2);  //2*nNodes++            // initialise the coefficients to zero+            Eigen::Array<Scalar, Eigen::Dynamic, 1> coefs = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_+1);+            Eigen::Array<Scalar, Eigen::Dynamic, 1> zeros = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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

@@ -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 <unsupported/Eigen/NumericalIntegration>
+ * \endcode
+ * at the start of your source file.
+ */
+}
+
+//#include <iomanip>
+#include <Eigen/Eigenvalues>
+//#include <unsupported/Eigen/MPRealSupport>
+
+#include "ComputeGaussKronrodNodesWeights.h"
+#include "GaussKronrodNodesWeights.h"
+#include "Integrator.h"
+
+#endif // EIGEN_NUMERICAL_INTEGRATION_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 <typename Scalar>+    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<Scalar, Eigen::Dynamic, 1>& abscGaussKronrod,+                            Eigen::Array<Scalar, Eigen::Dynamic, 1>& weightGaussKronrod,+                            Eigen::Array<Scalar, Eigen::Dynamic, 1>& weightGauss)+        {+            unsigned int arraySize = nNodes + 1;+            abscGaussKronrod = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize);+            weightGaussKronrod = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize);+            weightGauss = Eigen::Array<Scalar, Eigen::Dynamic, 1>::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<Scalar, Eigen::Dynamic, 1> tau(m);+            tau(0) = (aN + Scalar(2.0)) / (Scalar(2) * aN + Scalar(3.0));++            Eigen::Array<Scalar, Eigen::Dynamic, 1> 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<RealScalar>::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<Scalar,+                                      Eigen::Dynamic, 1> 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."<<std::endl;+                    return;+                }+            }++            // Computation of the weight.+            d0 = Scalar(1.);+            d1 = abscGaussKronrod;+            ai = Scalar(0.);++            for (size_t k = 0; k < nNodes - 1; ++k)+            {+                ai = ai + Scalar(1.);+                d2 = ((ai + ai + Scalar(1.)) * (abscGaussKronrod * d1) - (ai * d0)) / (ai + Scalar(1.));+                d0 = d1;+                d1 = d2;+            }++            weightGaussKronrod = chebCoeff / (fd * d2);+            return;+        }++        /**+        * \brief abscWeightGauss calculates a Gaussian abscissa and two weights.+        *+        *   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] nNodes The order of the Gauss rule.+        * \param[in] m The value of ( nNodes + 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.+        * \param[in,out] weightGauss The Gauss weight.+        */+        static void abscWeightGauss(unsigned int nNodes,+                                    unsigned int m,+                                    bool even,+                                    const Scalar& chebCoeff,+                                    Eigen::Array<Scalar, Eigen::Dynamic, 1> 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."<<std::endl;;+                    return;+                }+            }++            // Computation of the Gauss weight.+            Scalar aN = nNodes;+            weightGauss = Scalar(2.0) / (aN * pd2 * p0);++            // Initialize+            p1 = Scalar(0.0);+            p2 = betaCoeffs(m);+            yy = Scalar(4.0) * (abscGaussKronrod) * (abscGaussKronrod) - Scalar(2.0);++            for (size_t k = 1; k <= m; ++k)+            {+                p0 = p1;+                p1 = p2;+                p2 = yy * p1 - p0 + betaCoeffs(m - k);+            }++            if (even)+            {+                weightGaussKronrod = weightGauss + chebCoeff / (abscGaussKronrod * pd2 * (p2 - p1));+            }+            else+            {+                weightGaussKronrod = weightGauss + (Scalar(2.0) * chebCoeff) / (pd2 * (p2 - p0));+            }++            return;+        }++        /**+        * \brief machineEpsilon returns the machine precision roundoff error.+        * \returns Returns the machine precision round-off error for double precision+        */+        static Scalar machineEpsilon()+        {+            // ISO standard machine precision values+            // half precision:      2^-10   (9.76563e-04)+            // single precision:    2^-23   (1.1920928955078125e-07)+            // double precision:    2^-52   (2.22044604925031308e-16)+            // extended precision:  2^-63   (1.08420217248550443e-19)+            // quad precision:      2^-112  (1.9259299e-34)++            return NumTraits<Scalar>::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<Scalar, Eigen::Dynamic, 1> &abscGaussKronrod,+            Eigen::Array<Scalar, Eigen::Dynamic, 1> &weightGaussKronrod,+            Eigen::Array<Scalar, Eigen::Dynamic, 1> abscGauss,+            Eigen::Array<Scalar, Eigen::Dynamic, 1> &weightGauss)+        {+            Piessens::kronrod(nNodes,abscGaussKronrod,weightGaussKronrod,weightGauss);+            abscGauss = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(abscGauss.rows());+        }+    };++}//namespace Eigen+#endif //EIGEN_PIESSENS_H

@@ -0,0 +1,46 @@+#define _USE_MATH_DEFINES
+
+#include <iostream> 
+#include <cmath>
+#include <Eigen/Eigenvalues>
+#include <ComputeGaussKronrodNodesWeights.h>
+#include <GaussKronrodNodesWeights.h>
+#include <Integrator.h>
+#include <functional>
+
+extern "C" {
+
+class Integrand {
+ private:
+  std::function<double(double)> f;
+
+ public:
+  Integrand(std::function<double(double)>& 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<double(double)> f_ = [&](double x) { return f(x); };
+  Integrand integrand(f_);
+  // Define the integrator.
+  Eigen::Integrator<double> integrator(subdiv);
+  // Define a quadrature rule.
+  Eigen::Integrator<double>::QuadratureRule rule =
+      Eigen::Integrator<double>::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;
+}
+}

@@ -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

@@ -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