Evaluation of real-valued confluent hypergeometric function with specific complex arguments

Question

In my C++ code, I need to evaluate the confluent hypergeometric function ${}_1F_1(a,b;z)$ with complex arguments in a special case. More precisely, I have to compute $$ e^{-i\phi}{}_1F_1(\ell+1+iZ,2(\ell+1),2i\phi)\, , $$ where $Z,\phi$ are real parameters and $\ell\geq 0$ is an integer. I know that this expression is purely real (but I cannot for the love of god find the reference for this anymore...).

My problem is that the implementations of ${}_1F_1(a,b;z)$ in boost and gls only take real arguments.

• Do you have any ideas on how to evaluate the above expression in C++?
• Since we know that the above expression evaluates to real values, is there maybe a trick, that allows using the boost or gls implementations?

I tried using the integral represenation of ${}_1F_1(a,b;z)$ and perform numerical integration, but this only works for small arguments.

Thanks in advance for any help!

Answer 1

In the meantime, I found a solution myself, and in case someone else stumbles upon this question, I want to share it.

Shortly after posting this question, I found out that the type of expression I want to evaluate goes under the name of Coulomb wave functions. As so often, finding the proper terminology for the problem can be half the solution.

Coulomb wavefunctions are implemented in GSL, but the arguments are limited to certain domains. The C library arb also has an implementation and is applicable over a wider range of arguments. But it is also slower. So I ended up with the following combination of GSL and 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;
}

Thanks to everyone for your replies and inputs!

Answer 2

If you go by the power series route, Kummer's function can be expressed as:

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

(see https://dlmf.nist.gov/13.2)

The complex variable appears only in $z^s$ which can be evaluated with std::pow(std::complex) : https://en.cppreference.com/w/cpp/numeric/complex/pow.

There are numerous papers on less naive ways of evaluating these functions. But the basic approach should give you a good start.