CoCalc -- extract-dbl.c

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Path: gap4r8 / extern / gmp-6.0.0a / extract-dbl.c
Views: ⁴¹⁸³⁴⁶
1
/* __gmp_extract_double -- convert from double to array of mp_limb_t.
2

3
Copyright 1996, 1999-2002, 2006, 2012 Free Software Foundation, Inc.
4

5
This file is part of the GNU MP Library.
6

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

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

14
or
15

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

20
or both in parallel, as here.
21

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

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

31
#include "gmp.h"
32
#include "gmp-impl.h"
33

34
#ifdef XDEBUG
35
#undef _GMP_IEEE_FLOATS
36
#endif
37

38
#ifndef _GMP_IEEE_FLOATS
39
#define _GMP_IEEE_FLOATS 0
40
#endif
41

42
/* Extract a non-negative double in d.  */
43

44
int
45
__gmp_extract_double (mp_ptr rp, double d)
46
{
47
  long exp;
48
  unsigned sc;
49
#ifdef _LONG_LONG_LIMB
50
#define BITS_PER_PART 64	/* somewhat bogus */
51
  unsigned long long int manl;
52
#else
53
#define BITS_PER_PART GMP_LIMB_BITS
54
  unsigned long int manh, manl;
55
#endif
56

57
  /* BUGS
58

59
     1. Should handle Inf and NaN in IEEE specific code.
60
     2. Handle Inf and NaN also in default code, to avoid hangs.
61
     3. Generalize to handle all GMP_LIMB_BITS >= 32.
62
     4. This lits is incomplete and misspelled.
63
   */
64

65
  ASSERT (d >= 0.0);
66

67
  if (d == 0.0)
68
    {
69
      MPN_ZERO (rp, LIMBS_PER_DOUBLE);
70
      return 0;
71
    }
72

73
#if _GMP_IEEE_FLOATS
74
  {
75
#if defined (__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 8
76
    /* Work around alpha-specific bug in GCC 2.8.x.  */
77
    volatile
78
#endif
79
    union ieee_double_extract x;
80
    x.d = d;
81
    exp = x.s.exp;
82
#if BITS_PER_PART == 64		/* generalize this to BITS_PER_PART > BITS_IN_MANTISSA */
83
    manl = (((mp_limb_t) 1 << 63)
84
	    | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
85
    if (exp == 0)
86
      {
87
	/* Denormalized number.  Don't try to be clever about this,
88
	   since it is not an important case to make fast.  */
89
	exp = 1;
90
	do
91
	  {
92
	    manl = manl << 1;
93
	    exp--;
94
	  }
95
	while ((manl & GMP_LIMB_HIGHBIT) == 0);
96
      }
97
#endif
98
#if BITS_PER_PART == 32
99
    manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
100
    manl = x.s.manl << 11;
101
    if (exp == 0)
102
      {
103
	/* Denormalized number.  Don't try to be clever about this,
104
	   since it is not an important case to make fast.  */
105
	exp = 1;
106
	do
107
	  {
108
	    manh = (manh << 1) | (manl >> 31);
109
	    manl = manl << 1;
110
	    exp--;
111
	  }
112
	while ((manh & GMP_LIMB_HIGHBIT) == 0);
113
      }
114
#endif
115
#if BITS_PER_PART != 32 && BITS_PER_PART != 64
116
  You need to generalize the code above to handle this.
117
#endif
118
    exp -= 1022;		/* Remove IEEE bias.  */
119
  }
120
#else
121
  {
122
    /* Unknown (or known to be non-IEEE) double format.  */
123
    exp = 0;
124
    if (d >= 1.0)
125
      {
126
	ASSERT_ALWAYS (d * 0.5 != d);
127

128
	while (d >= 32768.0)
129
	  {
130
	    d *= (1.0 / 65536.0);
131
	    exp += 16;
132
	  }
133
	while (d >= 1.0)
134
	  {
135
	    d *= 0.5;
136
	    exp += 1;
137
	  }
138
      }
139
    else if (d < 0.5)
140
      {
141
	while (d < (1.0 / 65536.0))
142
	  {
143
	    d *=  65536.0;
144
	    exp -= 16;
145
	  }
146
	while (d < 0.5)
147
	  {
148
	    d *= 2.0;
149
	    exp -= 1;
150
	  }
151
      }
152

153
    d *= (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
154
#if BITS_PER_PART == 64
155
    manl = d;
156
#endif
157
#if BITS_PER_PART == 32
158
    manh = d;
159
    manl = (d - manh) * (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
160
#endif
161
  }
162
#endif /* IEEE */
163

164
  sc = (unsigned) (exp + 64 * GMP_NUMB_BITS) % GMP_NUMB_BITS;
165

166
  /* We add something here to get rounding right.  */
167
  exp = (exp + 64 * GMP_NUMB_BITS) / GMP_NUMB_BITS - 64 * GMP_NUMB_BITS / GMP_NUMB_BITS + 1;
168

169
#if BITS_PER_PART == 64 && LIMBS_PER_DOUBLE == 2
170
#if GMP_NAIL_BITS == 0
171
  if (sc != 0)
172
    {
173
      rp[1] = manl >> (GMP_LIMB_BITS - sc);
174
      rp[0] = manl << sc;
175
    }
176
  else
177
    {
178
      rp[1] = manl;
179
      rp[0] = 0;
180
      exp--;
181
    }
182
#else
183
  if (sc > GMP_NAIL_BITS)
184
    {
185
      rp[1] = manl >> (GMP_LIMB_BITS - sc);
186
      rp[0] = (manl << (sc - GMP_NAIL_BITS)) & GMP_NUMB_MASK;
187
    }
188
  else
189
    {
190
      if (sc == 0)
191
	{
192
	  rp[1] = manl >> GMP_NAIL_BITS;
193
	  rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
194
	  exp--;
195
	}
196
      else
197
	{
198
	  rp[1] = manl >> (GMP_LIMB_BITS - sc);
199
	  rp[0] = (manl >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
200
	}
201
    }
202
#endif
203
#endif
204

205
#if BITS_PER_PART == 64 && LIMBS_PER_DOUBLE == 3
206
  if (sc > GMP_NAIL_BITS)
207
    {
208
      rp[2] = manl >> (GMP_LIMB_BITS - sc);
209
      rp[1] = (manl << sc - GMP_NAIL_BITS) & GMP_NUMB_MASK;
210
      if (sc >= 2 * GMP_NAIL_BITS)
211
	rp[0] = 0;
212
      else
213
	rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS + sc) & GMP_NUMB_MASK;
214
    }
215
  else
216
    {
217
      if (sc == 0)
218
	{
219
	  rp[2] = manl >> GMP_NAIL_BITS;
220
	  rp[1] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
221
	  rp[0] = 0;
222
	  exp--;
223
	}
224
      else
225
	{
226
	  rp[2] = manl >> (GMP_LIMB_BITS - sc);
227
	  rp[1] = (manl >> GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
228
	  rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS + sc) & GMP_NUMB_MASK;
229
	}
230
    }
231
#endif
232

233
#if BITS_PER_PART == 32 && LIMBS_PER_DOUBLE == 3
234
#if GMP_NAIL_BITS == 0
235
  if (sc != 0)
236
    {
237
      rp[2] = manh >> (GMP_LIMB_BITS - sc);
238
      rp[1] = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
239
      rp[0] = manl << sc;
240
    }
241
  else
242
    {
243
      rp[2] = manh;
244
      rp[1] = manl;
245
      rp[0] = 0;
246
      exp--;
247
    }
248
#else
249
  if (sc > GMP_NAIL_BITS)
250
    {
251
      rp[2] = (manh >> (GMP_LIMB_BITS - sc));
252
      rp[1] = ((manh << (sc - GMP_NAIL_BITS)) |
253
	       (manl >> (GMP_LIMB_BITS - sc + GMP_NAIL_BITS))) & GMP_NUMB_MASK;
254
      if (sc >= 2 * GMP_NAIL_BITS)
255
	rp[0] = (manl << sc - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
256
      else
257
	rp[0] = manl >> (2 * GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
258
    }
259
  else
260
    {
261
      if (sc == 0)
262
	{
263
	  rp[2] = manh >> GMP_NAIL_BITS;
264
	  rp[1] = ((manh << GMP_NUMB_BITS - GMP_NAIL_BITS) | (manl >> 2 * GMP_NAIL_BITS)) & GMP_NUMB_MASK;
265
	  rp[0] = (manl << GMP_NUMB_BITS - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
266
	  exp--;
267
	}
268
      else
269
	{
270
	  rp[2] = (manh >> (GMP_LIMB_BITS - sc));
271
	  rp[1] = (manh >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
272
	  rp[0] = ((manh << (GMP_NUMB_BITS - GMP_NAIL_BITS + sc))
273
		   | (manl >> (GMP_LIMB_BITS - (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)))) & GMP_NUMB_MASK;
274
	}
275
    }
276
#endif
277
#endif
278

279
#if BITS_PER_PART == 32 && LIMBS_PER_DOUBLE > 3
280
  if (sc == 0)
281
    {
282
      int i;
283

284
      for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--)
285
	{
286
	  rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS);
287
	  manh = ((manh << GMP_NUMB_BITS)
288
		  | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS)));
289
	  manl = manl << GMP_NUMB_BITS;
290
	}
291
      exp--;
292
    }
293
  else
294
    {
295
      int i;
296

297
      rp[LIMBS_PER_DOUBLE - 1] = (manh >> (GMP_LIMB_BITS - sc));
298
      manh = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
299
      manl = (manl << sc);
300
      for (i = LIMBS_PER_DOUBLE - 2; i >= 0; i--)
301
	{
302
	  rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS);
303
	  manh = ((manh << GMP_NUMB_BITS)
304
		  | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS)));
305
	  manl = manl << GMP_NUMB_BITS;
306
	}
307
  }
308
#endif
309

310
  return exp;
311
}
312

313
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

Product

Resources

Company

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more, all in one place.

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
