numerical-integration 0.1.1.0 → 0.1.2.0
raw patch · 4 files changed
+473/−360 lines, 4 filesPVP: major bump suggested
API removals or changes: PVP suggests a major version bump
API changes (from Hackage documentation)
- Numerical.Integration: IntegralResult :: CDouble -> CDouble -> Int -> IntegralResult
+ Numerical.Integration: IntegralResult :: Double -> Double -> Int -> IntegralResult
- Numerical.Integration: [_error] :: IntegralResult -> CDouble
+ Numerical.Integration: [_error] :: IntegralResult -> Double
- Numerical.Integration: [_value] :: IntegralResult -> CDouble
+ Numerical.Integration: [_value] :: IntegralResult -> Double
- Numerical.Integration: integration :: (CDouble -> CDouble) -> CDouble -> CDouble -> CDouble -> CInt -> IO IntegralResult
+ Numerical.Integration: integration :: (CDouble -> CDouble) -> Double -> Double -> Double -> Int -> IO IntegralResult
Files
- CHANGELOG.md +20/−12
- README.md +91/−17
- numerical-integration.cabal +291/−289
- src/Numerical/Integration.hs +71/−42
CHANGELOG.md view
@@ -1,12 +1,20 @@-# Changelog for `numerical-integration`--## 0.1.0.0 - 2023-08-30--First release.---## 0.1.1.0 - 2023-09-18--- New data type `IntegralResult` to store the result (value, error, code), instead of in a triplet.--- Replaced `Double` with `CDouble` and `Int` with `CInt`.+# Changelog for `numerical-integration` + +## 0.1.0.0 - 2023-08-30 + +First release. + + +## 0.1.1.0 - 2023-09-18 + +- New data type `IntegralResult` to store the result (value, error, code), instead of in a triplet. + +- Replaced `Double` with `CDouble` and `Int` with `CInt`. + + +## 0.1.2.0 - 2023-09-19 + +- Avoided `CDouble`. + +- Completed README. +
README.md view
@@ -1,17 +1,91 @@-# numerical-integration--One-dimensional numerical integration using the -['NumericalIntegration'](https://github.com/tbs1980/NumericalIntegration) C++ library.--___--***Example.*** Integrate x² between 0 and 1 with desired relative error 1e-5 and -using 200 subdivisions. Exact value: 1/3.--```haskell-example :: IO IntegralResult -- value, error estimate, error code-example = integration (\x -> x*x) 0 1 1e-5 200--- IntegralResult {_value = 0.3333333333333334, _error = 3.7007434154171895e-15, _code = 0}-```--The error code 0 indicates the success.+# numerical-integration + +One-dimensional numerical integration using the +['NumericalIntegration'](https://github.com/tbs1980/NumericalIntegration) C++ library. + +___ + +***Example.*** Integrate x² between 0 and 1 with desired relative error 1e-5 and +using 200 subdivisions. Exact value: 1/3. + +```haskell +example :: IO IntegralResult -- value, error estimate, error code +example = integration (\x -> x*x) 0 1 1e-5 200 +-- IntegralResult { +-- _value = 0.3333333333333334, +-- _error = 3.7007434154171895e-15, +-- _code = 0 +-- } +``` + +The error code 0 indicates the success. + +___ + +***A highly oscillatory function.*** The function shown below is highly +oscillatory. It is know that the exact value of its integral from `0` to `1` +is `π exp(-10) / 2 ≈ 7.131404e-05`. + + + +Let's try to evaluate it with R with 200000 subdivisions. + +```r +f <- function(x) { + 5*cos(2*x/(1-x)) / (25*(1-x)**2 + x**2) +} +integrate(f, 0, 1, subdivisions = 200000) +# 7.76249e-05 with absolute error < 3.7e-05 +``` +R does not complain, however the result is not very good. + +Now let's try with the present library, with 100000 subdivisions. + +```haskell +intgr :: IO IntegralResult +intgr = integration (\t -> 5 * cos(2*t/(1-t)) / (25*(1-t)**2 + t**2)) 0 1 1e-4 100000 +-- IntegralResult { +-- _value = 7.131328051415349e-5, +-- _error = 4.991435083852171e-7, +-- _code = 2 +-- } +``` +As compared to R, the computation is very slow. But the result is quite better. +Note that the error code is 2, thereby indicating a failure of convergence. +So let's try 250000 subdivisions. This will take a while. + +```haskell +intgr :: IO IntegralResult +intgr = integration (\t -> 5 * cos(2*t/(1-t)) / (25*(1-t)**2 + t**2)) 0 1 1e-4 250000 +-- IntegralResult { +-- _value = 7.131328051415349e-5, +-- _error = 4.991435083852171e-7, +-- _code = 2 +-- } +``` + +The result is the same! + +Note that in the C++ code of 'numerical-integration', I fixed the desired +absolute error to 0. I'm not sure it is a good idea. + +___ + +***The tanh-sinh procedure.*** +I don't master the Haskell library 'integration' but I give it a try below. +It implements the tanh-sinh quadrature. + +```haskell +import Numeric.Integration.TanhSinh + +tanhsinh :: Result +tanhsinh = absolute 1e-6 $ parTrap (\t -> 5 * cos(2*t/(1-t)) / (25*(1-t)**2 + t**2)) 0 1 +-- Result { +-- result = -6.463872646093162e-3, +-- errorEstimate = 2.577460946077898e-2, +-- evaluations = 769 +-- } +``` +The result is totally wrong, which is not surprising in view of the weak +number of evaluations. But again, I don't master this library so I will +not conclude anything from this result.
numerical-integration.cabal view
@@ -1,289 +1,291 @@-cabal-version: 2.2-name: numerical-integration-version: 0.1.1.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/Cholesky- , cpp/Eigen/Core- , cpp/Eigen/Eigenvalues- , cpp/Eigen/Geometry- , cpp/Eigen/Householder- , cpp/Eigen/Jacobi- , cpp/Eigen/LU- , cpp/Eigen/QR- , cpp/Eigen/src/Cholesky/LDLT.h- , cpp/Eigen/src/Cholesky/LLT.h- , cpp/Eigen/src/Cholesky/LLT_LAPACKE.h- , cpp/Eigen/src/Core/arch/AltiVec/Complex.h- , cpp/Eigen/src/Core/arch/AltiVec/MathFunctions.h- , cpp/Eigen/src/Core/arch/AltiVec/MatrixProduct.h- , cpp/Eigen/src/Core/arch/AltiVec/MatrixProductCommon.h- , cpp/Eigen/src/Core/arch/AltiVec/MatrixProductMMA.h- , cpp/Eigen/src/Core/arch/AltiVec/PacketMath.h- , cpp/Eigen/src/Core/arch/AVX/Complex.h- , cpp/Eigen/src/Core/arch/AVX/MathFunctions.h- , cpp/Eigen/src/Core/arch/AVX/PacketMath.h- , cpp/Eigen/src/Core/arch/AVX/TypeCasting.h- , cpp/Eigen/src/Core/arch/AVX512/Complex.h- , cpp/Eigen/src/Core/arch/AVX512/MathFunctions.h- , cpp/Eigen/src/Core/arch/AVX512/PacketMath.h- , cpp/Eigen/src/Core/arch/AVX512/TypeCasting.h- , cpp/Eigen/src/Core/arch/CUDA/Complex.h- , cpp/Eigen/src/Core/arch/Default/BFloat16.h- , cpp/Eigen/src/Core/arch/Default/ConjHelper.h- , cpp/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h- , cpp/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h- , cpp/Eigen/src/Core/arch/Default/Half.h- , cpp/Eigen/src/Core/arch/Default/Settings.h- , cpp/Eigen/src/Core/arch/Default/TypeCasting.h- , cpp/Eigen/src/Core/arch/GPU/MathFunctions.h- , cpp/Eigen/src/Core/arch/GPU/PacketMath.h- , cpp/Eigen/src/Core/arch/GPU/TypeCasting.h- , cpp/Eigen/src/Core/arch/HIP/hcc/math_constants.h- , cpp/Eigen/src/Core/arch/MSA/Complex.h- , cpp/Eigen/src/Core/arch/MSA/MathFunctions.h- , cpp/Eigen/src/Core/arch/MSA/PacketMath.h- , cpp/Eigen/src/Core/arch/NEON/Complex.h- , cpp/Eigen/src/Core/arch/NEON/GeneralBlockPanelKernel.h- , cpp/Eigen/src/Core/arch/NEON/MathFunctions.h- , cpp/Eigen/src/Core/arch/NEON/PacketMath.h- , cpp/Eigen/src/Core/arch/NEON/TypeCasting.h- , cpp/Eigen/src/Core/arch/SSE/Complex.h- , cpp/Eigen/src/Core/arch/SSE/MathFunctions.h- , cpp/Eigen/src/Core/arch/SSE/PacketMath.h- , cpp/Eigen/src/Core/arch/SSE/TypeCasting.h- , cpp/Eigen/src/Core/arch/SVE/MathFunctions.h- , cpp/Eigen/src/Core/arch/SVE/PacketMath.h- , cpp/Eigen/src/Core/arch/SVE/TypeCasting.h- , cpp/Eigen/src/Core/arch/SYCL/InteropHeaders.h- , cpp/Eigen/src/Core/arch/SYCL/MathFunctions.h- , cpp/Eigen/src/Core/arch/SYCL/PacketMath.h- , cpp/Eigen/src/Core/arch/SYCL/SyclMemoryModel.h- , cpp/Eigen/src/Core/arch/SYCL/TypeCasting.h- , cpp/Eigen/src/Core/arch/ZVector/Complex.h- , cpp/Eigen/src/Core/arch/ZVector/MathFunctions.h- , cpp/Eigen/src/Core/arch/ZVector/PacketMath.h- , cpp/Eigen/src/Core/ArithmeticSequence.h- , cpp/Eigen/src/Core/Array.h- , cpp/Eigen/src/Core/ArrayBase.h- , cpp/Eigen/src/Core/ArrayWrapper.h- , cpp/Eigen/src/Core/Assign.h- , cpp/Eigen/src/Core/Assign_MKL.h- , cpp/Eigen/src/Core/AssignEvaluator.h- , cpp/Eigen/src/Core/BandMatrix.h- , cpp/Eigen/src/Core/Block.h- , cpp/Eigen/src/Core/BooleanRedux.h- , cpp/Eigen/src/Core/CommaInitializer.h- , cpp/Eigen/src/Core/ConditionEstimator.h- , cpp/Eigen/src/Core/CoreEvaluators.h- , cpp/Eigen/src/Core/CoreIterators.h- , cpp/Eigen/src/Core/CwiseBinaryOp.h- , cpp/Eigen/src/Core/CwiseNullaryOp.h- , cpp/Eigen/src/Core/CwiseTernaryOp.h- , cpp/Eigen/src/Core/CwiseUnaryOp.h- , cpp/Eigen/src/Core/CwiseUnaryView.h- , cpp/Eigen/src/Core/DenseBase.h- , cpp/Eigen/src/Core/DenseCoeffsBase.h- , cpp/Eigen/src/Core/DenseStorage.h- , cpp/Eigen/src/Core/Diagonal.h- , cpp/Eigen/src/Core/DiagonalMatrix.h- , cpp/Eigen/src/Core/DiagonalProduct.h- , cpp/Eigen/src/Core/Dot.h- , cpp/Eigen/src/Core/EigenBase.h- , cpp/Eigen/src/Core/ForceAlignedAccess.h- , cpp/Eigen/src/Core/functors/AssignmentFunctors.h- , cpp/Eigen/src/Core/functors/BinaryFunctors.h- , cpp/Eigen/src/Core/functors/NullaryFunctors.h- , cpp/Eigen/src/Core/functors/StlFunctors.h- , cpp/Eigen/src/Core/functors/TernaryFunctors.h- , cpp/Eigen/src/Core/functors/UnaryFunctors.h- , cpp/Eigen/src/Core/Fuzzy.h- , cpp/Eigen/src/Core/GeneralProduct.h- , cpp/Eigen/src/Core/GenericPacketMath.h- , cpp/Eigen/src/Core/GlobalFunctions.h- , cpp/Eigen/src/Core/IndexedView.h- , cpp/Eigen/src/Core/Inverse.h- , cpp/Eigen/src/Core/IO.h- , cpp/Eigen/src/Core/Map.h- , cpp/Eigen/src/Core/MapBase.h- , cpp/Eigen/src/Core/MathFunctions.h- , cpp/Eigen/src/Core/MathFunctionsImpl.h- , cpp/Eigen/src/Core/Matrix.h- , cpp/Eigen/src/Core/MatrixBase.h- , cpp/Eigen/src/Core/NestByValue.h- , cpp/Eigen/src/Core/NoAlias.h- , cpp/Eigen/src/Core/NumTraits.h- , cpp/Eigen/src/Core/PartialReduxEvaluator.h- , cpp/Eigen/src/Core/PermutationMatrix.h- , cpp/Eigen/src/Core/PlainObjectBase.h- , cpp/Eigen/src/Core/Product.h- , cpp/Eigen/src/Core/ProductEvaluators.h- , cpp/Eigen/src/Core/products/GeneralBlockPanelKernel.h- , cpp/Eigen/src/Core/products/GeneralMatrixMatrix.h- , cpp/Eigen/src/Core/products/GeneralMatrixMatrix_BLAS.h- , cpp/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h- , cpp/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h- , cpp/Eigen/src/Core/products/GeneralMatrixVector.h- , cpp/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h- , cpp/Eigen/src/Core/products/Parallelizer.h- , cpp/Eigen/src/Core/products/SelfadjointMatrixMatrix.h- , cpp/Eigen/src/Core/products/SelfadjointMatrixMatrix_BLAS.h- , cpp/Eigen/src/Core/products/SelfadjointMatrixVector.h- , cpp/Eigen/src/Core/products/SelfadjointMatrixVector_BLAS.h- , cpp/Eigen/src/Core/products/SelfadjointProduct.h- , cpp/Eigen/src/Core/products/SelfadjointRank2Update.h- , cpp/Eigen/src/Core/products/TriangularMatrixMatrix.h- , cpp/Eigen/src/Core/products/TriangularMatrixMatrix_BLAS.h- , cpp/Eigen/src/Core/products/TriangularMatrixVector.h- , cpp/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h- , cpp/Eigen/src/Core/products/TriangularSolverMatrix.h- , cpp/Eigen/src/Core/products/TriangularSolverMatrix_BLAS.h- , cpp/Eigen/src/Core/products/TriangularSolverVector.h- , cpp/Eigen/src/Core/Random.h- , cpp/Eigen/src/Core/Redux.h- , cpp/Eigen/src/Core/Ref.h- , cpp/Eigen/src/Core/Replicate.h- , cpp/Eigen/src/Core/Reshaped.h- , cpp/Eigen/src/Core/ReturnByValue.h- , cpp/Eigen/src/Core/Reverse.h- , cpp/Eigen/src/Core/Select.h- , cpp/Eigen/src/Core/SelfAdjointView.h- , cpp/Eigen/src/Core/SelfCwiseBinaryOp.h- , cpp/Eigen/src/Core/Solve.h- , cpp/Eigen/src/Core/SolverBase.h- , cpp/Eigen/src/Core/SolveTriangular.h- , cpp/Eigen/src/Core/StableNorm.h- , cpp/Eigen/src/Core/StlIterators.h- , cpp/Eigen/src/Core/Stride.h- , cpp/Eigen/src/Core/Swap.h- , cpp/Eigen/src/Core/Transpose.h- , cpp/Eigen/src/Core/Transpositions.h- , cpp/Eigen/src/Core/TriangularMatrix.h- , cpp/Eigen/src/Core/util/BlasUtil.h- , cpp/Eigen/src/Core/util/ConfigureVectorization.h- , cpp/Eigen/src/Core/util/Constants.h- , cpp/Eigen/src/Core/util/DisableStupidWarnings.h- , cpp/Eigen/src/Core/util/ForwardDeclarations.h- , cpp/Eigen/src/Core/util/IndexedViewHelper.h- , cpp/Eigen/src/Core/util/IntegralConstant.h- , cpp/Eigen/src/Core/util/Macros.h- , cpp/Eigen/src/Core/util/Memory.h- , cpp/Eigen/src/Core/util/Meta.h- , cpp/Eigen/src/Core/util/MKL_support.h- , cpp/Eigen/src/Core/util/NonMPL2.h- , cpp/Eigen/src/Core/util/ReenableStupidWarnings.h- , cpp/Eigen/src/Core/util/ReshapedHelper.h- , cpp/Eigen/src/Core/util/StaticAssert.h- , cpp/Eigen/src/Core/util/SymbolicIndex.h- , cpp/Eigen/src/Core/util/XprHelper.h- , cpp/Eigen/src/Core/VectorBlock.h- , cpp/Eigen/src/Core/VectorwiseOp.h- , cpp/Eigen/src/Core/Visitor.h- , cpp/Eigen/src/Eigenvalues/ComplexEigenSolver.h- , cpp/Eigen/src/Eigenvalues/ComplexSchur.h- , cpp/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h- , cpp/Eigen/src/Eigenvalues/EigenSolver.h- , cpp/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h- , cpp/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h- , cpp/Eigen/src/Eigenvalues/HessenbergDecomposition.h- , cpp/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h- , cpp/Eigen/src/Eigenvalues/RealQZ.h- , cpp/Eigen/src/Eigenvalues/RealSchur.h- , cpp/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h- , cpp/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h- , cpp/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h- , cpp/Eigen/src/Eigenvalues/Tridiagonalization.h- , cpp/Eigen/src/Geometry/AlignedBox.h- , cpp/Eigen/src/Geometry/AngleAxis.h- , cpp/Eigen/src/Geometry/arch/Geometry_SIMD.h- , cpp/Eigen/src/Geometry/EulerAngles.h- , cpp/Eigen/src/Geometry/Homogeneous.h- , cpp/Eigen/src/Geometry/Hyperplane.h- , cpp/Eigen/src/Geometry/OrthoMethods.h- , cpp/Eigen/src/Geometry/ParametrizedLine.h- , cpp/Eigen/src/Geometry/Quaternion.h- , cpp/Eigen/src/Geometry/Rotation2D.h- , cpp/Eigen/src/Geometry/RotationBase.h- , cpp/Eigen/src/Geometry/Scaling.h- , cpp/Eigen/src/Geometry/Transform.h- , cpp/Eigen/src/Geometry/Translation.h- , cpp/Eigen/src/Geometry/Umeyama.h- , cpp/Eigen/src/Householder/BlockHouseholder.h- , cpp/Eigen/src/Householder/Householder.h- , cpp/Eigen/src/Householder/HouseholderSequence.h- , cpp/Eigen/src/Jacobi/Jacobi.h- , cpp/Eigen/src/LU/arch/InverseSize4.h- , cpp/Eigen/src/LU/Determinant.h- , cpp/Eigen/src/LU/FullPivLU.h- , cpp/Eigen/src/LU/InverseImpl.h- , cpp/Eigen/src/LU/PartialPivLU.h- , cpp/Eigen/src/LU/PartialPivLU_LAPACKE.h- , cpp/Eigen/src/misc/Image.h- , cpp/Eigen/src/misc/Kernel.h- , cpp/Eigen/src/misc/lapacke.h- , cpp/Eigen/src/misc/RealSvd2x2.h- , cpp/Eigen/src/plugins/ArrayCwiseBinaryOps.h- , cpp/Eigen/src/plugins/ArrayCwiseUnaryOps.h- , cpp/Eigen/src/plugins/BlockMethods.h- , cpp/Eigen/src/plugins/CommonCwiseBinaryOps.h- , cpp/Eigen/src/plugins/CommonCwiseUnaryOps.h- , cpp/Eigen/src/plugins/IndexedViewMethods.h- , cpp/Eigen/src/plugins/MatrixCwiseBinaryOps.h- , cpp/Eigen/src/plugins/MatrixCwiseUnaryOps.h- , cpp/Eigen/src/plugins/ReshapedMethods.h- , cpp/Eigen/src/QR/ColPivHouseholderQR.h- , cpp/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h- , cpp/Eigen/src/QR/CompleteOrthogonalDecomposition.h- , cpp/Eigen/src/QR/FullPivHouseholderQR.h- , cpp/Eigen/src/QR/HouseholderQR.h- , cpp/Eigen/src/QR/HouseholderQR_LAPACKE.h- , cpp/Eigen/src/SVD/BDCSVD.h- , cpp/Eigen/src/SVD/JacobiSVD.h- , cpp/Eigen/src/SVD/JacobiSVD_LAPACKE.h- , cpp/Eigen/src/SVD/SVDBase.h- , cpp/Eigen/src/SVD/UpperBidiagonalization.h- , cpp/Eigen/SVD- , 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+cabal-version: 2.2 +name: numerical-integration +version: 0.1.2.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/Cholesky + , cpp/Eigen/Core + , cpp/Eigen/Eigenvalues + , cpp/Eigen/Geometry + , cpp/Eigen/Householder + , cpp/Eigen/Jacobi + , cpp/Eigen/LU + , cpp/Eigen/QR + , cpp/Eigen/src/Cholesky/LDLT.h + , cpp/Eigen/src/Cholesky/LLT.h + , cpp/Eigen/src/Cholesky/LLT_LAPACKE.h + , cpp/Eigen/src/Core/arch/AltiVec/Complex.h + , cpp/Eigen/src/Core/arch/AltiVec/MathFunctions.h + , cpp/Eigen/src/Core/arch/AltiVec/MatrixProduct.h + , cpp/Eigen/src/Core/arch/AltiVec/MatrixProductCommon.h + , cpp/Eigen/src/Core/arch/AltiVec/MatrixProductMMA.h + , cpp/Eigen/src/Core/arch/AltiVec/PacketMath.h + , cpp/Eigen/src/Core/arch/AVX/Complex.h + , cpp/Eigen/src/Core/arch/AVX/MathFunctions.h + , cpp/Eigen/src/Core/arch/AVX/PacketMath.h + , cpp/Eigen/src/Core/arch/AVX/TypeCasting.h + , cpp/Eigen/src/Core/arch/AVX512/Complex.h + , cpp/Eigen/src/Core/arch/AVX512/MathFunctions.h + , cpp/Eigen/src/Core/arch/AVX512/PacketMath.h + , cpp/Eigen/src/Core/arch/AVX512/TypeCasting.h + , cpp/Eigen/src/Core/arch/CUDA/Complex.h + , cpp/Eigen/src/Core/arch/Default/BFloat16.h + , cpp/Eigen/src/Core/arch/Default/ConjHelper.h + , cpp/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h + , cpp/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h + , cpp/Eigen/src/Core/arch/Default/Half.h + , cpp/Eigen/src/Core/arch/Default/Settings.h + , cpp/Eigen/src/Core/arch/Default/TypeCasting.h + , cpp/Eigen/src/Core/arch/GPU/MathFunctions.h + , cpp/Eigen/src/Core/arch/GPU/PacketMath.h + , cpp/Eigen/src/Core/arch/GPU/TypeCasting.h + , cpp/Eigen/src/Core/arch/HIP/hcc/math_constants.h + , cpp/Eigen/src/Core/arch/MSA/Complex.h + , cpp/Eigen/src/Core/arch/MSA/MathFunctions.h + , cpp/Eigen/src/Core/arch/MSA/PacketMath.h + , cpp/Eigen/src/Core/arch/NEON/Complex.h + , cpp/Eigen/src/Core/arch/NEON/GeneralBlockPanelKernel.h + , cpp/Eigen/src/Core/arch/NEON/MathFunctions.h + , cpp/Eigen/src/Core/arch/NEON/PacketMath.h + , cpp/Eigen/src/Core/arch/NEON/TypeCasting.h + , cpp/Eigen/src/Core/arch/SSE/Complex.h + , cpp/Eigen/src/Core/arch/SSE/MathFunctions.h + , cpp/Eigen/src/Core/arch/SSE/PacketMath.h + , cpp/Eigen/src/Core/arch/SSE/TypeCasting.h + , cpp/Eigen/src/Core/arch/SVE/MathFunctions.h + , cpp/Eigen/src/Core/arch/SVE/PacketMath.h + , cpp/Eigen/src/Core/arch/SVE/TypeCasting.h + , cpp/Eigen/src/Core/arch/SYCL/InteropHeaders.h + , cpp/Eigen/src/Core/arch/SYCL/MathFunctions.h + , cpp/Eigen/src/Core/arch/SYCL/PacketMath.h + , cpp/Eigen/src/Core/arch/SYCL/SyclMemoryModel.h + , cpp/Eigen/src/Core/arch/SYCL/TypeCasting.h + , cpp/Eigen/src/Core/arch/ZVector/Complex.h + , cpp/Eigen/src/Core/arch/ZVector/MathFunctions.h + , cpp/Eigen/src/Core/arch/ZVector/PacketMath.h + , cpp/Eigen/src/Core/ArithmeticSequence.h + , cpp/Eigen/src/Core/Array.h + , cpp/Eigen/src/Core/ArrayBase.h + , cpp/Eigen/src/Core/ArrayWrapper.h + , cpp/Eigen/src/Core/Assign.h + , cpp/Eigen/src/Core/Assign_MKL.h + , cpp/Eigen/src/Core/AssignEvaluator.h + , cpp/Eigen/src/Core/BandMatrix.h + , cpp/Eigen/src/Core/Block.h + , cpp/Eigen/src/Core/BooleanRedux.h + , cpp/Eigen/src/Core/CommaInitializer.h + , cpp/Eigen/src/Core/ConditionEstimator.h + , cpp/Eigen/src/Core/CoreEvaluators.h + , cpp/Eigen/src/Core/CoreIterators.h + , cpp/Eigen/src/Core/CwiseBinaryOp.h + , cpp/Eigen/src/Core/CwiseNullaryOp.h + , cpp/Eigen/src/Core/CwiseTernaryOp.h + , cpp/Eigen/src/Core/CwiseUnaryOp.h + , cpp/Eigen/src/Core/CwiseUnaryView.h + , cpp/Eigen/src/Core/DenseBase.h + , cpp/Eigen/src/Core/DenseCoeffsBase.h + , cpp/Eigen/src/Core/DenseStorage.h + , cpp/Eigen/src/Core/Diagonal.h + , cpp/Eigen/src/Core/DiagonalMatrix.h + , cpp/Eigen/src/Core/DiagonalProduct.h + , cpp/Eigen/src/Core/Dot.h + , cpp/Eigen/src/Core/EigenBase.h + , cpp/Eigen/src/Core/ForceAlignedAccess.h + , cpp/Eigen/src/Core/functors/AssignmentFunctors.h + , cpp/Eigen/src/Core/functors/BinaryFunctors.h + , cpp/Eigen/src/Core/functors/NullaryFunctors.h + , cpp/Eigen/src/Core/functors/StlFunctors.h + , cpp/Eigen/src/Core/functors/TernaryFunctors.h + , cpp/Eigen/src/Core/functors/UnaryFunctors.h + , cpp/Eigen/src/Core/Fuzzy.h + , cpp/Eigen/src/Core/GeneralProduct.h + , cpp/Eigen/src/Core/GenericPacketMath.h + , cpp/Eigen/src/Core/GlobalFunctions.h + , cpp/Eigen/src/Core/IndexedView.h + , cpp/Eigen/src/Core/Inverse.h + , cpp/Eigen/src/Core/IO.h + , cpp/Eigen/src/Core/Map.h + , cpp/Eigen/src/Core/MapBase.h + , cpp/Eigen/src/Core/MathFunctions.h + , cpp/Eigen/src/Core/MathFunctionsImpl.h + , cpp/Eigen/src/Core/Matrix.h + , cpp/Eigen/src/Core/MatrixBase.h + , cpp/Eigen/src/Core/NestByValue.h + , cpp/Eigen/src/Core/NoAlias.h + , cpp/Eigen/src/Core/NumTraits.h + , cpp/Eigen/src/Core/PartialReduxEvaluator.h + , cpp/Eigen/src/Core/PermutationMatrix.h + , cpp/Eigen/src/Core/PlainObjectBase.h + , cpp/Eigen/src/Core/Product.h + , cpp/Eigen/src/Core/ProductEvaluators.h + , cpp/Eigen/src/Core/products/GeneralBlockPanelKernel.h + , cpp/Eigen/src/Core/products/GeneralMatrixMatrix.h + , cpp/Eigen/src/Core/products/GeneralMatrixMatrix_BLAS.h + , cpp/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h + , cpp/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h + , cpp/Eigen/src/Core/products/GeneralMatrixVector.h + , cpp/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h + , cpp/Eigen/src/Core/products/Parallelizer.h + , cpp/Eigen/src/Core/products/SelfadjointMatrixMatrix.h + , cpp/Eigen/src/Core/products/SelfadjointMatrixMatrix_BLAS.h + , cpp/Eigen/src/Core/products/SelfadjointMatrixVector.h + , cpp/Eigen/src/Core/products/SelfadjointMatrixVector_BLAS.h + , cpp/Eigen/src/Core/products/SelfadjointProduct.h + , cpp/Eigen/src/Core/products/SelfadjointRank2Update.h + , cpp/Eigen/src/Core/products/TriangularMatrixMatrix.h + , cpp/Eigen/src/Core/products/TriangularMatrixMatrix_BLAS.h + , cpp/Eigen/src/Core/products/TriangularMatrixVector.h + , cpp/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h + , cpp/Eigen/src/Core/products/TriangularSolverMatrix.h + , cpp/Eigen/src/Core/products/TriangularSolverMatrix_BLAS.h + , cpp/Eigen/src/Core/products/TriangularSolverVector.h + , cpp/Eigen/src/Core/Random.h + , cpp/Eigen/src/Core/Redux.h + , cpp/Eigen/src/Core/Ref.h + , cpp/Eigen/src/Core/Replicate.h + , cpp/Eigen/src/Core/Reshaped.h + , cpp/Eigen/src/Core/ReturnByValue.h + , cpp/Eigen/src/Core/Reverse.h + , cpp/Eigen/src/Core/Select.h + , cpp/Eigen/src/Core/SelfAdjointView.h + , cpp/Eigen/src/Core/SelfCwiseBinaryOp.h + , cpp/Eigen/src/Core/Solve.h + , cpp/Eigen/src/Core/SolverBase.h + , cpp/Eigen/src/Core/SolveTriangular.h + , cpp/Eigen/src/Core/StableNorm.h + , cpp/Eigen/src/Core/StlIterators.h + , cpp/Eigen/src/Core/Stride.h + , cpp/Eigen/src/Core/Swap.h + , cpp/Eigen/src/Core/Transpose.h + , cpp/Eigen/src/Core/Transpositions.h + , cpp/Eigen/src/Core/TriangularMatrix.h + , cpp/Eigen/src/Core/util/BlasUtil.h + , cpp/Eigen/src/Core/util/ConfigureVectorization.h + , cpp/Eigen/src/Core/util/Constants.h + , cpp/Eigen/src/Core/util/DisableStupidWarnings.h + , cpp/Eigen/src/Core/util/ForwardDeclarations.h + , cpp/Eigen/src/Core/util/IndexedViewHelper.h + , cpp/Eigen/src/Core/util/IntegralConstant.h + , cpp/Eigen/src/Core/util/Macros.h + , cpp/Eigen/src/Core/util/Memory.h + , cpp/Eigen/src/Core/util/Meta.h + , cpp/Eigen/src/Core/util/MKL_support.h + , cpp/Eigen/src/Core/util/NonMPL2.h + , cpp/Eigen/src/Core/util/ReenableStupidWarnings.h + , cpp/Eigen/src/Core/util/ReshapedHelper.h + , cpp/Eigen/src/Core/util/StaticAssert.h + , cpp/Eigen/src/Core/util/SymbolicIndex.h + , cpp/Eigen/src/Core/util/XprHelper.h + , cpp/Eigen/src/Core/VectorBlock.h + , cpp/Eigen/src/Core/VectorwiseOp.h + , cpp/Eigen/src/Core/Visitor.h + , cpp/Eigen/src/Eigenvalues/ComplexEigenSolver.h + , cpp/Eigen/src/Eigenvalues/ComplexSchur.h + , cpp/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h + , cpp/Eigen/src/Eigenvalues/EigenSolver.h + , cpp/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h + , cpp/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h + , cpp/Eigen/src/Eigenvalues/HessenbergDecomposition.h + , cpp/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h + , cpp/Eigen/src/Eigenvalues/RealQZ.h + , cpp/Eigen/src/Eigenvalues/RealSchur.h + , cpp/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h + , cpp/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h + , cpp/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h + , cpp/Eigen/src/Eigenvalues/Tridiagonalization.h + , cpp/Eigen/src/Geometry/AlignedBox.h + , cpp/Eigen/src/Geometry/AngleAxis.h + , cpp/Eigen/src/Geometry/arch/Geometry_SIMD.h + , cpp/Eigen/src/Geometry/EulerAngles.h + , cpp/Eigen/src/Geometry/Homogeneous.h + , cpp/Eigen/src/Geometry/Hyperplane.h + , cpp/Eigen/src/Geometry/OrthoMethods.h + , cpp/Eigen/src/Geometry/ParametrizedLine.h + , cpp/Eigen/src/Geometry/Quaternion.h + , cpp/Eigen/src/Geometry/Rotation2D.h + , cpp/Eigen/src/Geometry/RotationBase.h + , cpp/Eigen/src/Geometry/Scaling.h + , cpp/Eigen/src/Geometry/Transform.h + , cpp/Eigen/src/Geometry/Translation.h + , cpp/Eigen/src/Geometry/Umeyama.h + , cpp/Eigen/src/Householder/BlockHouseholder.h + , cpp/Eigen/src/Householder/Householder.h + , cpp/Eigen/src/Householder/HouseholderSequence.h + , cpp/Eigen/src/Jacobi/Jacobi.h + , cpp/Eigen/src/LU/arch/InverseSize4.h + , cpp/Eigen/src/LU/Determinant.h + , cpp/Eigen/src/LU/FullPivLU.h + , cpp/Eigen/src/LU/InverseImpl.h + , cpp/Eigen/src/LU/PartialPivLU.h + , cpp/Eigen/src/LU/PartialPivLU_LAPACKE.h + , cpp/Eigen/src/misc/Image.h + , cpp/Eigen/src/misc/Kernel.h + , cpp/Eigen/src/misc/lapacke.h + , cpp/Eigen/src/misc/RealSvd2x2.h + , cpp/Eigen/src/plugins/ArrayCwiseBinaryOps.h + , cpp/Eigen/src/plugins/ArrayCwiseUnaryOps.h + , cpp/Eigen/src/plugins/BlockMethods.h + , cpp/Eigen/src/plugins/CommonCwiseBinaryOps.h + , cpp/Eigen/src/plugins/CommonCwiseUnaryOps.h + , cpp/Eigen/src/plugins/IndexedViewMethods.h + , cpp/Eigen/src/plugins/MatrixCwiseBinaryOps.h + , cpp/Eigen/src/plugins/MatrixCwiseUnaryOps.h + , cpp/Eigen/src/plugins/ReshapedMethods.h + , cpp/Eigen/src/QR/ColPivHouseholderQR.h + , cpp/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h + , cpp/Eigen/src/QR/CompleteOrthogonalDecomposition.h + , cpp/Eigen/src/QR/FullPivHouseholderQR.h + , cpp/Eigen/src/QR/HouseholderQR.h + , cpp/Eigen/src/QR/HouseholderQR_LAPACKE.h + , cpp/Eigen/src/SVD/BDCSVD.h + , cpp/Eigen/src/SVD/JacobiSVD.h + , cpp/Eigen/src/SVD/JacobiSVD_LAPACKE.h + , cpp/Eigen/src/SVD/SVDBase.h + , cpp/Eigen/src/SVD/UpperBidiagonalization.h + , cpp/Eigen/SVD + , 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 + if os(darwin) + cxx-options: -fPIC -std=gnu++11 + else + cxx-options: -fPIC -std=c++11 + +source-repository head + type: git + location: https://github.com/stla/numerical-integration
src/Numerical/Integration.hs view
@@ -1,42 +1,71 @@-{-# LANGUAGE ForeignFunctionInterface #-}-module Numerical.Integration- (integration, IntegralResult(..))- where-import Foreign.Marshal.Alloc (free, mallocBytes)-import Foreign.Ptr (FunPtr, Ptr, freeHaskellFunPtr)-import Foreign.Storable (peek, sizeOf)-import Foreign.C (CDouble(..), CInt(..)) --data IntegralResult = IntegralResult {- _value :: CDouble,- _error :: CDouble,- _code :: Int-} deriving Show--foreign import ccall safe "wrapper" funPtr- :: (CDouble -> CDouble) -> IO (FunPtr (CDouble -> CDouble))--foreign import ccall safe "integration" c_integration- :: FunPtr (CDouble -> CDouble) -> CDouble -> CDouble -> CDouble -> CInt- -> Ptr CDouble -> Ptr CInt -> IO CDouble---- | Numerical integration.-integration :: (CDouble -> CDouble) -- ^ integrand- -> CDouble -- ^ lower bound- -> CDouble -- ^ upper bound- -> CDouble -- ^ desired relative error- -> CInt -- ^ number of subdivisions- -> IO IntegralResult -- ^ value, error estimate, error code-integration f lower upper relError subdiv = do- errorEstimatePtr <- mallocBytes (sizeOf (0 :: CDouble))- errorCodePtr <- mallocBytes (sizeOf (0 :: CInt))- fPtr <- funPtr f- result <-- c_integration fPtr lower upper relError subdiv errorEstimatePtr errorCodePtr- errorEstimate <- peek errorEstimatePtr- errorCode <- fromIntegral <$> peek errorCodePtr- let out = IntegralResult {_value = result, _error = errorEstimate, _code = errorCode}- free errorEstimatePtr- free errorCodePtr- freeHaskellFunPtr fPtr- return out+{-# LANGUAGE ForeignFunctionInterface #-} +module Numerical.Integration + (integration, IntegralResult(..)) + where +import Foreign.Marshal.Alloc (free, mallocBytes) +import Foreign.Ptr (FunPtr, Ptr, freeHaskellFunPtr) +import Foreign.Storable (peek, sizeOf) +import Foreign.C (CDouble(..), CInt(..)) + +zeroDouble :: Double +zeroDouble = 0.0 + +-- zeroCDouble :: CDouble +-- zeroCDouble = 0.0 + +nanDouble :: Double +nanDouble = zeroDouble / zeroDouble + +-- nanCDouble :: CDouble +-- nanCDouble = zeroCDouble / zeroCDouble + +double2Cdouble :: Double -> CDouble +double2Cdouble = realToFrac + +-- cdouble2double :: CDouble -> Double +-- cdouble2double = realToFrac + + +data IntegralResult = IntegralResult { + _value :: Double, + _error :: Double, + _code :: Int +} deriving Show + +foreign import ccall safe "wrapper" funPtr + :: (CDouble -> CDouble) -> IO (FunPtr (CDouble -> CDouble)) + +foreign import ccall safe "integration" c_integration + :: FunPtr (CDouble -> CDouble) -> CDouble -> CDouble -> CDouble -> CInt + -> Ptr CDouble -> Ptr CInt -> IO CDouble + +-- | Numerical integration. +integration :: (CDouble -> CDouble) -- ^ integrand + -> Double -- ^ lower bound + -> Double -- ^ upper bound + -> Double -- ^ desired relative error + -> Int -- ^ number of subdivisions + -> IO IntegralResult -- ^ value, error estimate, error code +integration f lower upper relError subdiv = do + errorEstimatePtr <- mallocBytes (sizeOf (0 :: CDouble)) + errorCodePtr <- mallocBytes (sizeOf (0 :: CInt)) + fPtr <- funPtr f + let lower' = double2Cdouble lower + upper' = double2Cdouble upper + relError' = double2Cdouble relError + subdiv' = fromIntegral subdiv + result <- + c_integration fPtr lower' upper' relError' subdiv' errorEstimatePtr errorCodePtr + let result' = if isNaN result + then nanDouble + else realToFrac result + errorEstimate <- peek errorEstimatePtr + let errorEstimate' = if isNaN errorEstimate + then nanDouble + else realToFrac errorEstimate + errorCode <- fromIntegral <$> peek errorCodePtr + let out = IntegralResult {_value = result', _error = errorEstimate', _code = errorCode} + free errorEstimatePtr + free errorCodePtr + freeHaskellFunPtr fPtr + return out