/[guile]/guile/guile-core/libguile/numbers.c
ViewVC logotype

Diff of /guile/guile-core/libguile/numbers.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.205 by kryde, Sun Oct 19 00:27:00 2003 UTC revision 1.206 by kryde, Sun Oct 19 00:32:25 2003 UTC
# Line 167  scm_i_dbl2big (double d) Line 167  scm_i_dbl2big (double d)
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  }  }

Legend:
Removed from v.1.205  
changed lines
  Added in v.1.206

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26