我在 Julia 标准库 ( https://github.com/JuliaLang/julia/pull/7125 ) 中实现了一个复杂的 polygamma 函数。可以通过链接到 libjulia library从 C++ 调用它。

最简单的方法是调用 Julia 的@cfunction构造来获取指向已编译函数的 C/C++ 函数指针。然后你可以像任何其他 C++ 函数一样正常调用它。

特别是（1）下载并安装 Julia，（2）启动 Julia 并使用其包管理器（安装现在包含我的代码add SpecialFunctions的 Julia SpecialFunctions.jl 包polygamma），（3）编译与libjulia. 例如，以下 C++ 程序加载 Julia 并获取 C++ 函数指针polygamma以在 上调用 Julia polygamma 函数std::complex<double>：

#include <julia.h>
JULIA_DEFINE_FAST_TLS() // only define this once, in an executable (not in a shared library) if you want fast code.

#include <complex>
typedef std::complex<double> (*polygamma_ptr)(int m, std::complex<double> z);
polygamma_ptr polygamma = NULL;

#include <iostream>

int main(int argc, char *argv[])
{
    jl_init();

    // get a C function pointer to the Julia polygamma function compiled for
    // (int, complex<double>) arguments:
    jl_eval_string("import SpecialFunctions");
    jl_value_t *ret = jl_eval_string("@cfunction(SpecialFunctions.polygamma, ComplexF64, (Cint, ComplexF64))");
    polygamma = (polygamma_ptr) jl_unbox_voidpointer(ret);

    // example: call polygamma(3, 1.7+4.2i) and output it.
    std::cout << "polygamma(3, 1.7+4.2i) = " << polygamma(3, std::complex<double>(1.7, 4.2)) << "\n";

    jl_atexit_hook(0);
    return 0;
}

输出是polygamma(3, 1.7+4.2i) = (-0.0184344,0.0162077)。

查看 python scipy 库对其特殊功能所做的工作，如果请求零阶导数，则通过返回 digamma 来找到多伽马，否则返回其中是 gamma 函数，是两个参数 Riemann zeta 函数。假设这个恒等式适用于复数，它可以解释为什么你找不到任何明确的代码。$(-1)^{n+1}\Gamma(n+1)\zeta(n+1,z)$$\Gamma$$\zeta$

其实加个引用，这也是libstdc++的gnu实现中使用的方法，看这里