aboutsummaryrefslogtreecommitdiff
path: root/vendor/gmp-6.3.0/mpn/generic/get_d.c
blob: 8bef128108ae3d2d3b1bcc0a4e7a1f825709b4c0 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
/* mpn_get_d -- limbs to double conversion.

   THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
   CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
   FUTURE GNU MP RELEASES.

Copyright 2003, 2004, 2007, 2009, 2010, 2012, 2018 Free Software
Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 2 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library.  If not,
see https://www.gnu.org/licenses/.  */

#include "config.h"

#if HAVE_FLOAT_H
#include <float.h>  /* for DBL_MANT_DIG and FLT_RADIX */
#endif

#include "gmp-impl.h"
#include "longlong.h"

#ifndef _GMP_IEEE_FLOATS
#define _GMP_IEEE_FLOATS 0
#endif

/* To force use of the generic C code for testing, put
   "#define _GMP_IEEE_FLOATS 0" at this point.  */


/* In alpha gcc prior to 3.4, signed DI comparisons involving constants are
   rearranged from "x < n" to "x+(-n) < 0", which is of course hopelessly
   wrong if that addition overflows.

   The workaround here avoids this bug by ensuring n is not a literal constant.
   Note that this is alpha specific.  The offending transformation is/was in
   alpha.c alpha_emit_conditional_branch() under "We want to use cmpcc/bcc".

   Bizarrely, this happens also with Cray cc on alphaev5-cray-unicosmk2.0.6.X,
   and has the same solution.  Don't know why or how.  */

#if HAVE_HOST_CPU_FAMILY_alpha				\
  && ((defined (__GNUC__) && ! __GMP_GNUC_PREREQ(3,4))	\
      || defined (_CRAY))
static volatile const long CONST_1024 = 1024;
static volatile const long CONST_NEG_1023 = -1023;
static volatile const long CONST_NEG_1022_SUB_53 = -1022 - 53;
#else
#define CONST_1024	      (1024)
#define CONST_NEG_1023	      (-1023)
#define CONST_NEG_1022_SUB_53 (-1022 - 53)
#endif


/* Return the value {ptr,size}*2^exp, and negative if sign<0.  Must have
   size>=1, and a non-zero high limb ptr[size-1].

   When we know the fp format, the result is truncated towards zero.  This is
   consistent with other gmp conversions, like mpz_set_f or mpz_set_q, and is
   easy to implement and test.

   When we do not know the format, such truncation seems much harder.  One
   would need to defeat any rounding mode, including round-up.

   It's felt that GMP is not primarily concerned with hardware floats, and
   really isn't enhanced by getting involved with hardware rounding modes
   (which could even be some weird unknown style), so something unambiguous and
   straightforward is best.


   The IEEE code below is the usual case, it knows either a 32-bit or 64-bit
   limb and is done with shifts and masks.  The 64-bit case in particular
   should come out nice and compact.

   The generic code used to work one bit at a time, which was not only slow,
   but implicitly relied upon denorms for intermediates, since the lowest bits'
   weight of a perfectly valid fp number underflows in non-denorm.  Therefore,
   the generic code now works limb-per-limb, initially creating a number x such
   that 1 <= x <= BASE.  (BASE is reached only as result of rounding.)  Then
   x's exponent is scaled with explicit code (not ldexp to avoid libm
   dependency).  It is a tap-dance to avoid underflow or overflow, beware!


   Traps:

   Hardware traps for overflow to infinity, underflow to zero, or unsupported
   denorms may or may not be taken.  The IEEE code works bitwise and so
   probably won't trigger them, the generic code works by float operations and
   so probably will.  This difference might be thought less than ideal, but
   again its felt straightforward code is better than trying to get intimate
   with hardware exceptions (of perhaps unknown nature).


   Not done:

   mpz_get_d in the past handled size==1 with a cast limb->double.  This might
   still be worthwhile there (for up to the mantissa many bits), but for
   mpn_get_d here, the cost of applying "exp" to the resulting exponent would
   probably use up any benefit a cast may have over bit twiddling.  Also, if
   the exponent is pushed into denorm range then bit twiddling is the only
   option, to ensure the desired truncation is obtained.


   Other:

   For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL
   to the kernel for values >= 2^63.  This makes it slow, and worse the kernel
   Linux (what versions?) apparently uses untested code in its trap handling
   routines, and gets the sign wrong.  We don't use such a limb-to-double
   cast, neither in the IEEE or generic code.  */



