2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Have you tried using something like DLMF 13.4.1? It reduces the problem to an oscillatory integral over $[0,1]$, which can be done with, e.g., gsl's own integration routines. It might not work for very large arguments, in which case DLMF 13.7 or 13.11 might be useful. $\endgroup$ – Kirill Oct 5 '15 at 16:58
  • $\begingroup$ @Kirill That would be very slow. $\endgroup$ – becko Oct 5 '15 at 18:02
  • 2
    $\begingroup$ Slow relative to what? Special functions, in general, tend to be slow to evaluate, and reducing the problem to an integral, continued fraction, recurrence, infinite series, or an ODE, is just about the most common technique. Obviously it is better to have a more direct approximation. $\endgroup$ – Kirill Oct 5 '15 at 18:15
  • $\begingroup$ @Kirill The problem is that this was my initial approach, to compute Kummer's confluent as the Fourier transform of the Beta distribution. But it was too slow for my purposes, so that's why I decided to ask here. I am looking for something faster than this, if possible. $\endgroup$ – becko Dec 30 '15 at 16:18
  • $\begingroup$ I don't know what exactly you tried, but why Fourier transform? A Fourier transform would compute all the values on a grid of $x$'s, but your question is about a single $x$. There are much better ways to evaluate 1d oscillatory integrals than that. You can look at Fredrik Johansson's answer, and see what actual algorithms would be used (fredrikj.net/arb/acb_hypgeom.html) by his code, I think it's mostly power series and asymptotic expansions. $\endgroup$ – Kirill Dec 30 '15 at 23:20
3
$\begingroup$

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));
}
$\endgroup$
  • $\begingroup$ I am getting the following compilation error: error: ‘acb_set_d’ was not declared in this scope $\endgroup$ – becko Dec 24 '15 at 15:43
  • $\begingroup$ I installed arb 2.7.0 on my system. $\endgroup$ – becko Dec 24 '15 at 15:44
  • $\begingroup$ It used to require the git version. I have just released version 2.8.0 which is compatible with the code posted here. $\endgroup$ – Fredrik Johansson Dec 29 '15 at 18:42
2
$\begingroup$

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);

$\endgroup$
  • $\begingroup$ This does not seem to agree with Mathetatica's input: e.g. 1F1(1,1,1+2j) gives incorrect answers when inputted with hy_ge(1,1,1,1+2j) $\endgroup$ – Otto Jan 12 '17 at 6:59

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.