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 "arb.h"
13
14 void
arb_div_2expm1_ui(arb_t y,const arb_t x,ulong n,slong prec)15 arb_div_2expm1_ui(arb_t y, const arb_t x, ulong n, slong prec)
16 {
17 if (n < FLINT_BITS)
18 {
19 arb_div_ui(y, x, (UWORD(1) << n) - 1, prec);
20 }
21 else if (n < 1024 + prec / 32 || n > WORD_MAX / 4)
22 {
23 arb_t t;
24 fmpz_t e;
25
26 arb_init(t);
27 fmpz_init_set_ui(e, n);
28
29 arb_one(t);
30 arb_mul_2exp_fmpz(t, t, e);
31 arb_sub_ui(t, t, 1, prec);
32 arb_div(y, x, t, prec);
33
34 arb_clear(t);
35 fmpz_clear(e);
36 }
37 else
38 {
39 arb_t s, t;
40 slong i, b;
41
42 arb_init(s);
43 arb_init(t);
44
45 /* x / (2^n - 1) = sum_{k>=1} x * 2^(-k*n)*/
46 arb_mul_2exp_si(s, x, -n);
47 arb_set(t, s);
48 b = 1;
49
50 for (i = 2; i <= prec / n + 1; i++)
51 {
52 arb_mul_2exp_si(t, t, -n);
53 arb_add(s, s, t, prec);
54 b = i;
55 }
56
57 /* error bound: sum_{k>b} x * 2^(-k*n) <= x * 2^(-b*n - (n-1)) */
58 arb_mul_2exp_si(t, x, -b*n - (n-1));
59 arb_abs(t, t);
60 arb_add_error(s, t);
61
62 arb_set(y, s);
63
64 arb_clear(s);
65 arb_clear(t);
66 }
67 }
68
69