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 static void
arb_log1p_tiny(arb_t r,const arb_t z,slong prec)15 arb_log1p_tiny(arb_t r, const arb_t z, slong prec)
16 {
17 mag_t b, c;
18 arb_t t;
19
20 mag_init(b);
21 mag_init(c);
22 arb_init(t);
23
24 /* if |z| < 1, then |log(1+z) - [z - z^2/2]| <= |z|^3/(1-|z|) */
25 arb_get_mag(b, z);
26 mag_one(c);
27 mag_sub_lower(c, c, b);
28 mag_pow_ui(b, b, 3);
29 mag_div(b, b, c);
30
31 arb_mul(t, z, z, prec);
32 arb_mul_2exp_si(t, t, -1);
33 arb_sub(r, z, t, prec);
34
35 if (mag_is_finite(b))
36 arb_add_error_mag(r, b);
37 else
38 arb_indeterminate(r);
39
40 mag_clear(b);
41 mag_clear(c);
42 arb_clear(t);
43 }
44
45 void
arb_log1p(arb_t r,const arb_t z,slong prec)46 arb_log1p(arb_t r, const arb_t z, slong prec)
47 {
48 slong magz;
49
50 if (arb_is_zero(z))
51 {
52 arb_zero(r);
53 return;
54 }
55
56 magz = arf_abs_bound_lt_2exp_si(arb_midref(z));
57
58 if (magz < -prec)
59 {
60 arb_log1p_tiny(r, z, prec);
61 }
62 else
63 {
64 if (magz < 0)
65 arb_add_ui(r, z, 1, prec + (-magz) + 4);
66 else
67 arb_add_ui(r, z, 1, prec + 4);
68
69 arb_log(r, r, prec);
70 }
71 }
72
73