A question about using .dll in Mathematica with LibraryFunctionLoad

Question

I am trying to use a custom trilinear interpolation written in C in Mathematica (Windows).

This is the code I start with:

#include "pch.h"
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h> // Include OpenMP library

// Bilinear interpolation using its own grids for x and y
double bilinear_interpolation(double x, double y, double* grid_x, double* grid_y, double* zmax_data, int len_x, int len_y) {
    // Find indices for x
    int idx_x1 = 0;
    while (idx_x1 < len_x && x > grid_x[idx_x1]) idx_x1++;
    idx_x1--;

    idx_x1 = idx_x1 < 0 ? 0 : idx_x1;
    idx_x1 = idx_x1 >= len_x - 1 ? len_x - 2 : idx_x1;
    int idx_x2 = idx_x1 + 1;

    double x1 = grid_x[idx_x1], x2 = grid_x[idx_x2];
    double xd = (x - x1) / (x2 - x1);

    // Find indices for y
    int idx_y1 = 0;
    while (idx_y1 < len_y && y > grid_y[idx_y1]) idx_y1++;
    idx_y1--;

    idx_y1 = idx_y1 < 0 ? 0 : idx_y1;
    idx_y1 = idx_y1 >= len_y - 1 ? len_y - 2 : idx_y1;
    int idx_y2 = idx_y1 + 1;

    double y1 = grid_y[idx_y1], y2 = grid_y[idx_y2];
    double yd = (y - y1) / (y2 - y1);

    // Bilinear interpolation for z_max(x, y)
    double z11 = zmax_data[idx_x1 * len_y + idx_y1];
    double z21 = zmax_data[idx_x2 * len_y + idx_y1];
    double z12 = zmax_data[idx_x1 * len_y + idx_y2];
    double z22 = zmax_data[idx_x2 * len_y + idx_y2];

    double z1 = z11 * (1 - xd) + z21 * xd;
    double z2 = z12 * (1 - xd) + z22 * xd;

    return z1 * (1 - yd) + z2 * yd;
}

__declspec(dllexport) void trilinear_interpolation(int num_points, double fixed_x, double y_min, double y_max, double z_min,
    double* grid_x, double* grid_y, double* grid_z, double* grid_x_bilinear, double* grid_y_bilinear, double* zmax_data, double* distr,
    double* results, int len_x, int len_y, int len_z, int len_x_bilinear, int len_y_bilinear) {
    srand(time(NULL)); // Seed for random number generation

    // Fixed x index
    int idx_x1 = 0;
    while (idx_x1 < len_x && fixed_x > grid_x[idx_x1]) idx_x1++;
    idx_x1--;

    double x1 = grid_x[idx_x1], x2 = grid_x[idx_x1 + 1];
    double xd = (fixed_x - x1) / (x2 - x1);

#pragma omp parallel for
    for (int i = 0; i < num_points; i++) {
        double y = y_min + (y_max - y_min) * rand() / (double)RAND_MAX;
        double z_max = bilinear_interpolation(fixed_x, y, grid_x_bilinear, grid_y_bilinear, zmax_data, len_x_bilinear, len_y_bilinear);
        double z = z_min + (z_max - z_min) * rand() / (double)RAND_MAX;

        // Find indices for y
        int idx_y1 = 0;
        while (idx_y1 < len_y && y > grid_y[idx_y1]) idx_y1++;
        idx_y1--;

        idx_y1 = idx_y1 < 0 ? 0 : idx_y1;
        idx_y1 = idx_y1 >= len_y - 1 ? len_y - 2 : idx_y1;
        int idx_y2 = idx_y1 + 1;

        double y1 = grid_y[idx_y1], y2 = grid_y[idx_y2];
        double yd = (y - y1) / (y2 - y1);

        // Find indices for z
        int idx_z1 = 0;
        while (idx_z1 < len_z && z > grid_z[idx_z1]) idx_z1++;
        idx_z1--;

        idx_z1 = idx_z1 < 0 ? 0 : idx_z1;
        idx_z1 = idx_z1 >= len_z - 1 ? len_z - 2 : idx_z1;
        int idx_z2 = idx_z1 + 1;

        double z1 = grid_z[idx_z1], z2 = grid_z[idx_z2];
        double zd = (z - z1) / (z2 - z1);

        // Trilinear interpolation formula
        double z111 = distr[idx_x1 * len_y * len_z + idx_y1 * len_z + idx_z1];
        double z211 = distr[(idx_x1 + 1) * len_y * len_z + idx_y1 * len_z + idx_z1];
        double z121 = distr[idx_x1 * len_y * len_z + idx_y2 * len_z + idx_z1];
        double z221 = distr[(idx_x1 + 1) * len_y * len_z + idx_y2 * len_z + idx_z1];
        double z112 = distr[idx_x1 * len_y * len_z + idx_y1 * len_z + idx_z2];
        double z212 = distr[(idx_x1 + 1) * len_y * len_z + idx_y1 * len_z + idx_z2];
        double z122 = distr[idx_x1 * len_y * len_z + idx_y2 * len_z + idx_z2];
        double z222 = distr[(idx_x1 + 1) * len_y * len_z + idx_y2 * len_z + idx_z2];

        double c00 = z111 * (1 - xd) + z211 * xd;
        double c01 = z112 * (1 - xd) + z212 * xd;
        double c10 = z121 * (1 - xd) + z221 * xd;
        double c11 = z122 * (1 - xd) + z222 * xd;

        double c0 = c00 * (1 - yd) + c10 * yd;
        double c1 = c01 * (1 - yd) + c11 * yd;

        double final_z = exp(c0 * (1 - zd) + c1 * zd) * (exp(z_max) - z_min);
        results[i] = final_z;
    }
}

And then I export it as .dll.

When trying to import in Mathematica, I use

(*Load and call the C library function*)
Needs["CCompilerDriver`"]
lib = "libtrilinear_interpolation.dll";
trilinearInterpolation = 
  LibraryFunctionLoad[lib, 
   "trilinear_interpolation", {Integer,(*num_points*)Real,(*fixed_x*)
    Real,(*y_min*)Real,(*y_max*)
    Real,(*z_min*){Real, 1},(*grid_x*){Real, 1},(*grid_y*){Real, 
     1},(*grid_z*){Real, 1},(*bilinearGridX*){Real, 
     1},(*bilinearGridY*){Real, 1},(*zmax_data*){Real, 
     1},(*distr*){Real, 1},(*results*)Integer,(*len_x*)
    Integer,(*len_y*)Integer,(*len_z*)Integer,(*len_x_bilinear*)
    Integer  (*len_y_bilinear*)}, "Void"];

However, I get a simple error message

LibraryFunction::libload: The function trilinear_interpolation was not loaded from the file C:\Users\miksi\Dropbox\libtrilinear_interpolation.dll.

I probably miss some very trivial things and do it in a very wrong way. Could you please tell me where the problem can be?