如果采用幂级数路线，Kummer 函数可以表示为：

$M(a, b, z) = \Gamma(b) \sum_{s = 0}^\infty \cfrac{(a)_s}{\Gamma(b+s) s!} \, z^s$

（见https://dlmf.nist.gov/13.2）

复变量仅出现在$z^s$可以使用 std::pow(std::complex) 进行评估：https ://en.cppreference.com/w/cpp/numeric/complex/pow 。

关于评估这些功能的不那么天真的方法有很多论文。但是基本的方法应该给你一个好的开始。

同时，我自己找到了解决方案，以防其他人偶然发现这个问题，我想分享一下。

发布此问题后不久，我发现我要评估的表达式类型以Coulomb wave functions的名称命名。通常，为问题找到合适的术语可能是解决问题的一半。

库仑波函数在GSL中实现，但参数仅限于某些域。C 库arb也有一个实现，适用于更广泛的参数。但它也较慢。所以我最终得到了 GSL 和 arb 的以下组合。

#include "arb_hypgeom.h"
#include <gsl/gsl_sf_coulomb.h>

double Coulomb_Wave_ARB(int L, double eta, double rho)
{
    double result;
    slong prec;
    arb_t F, l, eta_2, rho_2;
    arb_init(F);
    arb_init(l);
    arb_init(eta_2);
    arb_init(rho_2);
    arb_set_d(l, L);
    arb_set_d(eta_2, eta);
    arb_set_d(rho_2, rho);
    for(prec = 80;; prec *= 2)
    {
        arb_hypgeom_coulomb(F, NULL, l, eta_2, rho_2, prec);
        if(arb_rel_accuracy_bits(F) >= 53)
        {
            result = arf_get_d(arb_midref(F), ARF_RND_NEAR);
            break;
        }
        else if(prec > 10000)
        {
            result = NAN;
            break;
        }
    }
    arb_clear(F);
    arb_clear(eta_2);
    arb_clear(rho_2);
    arb_clear(l);
    return result;
}

double Coulomb_Wave_GSL(int L, double eta, double rho, int& status)
{
    double fc_array[1];
    double F_exponent[1];
    status = gsl_sf_coulomb_wave_F_array(L, 0, eta, rho, fc_array, F_exponent);
    return fc_array[0];
}

double Coulomb_Wave(int L, double eta, double rho)
{
    int status;
    double cw = Coulomb_Wave_GSL(L, eta, rho, status);
    if(status != 0 || std::isnan(cw))
        cw = Coulomb_Wave_ARB(L, eta, rho);
    return cw;
}

感谢大家的回复和投入！