Raku调用C实现linspace的精度差异及long double适配问询

问题背景 #

我用C实现了NumPy风格的linspace函数（利用编译语言的循环效率优势），通过Raku NativeCall调用。功能正常，但结果和NumPy存在精度差异：Raku调用结果为1.0013333377777869，NumPy结果为1.0013333377777927。原以为double（对应Raku的num64）能保持15位十进制精度，测试后确认差异来自C的浮点误差，而非NativeCall本身。现需解决两个问题：

1. Raku NativeCall是否支持处理long double或long long？附示例
2. 如何让自定义linspace函数达到NumPy的精度水平

现有代码 #

C代码

//  C code
#include <stdio.h> 
#ifdef _WIN32 
#define DLLEXPORT __declspec(dllexport)
#else 
#define DLLEXPORT extern // if c++ code, requires extern "C"
#endif

DLLEXPORT void c_linspace(double start, double step, int num, double* vals) {
    for (int i = 0; i < num; i++)
    {
        vals[i] = start;
        start += step;
    }
}

Raku代码

// Raku code
sub c_linspace(num64, num64, int32, CArray[num64]) 
    is native('MYDYN') { * };

sub raku_linspace($start, $end, $num, :$endpoint = True, :$retstep = False) {
    my $step = $endpoint ?? ($end - $start)/($num - 1) !!  ($end-$start)/($num);
    my $vals = CArray[num64].allocate($num);
    c_linspace($start.Num, $step.Num, $num.Int, $vals);
    $retstep ?? ($vals.list, $step) !! $vals.list
}

系统信息：Windows 10 64位，gcc 13.2.0

解答 #

1. Raku NativeCall对long double和long long的支持 #

处理long long

Raku NativeCall直接提供int64类型对应C的long long，可以直接在子声明中使用：

C示例代码

DLLEXPORT long long add_long_long(long long a, long long b) {
    return a + b;
}

Raku调用代码

sub add_long_long(int64, int64) is native('your-dll-name') { * }
say add_long_long(9223372036854775806, 1); # 输出9223372036854775807

处理long double

Raku没有内置的long double对应类型，但可以通过Pointer和Blob结合nativecast来处理，需要注意平台差异（Windows下long double通常为80位，占16字节；部分Linux平台为128位）：

C示例代码

DLLEXPORT void set_long_double(long double* dest, long double input) {
    *dest = input;
}

DLLEXPORT long double get_long_double(long double val) {
    return val;
}

Raku调用代码

# 根据平台定义long double的字节大小，Windows下设为16
constant LONG_DOUBLE_SIZE = 16;

sub set_long_double(Pointer, num64) is native('your-dll-name') { * }
sub get_long_double(Pointer) returns num64 is native('your-dll-name') { * }

# 分配Blob存储long double的原始字节
my $ld-buffer = Blob.allocate(LONG_DOUBLE_SIZE);
my $ld-ptr = nativecast(Pointer, $ld-buffer);

# 设置long double值
set_long_double($ld-ptr, 3.14159265358979323846e0);

# 获取并输出值
my $result = get_long_double($ld-ptr);
say $result;

如果需要完整保留long double的精度，建议直接用Blob存储原始字节，避免转换为num64（double）时丢失精度。

2. 让linspace达到NumPy的精度水平 #

NumPy的linspace并未采用循环累加step的方式（这种方式会累积浮点误差），而是直接计算每个索引对应的数值，避免误差累积。修改你的C代码如下：

修改后的C代码

DLLEXPORT void c_linspace(double start, double end, int num, double* vals, int endpoint) {
    double step;
    if (endpoint && num > 1) {
        step = (end - start) / (num - 1);
    } else {
        step = (end - start) / num;
    }
    for (int i = 0; i < num; i++) {
        vals[i] = start + i * step; // 直接计算每个元素，避免累加误差
    }
}

修改后的Raku代码

sub c_linspace(num64, num64, int32, CArray[num64], int32) 
    is native('MYDYN') { * };

sub raku_linspace($start, $end, $num, :$endpoint = True, :$retstep = False) {
    my $step = $endpoint ?? ($end - $start)/($num - 1) !!  ($end-$start)/($num);
    my $vals = CArray[num64].allocate($num);
    c_linspace($start.Num, $end.Num, $num.Int, $vals, $endpoint ?? 1 !! 0);
    $retstep ?? ($vals.list, $step) !! $vals.list
}

这种方式每个元素都是独立计算的，不会累积循环中的浮点误差，结果将与NumPy完全对齐。

内容的提问来源于stack exchange，提问作者Suman