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 "mag.h"
13
14 void
mag_neg_log(mag_t z,const mag_t x)15 mag_neg_log(mag_t z, const mag_t x)
16 {
17 if (mag_is_special(x))
18 {
19 if (mag_is_zero(x))
20 mag_inf(z);
21 else
22 mag_zero(z);
23 }
24 else
25 {
26 double t;
27 fmpz exp = MAG_EXP(x);
28
29 if (!COEFF_IS_MPZ(exp))
30 {
31 if (exp >= 1)
32 {
33 mag_zero(z);
34 }
35 else if (exp > -(1000 - MAG_BITS))
36 {
37 t = ldexp(MAG_MAN(x), exp - MAG_BITS);
38 t = -mag_d_log_lower_bound(t);
39 mag_set_d(z, t);
40 }
41 else
42 {
43 /* -log(2^(exp-1) * (2*v)) = -exp*log(2) - log(2*v), 1 <= 2v < 2 */
44 t = MAG_MAN(x) * ldexp(1.0, 1 - MAG_BITS);
45 t = -mag_d_log_lower_bound(t);
46 t -= (exp - 1) * 0.69314718055994530942;
47 t *= (1 + 1e-13);
48 mag_set_d(z, t);
49 }
50 }
51 else if (fmpz_sgn(MAG_EXPREF(x)) > 0)
52 {
53 mag_zero(z);
54 }
55 else
56 {
57 mag_inv(z, x);
58 mag_log(z, z);
59 }
60 }
61 }
62
63 void
mag_neg_log_lower(mag_t z,const mag_t x)64 mag_neg_log_lower(mag_t z, const mag_t x)
65 {
66 if (mag_is_special(x))
67 {
68 if (mag_is_zero(x))
69 mag_inf(z);
70 else
71 mag_zero(z);
72 }
73 else
74 {
75 double t;
76 fmpz exp = MAG_EXP(x);
77
78 if (!COEFF_IS_MPZ(exp))
79 {
80 if (exp >= 1)
81 {
82 mag_zero(z);
83 }
84 else if (exp > -(1000 - MAG_BITS))
85 {
86 t = ldexp(MAG_MAN(x), exp - MAG_BITS);
87 t = -mag_d_log_upper_bound(t);
88 mag_set_d_lower(z, t);
89 }
90 else
91 {
92 /* -log(2^(exp-1) * (2*v)) = -exp*log(2) - log(2*v), 1 <= 2v < 2 */
93 t = MAG_MAN(x) * ldexp(1.0, 1 - MAG_BITS);
94 t = -mag_d_log_upper_bound(t);
95 t -= (exp - 1) * 0.69314718055994530942;
96 t *= (1 - 1e-13);
97 mag_set_d_lower(z, t);
98 }
99 }
100 else if (fmpz_sgn(MAG_EXPREF(x)) > 0)
101 {
102 mag_zero(z);
103 }
104 else
105 {
106 mag_inv_lower(z, x);
107 mag_log_lower(z, z);
108 }
109 }
110 }
111
112 void
mag_log(mag_t z,const mag_t x)113 mag_log(mag_t z, const mag_t x)
114 {
115 if (mag_is_special(x))
116 {
117 if (mag_is_zero(x))
118 mag_zero(z);
119 else
120 mag_inf(z);
121 }
122 else
123 {
124 double t;
125 fmpz exp = MAG_EXP(x);
126
127 if (!COEFF_IS_MPZ(exp))
128 {
129 if (exp <= 0 || (exp == 1 && MAG_MAN(x) == MAG_ONE_HALF))
130 {
131 mag_zero(z);
132 }
133 else if (exp < 1000) /* safely within the double range */
134 {
135 t = ldexp(MAG_MAN(x), exp - MAG_BITS);
136 t = mag_d_log_upper_bound(t);
137 mag_set_d(z, t);
138 }
139 else
140 {
141 /* log(2^(exp-1)*(2v)) = (exp-1)*log(2) + log(2v), 1 <= 2v < 2 */
142 t = MAG_MAN(x) * ldexp(1.0, 1 - MAG_BITS); /* exact */
143 t = mag_d_log_upper_bound(t);
144 t += (exp - 1.0) * 0.69314718055994530942;
145 t *= (1 + 1e-13); /* last step could have rounded down */
146 mag_set_d(z, t);
147 }
148 }
149 else if (fmpz_sgn(MAG_EXPREF(x)) < 0)
150 {
151 mag_zero(z);
152 }
153 else
154 {
155 /* log(2^exp) = exp*log(2), log(2) < 744261118/2^30 */
156 fmpz_t t;
157 fmpz_init(t);
158 fmpz_mul_ui(t, MAG_EXPREF(x), 744261118);
159 mag_set_fmpz(z, t);
160 mag_mul_2exp_si(z, z, -30);
161 fmpz_clear(t);
162 }
163 }
164 }
165
166 void
mag_log_lower(mag_t z,const mag_t x)167 mag_log_lower(mag_t z, const mag_t x)
168 {
169 if (mag_is_special(x))
170 {
171 if (mag_is_zero(x))
172 mag_zero(z);
173 else
174 mag_inf(z);
175 }
176 else
177 {
178 double t;
179 fmpz exp = MAG_EXP(x);
180
181 if (!COEFF_IS_MPZ(exp))
182 {
183 if (exp <= 0 || (exp == 1 && MAG_MAN(x) == MAG_ONE_HALF))
184 {
185 mag_zero(z);
186 }
187 else if (exp < 1000)
188 {
189 t = ldexp(MAG_MAN(x), exp - MAG_BITS);
190 t = mag_d_log_lower_bound(t);
191 mag_set_d_lower(z, t);
192 }
193 else
194 {
195 /* log(2^(exp-1) * (2*v)) = exp*log(2) + log(2*v), 1 <= 2v < 2 */
196 t = MAG_MAN(x) * ldexp(1.0, 1 - MAG_BITS);
197 t = mag_d_log_lower_bound(t);
198 t += (exp - 1) * 0.69314718055994530942;
199 t *= (1 - 1e-13);
200 mag_set_d_lower(z, t);
201 }
202 }
203 else if (fmpz_sgn(MAG_EXPREF(x)) < 0)
204 {
205 mag_zero(z);
206 }
207 else
208 {
209 /* log(2^exp) = exp*log(2), log(2) > 744261117/2^30 */
210 fmpz_t t;
211 fmpz_init(t);
212 fmpz_sub_ui(t, MAG_EXPREF(x), 1); /* lower bound for x */
213 fmpz_mul_ui(t, t, 744261117);
214 mag_set_fmpz_lower(z, t);
215 mag_mul_2exp_si(z, z, -30);
216 fmpz_clear(t);
217 }
218 }
219 }
220
221