Kummer's confluent hypergeometric for complex arguments in C/C++?

Question

I need to evaluate Kummer's confluent hypergeometric function for imaginary arguments:

$$_1F_1(a,b;ix)$$

where $i$ is the imaginary unit, $a,b,x$ are real, and $a,b>0$. Is there a routine available in C/C++? I already saw GSL, but it only computes $_1F_1$ for real arguments.

Answer 1

One possibility is to wrap the arbitrary precision interval implementation in Arb. This will not be as fast as a dedicated double precision implementation, but it might still be fast enough.

Note 1: this code requires Arb version 2.8.0 or later.

#include "acb_hypgeom.h"

void
hyp1f1ix(double * re, double * im, double a, double b, double x)
{
    long prec;
    acb_t aa, bb, xx, rr;
    acb_init(aa); acb_init(bb); acb_init(xx); acb_init(rr);

    acb_set_d(aa, a);
    acb_set_d(bb, b);
    acb_set_d(xx, x);
    acb_mul_onei(xx, xx);

    for (prec = 64; ; prec *= 2)
    {
        acb_hypgeom_m(rr, aa, bb, xx, 0, prec);
        if (acb_rel_accuracy_bits(rr) >= 53)
            break;
    }

    *re = arf_get_d(arb_midref(acb_realref(rr)), ARF_RND_DOWN);
    *im = arf_get_d(arb_midref(acb_imagref(rr)), ARF_RND_DOWN);

    acb_clear(aa); acb_clear(bb); acb_clear(xx); acb_clear(rr);
}

int main()
{
    double re, im;
    hyp1f1ix(&re, &im, 3.14, 2.78, 2015.1130);
    printf("%.15g %.15g\n", re, im);
}

Note 2: I have prepared a file arbcmath.h (https://github.com/fredrik-johansson/arbcmath) that wraps 1F1 and many other functions (2F1, Bessel, incomplete gamma, etc.) for use with C99 complex doubles, so you don't need to implement the wrapper code yourself.

#include "arbcmath.h"

int main()
{
    double complex w;
    w = ac_hyp1f1(3.14, 2.78, 2015.1130*I);
    printf("%.15g + %.15g*I\n", creal(w), cimag(w));
}

Answer 2

The following code:

// g++ hg.cc -o hg -L/usr/local/lib -lgsl -lgslcblas -lm && ./hg
#include <stdio.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

// ----------------------------------------------------
gsl_complex operator/(   gsl_complex  z1 , const gsl_complex & z2) { 
    return gsl_complex_div(z1, z2); 
}
gsl_complex operator*(   gsl_complex  z1 , const gsl_complex & z2) { 
    return gsl_complex_mul(z1, z2); 
}
gsl_complex operator+(   gsl_complex  z1 , const gsl_complex & z2) { 
    return gsl_complex_add(z1, z2); 
}
gsl_complex operator+(   gsl_complex z2  , const double  & x ) { 
    return gsl_complex_add( z2,  gsl_complex_rect(x,0)); 
}
gsl_complex operator/(   gsl_complex z2  , const double  & x ) { 
    return gsl_complex_div( z2,  gsl_complex_rect(x,0)); 
}
// ----------------------------------------------------

gsl_complex hy_ge( const gsl_complex & a 
,                  const gsl_complex & b 
,                  const gsl_complex & c 
,                  const gsl_complex & z
,           unsigned int n = 100                   // accuracy
,           unsigned int i = 0                     // itteration step
,            gsl_complex t = gsl_complex_rect(1,0) // coefficient t for internal calculations 
) {

    gsl_complex  t_next = (a+i)*(b+i)/(c+i)/(i+1) * t;

    return (i==n+1) ? gsl_complex_rect(0,0) : t*gsl_complex_pow_real(z,i)+hy_ge(a,b,c,z,n,i+1,t_next); 
}

// ----------------------------------------------------
int main (void)
{

    gsl_complex a = gsl_complex_rect(  1,  0);
    gsl_complex b = gsl_complex_rect(  2,  0);
    gsl_complex c = gsl_complex_rect(  2,  2);
    gsl_complex z = gsl_complex_rect(0.2,0.1);

    gsl_complex z_out = hy_ge(a,b,c,z);

    printf( " %.10f%+.10f i \n",  GSL_REAL( z_out ), GSL_IMAG( z_out ) ) ;

    return 0;
}
// ----------------------------------------------------

calculates: 1.1833113578-0.0657194591 i

in Your case, You may set a=1: hy_ge(1,b,c,z);