#undef FORMAT_RECOGNIZED

double
mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp)
{
  int lshift, nbits;
  mp_limb_t x, mhi, mlo;

  ASSERT (size >= 0);
  ASSERT_MPN (up, size);
  ASSERT (size == 0 || up[size-1] != 0);

  if (size == 0)
    return 0.0;

  /* Adjust exp to a radix point just above {up,size}, guarding against
     overflow.  After this exp can of course be reduced to anywhere within
     the {up,size} region without underflow.  */
  if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size)
		> ((unsigned long) LONG_MAX - exp)))
    {
#if _GMP_IEEE_FLOATS
      goto ieee_infinity;
#endif

      /* generic */
      exp = LONG_MAX;
    }
  else
    {
      exp += GMP_NUMB_BITS * size;
    }

#if _GMP_IEEE_FLOATS
    {
      union ieee_double_extract u;

      up += size;

#if GMP_LIMB_BITS == 64
      mlo = up[-1];
      count_leading_zeros (lshift, mlo);

      exp -= (lshift - GMP_NAIL_BITS) + 1;
      mlo <<= lshift;

      nbits = GMP_LIMB_BITS - lshift;

      if (nbits < 53 && size > 1)
	{
	  x = up[-2];
	  x <<= GMP_NAIL_BITS;
	  x >>= nbits;
	  mlo |= x;
	  nbits += GMP_NUMB_BITS;

	  if (LIMBS_PER_DOUBLE >= 3 && nbits < 53 && size > 2)
	    {
	      x = up[-3];
	      x <<= GMP_NAIL_BITS;
	      x >>= nbits;
	      mlo |= x;
	      nbits += GMP_NUMB_BITS;
	    }
	}
      mhi = mlo >> (32 + 11);
      mlo = mlo >> 11;		/* later implicitly truncated to 32 bits */
#endif
#if GMP_LIMB_BITS == 32
      x = *--up;
      count_leading_zeros (lshift, x);

      exp -= (lshift - GMP_NAIL_BITS) + 1;
      x <<= lshift;
      mhi = x >> 11;

      if (lshift < 11)		/* FIXME: never true if NUMB < 20 bits */
	{
	  /* All 20 bits in mhi */
	  mlo = x << 21;
	  /* >= 1 bit in mlo */
	  nbits = GMP_LIMB_BITS - lshift - 21;
	}
      else
	{
	  if (size > 1)
	    {
	      nbits = GMP_LIMB_BITS - lshift;

	      x = *--up, size--;
	      x <<= GMP_NAIL_BITS;
	      mhi |= x >> nbits >> 11;

	      mlo = x << (GMP_LIMB_BITS - nbits - 11);
	      nbits = nbits + 11 - GMP_NAIL_BITS;
	    }
	  else
	    {
	      mlo = 0;
	      goto done;
	    }
	}

      /* Now all needed bits in mhi have been accumulated.  Add bits to mlo.  */

      if (LIMBS_PER_DOUBLE >= 2 && nbits < 32 && size > 1)
	{
	  x = up[-1];
	  x <<= GMP_NAIL_BITS;
	  x >>= nbits;
	  mlo |= x;
	  nbits += GMP_NUMB_BITS;

	  if (LIMBS_PER_DOUBLE >= 3 && nbits < 32 && size > 2)
	    {
	      x = up[-2];
	      x <<= GMP_NAIL_BITS;
	      x >>= nbits;
	      mlo |= x;
	      nbits += GMP_NUMB_BITS;

	      if (LIMBS_PER_DOUBLE >= 4 && nbits < 32 && size > 3)
		{
		  x = up[-3];
		  x <<= GMP_NAIL_BITS;
		  x >>= nbits;
		  mlo |= x;
		  nbits += GMP_NUMB_BITS;
		}
	    }
	}

    done:;

#endif
      if (UNLIKELY (exp >= CONST_1024))
	{
	  /* overflow, return infinity */
	ieee_infinity:
	  mhi = 0;
	  mlo = 0;
	  exp = 1024;
	}
      else if (UNLIKELY (exp <= CONST_NEG_1023))
	{
	  int rshift;

	  if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
	    return 0.0;	 /* denorm underflows to zero */

	  rshift = -1022 - exp;
	  ASSERT (rshift > 0 && rshift < 53);
#if GMP_LIMB_BITS > 53
	  mlo >>= rshift;
	  mhi = mlo >> 32;
#else
	  if (rshift >= 32)
	    {
	      mlo = mhi;
	      mhi = 0;
	      rshift -= 32;
	    }
	  lshift = GMP_LIMB_BITS - rshift;
	  mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift);
	  mhi >>= rshift;
#endif
	  exp = -1023;
	}
      u.s.manh = mhi;
      u.s.manl = mlo;
      u.s.exp = exp + 1023;
      u.s.sig = (sign < 0);
      return u.d;
    }
