Kummer 的 C/C++ 中复杂参数的汇合超几何?

计算科学 C++ C 特殊功能
2021-12-20 18:10:38

我需要评估 Kummer 的融合超几何函数以获取虚构参数:

1F1(a,b;ix)

在哪里i是虚数单位,a,b,x是真实的,并且a,b>0. C/C++ 中是否有可用的例程?我已经看过 GSL,但它只计算1F1对于真正的论点。

2个回答

一种可能性是将任意精度间隔实现包装在Arb中。这不会像专用的双精度实现那样快,但它可能仍然足够快。

注 1:此代码需要 Arb 版本 2.8.0 或更高版本。

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

注意 2:我准备了一个文件 arbcmath.h ( https://github.com/fredrik-johansson/arbcmath ),其中包含 1F1 和许多其他函数(2F1、Bessel、不完全 gamma 等),用于 C99 复数双打,因此您不需要自己实现包装器代码。

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

以下代码:

// 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;
}
// ----------------------------------------------------

计算:1.1833113578-0.0657194591 i

在您的情况下,您可以设置 a=1: hy_ge(1,b,c,z);