FFI returns a huge integer value instead of 0

Question

I've wrapped the 'NumericalIntegration' C++ library in Haskell. Here is the latest version of the package (the version on Hackage is older).

Here is the main part of the C++ code :

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

And here is the main part of the Haskell code:

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

This works but there's a problem regarding the error code of the integration. When the integration is ok, the error code should be 0. Sometimes it is 0, as expected. But sometimes it is a huge integer number, nonsensical, though the integration is fine.

Would you have an idea about this issue? Why this strange error code? Is there something bad in my code? I'm not fluent in C++ (nor in Haskell). But apart this strange error code, the library seems to work very well.

Answer 1

On the x86_64 platform, a Haskell Int is 64 bits while a C int is 32 bits. (A C long is 64 bits.) In your code, you are picking up garbage in the upper bytes and getting an absurdly large 64-bit integer whose lowest 32-bits are zero, no doubt.

Anyway, there's a module Foreign.C that contains newtypes CInt, CDouble, etc. that are intended to match the corresponding C types on the target platform, and I think it's considered good practice to always use those:

import Foreign
import Foreign.C

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

Of course, since these are newtypes, there is the bother of wrapping and unwrapping, though usually fromIntegral takes care of that automatically, while also converting across bit sizes:

errorCode <- fromIntegral <$> peek errorCodePtr

But, as a less portable alternative, you could stick with Double and Int32 in your foreign import declarations, provided you are only targeting platforms where C ints are 32 bits.

Also note that, if you get the types right, then malloc is a type-safe alternative to mallocBytes:

errorCodePtr <- malloc

The type of the Ptr determines the correct number of bytes to allocate here.