numpy  2.0.0
src/npymath/npy_math_private.h
Go to the documentation of this file.
00001 /*
00002  *
00003  * ====================================================
00004  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00005  *
00006  * Developed at SunPro, a Sun Microsystems, Inc. business.
00007  * Permission to use, copy, modify, and distribute this
00008  * software is freely granted, provided that this notice
00009  * is preserved.
00010  * ====================================================
00011  */
00012 
00013 /*
00014  * from: @(#)fdlibm.h 5.1 93/09/24
00015  * $FreeBSD$
00016  */
00017 
00018 #ifndef _NPY_MATH_PRIVATE_H_
00019 #define _NPY_MATH_PRIVATE_H_
00020 
00021 #include <Python.h>
00022 #include <math.h>
00023 
00024 #include "npy_config.h"
00025 #include "npy_fpmath.h"
00026 
00027 #include "numpy/npy_math.h"
00028 #include "numpy/npy_cpu.h"
00029 #include "numpy/npy_endian.h"
00030 #include "numpy/npy_common.h"
00031 
00032 /*
00033  * The original fdlibm code used statements like:
00034  *      n0 = ((*(int*)&one)>>29)^1;             * index of high word *
00035  *      ix0 = *(n0+(int*)&x);                   * high word of x *
00036  *      ix1 = *((1-n0)+(int*)&x);               * low word of x *
00037  * to dig two 32 bit words out of the 64 bit IEEE floating point
00038  * value.  That is non-ANSI, and, moreover, the gcc instruction
00039  * scheduler gets it wrong.  We instead use the following macros.
00040  * Unlike the original code, we determine the endianness at compile
00041  * time, not at run time; I don't see much benefit to selecting
00042  * endianness at run time.
00043  */
00044 
00045 /*
00046  * A union which permits us to convert between a double and two 32 bit
00047  * ints.
00048  */
00049 
00050 /* XXX: not really, but we already make this assumption elsewhere. Will have to
00051  * fix this at some point */
00052 #define IEEE_WORD_ORDER NPY_BYTE_ORDER
00053 
00054 #if IEEE_WORD_ORDER == NPY_BIG_ENDIAN
00055 
00056 typedef union
00057 {
00058   double value;
00059   struct
00060   {
00061     npy_uint32 msw;
00062     npy_uint32 lsw;
00063   } parts;
00064 } ieee_double_shape_type;
00065 
00066 #endif
00067 
00068 #if IEEE_WORD_ORDER == NPY_LITTLE_ENDIAN
00069 
00070 typedef union
00071 {
00072   double value;
00073   struct
00074   {
00075     npy_uint32 lsw;
00076     npy_uint32 msw;
00077   } parts;
00078 } ieee_double_shape_type;
00079 
00080 #endif
00081 
00082 /* Get two 32 bit ints from a double.  */
00083 
00084 #define EXTRACT_WORDS(ix0,ix1,d)                                \
00085 do {                                                            \
00086   ieee_double_shape_type ew_u;                                  \
00087   ew_u.value = (d);                                             \
00088   (ix0) = ew_u.parts.msw;                                       \
00089   (ix1) = ew_u.parts.lsw;                                       \
00090 } while (0)
00091 
00092 /* Get the more significant 32 bit int from a double.  */
00093 
00094 #define GET_HIGH_WORD(i,d)                                      \
00095 do {                                                            \
00096   ieee_double_shape_type gh_u;                                  \
00097   gh_u.value = (d);                                             \
00098   (i) = gh_u.parts.msw;                                         \
00099 } while (0)
00100 
00101 /* Get the less significant 32 bit int from a double.  */
00102 
00103 #define GET_LOW_WORD(i,d)                                       \
00104 do {                                                            \
00105   ieee_double_shape_type gl_u;                                  \
00106   gl_u.value = (d);                                             \
00107   (i) = gl_u.parts.lsw;                                         \
00108 } while (0)
00109 
00110 /* Set the more significant 32 bits of a double from an int.  */
00111 
00112 #define SET_HIGH_WORD(d,v)                                      \
00113 do {                                                            \
00114   ieee_double_shape_type sh_u;                                  \
00115   sh_u.value = (d);                                             \
00116   sh_u.parts.msw = (v);                                         \
00117   (d) = sh_u.value;                                             \
00118 } while (0)
00119 
00120 /* Set the less significant 32 bits of a double from an int.  */
00121 
00122 #define SET_LOW_WORD(d,v)                                       \
00123 do {                                                            \
00124   ieee_double_shape_type sl_u;                                  \
00125   sl_u.value = (d);                                             \
00126   sl_u.parts.lsw = (v);                                         \
00127   (d) = sl_u.value;                                             \
00128 } while (0)
00129 
00130 /* Set a double from two 32 bit ints.  */
00131 
00132 #define INSERT_WORDS(d,ix0,ix1)                                 \
00133 do {                                                            \
00134   ieee_double_shape_type iw_u;                                  \
00135   iw_u.parts.msw = (ix0);                                       \
00136   iw_u.parts.lsw = (ix1);                                       \
00137   (d) = iw_u.value;                                             \
00138 } while (0)
00139 
00140 /*
00141  * A union which permits us to convert between a float and a 32 bit
00142  * int.
00143  */
00144 
00145 typedef union
00146 {
00147   float value;
00148   /* FIXME: Assumes 32 bit int.  */
00149   npy_uint32 word;
00150 } ieee_float_shape_type;
00151 
00152 /* Get a 32 bit int from a float.  */
00153 
00154 #define GET_FLOAT_WORD(i,d)                                     \
00155 do {                                                            \
00156   ieee_float_shape_type gf_u;                                   \
00157   gf_u.value = (d);                                             \
00158   (i) = gf_u.word;                                              \
00159 } while (0)
00160 
00161 /* Set a float from a 32 bit int.  */
00162 
00163 #define SET_FLOAT_WORD(d,i)                                     \
00164 do {                                                            \
00165   ieee_float_shape_type sf_u;                                   \
00166   sf_u.word = (i);                                              \
00167   (d) = sf_u.value;                                             \
00168 } while (0)
00169 
00170 #ifdef NPY_USE_C99_COMPLEX
00171 #include <complex.h>
00172 #endif
00173 
00174 /*
00175  * Long double support
00176  */
00177 #if defined(HAVE_LDOUBLE_INTEL_EXTENDED_12_BYTES_LE)
00178     /*
00179      * Intel extended 80 bits precision. Bit representation is
00180      *          |  junk  |     s  |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm|
00181      *          | 16 bits|  1 bit |    15 bits    |            64 bits            |
00182      *          |             a[2]                |     a[1]     |    a[0]        |
00183      *
00184      * 16 stronger bits of a[2] are junk
00185      */
00186     typedef npy_uint32 IEEEl2bitsrep_part;
00187 
00188 /* my machine */
00189 
00190     union IEEEl2bitsrep {
00191         npy_longdouble     e;
00192         IEEEl2bitsrep_part a[3];
00193     };
00194 
00195     #define LDBL_MANL_INDEX     0
00196     #define LDBL_MANL_MASK      0xFFFFFFFF
00197     #define LDBL_MANL_SHIFT     0
00198 
00199     #define LDBL_MANH_INDEX     1
00200     #define LDBL_MANH_MASK      0xFFFFFFFF
00201     #define LDBL_MANH_SHIFT     0
00202 
00203     #define LDBL_EXP_INDEX      2
00204     #define LDBL_EXP_MASK       0x7FFF
00205     #define LDBL_EXP_SHIFT      0
00206 
00207     #define LDBL_SIGN_INDEX     2
00208     #define LDBL_SIGN_MASK      0x8000
00209     #define LDBL_SIGN_SHIFT     15
00210 
00211     #define LDBL_NBIT           0x80000000
00212 
00213     typedef npy_uint32 ldouble_man_t;
00214     typedef npy_uint32 ldouble_exp_t;
00215     typedef npy_uint32 ldouble_sign_t;
00216 #elif defined(HAVE_LDOUBLE_INTEL_EXTENDED_16_BYTES_LE)
00217     /*
00218      * Intel extended 80 bits precision, 16 bytes alignment.. Bit representation is
00219      *          |  junk  |     s  |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm|
00220      *          | 16 bits|  1 bit |    15 bits    |            64 bits            |
00221      *          |             a[2]                |     a[1]     |    a[0]        |
00222      *
00223      * a[3] and 16 stronger bits of a[2] are junk
00224      */
00225     typedef npy_uint32 IEEEl2bitsrep_part;
00226 
00227     union IEEEl2bitsrep {
00228         npy_longdouble     e;
00229         IEEEl2bitsrep_part a[4];
00230     };
00231 
00232     #define LDBL_MANL_INDEX     0
00233     #define LDBL_MANL_MASK      0xFFFFFFFF
00234     #define LDBL_MANL_SHIFT     0
00235 
00236     #define LDBL_MANH_INDEX     1
00237     #define LDBL_MANH_MASK      0xFFFFFFFF
00238     #define LDBL_MANH_SHIFT     0
00239 
00240     #define LDBL_EXP_INDEX      2
00241     #define LDBL_EXP_MASK       0x7FFF
00242     #define LDBL_EXP_SHIFT      0
00243 
00244     #define LDBL_SIGN_INDEX     2
00245     #define LDBL_SIGN_MASK      0x8000
00246     #define LDBL_SIGN_SHIFT     15
00247 
00248     #define LDBL_NBIT           0x800000000
00249 
00250     typedef npy_uint32 ldouble_man_t;
00251     typedef npy_uint32 ldouble_exp_t;
00252     typedef npy_uint32 ldouble_sign_t;
00253 #elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_16_BYTES_BE) || \
00254       defined(HAVE_LDOUBLE_IEEE_DOUBLE_BE)
00255     /* 64 bits IEEE double precision aligned on 16 bytes: used by ppc arch on
00256      * Mac OS X */
00257 
00258     /*
00259      * IEEE double precision. Bit representation is
00260      *          |  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
00261      *          |1 bit|  11 bits  |            52 bits            |
00262      *          |          a[0]         |         a[1]            |
00263      */
00264     typedef npy_uint32 IEEEl2bitsrep_part;
00265 
00266     union IEEEl2bitsrep {
00267         npy_longdouble     e;
00268         IEEEl2bitsrep_part a[2];
00269     };
00270 
00271     #define LDBL_MANL_INDEX     1
00272     #define LDBL_MANL_MASK      0xFFFFFFFF
00273     #define LDBL_MANL_SHIFT     0
00274 
00275     #define LDBL_MANH_INDEX     0
00276     #define LDBL_MANH_MASK      0x000FFFFF
00277     #define LDBL_MANH_SHIFT     0
00278 
00279     #define LDBL_EXP_INDEX      0
00280     #define LDBL_EXP_MASK       0x7FF00000
00281     #define LDBL_EXP_SHIFT      20
00282 
00283     #define LDBL_SIGN_INDEX     0
00284     #define LDBL_SIGN_MASK      0x80000000
00285     #define LDBL_SIGN_SHIFT     31
00286 
00287     #define LDBL_NBIT           0
00288 
00289     typedef npy_uint32 ldouble_man_t;
00290     typedef npy_uint32 ldouble_exp_t;
00291     typedef npy_uint32 ldouble_sign_t;
00292 #elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_LE)
00293     /* 64 bits IEEE double precision, Little Endian. */
00294 
00295     /*
00296      * IEEE double precision. Bit representation is
00297      *          |  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
00298      *          |1 bit|  11 bits  |            52 bits            |
00299      *          |          a[1]         |         a[0]            |
00300      */
00301     typedef npy_uint32 IEEEl2bitsrep_part;
00302 
00303     union IEEEl2bitsrep {
00304         npy_longdouble     e;
00305         IEEEl2bitsrep_part a[2];
00306     };
00307 
00308     #define LDBL_MANL_INDEX     0
00309     #define LDBL_MANL_MASK      0xFFFFFFFF
00310     #define LDBL_MANL_SHIFT     0
00311 
00312     #define LDBL_MANH_INDEX     1
00313     #define LDBL_MANH_MASK      0x000FFFFF
00314     #define LDBL_MANH_SHIFT     0
00315 
00316     #define LDBL_EXP_INDEX      1
00317     #define LDBL_EXP_MASK       0x7FF00000
00318     #define LDBL_EXP_SHIFT      20
00319 
00320     #define LDBL_SIGN_INDEX     1
00321     #define LDBL_SIGN_MASK      0x80000000
00322     #define LDBL_SIGN_SHIFT     31
00323 
00324     #define LDBL_NBIT           0x00000080
00325 
00326     typedef npy_uint32 ldouble_man_t;
00327     typedef npy_uint32 ldouble_exp_t;
00328     typedef npy_uint32 ldouble_sign_t;
00329 #elif defined(HAVE_LDOUBLE_IEEE_QUAD_BE)
00330     /*
00331      * IEEE quad precision, Big Endian. Bit representation is
00332      *          |  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
00333      *          |1 bit|  15 bits  |            112 bits           |
00334      *          |          a[0]         |           a[1]          |
00335      */
00336     typedef npy_uint64 IEEEl2bitsrep_part;
00337 
00338     union IEEEl2bitsrep {
00339         npy_longdouble     e;
00340         IEEEl2bitsrep_part a[2];
00341     };
00342 
00343     #define LDBL_MANL_INDEX     1
00344     #define LDBL_MANL_MASK      0xFFFFFFFFFFFFFFFF
00345     #define LDBL_MANL_SHIFT     0
00346 
00347     #define LDBL_MANH_INDEX     0
00348     #define LDBL_MANH_MASK      0x0000FFFFFFFFFFFF
00349     #define LDBL_MANH_SHIFT     0
00350 
00351     #define LDBL_EXP_INDEX      0
00352     #define LDBL_EXP_MASK       0x7FFF000000000000
00353     #define LDBL_EXP_SHIFT      48
00354 
00355     #define LDBL_SIGN_INDEX     0
00356     #define LDBL_SIGN_MASK      0x8000000000000000
00357     #define LDBL_SIGN_SHIFT     63
00358 
00359     #define LDBL_NBIT           0
00360 
00361     typedef npy_uint64 ldouble_man_t;
00362     typedef npy_uint64 ldouble_exp_t;
00363     typedef npy_uint32 ldouble_sign_t;
00364 #elif defined(HAVE_LDOUBLE_IEEE_QUAD_LE)
00365     /*
00366      * IEEE quad precision, Little Endian. Bit representation is
00367      *          |  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
00368      *          |1 bit|  15 bits  |            112 bits           |
00369      *          |          a[1]         |           a[0]          |
00370      */
00371     typedef npy_uint64 IEEEl2bitsrep_part;
00372 
00373     union IEEEl2bitsrep {
00374         npy_longdouble     e;
00375         IEEEl2bitsrep_part a[2];
00376     };
00377 
00378     #define LDBL_MANL_INDEX     0
00379     #define LDBL_MANL_MASK      0xFFFFFFFFFFFFFFFF
00380     #define LDBL_MANL_SHIFT     0
00381 
00382     #define LDBL_MANH_INDEX     1
00383     #define LDBL_MANH_MASK      0x0000FFFFFFFFFFFF
00384     #define LDBL_MANH_SHIFT     0
00385 
00386     #define LDBL_EXP_INDEX      1
00387     #define LDBL_EXP_MASK       0x7FFF000000000000
00388     #define LDBL_EXP_SHIFT      48
00389 
00390     #define LDBL_SIGN_INDEX     1
00391     #define LDBL_SIGN_MASK      0x8000000000000000
00392     #define LDBL_SIGN_SHIFT     63
00393 
00394     #define LDBL_NBIT           0
00395 
00396     typedef npy_uint64 ldouble_man_t;
00397     typedef npy_uint64 ldouble_exp_t;
00398     typedef npy_uint32 ldouble_sign_t;
00399 #endif
00400 
00401 #ifndef HAVE_LDOUBLE_DOUBLE_DOUBLE_BE
00402 /* Get the sign bit of x. x should be of type IEEEl2bitsrep */
00403 #define GET_LDOUBLE_SIGN(x) \
00404     (((x).a[LDBL_SIGN_INDEX] & LDBL_SIGN_MASK) >> LDBL_SIGN_SHIFT)
00405 
00406 /* Set the sign bit of x to v. x should be of type IEEEl2bitsrep */
00407 #define SET_LDOUBLE_SIGN(x, v) \
00408   ((x).a[LDBL_SIGN_INDEX] =                                             \
00409    ((x).a[LDBL_SIGN_INDEX] & ~LDBL_SIGN_MASK) |                         \
00410    (((IEEEl2bitsrep_part)(v) << LDBL_SIGN_SHIFT) & LDBL_SIGN_MASK))
00411 
00412 /* Get the exp bits of x. x should be of type IEEEl2bitsrep */
00413 #define GET_LDOUBLE_EXP(x) \
00414     (((x).a[LDBL_EXP_INDEX] & LDBL_EXP_MASK) >> LDBL_EXP_SHIFT)
00415 
00416 /* Set the exp bit of x to v. x should be of type IEEEl2bitsrep */
00417 #define SET_LDOUBLE_EXP(x, v) \
00418   ((x).a[LDBL_EXP_INDEX] =                                              \
00419    ((x).a[LDBL_EXP_INDEX] & ~LDBL_EXP_MASK) |                           \
00420    (((IEEEl2bitsrep_part)(v) << LDBL_EXP_SHIFT) & LDBL_EXP_MASK))
00421 
00422 /* Get the manl bits of x. x should be of type IEEEl2bitsrep */
00423 #define GET_LDOUBLE_MANL(x) \
00424     (((x).a[LDBL_MANL_INDEX] & LDBL_MANL_MASK) >> LDBL_MANL_SHIFT)
00425 
00426 /* Set the manl bit of x to v. x should be of type IEEEl2bitsrep */
00427 #define SET_LDOUBLE_MANL(x, v) \
00428   ((x).a[LDBL_MANL_INDEX] =                                             \
00429    ((x).a[LDBL_MANL_INDEX] & ~LDBL_MANL_MASK) |                         \
00430    (((IEEEl2bitsrep_part)(v) << LDBL_MANL_SHIFT) & LDBL_MANL_MASK))
00431 
00432 /* Get the manh bits of x. x should be of type IEEEl2bitsrep */
00433 #define GET_LDOUBLE_MANH(x) \
00434     (((x).a[LDBL_MANH_INDEX] & LDBL_MANH_MASK) >> LDBL_MANH_SHIFT)
00435 
00436 /* Set the manh bit of x to v. x should be of type IEEEl2bitsrep */
00437 #define SET_LDOUBLE_MANH(x, v) \
00438     ((x).a[LDBL_MANH_INDEX] = \
00439      ((x).a[LDBL_MANH_INDEX] & ~LDBL_MANH_MASK) |                       \
00440      (((IEEEl2bitsrep_part)(v) << LDBL_MANH_SHIFT) & LDBL_MANH_MASK))
00441 
00442 #endif /* #ifndef HAVE_LDOUBLE_DOUBLE_DOUBLE_BE */
00443 
00444 /*
00445  * Those unions are used to convert a pointer of npy_cdouble to native C99
00446  * complex or our own complex type independently on whether C99 complex
00447  * support is available
00448  */
00449 #ifdef NPY_USE_C99_COMPLEX
00450 typedef union {
00451         npy_cdouble npy_z;
00452         complex double c99_z;
00453 } __npy_cdouble_to_c99_cast;
00454 
00455 typedef union {
00456         npy_cfloat npy_z;
00457         complex float c99_z;
00458 } __npy_cfloat_to_c99_cast;
00459 
00460 typedef union {
00461         npy_clongdouble npy_z;
00462         complex long double c99_z;
00463 } __npy_clongdouble_to_c99_cast;
00464 #else
00465 typedef union {
00466         npy_cdouble npy_z;
00467         npy_cdouble c99_z;
00468 } __npy_cdouble_to_c99_cast;
00469 
00470 typedef union {
00471         npy_cfloat npy_z;
00472         npy_cfloat c99_z;
00473 } __npy_cfloat_to_c99_cast;
00474 
00475 typedef union {
00476         npy_clongdouble npy_z;
00477         npy_clongdouble c99_z;
00478 } __npy_clongdouble_to_c99_cast;
00479 #endif
00480 
00481 #endif /* !_NPY_MATH_PRIVATE_H_ */