contrib/gmp/mpn/generic/get_d.c

   1 /* mpn_get_d -- limbs to double conversion.
   2
   3    THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
   4    CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
   5    FUTURE GNU MP RELEASES.
   6
   7 Copyright 2003, 2004 Free Software Foundation, Inc.
   8
   9 This file is part of the GNU MP Library.
  10
  11 The GNU MP Library is free software; you can redistribute it and/or modify
  12 it under the terms of the GNU Lesser General Public License as published by
  13 the Free Software Foundation; either version 3 of the License, or (at your
  14 option) any later version.
  15
  16 The GNU MP Library is distributed in the hope that it will be useful, but
  17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  18 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  19 License for more details.
  20
  21 You should have received a copy of the GNU Lesser General Public License
  22 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  23
  24 #include "gmp.h"
  25 #include "gmp-impl.h"
  26 #include "longlong.h"
  27
  28 #ifndef _GMP_IEEE_FLOATS
  29 #define _GMP_IEEE_FLOATS 0
  30 #endif
  31
  32 #if ! _GMP_IEEE_FLOATS
  33 /* dummy definition, just to let dead code compile */
  34 union ieee_double_extract {
  35   struct {
  36     int manh, manl, sig, exp;
  37   } s;
  38   double d;
  39 };
  40 #endif
  41
  42 /* To force use of the generic C code for testing, put
  43    "#define _GMP_IEEE_FLOATS 0" at this point.  */
  44
  45
  46
  47 /* In alpha gcc prior to 3.4, signed DI comparisons involving constants are
  48    rearranged from "x < n" to "x+(-n) < 0", which is of course hopelessly
  49    wrong if that addition overflows.
  50
  51    The workaround here avoids this bug by ensuring n is not a literal
  52    constant.  Note that this is alpha specific.  The offending transformation
  53    is/was in alpha.c alpha_emit_conditional_branch() under "We want to use
  54    cmpcc/bcc".
  55
  56    Bizarrely, it turns out this happens also with Cray cc on
  57    alphaev5-cray-unicosmk2.0.6.X, and has the same solution.  Don't know why
  58    or how.  */
  59
  60 #if HAVE_HOST_CPU_FAMILY_alpha                          \
  61   && ((defined (__GNUC__) && ! __GMP_GNUC_PREREQ(3,4))  \
  62       || defined (_CRAY))
  63 static volatile const long CONST_1024 = 1024;
  64 static volatile const long CONST_NEG_1023 = -1023;
  65 static volatile const long CONST_NEG_1022_SUB_53 = -1022 - 53;
  66 #else
  67 #define CONST_1024            (1024)
  68 #define CONST_NEG_1023        (-1023)
  69 #define CONST_NEG_1022_SUB_53 (-1022 - 53)
  70 #endif
  71
  72
  73
  74 /* Return the value {ptr,size}*2^exp, and negative if sign<0.
  75    Must have size>=1, and a non-zero high limb ptr[size-1].
  76
  77    {ptr,size} is truncated towards zero.  This is consistent with other gmp
  78    conversions, like mpz_set_f or mpz_set_q, and is easy to implement and
  79    test.
  80
  81    In the past conversions had attempted (imperfectly) to let the hardware
  82    float rounding mode take effect, but that gets tricky since multiple
  83    roundings need to be avoided, or taken into account, and denorms mean the
  84    effective precision of the mantissa is not constant.  (For reference,
  85    mpz_get_d on IEEE systems was ok, except it operated on the absolute
  86    value.  mpf_get_d and mpq_get_d suffered from multiple roundings and from
  87    not always using enough bits to get the rounding right.)
  88
  89    It's felt that GMP is not primarily concerned with hardware floats, and
  90    really isn't enhanced by getting involved with hardware rounding modes
  91    (which could even be some weird unknown style), so something unambiguous
  92    and straightforward is best.
  93
  94
  95    The IEEE code below is the usual case, it knows either a 32-bit or 64-bit
  96    limb and is done with shifts and masks.  The 64-bit case in particular
  97    should come out nice and compact.
  98
  99    The generic code works one bit at a time, which will be quite slow, but
 100    should support any binary-based "double" and be safe against any rounding
 101    mode.  Note in particular it works on IEEE systems too.
 102
 103
 104    Traps:
 105
 106    Hardware traps for overflow to infinity, underflow to zero, or
 107    unsupported denorms may or may not be taken.  The IEEE code works bitwise
 108    and so probably won't trigger them, the generic code works by float
 109    operations and so probably will.  This difference might be thought less
 110    than ideal, but again its felt straightforward code is better than trying
 111    to get intimate with hardware exceptions (of perhaps unknown nature).
 112
 113
 114    Not done:
 115
 116    mpz_get_d in the past handled size==1 with a cast limb->double.  This
 117    might still be worthwhile there (for up to the mantissa many bits), but
 118    for mpn_get_d here, the cost of applying "exp" to the resulting exponent
 119    would probably use up any benefit a cast may have over bit twiddling.
 120    Also, if the exponent is pushed into denorm range then bit twiddling is
 121    the only option, to ensure the desired truncation is obtained.
 122
 123
 124    Other:
 125
 126    For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL
 127    to the kernel for values >= 2^63.  This makes it slow, and worse the
 128    Linux kernel (what versions?) apparently uses untested code in its trap
 129    handling routines, and gets the sign wrong.  We don't use such a limb to
 130    double cast, neither in the IEEE or generic code.  */
 131
 132
 133 double
 134 mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp)
 135 {
 136   ASSERT (size >= 0);
 137   ASSERT_MPN (up, size);
 138   ASSERT (size == 0 || up[size-1] != 0);
 139
 140   if (size == 0)
 141     return 0.0;
 142
 143   /* Adjust exp to a radix point just above {up,size}, guarding against
 144      overflow.  After this exp can of course be reduced to anywhere within
 145      the {up,size} region without underflow.  */
 146   if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size)
 147                 > (unsigned long) (LONG_MAX - exp)))
 148     {
 149       if (_GMP_IEEE_FLOATS)
 150         goto ieee_infinity;
 151
 152       /* generic */
 153       exp = LONG_MAX;
 154     }
 155   else
 156     {
 157       exp += GMP_NUMB_BITS * size;
 158     }
 159
 160
 161 #if 1
 162 {
 163   int lshift, nbits;
 164   union ieee_double_extract u;
 165   mp_limb_t x, mhi, mlo;
 166 #if GMP_LIMB_BITS == 64
 167   mp_limb_t m;
 168   up += size;
 169   m = *--up;
 170   count_leading_zeros (lshift, m);
 171
 172   exp -= (lshift - GMP_NAIL_BITS) + 1;
 173   m <<= lshift;
 174
 175   nbits = GMP_LIMB_BITS - lshift;
 176
 177   if (nbits < 53 && size > 1)
 178     {
 179       x = *--up;
 180       x <<= GMP_NAIL_BITS;
 181       x >>= nbits;
 182       m |= x;
 183       nbits += GMP_NUMB_BITS;
 184
 185       if (LIMBS_PER_DOUBLE >= 3 && nbits < 53 && size > 2)
 186         {
 187           x = *--up;
 188           x <<= GMP_NAIL_BITS;
 189           x >>= nbits;
 190           m |= x;
 191           nbits += GMP_NUMB_BITS;
 192         }
 193     }
 194   mhi = m >> (32 + 11);
 195   mlo = m >> 11;
 196 #endif
 197 #if GMP_LIMB_BITS == 32
 198   up += size;
 199   x = *--up, size--;
 200   count_leading_zeros (lshift, x);
 201
 202   exp -= (lshift - GMP_NAIL_BITS) + 1;
 203   x <<= lshift;
 204   mhi = x >> 11;
 205
 206   if (lshift < 11)              /* FIXME: never true if NUMB < 20 bits */
 207     {
 208       /* All 20 bits in mhi */
 209       mlo = x << 21;
 210       /* >= 1 bit in mlo */
 211       nbits = GMP_LIMB_BITS - lshift - 21;
 212     }
 213   else
 214     {
 215       if (size != 0)
 216         {
 217           nbits = GMP_LIMB_BITS - lshift;
 218
 219           x = *--up, size--;
 220           x <<= GMP_NAIL_BITS;
 221           mhi |= x >> nbits >> 11;
 222
 223           mlo = x << GMP_LIMB_BITS - nbits - 11;
 224           nbits = nbits + 11 - GMP_NAIL_BITS;
 225         }
 226       else
 227         {
 228           mlo = 0;
 229           goto done;
 230         }
 231     }
 232
 233   if (LIMBS_PER_DOUBLE >= 2 && nbits < 32 && size != 0)
 234     {
 235       x = *--up, size--;
 236       x <<= GMP_NAIL_BITS;
 237       x >>= nbits;
 238       mlo |= x;
 239       nbits += GMP_NUMB_BITS;
 240
 241       if (LIMBS_PER_DOUBLE >= 3 && nbits < 32 && size != 0)
 242         {
 243           x = *--up, size--;
 244           x <<= GMP_NAIL_BITS;
 245           x >>= nbits;
 246           mlo |= x;
 247           nbits += GMP_NUMB_BITS;
 248
 249           if (LIMBS_PER_DOUBLE >= 4 && nbits < 32 && size != 0)
 250             {
 251               x = *--up;
 252               x <<= GMP_NAIL_BITS;
 253               x >>= nbits;
 254               mlo |= x;
 255               nbits += GMP_NUMB_BITS;
 256             }
 257         }
 258     }
 259
 260  done:;
 261
 262 #endif
 263   {
 264     if (UNLIKELY (exp >= CONST_1024))
 265       {
 266         /* overflow, return infinity */
 267       ieee_infinity:
 268         mhi = 0;
 269         mlo = 0;
 270         exp = 1024;
 271       }
 272     else if (UNLIKELY (exp <= CONST_NEG_1023))
 273       {
 274         int rshift;
 275
 276         if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
 277           return 0.0;    /* denorm underflows to zero */
 278
 279         rshift = -1022 - exp;
 280         ASSERT (rshift > 0 && rshift < 53);
 281 #if GMP_LIMB_BITS > 53
 282         mlo >>= rshift;
 283         mhi = mlo >> 32;
 284 #else
 285         if (rshift >= 32)
 286           {
 287             mlo = mhi;
 288             mhi = 0;
 289             rshift -= 32;
 290           }
 291         lshift = GMP_LIMB_BITS - rshift;
 292         mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift);
 293         mhi >>= rshift;
 294 #endif
 295         exp = -1023;
 296       }
 297   }
 298   u.s.manh = mhi;
 299   u.s.manl = mlo;
 300   u.s.exp = exp + 1023;
 301   u.s.sig = (sign < 0);
 302   return u.d;
 303 }
 304 #else
 305
 306
 307 #define ONE_LIMB    (GMP_LIMB_BITS == 64 && 2*GMP_NUMB_BITS >= 53)
 308 #define TWO_LIMBS   (GMP_LIMB_BITS == 32 && 3*GMP_NUMB_BITS >= 53)
 309
 310   if (_GMP_IEEE_FLOATS && (ONE_LIMB || TWO_LIMBS))
 311     {
 312       union ieee_double_extract  u;
 313       mp_limb_t  m0, m1, m2, rmask;
 314       int        lshift, rshift;
 315
 316       m0 = up[size-1];                      /* high limb */
 317       m1 = (size >= 2 ? up[size-2] : 0);   /* second highest limb */
 318       count_leading_zeros (lshift, m0);
 319
 320       /* relative to just under high non-zero bit */
 321       exp -= (lshift - GMP_NAIL_BITS) + 1;
 322
 323       if (ONE_LIMB)
 324         {
 325           /* lshift to have high of m0 non-zero, and collapse nails */
 326           rshift = GMP_LIMB_BITS - lshift;
 327           m1 <<= GMP_NAIL_BITS;
 328           rmask = GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX;
 329           m0 = (m0 << lshift) | ((m1 >> rshift) & rmask);
 330
 331           /* rshift back to have bit 53 of m0 the high non-zero */
 332           m0 >>= 11;
 333         }
 334       else /* TWO_LIMBS */
 335         {
 336           m2 = (size >= 3 ? up[size-3] : 0);  /* third highest limb */
 337
 338           /* collapse nails from m1 and m2 */
 339 #if GMP_NAIL_BITS != 0
 340           m1 = (m1 << GMP_NAIL_BITS) | (m2 >> (GMP_NUMB_BITS-GMP_NAIL_BITS));
 341           m2 <<= 2*GMP_NAIL_BITS;
 342 #endif
 343
 344           /* lshift to have high of m0:m1 non-zero, collapse nails from m0 */
 345           rshift = GMP_LIMB_BITS - lshift;
 346           rmask = (GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX);
 347           m0 = (m0 << lshift) | ((m1 >> rshift) & rmask);
 348           m1 = (m1 << lshift) | ((m2 >> rshift) & rmask);
 349
 350           /* rshift back to have bit 53 of m0:m1 the high non-zero */
 351           m1 = (m1 >> 11) | (m0 << (GMP_LIMB_BITS-11));
 352           m0 >>= 11;
 353         }
 354
 355       if (UNLIKELY (exp >= CONST_1024))
 356         {
 357           /* overflow, return infinity */
 358         ieee_infinity:
 359           m0 = 0;
 360           m1 = 0;
 361           exp = 1024;
 362         }
 363       else if (UNLIKELY (exp <= CONST_NEG_1023))
 364         {
 365           if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
 366             return 0.0;  /* denorm underflows to zero */
 367
 368           rshift = -1022 - exp;
 369           ASSERT (rshift > 0 && rshift < 53);
 370           if (ONE_LIMB)
 371             {
 372               m0 >>= rshift;
 373             }
 374           else /* TWO_LIMBS */
 375             {
 376               if (rshift >= 32)
 377                 {
 378                   m1 = m0;
 379                   m0 = 0;
 380                   rshift -= 32;
 381                 }
 382               lshift = GMP_LIMB_BITS - rshift;
 383               m1 = (m1 >> rshift) | (rshift == 0 ? 0 : m0 << lshift);
 384               m0 >>= rshift;
 385             }
 386           exp = -1023;
 387         }
 388
 389       if (ONE_LIMB)
 390         {
 391 #if GMP_LIMB_BITS > 32  /* avoid compiler warning about big shift */
 392           u.s.manh = m0 >> 32;
 393 #endif
 394           u.s.manl = m0;
 395         }
 396       else /* TWO_LIMBS */
 397         {
 398           u.s.manh = m0;
 399           u.s.manl = m1;
 400         }
 401
 402       u.s.exp = exp + 1023;
 403       u.s.sig = (sign < 0);
 404       return u.d;
 405     }
 406   else
 407     {
 408       /* Non-IEEE or strange limb size, do something generic. */
 409
 410       mp_size_t      i;
 411       mp_limb_t      limb, bit;
 412       int            shift;
 413       double         base, factor, prev_factor, d, new_d, diff;
 414
 415       /* "limb" is "up[i]" the limb being examined, "bit" is a mask for the
 416          bit being examined, initially the highest non-zero bit.  */
 417       i = size-1;
 418       limb = up[i];
 419       count_leading_zeros (shift, limb);
 420       bit = GMP_LIMB_HIGHBIT >> shift;
 421
 422       /* relative to just under high non-zero bit */
 423       exp -= (shift - GMP_NAIL_BITS) + 1;
 424
 425       /* Power up "factor" to 2^exp, being the value of the "bit" in "limb"
 426          being examined.  */
 427       base = (exp >= 0 ? 2.0 : 0.5);
 428       exp = ABS (exp);
 429       factor = 1.0;
 430       for (;;)
 431         {
 432           if (exp & 1)
 433             {
 434               prev_factor = factor;
 435               factor *= base;
 436               FORCE_DOUBLE (factor);
 437               if (factor == 0.0)
 438                 return 0.0;     /* underflow */
 439               if (factor == prev_factor)
 440                 {
 441                   d = factor;     /* overflow, apparent infinity */
 442                   goto generic_done;
 443                 }
 444             }
 445           exp >>= 1;
 446           if (exp == 0)
 447             break;
 448           base *= base;
 449         }
 450
 451       /* Add a "factor" for each non-zero bit, working from high to low.
 452          Stop if any rounding occurs, hence implementing a truncation.
 453
 454          Note no attention is paid to DBL_MANT_DIG, since the effective
 455          number of bits in the mantissa isn't constant when in denorm range.
 456          We also encountered an ARM system with apparently somewhat doubtful
 457          software floats where DBL_MANT_DIG claimed 53 bits but only 32
 458          actually worked.  */
 459
 460       d = factor;  /* high bit */
 461       for (;;)
 462         {
 463           factor *= 0.5;  /* next bit */
 464           bit >>= 1;
 465           if (bit == 0)
 466             {
 467               /* next limb, if any */
 468               i--;
 469               if (i < 0)
 470                 break;
 471               limb = up[i];
 472               bit = GMP_NUMB_HIGHBIT;
 473             }
 474
 475           if (bit & limb)
 476             {
 477               new_d = d + factor;
 478               FORCE_DOUBLE (new_d);
 479               diff = new_d - d;
 480               if (diff != factor)
 481                 break;   /* rounding occured, stop now */
 482               d = new_d;
 483             }
 484         }
 485
 486     generic_done:
 487       return (sign >= 0 ? d : -d);
 488     }
 489 #endif
 490 }