source/libs/gmp/gmp-src/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, 2007, 2009, 2010, 2012 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 either:
  13
  14   * the GNU Lesser General Public License as published by the Free
  15     Software Foundation; either version 3 of the License, or (at your
  16     option) any later version.
  17
  18 or
  19
  20   * the GNU General Public License as published by the Free Software
  21     Foundation; either version 2 of the License, or (at your option) any
  22     later version.
  23
  24 or both in parallel, as here.
  25
  26 The GNU MP Library is distributed in the hope that it will be useful, but
  27 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  28 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  29 for more details.
  30
  31 You should have received copies of the GNU General Public License and the
  32 GNU Lesser General Public License along with the GNU MP Library.  If not,
  33 see https://www.gnu.org/licenses/.  */
  34
  35 #include "gmp.h"
  36 #include "gmp-impl.h"
  37 #include "longlong.h"
  38
  39 #ifndef _GMP_IEEE_FLOATS
  40 #define _GMP_IEEE_FLOATS 0
  41 #endif
  42
  43 /* To force use of the generic C code for testing, put
  44    "#define _GMP_IEEE_FLOATS 0" at this point.  */
  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 constant.
  52    Note that this is alpha specific.  The offending transformation is/was in
  53    alpha.c alpha_emit_conditional_branch() under "We want to use cmpcc/bcc".
  54
  55    Bizarrely, this happens also with Cray cc on alphaev5-cray-unicosmk2.0.6.X,
  56    and has the same solution.  Don't know why or how.  */
  57
  58 #if HAVE_HOST_CPU_FAMILY_alpha                          \
  59   && ((defined (__GNUC__) && ! __GMP_GNUC_PREREQ(3,4))  \
  60       || defined (_CRAY))
  61 static volatile const long CONST_1024 = 1024;
  62 static volatile const long CONST_NEG_1023 = -1023;
  63 static volatile const long CONST_NEG_1022_SUB_53 = -1022 - 53;
  64 #else
  65 #define CONST_1024            (1024)
  66 #define CONST_NEG_1023        (-1023)
  67 #define CONST_NEG_1022_SUB_53 (-1022 - 53)
  68 #endif
  69
  70
  71 /* Return the value {ptr,size}*2^exp, and negative if sign<0.  Must have
  72    size>=1, and a non-zero high limb ptr[size-1].
  73
  74    When we know the fp format, the result is truncated towards zero.  This is
  75    consistent with other gmp conversions, like mpz_set_f or mpz_set_q, and is
  76    easy to implement and test.
  77
  78    When we do not know the format, such truncation seems much harder.  One
  79    would need to defeat any rounding mode, including round-up.
  80
  81    It's felt that GMP is not primarily concerned with hardware floats, and
  82    really isn't enhanced by getting involved with hardware rounding modes
  83    (which could even be some weird unknown style), so something unambiguous and
  84    straightforward is best.
  85
  86
  87    The IEEE code below is the usual case, it knows either a 32-bit or 64-bit
  88    limb and is done with shifts and masks.  The 64-bit case in particular
  89    should come out nice and compact.
  90
  91    The generic code used to work one bit at a time, which was not only slow,
  92    but implicitly relied upon denorms for intermediates, since the lowest bits'
  93    weight of a perfectly valid fp number underflows in non-denorm.  Therefore,
  94    the generic code now works limb-per-limb, initially creating a number x such
  95    that 1 <= x <= BASE.  (BASE is reached only as result of rounding.)  Then
  96    x's exponent is scaled with explicit code (not ldexp to avoid libm
  97    dependency).  It is a tap-dance to avoid underflow or overflow, beware!
  98
  99
 100    Traps:
 101
 102    Hardware traps for overflow to infinity, underflow to zero, or unsupported
 103    denorms may or may not be taken.  The IEEE code works bitwise and so
 104    probably won't trigger them, the generic code works by float operations and
 105    so probably will.  This difference might be thought less than ideal, but
 106    again its felt straightforward code is better than trying to get intimate
 107    with hardware exceptions (of perhaps unknown nature).
 108
 109
 110    Not done:
 111
 112    mpz_get_d in the past handled size==1 with a cast limb->double.  This might
 113    still be worthwhile there (for up to the mantissa many bits), but for
 114    mpn_get_d here, the cost of applying "exp" to the resulting exponent would
 115    probably use up any benefit a cast may have over bit twiddling.  Also, if
 116    the exponent is pushed into denorm range then bit twiddling is the only
 117    option, to ensure the desired truncation is obtained.
 118
 119
 120    Other:
 121
 122    For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL
 123    to the kernel for values >= 2^63.  This makes it slow, and worse the kernel
 124    Linux (what versions?) apparently uses untested code in its trap handling
 125    routines, and gets the sign wrong.  We don't use such a limb-to-double
 126    cast, neither in the IEEE or generic code.  */
 127
 128
 129
 130 #undef FORMAT_RECOGNIZED
 131
 132 double
 133 mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp)
 134 {
 135   int lshift, nbits;
 136   mp_limb_t x, mhi, mlo;
 137
 138   ASSERT (size >= 0);
 139   ASSERT_MPN (up, size);
 140   ASSERT (size == 0 || up[size-1] != 0);
 141
 142   if (size == 0)
 143     return 0.0;
 144
 145   /* Adjust exp to a radix point just above {up,size}, guarding against
 146      overflow.  After this exp can of course be reduced to anywhere within
 147      the {up,size} region without underflow.  */
 148   if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size)
 149                 > ((unsigned long) LONG_MAX - exp)))
 150     {
 151 #if _GMP_IEEE_FLOATS
 152       goto ieee_infinity;
 153 #endif
 154
 155       /* generic */
 156       exp = LONG_MAX;
 157     }
 158   else
 159     {
 160       exp += GMP_NUMB_BITS * size;
 161     }
 162
 163 #if _GMP_IEEE_FLOATS
 164     {
 165       union ieee_double_extract u;
 166
 167       up += size;
 168
 169 #if GMP_LIMB_BITS == 64
 170       mlo = up[-1];
 171       count_leading_zeros (lshift, mlo);
 172
 173       exp -= (lshift - GMP_NAIL_BITS) + 1;
 174       mlo <<= lshift;
 175
 176       nbits = GMP_LIMB_BITS - lshift;
 177
 178       if (nbits < 53 && size > 1)
 179         {
 180           x = up[-2];
 181           x <<= GMP_NAIL_BITS;
 182           x >>= nbits;
 183           mlo |= x;
 184           nbits += GMP_NUMB_BITS;
 185
 186           if (LIMBS_PER_DOUBLE >= 3 && nbits < 53 && size > 2)
 187             {
 188               x = up[-3];
 189               x <<= GMP_NAIL_BITS;
 190               x >>= nbits;
 191               mlo |= x;
 192               nbits += GMP_NUMB_BITS;
 193             }
 194         }
 195       mhi = mlo >> (32 + 11);
 196       mlo = mlo >> 11;          /* later implicitly truncated to 32 bits */
 197 #endif
 198 #if GMP_LIMB_BITS == 32
 199       x = *--up;
 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 > 1)
 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       /* Now all needed bits in mhi have been accumulated.  Add bits to mlo.  */
 234
 235       if (LIMBS_PER_DOUBLE >= 2 && nbits < 32 && size > 1)
 236         {
 237           x = up[-1];
 238           x <<= GMP_NAIL_BITS;
 239           x >>= nbits;
 240           mlo |= x;
 241           nbits += GMP_NUMB_BITS;
 242
 243           if (LIMBS_PER_DOUBLE >= 3 && nbits < 32 && size > 2)
 244             {
 245               x = up[-2];
 246               x <<= GMP_NAIL_BITS;
 247               x >>= nbits;
 248               mlo |= x;
 249               nbits += GMP_NUMB_BITS;
 250
 251               if (LIMBS_PER_DOUBLE >= 4 && nbits < 32 && size > 3)
 252                 {
 253                   x = up[-3];
 254                   x <<= GMP_NAIL_BITS;
 255                   x >>= nbits;
 256                   mlo |= x;
 257                   nbits += GMP_NUMB_BITS;
 258                 }
 259             }
 260         }
 261
 262     done:;
 263
 264 #endif
 265       if (UNLIKELY (exp >= CONST_1024))
 266         {
 267           /* overflow, return infinity */
 268         ieee_infinity:
 269           mhi = 0;
 270           mlo = 0;
 271           exp = 1024;
 272         }
 273       else if (UNLIKELY (exp <= CONST_NEG_1023))
 274         {
 275           int rshift;
 276
 277           if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
 278             return 0.0;  /* denorm underflows to zero */
 279
 280           rshift = -1022 - exp;
 281           ASSERT (rshift > 0 && rshift < 53);
 282 #if GMP_LIMB_BITS > 53
 283           mlo >>= rshift;
 284           mhi = mlo >> 32;
 285 #else
 286           if (rshift >= 32)
 287             {
 288               mlo = mhi;
 289               mhi = 0;
 290               rshift -= 32;
 291             }
 292           lshift = GMP_LIMB_BITS - rshift;
 293           mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift);
 294           mhi >>= rshift;
 295 #endif
 296           exp = -1023;
 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 #define FORMAT_RECOGNIZED 1
 305 #endif
 306
 307 #if HAVE_DOUBLE_VAX_D
 308     {
 309       union double_extract u;
 310
 311       up += size;
 312
 313       mhi = up[-1];
 314
 315       count_leading_zeros (lshift, mhi);
 316       exp -= lshift;
 317       mhi <<= lshift;
 318
 319       mlo = 0;
 320       if (size > 1)
 321         {
 322           mlo = up[-2];
 323           if (lshift != 0)
 324             mhi += mlo >> (GMP_LIMB_BITS - lshift);
 325           mlo <<= lshift;
 326
 327           if (size > 2 && lshift > 8)
 328             {
 329               x = up[-3];
 330               mlo += x >> (GMP_LIMB_BITS - lshift);
 331             }
 332         }
 333
 334       if (UNLIKELY (exp >= 128))
 335         {
 336           /* overflow, return maximum number */
 337           mhi = 0xffffffff;
 338           mlo = 0xffffffff;
 339           exp = 127;
 340         }
 341       else if (UNLIKELY (exp < -128))
 342         {
 343           return 0.0;    /* underflows to zero */
 344         }
 345
 346       u.s.man3 = mhi >> 24;     /* drop msb, since implicit */
 347       u.s.man2 = mhi >> 8;
 348       u.s.man1 = (mhi << 8) + (mlo >> 24);
 349       u.s.man0 = mlo >> 8;
 350       u.s.exp = exp + 128;
 351       u.s.sig = sign < 0;
 352       return u.d;
 353     }
 354 #define FORMAT_RECOGNIZED 1
 355 #endif
 356
 357 #if ! FORMAT_RECOGNIZED
 358     {      /* Non-IEEE or strange limb size, do something generic. */
 359       mp_size_t i;
 360       double d, weight;
 361       unsigned long uexp;
 362
 363       /* First generate an fp number disregarding exp, instead keeping things
 364          within the numb base factor from 1, which should prevent overflow and
 365          underflow even for the most exponent limited fp formats.  The
 366          termination criteria should be refined, since we now include too many
 367          limbs.  */
 368       weight = 1/MP_BASE_AS_DOUBLE;
 369       d = up[size - 1];
 370       for (i = size - 2; i >= 0; i--)
 371         {
 372           d += up[i] * weight;
 373           weight /= MP_BASE_AS_DOUBLE;
 374           if (weight == 0)
 375             break;
 376         }
 377
 378       /* Now apply exp.  */
 379       exp -= GMP_NUMB_BITS;
 380       if (exp > 0)
 381         {
 382           weight = 2.0;
 383           uexp = exp;
 384         }
 385       else
 386         {
 387           weight = 0.5;
 388           uexp = 1 - (unsigned long) (exp + 1);
 389         }
 390 #if 1
 391       /* Square-and-multiply exponentiation.  */
 392       if (uexp & 1)
 393         d *= weight;
 394       while (uexp >>= 1)
 395         {
 396           weight *= weight;
 397           if (uexp & 1)
 398             d *= weight;
 399         }
 400 #else
 401       /* Plain exponentiation.  */
 402       while (uexp > 0)
 403         {
 404           d *= weight;
 405           uexp--;
 406         }
 407 #endif
 408
 409       return sign >= 0 ? d : -d;
 410     }
 411 #endif
 412 }