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 "arb.h"
13
14 slong _arb_compute_bs_exponents(slong * tab, slong n);
15
16 slong _arb_get_exp_pos(const slong * tab, slong step);
17
18 static void
bsplit(fmpz_t T,fmpz_t Q,flint_bitcnt_t * Qexp,const slong * xexp,const fmpz * xpow,flint_bitcnt_t r,slong a,slong b)19 bsplit(fmpz_t T, fmpz_t Q, flint_bitcnt_t * Qexp,
20 const slong * xexp,
21 const fmpz * xpow, flint_bitcnt_t r, slong a, slong b)
22 {
23 if (b - a == 1)
24 {
25 fmpz_set(T, xpow);
26
27 if (a % 2 == 0)
28 fmpz_neg_ui(Q, 2 * a + 3);
29 else
30 fmpz_set_ui(Q, 2 * a + 3);
31
32 *Qexp = 2 * r;
33 }
34 else if (b - a == 2)
35 {
36 fmpz_mul_ui(T, xpow, 2 * a + 5);
37 fmpz_mul_2exp(T, T, 2 * r);
38 fmpz_submul_ui(T, xpow + 1, 2 * a + 3);
39
40 if (a % 2 == 1)
41 fmpz_neg(T, T);
42
43 fmpz_neg_ui(Q, 2 * a + 3);
44 fmpz_mul_ui(Q, Q, 2 * a + 5);
45 *Qexp = 4 * r;
46 }
47 else
48 {
49 slong step, m, i;
50 flint_bitcnt_t Q2exp[1];
51 fmpz_t Q2, T2;
52
53 step = (b - a) / 2;
54 m = a + step;
55
56 fmpz_init(Q2);
57 fmpz_init(T2);
58
59 bsplit(T, Q, Qexp, xexp, xpow, r, a, m);
60 bsplit(T2, Q2, Q2exp, xexp, xpow, r, m, b);
61
62 fmpz_mul(T, T, Q2);
63 fmpz_mul_2exp(T, T, *Q2exp);
64
65 /* find x^step in table */
66 i = _arb_get_exp_pos(xexp, step);
67 fmpz_mul(T2, T2, Q);
68 fmpz_addmul(T, xpow + i, T2);
69 fmpz_clear(T2);
70
71 fmpz_mul(Q, Q, Q2);
72 *Qexp = *Qexp + *Q2exp;
73 fmpz_clear(Q2);
74 }
75 }
76
77 void
_arb_atan_sum_bs_powtab(fmpz_t T,fmpz_t Q,flint_bitcnt_t * Qexp,const fmpz_t x,flint_bitcnt_t r,slong N)78 _arb_atan_sum_bs_powtab(fmpz_t T, fmpz_t Q, flint_bitcnt_t * Qexp,
79 const fmpz_t x, flint_bitcnt_t r, slong N)
80 {
81 slong * xexp;
82 slong length, i;
83 fmpz * xpow;
84
85 /* compute the powers of x^2 that will appear (at least x^2) */
86 xexp = flint_calloc(2 * FLINT_BITS, sizeof(slong));
87 length = _arb_compute_bs_exponents(xexp, N);
88
89 xpow = _fmpz_vec_init(length);
90 fmpz_mul(xpow, x, x);
91
92 /* build x^i table */
93 for (i = 1; i < length; i++)
94 {
95 if (xexp[i] == 2 * xexp[i-1])
96 {
97 fmpz_mul(xpow + i, xpow + i - 1, xpow + i - 1);
98 }
99 else if (xexp[i] == 2 * xexp[i-2]) /* prefer squaring if possible */
100 {
101 fmpz_mul(xpow + i, xpow + i - 2, xpow + i - 2);
102 }
103 else if (xexp[i] == 2 * xexp[i-1] + 1)
104 {
105 fmpz_mul(xpow + i, xpow + i - 1, xpow + i - 1);
106 fmpz_mul(xpow + i, xpow + i, xpow);
107 }
108 else if (xexp[i] == 2 * xexp[i-2] + 1)
109 {
110 fmpz_mul(xpow + i, xpow + i - 2, xpow + i - 2);
111 fmpz_mul(xpow + i, xpow + i, xpow);
112 }
113 else
114 {
115 flint_printf("power table has the wrong structure!\n");
116 flint_abort();
117 }
118 }
119
120 bsplit(T, Q, Qexp, xexp, xpow, r, 0, N);
121
122 _fmpz_vec_clear(xpow, length);
123 flint_free(xexp);
124 }
125
126