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!

  • 1
    $\begingroup$ Related: scicomp.stackexchange.com/questions/20927/… $\endgroup$
    – Tyberius
    Apr 20, 2021 at 18:43
  • 3
    $\begingroup$ If Fortran code is acceptable, have a look at: M. Nardin, W. F. Perger, and A. Bhalla, "Algorithm 707: CONHYP: A numerical evaluator of the confluent hypergeometric function for complex arguments of large magnitudes." ACM Trans. Math. Software 18 (3), pp. 345–349. $\endgroup$
    – njuffa
    Apr 21, 2021 at 0:07
  • $\begingroup$ I will check out the Fortran code, thank you. $\endgroup$
    – user157765
    Apr 21, 2021 at 9:29
  • 2
    $\begingroup$ See a C++ implementation of CONHYP here: github.com/spring01/libPSI/blob/… $\endgroup$ Apr 22, 2021 at 14:45
  • $\begingroup$ I will give that one a try $\endgroup$
    – user157765
    Apr 23, 2021 at 7:47

2 Answers 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.

  • $\begingroup$ Well, you did not touch the main problem here which is $a$ is complex here. Is there any guarantee that your series would work when $a$ is complex? $\endgroup$ Apr 20, 2021 at 19:52
  • $\begingroup$ @AloneProgrammer The summation can be split into a real part and an imaginary part (in $a$) and is not too problematic. $\endgroup$ Apr 20, 2021 at 22:43
  • 1
    $\begingroup$ @user157765: Which Python library are you using? The standard ones in scipy all assume real $a$. $\endgroup$ Apr 21, 2021 at 20:29
  • 1
    $\begingroup$ I used mpmath in python. $\endgroup$
    – user157765
    Apr 23, 2021 at 7:47
  • 1
    $\begingroup$ @BiswajitBanerjee It seems indeed mpmath is able to evaluate the $_{1}F_{1}$ for complex $a$ and $b$ values: mpmath.org/doc/current/functions/hypergeometric.html#hyp1f1 $\endgroup$ Apr 23, 2021 at 14:27

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_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);
        else if(prec > 10000)
            result = NAN;
    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!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.