1 /*
2 Copyright (C) 2015 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "acb_hypgeom.h"
13
14 void
acb_hypgeom_0f1_asymp(acb_t res,const acb_t a,const acb_t z,int regularized,slong prec)15 acb_hypgeom_0f1_asymp(acb_t res, const acb_t a, const acb_t z, int regularized, slong prec)
16 {
17 acb_t t, u, v;
18 int neg;
19
20 acb_init(t);
21 acb_init(u);
22 acb_init(v);
23
24 /* both expansions are correct, but we want the one that works better
25 on the real line */
26 neg = arf_sgn(arb_midref(acb_realref(z))) < 0;
27
28 if (neg)
29 acb_neg(t, z);
30 else
31 acb_set(t, z);
32
33 acb_sqrt(t, t, prec);
34 acb_mul_2exp_si(v, t, 1);
35 acb_sub_ui(u, a, 1, prec);
36
37 if (neg)
38 acb_hypgeom_bessel_j_asymp(v, u, v, prec);
39 else
40 acb_hypgeom_bessel_i_asymp(v, u, v, 0, prec);
41
42 acb_neg(u, u);
43 acb_pow(t, t, u, prec);
44 acb_mul(v, v, t, prec);
45
46 if (!regularized)
47 {
48 acb_gamma(t, a, prec);
49 acb_mul(v, v, t, prec);
50 }
51
52 acb_set(res, v);
53
54 acb_clear(t);
55 acb_clear(u);
56 acb_clear(v);
57 }
58
59 void
acb_hypgeom_0f1_direct(acb_t res,const acb_t a,const acb_t z,int regularized,slong prec)60 acb_hypgeom_0f1_direct(acb_t res, const acb_t a, const acb_t z, int regularized, slong prec)
61 {
62 if (regularized)
63 {
64 if (acb_is_int(a) && arf_sgn(arb_midref(acb_realref(a))) <= 0)
65 {
66 acb_t t, u;
67 acb_init(t);
68 acb_init(u);
69 acb_sub_ui(t, a, 2, prec);
70 acb_neg(t, t);
71 acb_sub_ui(u, a, 1, prec);
72 acb_neg(u, u);
73 acb_pow(u, z, u, prec);
74 /* this cannot recurse infinitely, because t will
75 either be an exact positive integer, or inexact */
76 acb_hypgeom_0f1_direct(res, t, z, regularized, prec);
77 acb_mul(res, res, u, prec);
78 acb_clear(t);
79 acb_clear(u);
80 }
81 else /* todo: could skip when a=1 or a=2 */
82 {
83 acb_t t;
84 acb_init(t);
85 acb_rgamma(t, a, prec);
86 acb_hypgeom_0f1_direct(res, a, z, 0, prec);
87 acb_mul(res, res, t, prec);
88 acb_clear(t);
89 }
90 }
91 else
92 {
93 acb_struct bb[2];
94 bb[0] = *a;
95 acb_init(bb + 1);
96 acb_one(bb + 1);
97 acb_hypgeom_pfq_direct(res, NULL, 0, bb, 2, z, -1, prec);
98 acb_clear(bb + 1);
99 }
100 }
101
102 int
acb_hypgeom_0f1_use_asymp(const acb_t z,slong prec)103 acb_hypgeom_0f1_use_asymp(const acb_t z, slong prec)
104 {
105 double x, y, c;
106
107 if ((arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 0) < 0 &&
108 arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 0) < 0))
109 {
110 return 0;
111 }
112
113 if ((arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 0) > 128 ||
114 arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 0) > 128))
115 {
116 return 1;
117 }
118
119 x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
120 y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
121
122 c = prec * 0.69314718055994530942;
123 c = c * c;
124 c = c * c;
125
126 return x * x + y * y > c;
127 }
128
129 void
acb_hypgeom_0f1(acb_t res,const acb_t a,const acb_t z,int regularized,slong prec)130 acb_hypgeom_0f1(acb_t res, const acb_t a, const acb_t z, int regularized, slong prec)
131 {
132 if (acb_hypgeom_0f1_use_asymp(z, prec))
133 acb_hypgeom_0f1_asymp(res, a, z, regularized, prec);
134 else
135 acb_hypgeom_0f1_direct(res, a, z, regularized, prec);
136 }
137
138