167 |
return z; |
return z; |
168 |
} |
} |
169 |
|
|
170 |
SCM_C_INLINE_KEYWORD double |
/* scm_i_big2dbl() rounds to the closest representable double, in accordance |
171 |
|
with R5RS exact->inexact. |
172 |
|
|
173 |
|
The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits |
174 |
|
(ie. it truncates towards zero), then adjust to get the closest double by |
175 |
|
examining the next lower bit and adding 1 if necessary. |
176 |
|
|
177 |
|
Note that bignums exactly half way between representable doubles are |
178 |
|
rounded to the next higher absolute value (ie. away from zero). This |
179 |
|
seems like an adequate interpretation of R5RS "numerically closest", and |
180 |
|
it's easier and faster than a full "nearest-even" style. |
181 |
|
|
182 |
|
The bit test is done on the absolute value of the mpz_t, which means we |
183 |
|
must use mpz_getlimbn. mpz_tstbit is not right, it treats negatives as |
184 |
|
twos complement. |
185 |
|
|
186 |
|
Prior to GMP 4.2, the rounding done by mpz_get_d was unspecified. It |
187 |
|
happened to follow the hardware rounding mode, but on the absolute value |
188 |
|
of its operand. This is not what we want, so we put the high |
189 |
|
DBL_MANT_DIG bits into a temporary. This extra init/clear is a slowdown, |
190 |
|
but doesn't matter too much since it's only for older GMP. */ |
191 |
|
|
192 |
|
double |
193 |
scm_i_big2dbl (SCM b) |
scm_i_big2dbl (SCM b) |
194 |
{ |
{ |
195 |
double result = mpz_get_d (SCM_I_BIG_MPZ (b)); |
double result; |
196 |
|
size_t bits; |
197 |
|
|
198 |
|
bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2); |
199 |
|
|
200 |
|
#if __GNU_MP_VERSION < 4 \ |
201 |
|
|| (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2) |
202 |
|
{ |
203 |
|
/* GMP prior to 4.2, force truncate towards zero */ |
204 |
|
mpz_t tmp; |
205 |
|
if (bits > DBL_MANT_DIG) |
206 |
|
{ |
207 |
|
size_t shift = bits - DBL_MANT_DIG; |
208 |
|
mpz_init2 (tmp, DBL_MANT_DIG); |
209 |
|
mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift); |
210 |
|
result = ldexp (mpz_get_d (tmp), shift); |
211 |
|
mpz_clear (tmp); |
212 |
|
} |
213 |
|
else |
214 |
|
{ |
215 |
|
result = mpz_get_d (SCM_I_BIG_MPZ (b)); |
216 |
|
} |
217 |
|
} |
218 |
|
#else |
219 |
|
/* GMP 4.2 and up */ |
220 |
|
result = mpz_get_d (SCM_I_BIG_MPZ (b)); |
221 |
|
#endif |
222 |
|
|
223 |
|
if (bits > DBL_MANT_DIG) |
224 |
|
{ |
225 |
|
unsigned long pos = bits - DBL_MANT_DIG - 1; |
226 |
|
/* test bit number "pos" in absolute value */ |
227 |
|
if (mpz_getlimbn (SCM_I_BIG_MPZ (b), pos / GMP_NUMB_BITS) |
228 |
|
& ((mp_limb_t) 1 << (pos % GMP_NUMB_BITS))) |
229 |
|
{ |
230 |
|
result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1); |
231 |
|
} |
232 |
|
} |
233 |
|
|
234 |
scm_remember_upto_here_1 (b); |
scm_remember_upto_here_1 (b); |
235 |
return result; |
return result; |
236 |
} |
} |