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 void
_arb_atan_taylor_naive(mp_ptr y,mp_limb_t * error,mp_srcptr x,mp_size_t xn,ulong N,int alternating)15 _arb_atan_taylor_naive(mp_ptr y, mp_limb_t * error,
16 mp_srcptr x, mp_size_t xn, ulong N, int alternating)
17 {
18 ulong k;
19 mp_ptr s, t, x1, x2, u;
20 mp_size_t nn = xn + 1;
21
22 if (N == 0)
23 {
24 flint_mpn_zero(y, xn);
25 error[0] = 0;
26 return;
27 }
28
29 if (N == 1)
30 {
31 flint_mpn_copyi(y, x, xn);
32 error[0] = 0;
33 }
34
35 s = flint_malloc(sizeof(mp_limb_t) * nn);
36 t = flint_malloc(sizeof(mp_limb_t) * nn);
37 u = flint_malloc(sizeof(mp_limb_t) * 2 * nn);
38 x1 = flint_malloc(sizeof(mp_limb_t) * nn);
39 x2 = flint_malloc(sizeof(mp_limb_t) * nn);
40
41 flint_mpn_zero(s, nn);
42 flint_mpn_zero(t, nn);
43 flint_mpn_zero(u, 2 * nn);
44 flint_mpn_zero(x1, nn);
45 flint_mpn_zero(x2, nn);
46
47 /* x1 = x */
48 flint_mpn_copyi(x1 + 1, x, xn);
49
50 /* x2 = x * x */
51 mpn_mul_n(u, x1, x1, nn);
52 flint_mpn_copyi(x2, u + nn, nn);
53
54 /* s = t = x */
55 flint_mpn_copyi(s, x1, nn);
56 flint_mpn_copyi(t, x1, nn);
57
58 for (k = 1; k < N; k++)
59 {
60 /* t = t * x2 */
61 mpn_mul_n(u, t, x2, nn);
62 flint_mpn_copyi(t, u + nn, nn);
63
64 /* u = t / (2k+1) */
65 mpn_divrem_1(u, 0, t, nn, 2 * k + 1);
66
67 if (alternating & k)
68 mpn_sub_n(s, s, u, nn);
69 else
70 mpn_add_n(s, s, u, nn);
71 }
72
73 flint_mpn_copyi(y, s + 1, xn);
74 error[0] = 2;
75
76 flint_free(s);
77 flint_free(t);
78 flint_free(u);
79 flint_free(x1);
80 flint_free(x2);
81 }
82
83