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$
5
  • $\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, 2015 at 16:58
  • $\begingroup$ @Kirill That would be very slow. $\endgroup$
    – a06e
    Oct 5, 2015 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, 2015 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$
    – a06e
    Dec 30, 2015 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, 2015 at 23:20

2 Answers 2

5
$\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$
3
  • $\begingroup$ I am getting the following compilation error: error: ‘acb_set_d’ was not declared in this scope $\endgroup$
    – a06e
    Dec 24, 2015 at 15:43
  • $\begingroup$ I installed arb 2.7.0 on my system. $\endgroup$
    – a06e
    Dec 24, 2015 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$ Dec 29, 2015 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$
1
  • $\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$
    – OTH
    Jan 12, 2017 at 6:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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