使用mpfr::mpreal获取复变量不完全伽马函数的高精度值

问题背景 #

我需要计算复自变量的不完全伽马函数值，具体为实参数a、b对应的Gamma[0, a + i*b]。代码其余部分均使用mpfr::mpreal类型，因此希望伽马函数的输出也保持该类型，将计算结果的实部和虚部分别存储在两个mpfr::mpreal变量中。

目前我通过接收并返回double类型的arb包装器实现计算，代码如下：

complex_double zero;
zero.real = 0;
zero.imag = 0;

complex_double arg1;
arg1.real = a;
arg1.imag = b;

arb_fpwrap_cdouble_gamma_upper(&gamma1, zero, arg1, regularized, flags); // regularized和flags参数与示例无关
// 实部为gamma1.real，虚部为gamma1.imag

但该方法因使用double类型会损失精度，需要找到全程保持精度、输出实部虚部为mpfr::mpreal的方案。

可行解决方案 #

1. 直接使用Arb库的高精度接口 #

Arb本身支持任意精度复数运算，无需依赖double包装器。可以直接通过Arb的acb（高精度复数）类型完成计算，再将结果转换为mpfr::mpreal：

• 将mpfr::mpreal类型的a、b转换为Arb的arb类型（通过mpfr_ptr获取底层MPFR指针，调用arb_set_mpfr赋值）。
• 构造acb类型的自变量a + i*b。
• 调用acb_gamma_upper计算上不完全伽马函数（对应Gamma[0, z]）。
• 提取结果的实部、虚部（arb类型），再转换回mpfr::mpreal（通过arb_get_mpfr获取MPFR指针，赋值给mpfr::mpreal的底层结构）。

示例代码片段：

// 假设a、b为mpfr::mpreal类型
arb_t ar, ai;
arb_init(ar);
arb_init(ai);
acb_t z, res;
acb_init(z);
acb_init(res);

// 将mpfr::mpreal转换为arb类型
arb_set_mpfr(ar, a.mpfr_ptr());
arb_set_mpfr(ai, b.mpfr_ptr());
acb_set_arb_arb(z, ar, ai);

// 计算上不完全伽马函数Gamma(0, z)
acb_gamma_upper(res, NULL, z, 0, MPFR_RNDN); // NULL表示第一个参数（a）为0，最后为舍入模式

// 将结果实部、虚部提取到mpfr::mpreal变量
mpfr::mpreal real_part, imag_part;
arb_get_mpfr(real_part.mpfr_ptr(), acb_realref(res), MPFR_RNDN);
arb_get_mpfr(imag_part.mpfr_ptr(), acb_imagref(res), MPFR_RNDN);

// 清理资源
arb_clear(ar);
arb_clear(ai);
acb_clear(z);
acb_clear(res);

2. 基于MPFR兼容库实现转换 #

如果不想直接操作Arb底层类型，可选择以下方向：

• 使用MPMath-C++：部分接口支持直接传入mpfr::mpreal类型的复数（实部+虚部）计算不完全伽马函数，返回结果可拆分存储。
• 基于mpc库（MPFR的复数扩展）：先将mpfr::mpreal转换为mpc类型，完成计算后再转回mpfr::mpreal的实部和虚部。

注意事项 #

• 确保Arb库的精度设置与mpfr::mpreal的精度一致，避免精度不匹配导致的损失。
• 转换过程中需指定正确的舍入模式（如MPFR_RNDN），保证数值精度一致性。

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