Haskell封装C++数值积分库：FFI返回异常超大错误码而非0

我用Haskell封装了NumericalIntegration C++库，Hackage上的版本较旧，最新版本可在GitHub获取。

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;
}

Haskell核心代码 #

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 IntegralResult        -- ^ 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 = IntegralResult {_value = result, _error = errorEstimate, _code = errorCode}
  free errorEstimatePtr
  free errorCodePtr
  freeHaskellFunPtr fPtr
  return out

问题描述 #

封装功能整体正常，但存在错误码异常问题：当积分正常完成时，错误码应为0，有时符合预期，但偶尔会返回无意义的超大整数，而积分结果本身是正确的。

问题原因及解决方案 #

核心原因：Haskell与C++的类型不匹配 #

Haskell的Int类型是平台相关的（通常等于机器字长，64位系统下为64位），而C的int类型在多数平台下是32位。当C向int*指针写入32位的错误码时，Haskell用Ptr Int去读取64位数据，会把内存中额外的32位垃圾数据一起读入，导致出现无意义的超大整数。

修复步骤 #

1. 引入Foreign.C.Types模块，使用CInt类型对应C++的int：

   import Foreign.C.Types (CInt)
   

2. 修正foreign import的类型签名，将Int替换为CInt：

   foreign import ccall safe "integration" c_integration
       :: FunPtr (Double -> Double) -> Double -> Double -> Double -> CInt
       -> Ptr Double -> Ptr CInt -> IO Double
   

3. 调整integration函数的参数和内存操作：

   • 将subdiv参数类型改为CInt
   • 用Ptr CInt替代Ptr Int处理错误码
   • 确保IntegralResult的_code字段类型为CInt（如果定义了该数据类型）

   修改后的integration函数示例：

   -- | Numerical integration.
   integration :: (Double -> Double)       -- ^ integrand
               -> Double                   -- ^ lower bound
               -> Double                   -- ^ upper bound
               -> Double                   -- ^ 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 :: Double))
     errorCodePtr <- mallocBytes (sizeOf (0 :: CInt))
     fPtr <- funPtr f
     result <-
       c_integration fPtr lower upper relError subdiv errorEstimatePtr errorCodePtr
     errorEstimate <- peek errorEstimatePtr
     errorCode <- peek errorCodePtr
     let out = IntegralResult {_value = result, _error = errorEstimate, _code = errorCode}
     free errorEstimatePtr
     free errorCodePtr
     freeHaskellFunPtr fPtr
     return out
   

4. （可选）用alloca替代mallocBytes：
   alloca会自动处理内存对齐问题，避免因内存未对齐导致的未定义行为，代码更简洁安全：

   integration f lower upper relError subdiv = alloca $ \errorEstimatePtr ->
     alloca $ \errorCodePtr -> do
       fPtr <- funPtr f
       result <- c_integration fPtr lower upper relError subdiv errorEstimatePtr errorCodePtr
       errorEstimate <- peek errorEstimatePtr
       errorCode <- peek errorCodePtr
       freeHaskellFunPtr fPtr
       return $ IntegralResult result errorEstimate errorCode
   

额外检查点 #

确认Eigen库的Integrator::errorCode()返回类型确实是int，如果返回的是其他整数类型（如long），需要对应使用Haskell的CLong类型。

内容的提问来源于stack exchange，提问作者Stéphane Laurent