1 /*
2 Copyright (C) 2014 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.h"
13
14 static void
acb_rising_get_mag2_right(mag_t bound,const arb_t a,const arb_t b,ulong n)15 acb_rising_get_mag2_right(mag_t bound, const arb_t a, const arb_t b, ulong n)
16 {
17 mag_t t, u;
18 ulong k;
19
20 mag_init(t);
21 mag_init(u);
22
23 arb_get_mag(t, a);
24 arb_get_mag(u, b);
25
26 mag_mul(bound, t, t);
27 mag_addmul(bound, u, u);
28 mag_set(u, bound);
29 mag_mul_2exp_si(t, t, 1);
30
31 for (k = 1; k < n; k++)
32 {
33 mag_add_ui_2exp_si(u, u, 2 * k - 1, 0);
34 mag_add(u, u, t);
35 mag_mul(bound, bound, u);
36 }
37
38 mag_clear(t);
39 mag_clear(u);
40 }
41
42 void
acb_rising_ui_get_mag(mag_t bound,const acb_t s,ulong n)43 acb_rising_ui_get_mag(mag_t bound, const acb_t s, ulong n)
44 {
45 if (n == 0)
46 {
47 mag_one(bound);
48 return;
49 }
50
51 if (n == 1)
52 {
53 acb_get_mag(bound, s);
54 return;
55 }
56
57 if (!acb_is_finite(s))
58 {
59 mag_inf(bound);
60 return;
61 }
62
63 if (arf_sgn(arb_midref(acb_realref(s))) >= 0)
64 {
65 acb_rising_get_mag2_right(bound, acb_realref(s), acb_imagref(s), n);
66 }
67 else
68 {
69 arb_t a;
70 slong k;
71 mag_t bound2, t, u;
72
73 arb_init(a);
74 mag_init(bound2);
75 mag_init(t);
76 mag_init(u);
77
78 arb_get_mag(u, acb_imagref(s));
79 mag_mul(u, u, u);
80 mag_one(bound);
81
82 for (k = 0; k < n; k++)
83 {
84 arb_add_ui(a, acb_realref(s), k, MAG_BITS);
85
86 if (arf_sgn(arb_midref(a)) >= 0)
87 {
88 acb_rising_get_mag2_right(bound2, a, acb_imagref(s), n - k);
89 mag_mul(bound, bound, bound2);
90 break;
91 }
92 else
93 {
94 arb_get_mag(t, a);
95 mag_mul(t, t, t);
96 mag_add(t, t, u);
97 mag_mul(bound, bound, t);
98 }
99 }
100
101 arb_clear(a);
102 mag_clear(bound2);
103 mag_clear(t);
104 mag_clear(u);
105 }
106
107 mag_sqrt(bound, bound);
108 }
109
110