#define FORMAT_RECOGNIZED 1
#endif

#if HAVE_DOUBLE_VAX_D
    {
      union double_extract u;

      up += size;

      mhi = up[-1];

      count_leading_zeros (lshift, mhi);
      exp -= lshift;
      mhi <<= lshift;

      mlo = 0;
      if (size > 1)
	{
	  mlo = up[-2];
	  if (lshift != 0)
	    mhi += mlo >> (GMP_LIMB_BITS - lshift);
	  mlo <<= lshift;

	  if (size > 2 && lshift > 8)
	    {
	      x = up[-3];
	      mlo += x >> (GMP_LIMB_BITS - lshift);
	    }
	}

      if (UNLIKELY (exp >= 128))
	{
	  /* overflow, return maximum number */
	  mhi = 0xffffffff;
	  mlo = 0xffffffff;
	  exp = 127;
	}
      else if (UNLIKELY (exp < -128))
	{
	  return 0.0;	 /* underflows to zero */
	}

      u.s.man3 = mhi >> 24;	/* drop msb, since implicit */
      u.s.man2 = mhi >> 8;
      u.s.man1 = (mhi << 8) + (mlo >> 24);
      u.s.man0 = mlo >> 8;
      u.s.exp = exp + 128;
      u.s.sig = sign < 0;
      return u.d;
    }
#define FORMAT_RECOGNIZED 1
#endif

#if ! FORMAT_RECOGNIZED

#if !defined(GMP_DBL_MANT_BITS)
#if defined(DBL_MANT_DIG) && FLT_RADIX == 2
#define GMP_DBL_MANT_BITS DBL_MANT_DIG
#else
/* FIXME: Chose a smarter default value. */
#define GMP_DBL_MANT_BITS (16 * sizeof (double))
#endif
#endif

    { /* Non-IEEE or strange limb size, generically convert
	 GMP_DBL_MANT_BITS bits. */
      mp_limb_t l;
      int m;
      mp_size_t i;
      double d, weight;
      unsigned long uexp;

      /* First generate an fp number disregarding exp, instead keeping things
	 within the numb base factor from 1, which should prevent overflow and
	 underflow even for the most exponent limited fp formats.  */
      i = size - 1;
      l = up[i];
      count_leading_zeros (m, l);
      m = m + GMP_DBL_MANT_BITS - GMP_LIMB_BITS;
      if (m < 0)
	l &= GMP_NUMB_MAX << -m;
      d = l;
      for (weight = 1/MP_BASE_AS_DOUBLE; m > 0 && --i >= 0;)
	{
	  l = up[i];
	  m -= GMP_NUMB_BITS;
	  if (m < 0)
	    l &= GMP_NUMB_MAX << -m;
	  d += l * weight;
	  weight /= MP_BASE_AS_DOUBLE;
	  if (weight == 0)
	    break;
	}

      /* Now apply exp.  */
      exp -= GMP_NUMB_BITS;
      if (exp > 0)
	{
	  weight = 2.0;
	  uexp = exp;
	}
      else
	{
	  weight = 0.5;
	  uexp = NEG_CAST (unsigned long, exp);
	}
#if 1
      /* Square-and-multiply exponentiation.  */
      if (uexp & 1)
	d *= weight;
      while (uexp >>= 1)
	{
	  weight *= weight;
	  if (uexp & 1)
	    d *= weight;
	}
#else
      /* Plain exponentiation.  */
      while (uexp > 0)
	{
	  d *= weight;
	  uexp--;
	}
#endif

      return sign >= 0 ? d : -d;
    }
#endif
}