3
* ====================================================
4
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6
* Developed at SunPro, a Sun Microsystems, Inc. business.
7
* Permission to use, copy, modify, and distribute this
8
* software is freely granted, provided that this notice
10
* ====================================================
14
* from: @(#)fdlibm.h 5.1 93/09/24
18
#ifndef _NPY_MATH_PRIVATE_H_
19
#define _NPY_MATH_PRIVATE_H_
24
#include "npy_config.h"
25
#include "npy_fpmath.h"
27
#include "numpy/npy_math.h"
28
#include "numpy/npy_cpu.h"
29
#include "numpy/npy_endian.h"
30
#include "numpy/npy_common.h"
33
* The original fdlibm code used statements like:
34
* n0 = ((*(int*)&one)>>29)^1; * index of high word *
35
* ix0 = *(n0+(int*)&x); * high word of x *
36
* ix1 = *((1-n0)+(int*)&x); * low word of x *
37
* to dig two 32 bit words out of the 64 bit IEEE floating point
38
* value. That is non-ANSI, and, moreover, the gcc instruction
39
* scheduler gets it wrong. We instead use the following macros.
40
* Unlike the original code, we determine the endianness at compile
41
* time, not at run time; I don't see much benefit to selecting
42
* endianness at run time.
46
* A union which permits us to convert between a double and two 32 bit
50
/* XXX: not really, but we already make this assumption elsewhere. Will have to
51
* fix this at some point */
52
#define IEEE_WORD_ORDER NPY_BYTE_ORDER
54
#if IEEE_WORD_ORDER == NPY_BIG_ENDIAN
64
} ieee_double_shape_type;
68
#if IEEE_WORD_ORDER == NPY_LITTLE_ENDIAN
78
} ieee_double_shape_type;
82
/* Get two 32 bit ints from a double. */
84
#define EXTRACT_WORDS(ix0,ix1,d) \
86
ieee_double_shape_type ew_u; \
88
(ix0) = ew_u.parts.msw; \
89
(ix1) = ew_u.parts.lsw; \
92
/* Get the more significant 32 bit int from a double. */
94
#define GET_HIGH_WORD(i,d) \
96
ieee_double_shape_type gh_u; \
98
(i) = gh_u.parts.msw; \
101
/* Get the less significant 32 bit int from a double. */
103
#define GET_LOW_WORD(i,d) \
105
ieee_double_shape_type gl_u; \
107
(i) = gl_u.parts.lsw; \
110
/* Set the more significant 32 bits of a double from an int. */
112
#define SET_HIGH_WORD(d,v) \
114
ieee_double_shape_type sh_u; \
116
sh_u.parts.msw = (v); \
120
/* Set the less significant 32 bits of a double from an int. */
122
#define SET_LOW_WORD(d,v) \
124
ieee_double_shape_type sl_u; \
126
sl_u.parts.lsw = (v); \
130
/* Set a double from two 32 bit ints. */
132
#define INSERT_WORDS(d,ix0,ix1) \
134
ieee_double_shape_type iw_u; \
135
iw_u.parts.msw = (ix0); \
136
iw_u.parts.lsw = (ix1); \
141
* A union which permits us to convert between a float and a 32 bit
148
/* FIXME: Assumes 32 bit int. */
150
} ieee_float_shape_type;
152
/* Get a 32 bit int from a float. */
154
#define GET_FLOAT_WORD(i,d) \
156
ieee_float_shape_type gf_u; \
161
/* Set a float from a 32 bit int. */
163
#define SET_FLOAT_WORD(d,i) \
165
ieee_float_shape_type sf_u; \
170
#ifdef NPY_USE_C99_COMPLEX
175
* Long double support
177
#if defined(HAVE_LDOUBLE_INTEL_EXTENDED_12_BYTES_LE)
179
* Intel extended 80 bits precision. Bit representation is
180
* | junk | s |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm|
181
* | 16 bits| 1 bit | 15 bits | 64 bits |
182
* | a[2] | a[1] | a[0] |
184
* 16 stronger bits of a[2] are junk
186
typedef npy_uint32 IEEEl2bitsrep_part;
190
union IEEEl2bitsrep {
192
IEEEl2bitsrep_part a[3];
195
#define LDBL_MANL_INDEX 0
196
#define LDBL_MANL_MASK 0xFFFFFFFF
197
#define LDBL_MANL_SHIFT 0
199
#define LDBL_MANH_INDEX 1
200
#define LDBL_MANH_MASK 0xFFFFFFFF
201
#define LDBL_MANH_SHIFT 0
203
#define LDBL_EXP_INDEX 2
204
#define LDBL_EXP_MASK 0x7FFF
205
#define LDBL_EXP_SHIFT 0
207
#define LDBL_SIGN_INDEX 2
208
#define LDBL_SIGN_MASK 0x8000
209
#define LDBL_SIGN_SHIFT 15
211
#define LDBL_NBIT 0x80000000
213
typedef npy_uint32 ldouble_man_t;
214
typedef npy_uint32 ldouble_exp_t;
215
typedef npy_uint32 ldouble_sign_t;
216
#elif defined(HAVE_LDOUBLE_INTEL_EXTENDED_16_BYTES_LE)
218
* Intel extended 80 bits precision, 16 bytes alignment.. Bit representation is
219
* | junk | s |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm|
220
* | 16 bits| 1 bit | 15 bits | 64 bits |
221
* | a[2] | a[1] | a[0] |
223
* a[3] and 16 stronger bits of a[2] are junk
225
typedef npy_uint32 IEEEl2bitsrep_part;
227
union IEEEl2bitsrep {
229
IEEEl2bitsrep_part a[4];
232
#define LDBL_MANL_INDEX 0
233
#define LDBL_MANL_MASK 0xFFFFFFFF
234
#define LDBL_MANL_SHIFT 0
236
#define LDBL_MANH_INDEX 1
237
#define LDBL_MANH_MASK 0xFFFFFFFF
238
#define LDBL_MANH_SHIFT 0
240
#define LDBL_EXP_INDEX 2
241
#define LDBL_EXP_MASK 0x7FFF
242
#define LDBL_EXP_SHIFT 0
244
#define LDBL_SIGN_INDEX 2
245
#define LDBL_SIGN_MASK 0x8000
246
#define LDBL_SIGN_SHIFT 15
248
#define LDBL_NBIT 0x800000000
250
typedef npy_uint32 ldouble_man_t;
251
typedef npy_uint32 ldouble_exp_t;
252
typedef npy_uint32 ldouble_sign_t;
253
#elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_16_BYTES_BE) || \
254
defined(HAVE_LDOUBLE_IEEE_DOUBLE_BE)
255
/* 64 bits IEEE double precision aligned on 16 bytes: used by ppc arch on
259
* IEEE double precision. Bit representation is
260
* | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
261
* |1 bit| 11 bits | 52 bits |
264
typedef npy_uint32 IEEEl2bitsrep_part;
266
union IEEEl2bitsrep {
268
IEEEl2bitsrep_part a[2];
271
#define LDBL_MANL_INDEX 1
272
#define LDBL_MANL_MASK 0xFFFFFFFF
273
#define LDBL_MANL_SHIFT 0
275
#define LDBL_MANH_INDEX 0
276
#define LDBL_MANH_MASK 0x000FFFFF
277
#define LDBL_MANH_SHIFT 0
279
#define LDBL_EXP_INDEX 0
280
#define LDBL_EXP_MASK 0x7FF00000
281
#define LDBL_EXP_SHIFT 20
283
#define LDBL_SIGN_INDEX 0
284
#define LDBL_SIGN_MASK 0x80000000
285
#define LDBL_SIGN_SHIFT 31
289
typedef npy_uint32 ldouble_man_t;
290
typedef npy_uint32 ldouble_exp_t;
291
typedef npy_uint32 ldouble_sign_t;
292
#elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_LE)
293
/* 64 bits IEEE double precision, Little Endian. */
296
* IEEE double precision. Bit representation is
297
* | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
298
* |1 bit| 11 bits | 52 bits |
301
typedef npy_uint32 IEEEl2bitsrep_part;
303
union IEEEl2bitsrep {
305
IEEEl2bitsrep_part a[2];
308
#define LDBL_MANL_INDEX 0
309
#define LDBL_MANL_MASK 0xFFFFFFFF
310
#define LDBL_MANL_SHIFT 0
312
#define LDBL_MANH_INDEX 1
313
#define LDBL_MANH_MASK 0x000FFFFF
314
#define LDBL_MANH_SHIFT 0
316
#define LDBL_EXP_INDEX 1
317
#define LDBL_EXP_MASK 0x7FF00000
318
#define LDBL_EXP_SHIFT 20
320
#define LDBL_SIGN_INDEX 1
321
#define LDBL_SIGN_MASK 0x80000000
322
#define LDBL_SIGN_SHIFT 31
324
#define LDBL_NBIT 0x00000080
326
typedef npy_uint32 ldouble_man_t;
327
typedef npy_uint32 ldouble_exp_t;
328
typedef npy_uint32 ldouble_sign_t;
329
#elif defined(HAVE_LDOUBLE_IEEE_QUAD_BE)
331
* IEEE quad precision, Big Endian. Bit representation is
332
* | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
333
* |1 bit| 15 bits | 112 bits |
336
typedef npy_uint64 IEEEl2bitsrep_part;
338
union IEEEl2bitsrep {
340
IEEEl2bitsrep_part a[2];
343
#define LDBL_MANL_INDEX 1
344
#define LDBL_MANL_MASK 0xFFFFFFFFFFFFFFFF
345
#define LDBL_MANL_SHIFT 0
347
#define LDBL_MANH_INDEX 0
348
#define LDBL_MANH_MASK 0x0000FFFFFFFFFFFF
349
#define LDBL_MANH_SHIFT 0
351
#define LDBL_EXP_INDEX 0
352
#define LDBL_EXP_MASK 0x7FFF000000000000
353
#define LDBL_EXP_SHIFT 48
355
#define LDBL_SIGN_INDEX 0
356
#define LDBL_SIGN_MASK 0x8000000000000000
357
#define LDBL_SIGN_SHIFT 63
361
typedef npy_uint64 ldouble_man_t;
362
typedef npy_uint64 ldouble_exp_t;
363
typedef npy_uint32 ldouble_sign_t;
366
/* Get the sign bit of x. x should be of type IEEEl2bitsrep */
367
#define GET_LDOUBLE_SIGN(x) \
368
(((x).a[LDBL_SIGN_INDEX] & LDBL_SIGN_MASK) >> LDBL_SIGN_SHIFT)
370
/* Set the sign bit of x to v. x should be of type IEEEl2bitsrep */
371
#define SET_LDOUBLE_SIGN(x, v) \
372
((x).a[LDBL_SIGN_INDEX] = \
373
((x).a[LDBL_SIGN_INDEX] & ~LDBL_SIGN_MASK) | \
374
(((IEEEl2bitsrep_part)(v) << LDBL_SIGN_SHIFT) & LDBL_SIGN_MASK))
376
/* Get the exp bits of x. x should be of type IEEEl2bitsrep */
377
#define GET_LDOUBLE_EXP(x) \
378
(((x).a[LDBL_EXP_INDEX] & LDBL_EXP_MASK) >> LDBL_EXP_SHIFT)
380
/* Set the exp bit of x to v. x should be of type IEEEl2bitsrep */
381
#define SET_LDOUBLE_EXP(x, v) \
382
((x).a[LDBL_EXP_INDEX] = \
383
((x).a[LDBL_EXP_INDEX] & ~LDBL_EXP_MASK) | \
384
(((IEEEl2bitsrep_part)(v) << LDBL_EXP_SHIFT) & LDBL_EXP_MASK))
386
/* Get the manl bits of x. x should be of type IEEEl2bitsrep */
387
#define GET_LDOUBLE_MANL(x) \
388
(((x).a[LDBL_MANL_INDEX] & LDBL_MANL_MASK) >> LDBL_MANL_SHIFT)
390
/* Set the manl bit of x to v. x should be of type IEEEl2bitsrep */
391
#define SET_LDOUBLE_MANL(x, v) \
392
((x).a[LDBL_MANL_INDEX] = \
393
((x).a[LDBL_MANL_INDEX] & ~LDBL_MANL_MASK) | \
394
(((IEEEl2bitsrep_part)(v) << LDBL_MANL_SHIFT) & LDBL_MANL_MASK))
396
/* Get the manh bits of x. x should be of type IEEEl2bitsrep */
397
#define GET_LDOUBLE_MANH(x) \
398
(((x).a[LDBL_MANH_INDEX] & LDBL_MANH_MASK) >> LDBL_MANH_SHIFT)
400
/* Set the manh bit of x to v. x should be of type IEEEl2bitsrep */
401
#define SET_LDOUBLE_MANH(x, v) \
402
((x).a[LDBL_MANH_INDEX] = \
403
((x).a[LDBL_MANH_INDEX] & ~LDBL_MANH_MASK) | \
404
(((IEEEl2bitsrep_part)(v) << LDBL_MANH_SHIFT) & LDBL_MANH_MASK))
407
* Those unions are used to convert a pointer of npy_cdouble to native C99
408
* complex or our own complex type independently on whether C99 complex
409
* support is available
411
#ifdef NPY_USE_C99_COMPLEX
414
complex double c99_z;
415
} __npy_cdouble_to_c99_cast;
420
} __npy_cfloat_to_c99_cast;
423
npy_clongdouble npy_z;
424
complex long double c99_z;
425
} __npy_clongdouble_to_c99_cast;
430
} __npy_cdouble_to_c99_cast;
435
} __npy_cfloat_to_c99_cast;
438
npy_clongdouble npy_z;
439
npy_clongdouble c99_z;
440
} __npy_clongdouble_to_c99_cast;
443
#endif /* !_NPY_MATH_PRIVATE_H_ */