I am searching for a C++ code implementing the complex polygamma function

Looking at what the python scipy library does for its special functions, the polygamma is found by returning the digamma if the zeroth derivative is requested, otherwise return $(-1)^{n+1}Gamma(n+1)zeta(n+1,z)$ where $Gamma$ is the gamma function and $zeta$ the two argument Riemann zeta function. Assuming that this identity holds for the complex numbers, it may explain why you can't find any explicit code.

In fact to add a reference, this is also the method used in the gnu implementation of libstdc++, see here

I implemented a complex polygamma function in the Julia standard library (https://github.com/JuliaLang/julia/pull/7125). It is possible to call this from C++ by linking to the libjulia library.

The easiest way to do this is to call Julia's @cfunction construct to get a C/C++ function pointer to the compiled function. Then you can call it normally like any other C++ function.

In particular (1) download and install Julia, (2) launch Julia and use its package manager to add SpecialFunctions (to install the Julia SpecialFunctions.jl package that now contains my polygamma code), (3) compile your C++ code linked with libjulia. For example, the following C++ program loads Julia and gets a C++ function pointer polygamma to call the Julia polygamma function on 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;
}

The output is polygamma(3, 1.7+4.2i) = (-0.0184344,0.0162077).