1 /* eta -- Functions for computing the Dedekind eta function
2
3 Copyright (C) 2022 INRIA
4
5 This file is part of GNU MPC.
6
7 GNU MPC is free software; you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with this program. If not, see http://www.gnu.org/licenses/ .
19 */
20
21 #include <math.h>
22 #include <limits.h> /* for CHAR_BIT */
23 #include "mpc-impl.h"
24
25 static void
eta_series(mpcb_ptr eta,mpcb_srcptr q,mpfr_exp_t expq,int N)26 eta_series (mpcb_ptr eta, mpcb_srcptr q, mpfr_exp_t expq, int N)
27 /* Evaluate 2N+1 terms of the Dedekind eta function without the q^(1/24)
28 factor (where internally N is taken to be at least 1).
29 expq is an upper bound on the exponent of |q|, valid everywhere
30 inside the ball; for the error analysis to hold the function assumes
31 that expq < -1, which implies |q| < 1/4. */
32 {
33 const mpfr_prec_t p = mpcb_get_prec (q);
34 mpcb_t q2, qn, q2n1, q3nm1, q3np1;
35 int M, n;
36 mpcr_t r, r2;
37
38 mpcb_init (q2);
39 mpcb_init (qn);
40 mpcb_init (q2n1);
41 mpcb_init (q3nm1);
42 mpcb_init (q3np1);
43
44 mpcb_sqr (q2, q);
45
46 /* n = 0 */
47 mpcb_set_ui_ui (eta, 1, 0, p);
48
49 /* n = 1 */
50 mpcb_set (qn, q); /* q^n */
51 mpcb_neg (q2n1, q); /* -q^(2n-1) */
52 mpcb_neg (q3nm1, q); /* +- q^((3n-1)*n/2) */
53 mpcb_neg (q3np1, q2); /* +- q^(3n+1)*n/2) */
54 mpcb_add (eta, eta, q3nm1);
55 mpcb_add (eta, eta, q3np1);
56
57 N = MPC_MAX (1, N);
58 for (n = 2; n <= N; n++) {
59 mpcb_mul (qn, qn, q);
60 mpcb_mul (q2n1, q2n1, q2);
61 mpcb_mul (q3nm1, q3np1, q2n1);
62 mpcb_mul (q3np1, q3nm1, qn);
63 mpcb_add (eta, eta, q3nm1);
64 mpcb_add (eta, eta, q3np1);
65 }
66
67 /* Compute the relative error due to the truncation of the series
68 as explained in algorithms.tex. */
69 M = (3 * (N+1) - 1) * (N+1) / 2;
70 mpcr_set_one (r);
71 mpcr_div_2ui (r, r, (unsigned long int) (- (M * expq + 1)));
72
73 /* Compose the two relative errors. */
74 mpcr_mul (r2, r, eta->r);
75 mpcr_add (eta->r, eta->r, r);
76 mpcr_add (eta->r, eta->r, r2);
77
78 mpcb_clear (q2);
79 mpcb_clear (qn);
80 mpcb_clear (q2n1);
81 mpcb_clear (q3nm1);
82 mpcb_clear (q3np1);
83 }
84
85
86 static void
mpcb_eta_q24(mpcb_ptr eta,mpcb_srcptr q24)87 mpcb_eta_q24 (mpcb_ptr eta, mpcb_srcptr q24)
88 /* Assuming that q24 is a ball containing
89 q^{1/24} = exp (2 * pi * i * z / 24) for z in the fundamental domain,
90 the function computes eta (z).
91 In fact it works on a larger domain and checks that |q|=|q24^24| < 1/4
92 in the ball; otherwise or if in doubt it returns infinity. */
93 {
94 mpcb_t q;
95 mpfr_exp_t expq;
96 int N;
97
98 mpcb_init (q);
99
100 mpcb_pow_ui (q, q24, 24);
101
102 /* We need an upper bound on the exponent of |q|. Writing q as having
103 the centre x+i*y and the radius r, we have
104 |q| = sqrt (x^2+y^2) |1+\theta| with |theta| <= r
105 <= (1 + r) \sqrt 2 max (|x|, |y|)
106 < 2^{max (Exp x, Exp y) + 1}
107 assuming that r < sqrt 2 - 1, which is the case for r < 1/4
108 or Exp r < -1.
109 Then Exp (|q|) <= max (Exp x, Exp y) + 1. */
110 if (mpcr_inf_p (q->r) || mpcr_get_exp (q->r) >= -1)
111 mpcb_set_inf (eta);
112 else {
113 expq = MPC_MAX (mpfr_get_exp (mpc_realref (q->c)),
114 mpfr_get_exp (mpc_imagref (q->c))) + 1;
115 if (expq >= -1)
116 mpcb_set_inf (eta);
117 else {
118 /* Compute an approximate N such that
119 (3*N+1)*N/2 * |expq| > prec. */
120 N = (int) sqrt (2 * mpcb_get_prec (q24) / 3.0 / (-expq)) + 1;
121 eta_series (eta, q, expq, N);
122 mpcb_mul (eta, eta, q24);
123 }
124 }
125
126 mpcb_clear (q);
127 }
128
129
130 static void
q24_from_z(mpcb_ptr q24,mpc_srcptr z,unsigned long int err_re,unsigned long int err_im)131 q24_from_z (mpcb_ptr q24, mpc_srcptr z, unsigned long int err_re,
132 unsigned long int err_im)
133 /* Given z=x+i*y, compute q24 = exp (pi*i*z/12).
134 err_re and err_im are a priori errors (in 1/2 ulp) of x and y,
135 respectively; they can be 0 if a part is exact. In particular we
136 need err_re=0 when x=0.
137 The function requires and checks that |x|<=5/8 and y>=1/2.
138 Moreover if err_im != 0, it assumes (but cannot check, so this must be
139 assured by the caller) that y is a lower bound on the correct value.
140 The algorithm is taken from algorithms.tex.
141 The precision of q24 is computed from z with a little extra so that
142 the series has a good chance of being rounded to the precision of z. */
143 {
144 const mpfr_prec_t pz = MPC_MAX_PREC (z);
145 int xzero;
146 long int Y, err_a, err_b;
147 mpfr_prec_t p, min;
148 mpfr_t pi, u, v, t, r, s;
149 mpc_t q24c;
150
151 xzero = mpfr_zero_p (mpc_realref (z));
152 if ( mpfr_cmp_d (mpc_realref (z), 0.625) > 0
153 || mpfr_cmp_d (mpc_realref (z), -0.625) < 0
154 || mpfr_cmp_d (mpc_imagref (z), 0.5) < 0
155 || (xzero && err_re > 0))
156 mpcb_set_inf (q24);
157 else {
158 /* Experiments seem to imply that it is enough to add 20 bits to the
159 target precision; to be on the safe side, we also add 1%. */
160 p = pz * 101 / 100 + 20;
161 /* We need 2^p >= 240 + 66 k_x = 240 + 33 err_re. */
162 if (p < (mpfr_prec_t) (CHAR_BIT * sizeof (mpfr_prec_t))) {
163 min = (240 + 33 * err_re) >> p;
164 while (min > 0) {
165 p++;
166 min >>= 1;
167 }
168 }
169
170 mpfr_init2 (pi, p);
171 mpfr_init2 (u, p);
172 mpfr_init2 (v, p);
173 mpfr_init2 (t, p);
174 mpfr_init2 (r, p);
175 mpfr_init2 (s, p);
176 mpc_init2 (q24c, p);
177
178 mpfr_const_pi (pi, MPFR_RNDD);
179 mpfr_div_ui (pi, pi, 12, MPFR_RNDD);
180 mpfr_mul (u, mpc_imagref (z), pi, MPFR_RNDD);
181 mpfr_neg (u, u, MPFR_RNDU);
182 mpfr_mul (v, mpc_realref (z), pi, MPFR_RNDN);
183 mpfr_exp (t, u, MPFR_RNDU);
184 if (xzero) {
185 mpfr_set (mpc_realref (q24c), t, MPFR_RNDN);
186 mpfr_set_ui (mpc_imagref (q24c), 0, MPFR_RNDN);
187 }
188 else {
189 /* Unfortunately we cannot round in two different directions with
190 mpfr_sin_cos. */
191 mpfr_cos (r, v, MPFR_RNDZ);
192 mpfr_sin (s, v, MPFR_RNDA);
193 mpfr_mul (mpc_realref (q24c), t, r, MPFR_RNDN);
194 mpfr_mul (mpc_imagref (q24c), t, s, MPFR_RNDN);
195 }
196 Y = mpfr_get_exp (mpc_imagref (z));
197 if (xzero) {
198 Y = (224 + 33 * err_im + 63) / 64 << Y;
199 err_a = Y + 1;
200 err_b = 0;
201 }
202 else {
203 if (Y >= 2)
204 Y = (32 + 5 * err_im) << (Y - 2);
205 else if (Y == 1)
206 Y = 16 + (5 * err_im + 1) / 2;
207 else /* Y == 0 */
208 Y = 8 + (5 * err_im + 3) / 4;
209 err_a = Y + 9 + err_re;
210 err_b = Y + (67 + 9 * err_re + 1) / 2;
211 }
212 mpcb_set_c (q24, q24c, p, err_a, err_b);
213
214 mpfr_clear (pi);
215 mpfr_clear (u);
216 mpfr_clear (v);
217 mpfr_clear (t);
218 mpfr_clear (r);
219 mpfr_clear (s);
220 mpc_clear (q24c);
221 }
222 }
223
224
225 void
mpcb_eta_err(mpcb_ptr eta,mpc_srcptr z,unsigned long int err_re,unsigned long int err_im)226 mpcb_eta_err (mpcb_ptr eta, mpc_srcptr z, unsigned long int err_re,
227 unsigned long int err_im)
228 /* Given z=x+i*y in the fundamental domain, compute eta (z).
229 err_re and err_im are a priori errors (in 1/2 ulp) of x and y,
230 respectively; they can be 0 if a part is exact. In particular we
231 need err_re=0 when x=0.
232 The function requires (and checks through the call to q24_from_z)
233 that |x|<=5/8 and y>=1/2.
234 Moreover if err_im != 0, it assumes (but cannot check, so this must
235 be assured by the caller) that y is a lower bound on the correct
236 value. */
237 {
238 mpcb_t q24;
239
240 mpcb_init (q24);
241
242 q24_from_z (q24, z, err_re, err_im);
243 mpcb_eta_q24 (eta, q24);
244
245 mpcb_clear (q24);
246 }
247
248
249 int
mpc_eta_fund(mpc_ptr rop,mpc_srcptr z,mpc_rnd_t rnd)250 mpc_eta_fund (mpc_ptr rop, mpc_srcptr z, mpc_rnd_t rnd)
251 /* Given z in the fundamental domain for Sl_2 (Z), that is,
252 |Re z| <= 1/2 and |z| >= 1, compute Dedekind eta (z).
253 Outside the fundamental domain, the function may loop
254 indefinitely. */
255 {
256 mpfr_prec_t prec;
257 mpc_t zl;
258 mpcb_t eta;
259 int xzero, ok, inex;
260
261 mpc_init2 (zl, 2);
262 mpcb_init (eta);
263
264 xzero = mpfr_zero_p (mpc_realref (z));
265 prec = MPC_MAX (MPC_MAX_PREC (rop), MPC_MAX_PREC (z));
266 do {
267 mpc_set_prec (zl, prec);
268 mpc_set (zl, z, MPC_RNDNN); /* exact */
269 mpcb_eta_err (eta, zl, 0, 0);
270
271 if (!xzero)
272 ok = mpcb_can_round (eta, MPC_PREC_RE (rop), MPC_PREC_IM (rop),
273 rnd);
274 else {
275 /* TODO: The result is real, so the ball contains part of the
276 imaginary axis, and rounding to a complex number is impossible
277 independently of the precision.
278 It would be best to project to a real interval and to decide
279 whether we can round. Lacking such functionality, we add
280 the non-representable number 0.1*I (in ball arithmetic) and
281 check whether rounding is possible then. */
282 mpc_t fuzz;
283 mpcb_t fuzzb;
284 mpc_init2 (fuzz, prec);
285 mpcb_init (fuzzb);
286 mpc_set_ui_ui (fuzz, 0, 1, MPC_RNDNN);
287 mpc_div_ui (fuzz, fuzz, 10, MPC_RNDNN);
288 mpcb_set_c (fuzzb, fuzz, prec, 0, 1);
289 ok = mpfr_zero_p (mpc_imagref (eta->c));
290 mpcb_add (eta, eta, fuzzb);
291 ok &= mpcb_can_round (eta, MPC_PREC_RE (rop), 2, rnd);
292 mpc_clear (fuzz);
293 mpcb_clear (fuzzb);
294 }
295
296 prec += 20;
297 } while (!ok);
298
299 if (!xzero)
300 inex = mpcb_round (rop, eta, rnd);
301 else
302 inex = MPC_INEX (mpfr_set (mpc_realref (rop), mpc_realref (eta->c),
303 MPC_RND_RE (rnd)),
304 mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN));
305
306 mpc_clear (zl);
307 mpcb_clear (eta);
308
309 return inex;
310 }
311
312