1 /*
2 Copyright (C) 2017 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_log(arb_t res,const arb_t x,slong prec)15 arb_log(arb_t res, const arb_t x, slong prec)
16 {
17 if (arb_is_exact(x))
18 {
19 arb_log_arf(res, arb_midref(x), prec);
20 }
21 else
22 {
23 /*
24 Let the input be [a-b, a+b]. We require a > b >= 0 (otherwise the
25 interval contains zero or a negative number and the logarithm is not
26 defined). The error is largest at a-b, and we have
27
28 log(a) - log(a-b) = log(1 + b/(a-b)).
29 */
30 mag_t t;
31 mag_init(t);
32
33 arb_get_mag_lower_nonnegative(t, x);
34
35 if (mag_is_zero(t))
36 {
37 arf_nan(arb_midref(res));
38 mag_inf(arb_radref(res));
39 }
40 else if (mag_is_inf(t))
41 {
42 arf_pos_inf(arb_midref(res));
43 mag_zero(arb_radref(res));
44 }
45 else
46 {
47 slong acc;
48
49 acc = arb_rel_accuracy_bits(x);
50 acc = FLINT_MIN(acc, prec);
51 acc += fmpz_bits(MAG_EXPREF(t));
52
53 if (acc < 20)
54 {
55 mag_t u;
56 mag_init(u);
57 arb_get_mag(u, x);
58
59 if (mag_cmp_2exp_si(t, 0) >= 0)
60 {
61 mag_log_lower(t, t);
62 mag_log(u, u);
63 arb_set_interval_mag(res, t, u, prec);
64 }
65 else if (mag_cmp_2exp_si(u, 0) <= 0)
66 {
67 mag_neg_log_lower(u, u);
68 mag_neg_log(t, t);
69 arb_set_interval_mag(res, u, t, prec);
70 arb_neg(res, res);
71 }
72 else
73 {
74 mag_neg_log(t, t);
75 mag_log(u, u);
76 arb_set_interval_neg_pos_mag(res, t, u, prec);
77 }
78
79 mag_clear(u);
80 }
81 else
82 {
83 acc = FLINT_MAX(acc, 0);
84 acc = FLINT_MIN(acc, prec);
85 prec = FLINT_MIN(prec, acc + MAG_BITS);
86
87 mag_div(t, arb_radref(x), t);
88 mag_log1p(t, t);
89 arb_log_arf(res, arb_midref(x), prec);
90 mag_add(arb_radref(res), arb_radref(res), t);
91 }
92 }
93
94 mag_clear(t);
95 }
96 }
97
98 void
arb_log_ui(arb_t z,ulong x,slong prec)99 arb_log_ui(arb_t z, ulong x, slong prec)
100 {
101 if (x == 2)
102 {
103 arb_const_log2(z, prec);
104 }
105 else if (x == 10)
106 {
107 arb_const_log10(z, prec);
108 }
109 else
110 {
111 arf_t t;
112 arf_init(t);
113 arf_set_ui(t, x);
114 arb_log_arf(z, t, prec);
115 arf_clear(t);
116 }
117 }
118
119 void
arb_log_fmpz(arb_t z,const fmpz_t x,slong prec)120 arb_log_fmpz(arb_t z, const fmpz_t x, slong prec)
121 {
122 arf_t t;
123 arf_init(t);
124 arf_set_fmpz(t, x);
125 arb_log_arf(z, t, prec);
126 arf_clear(t);
127 }
128
129