我需要评估 Kummer 的融合超几何函数以获取虚构参数:
在哪里是虚数单位,是真实的,并且. C/C++ 中是否有可用的例程?我已经看过 GSL,但它只计算对于真正的论点。
我需要评估 Kummer 的融合超几何函数以获取虚构参数:
在哪里是虚数单位,是真实的,并且. C/C++ 中是否有可用的例程?我已经看过 GSL,但它只计算对于真正的论点。
一种可能性是将任意精度间隔实现包装在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);