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 
00189     union IEEEl2bitsrep {
00190         npy_longdouble     e;
00191         IEEEl2bitsrep_part a[3];
00192     };
00193 
00194     #define LDBL_MANL_INDEX     0
00195     #define LDBL_MANL_MASK      0xFFFFFFFF
00196     #define LDBL_MANL_SHIFT     0
00197 
00198     #define LDBL_MANH_INDEX     1
00199     #define LDBL_MANH_MASK      0xFFFFFFFF
00200     #define LDBL_MANH_SHIFT     0
00201 
00202     #define LDBL_EXP_INDEX      2
00203     #define LDBL_EXP_MASK       0x7FFF
00204     #define LDBL_EXP_SHIFT      0
00205 
00206     #define LDBL_SIGN_INDEX     2
00207     #define LDBL_SIGN_MASK      0x8000
00208     #define LDBL_SIGN_SHIFT     15
00209 
00210     #define LDBL_NBIT           0x80000000
00211 
00212     typedef npy_uint32 ldouble_man_t;
00213     typedef npy_uint32 ldouble_exp_t;
00214     typedef npy_uint32 ldouble_sign_t;
00215 #elif defined(HAVE_LDOUBLE_INTEL_EXTENDED_16_BYTES_LE)
00216     /*
00217      * Intel extended 80 bits precision, 16 bytes alignment.. Bit representation is
00218      *          |  junk  |     s  |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm|
00219      *          | 16 bits|  1 bit |    15 bits    |            64 bits            |
00220      *          |             a[2]                |     a[1]     |    a[0]        |
00221      *
00222      * a[3] and 16 stronger bits of a[2] are junk
00223      */
00224     typedef npy_uint32 IEEEl2bitsrep_part;
00225 
00226     union IEEEl2bitsrep {
00227         npy_longdouble     e;
00228         IEEEl2bitsrep_part a[4];
00229     };
00230 
00231     #define LDBL_MANL_INDEX     0
00232     #define LDBL_MANL_MASK      0xFFFFFFFF
00233     #define LDBL_MANL_SHIFT     0
00234 
00235     #define LDBL_MANH_INDEX     1
00236     #define LDBL_MANH_MASK      0xFFFFFFFF
00237     #define LDBL_MANH_SHIFT     0
00238 
00239     #define LDBL_EXP_INDEX      2
00240     #define LDBL_EXP_MASK       0x7FFF
00241     #define LDBL_EXP_SHIFT      0
00242 
00243     #define LDBL_SIGN_INDEX     2
00244     #define LDBL_SIGN_MASK      0x8000
00245     #define LDBL_SIGN_SHIFT     15
00246 
00247     #define LDBL_NBIT           0x800000000
00248 
00249     typedef npy_uint32 ldouble_man_t;
00250     typedef npy_uint32 ldouble_exp_t;
00251     typedef npy_uint32 ldouble_sign_t;
00252 #elif defined(HAVE_LDOUBLE_MOTOROLA_EXTENDED_12_BYTES_BE)
00253     /*
00254      * Motorola extended 80 bits precision. Bit representation is
00255      *          |     s  |eeeeeeeeeeeeeee|  junk  |mmmmmmmm................mmmmmmm|
00256      *          |  1 bit |    15 bits    | 16 bits|            64 bits            |
00257      *          |             a[0]                |     a[1]     |    a[2]        |
00258      *
00259      * 16 low bits of a[0] are junk
00260      */
00261     typedef npy_uint32 IEEEl2bitsrep_part;
00262 
00263 
00264     union IEEEl2bitsrep {
00265         npy_longdouble     e;
00266         IEEEl2bitsrep_part a[3];
00267     };
00268 
00269     #define LDBL_MANL_INDEX     2
00270     #define LDBL_MANL_MASK      0xFFFFFFFF
00271     #define LDBL_MANL_SHIFT     0
00272 
00273     #define LDBL_MANH_INDEX     1
00274     #define LDBL_MANH_MASK      0xFFFFFFFF
00275     #define LDBL_MANH_SHIFT     0
00276 
00277     #define LDBL_EXP_INDEX      0
00278     #define LDBL_EXP_MASK       0x7FFF0000
00279     #define LDBL_EXP_SHIFT      16
00280 
00281     #define LDBL_SIGN_INDEX     0
00282     #define LDBL_SIGN_MASK      0x80000000
00283     #define LDBL_SIGN_SHIFT     31
00284 
00285     #define LDBL_NBIT           0x80000000
00286 
00287     typedef npy_uint32 ldouble_man_t;
00288     typedef npy_uint32 ldouble_exp_t;
00289     typedef npy_uint32 ldouble_sign_t;
00290 #elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_16_BYTES_BE) || \
00291       defined(HAVE_LDOUBLE_IEEE_DOUBLE_BE)
00292     /* 64 bits IEEE double precision aligned on 16 bytes: used by ppc arch on
00293      * Mac OS X */
00294 
00295     /*
00296      * IEEE double precision. Bit representation is
00297      *          |  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
00298      *          |1 bit|  11 bits  |            52 bits            |
00299      *          |          a[0]         |         a[1]            |
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     1
00309     #define LDBL_MANL_MASK      0xFFFFFFFF
00310     #define LDBL_MANL_SHIFT     0
00311 
00312     #define LDBL_MANH_INDEX     0
00313     #define LDBL_MANH_MASK      0x000FFFFF
00314     #define LDBL_MANH_SHIFT     0
00315 
00316     #define LDBL_EXP_INDEX      0
00317     #define LDBL_EXP_MASK       0x7FF00000
00318     #define LDBL_EXP_SHIFT      20
00319 
00320     #define LDBL_SIGN_INDEX     0
00321     #define LDBL_SIGN_MASK      0x80000000
00322     #define LDBL_SIGN_SHIFT     31
00323 
00324     #define LDBL_NBIT           0
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_DOUBLE_LE)
00330     /* 64 bits IEEE double precision, Little Endian. */
00331 
00332     /*
00333      * IEEE double precision. Bit representation is
00334      *          |  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
00335      *          |1 bit|  11 bits  |            52 bits            |
00336      *          |          a[1]         |         a[0]            |
00337      */
00338     typedef npy_uint32 IEEEl2bitsrep_part;
00339 
00340     union IEEEl2bitsrep {
00341         npy_longdouble     e;
00342         IEEEl2bitsrep_part a[2];
00343     };
00344 
00345     #define LDBL_MANL_INDEX     0
00346     #define LDBL_MANL_MASK      0xFFFFFFFF
00347     #define LDBL_MANL_SHIFT     0
00348 
00349     #define LDBL_MANH_INDEX     1
00350     #define LDBL_MANH_MASK      0x000FFFFF
00351     #define LDBL_MANH_SHIFT     0
00352 
00353     #define LDBL_EXP_INDEX      1
00354     #define LDBL_EXP_MASK       0x7FF00000
00355     #define LDBL_EXP_SHIFT      20
00356 
00357     #define LDBL_SIGN_INDEX     1
00358     #define LDBL_SIGN_MASK      0x80000000
00359     #define LDBL_SIGN_SHIFT     31
00360 
00361     #define LDBL_NBIT           0x00000080
00362 
00363     typedef npy_uint32 ldouble_man_t;
00364     typedef npy_uint32 ldouble_exp_t;
00365     typedef npy_uint32 ldouble_sign_t;
00366 #elif defined(HAVE_LDOUBLE_IEEE_QUAD_BE)
00367     /*
00368      * IEEE quad precision, Big Endian. Bit representation is
00369      *          |  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
00370      *          |1 bit|  15 bits  |            112 bits           |
00371      *          |          a[0]         |           a[1]          |
00372      */
00373     typedef npy_uint64 IEEEl2bitsrep_part;
00374 
00375     union IEEEl2bitsrep {
00376         npy_longdouble     e;
00377         IEEEl2bitsrep_part a[2];
00378     };
00379 
00380     #define LDBL_MANL_INDEX     1
00381     #define LDBL_MANL_MASK      0xFFFFFFFFFFFFFFFF
00382     #define LDBL_MANL_SHIFT     0
00383 
00384     #define LDBL_MANH_INDEX     0
00385     #define LDBL_MANH_MASK      0x0000FFFFFFFFFFFF
00386     #define LDBL_MANH_SHIFT     0
00387 
00388     #define LDBL_EXP_INDEX      0
00389     #define LDBL_EXP_MASK       0x7FFF000000000000
00390     #define LDBL_EXP_SHIFT      48
00391 
00392     #define LDBL_SIGN_INDEX     0
00393     #define LDBL_SIGN_MASK      0x8000000000000000
00394     #define LDBL_SIGN_SHIFT     63
00395 
00396     #define LDBL_NBIT           0
00397 
00398     typedef npy_uint64 ldouble_man_t;
00399     typedef npy_uint64 ldouble_exp_t;
00400     typedef npy_uint32 ldouble_sign_t;
00401 #elif defined(HAVE_LDOUBLE_IEEE_QUAD_LE)
00402     /*
00403      * IEEE quad precision, Little Endian. Bit representation is
00404      *          |  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
00405      *          |1 bit|  15 bits  |            112 bits           |
00406      *          |          a[1]         |           a[0]          |
00407      */
00408     typedef npy_uint64 IEEEl2bitsrep_part;
00409 
00410     union IEEEl2bitsrep {
00411         npy_longdouble     e;
00412         IEEEl2bitsrep_part a[2];
00413     };
00414 
00415     #define LDBL_MANL_INDEX     0
00416     #define LDBL_MANL_MASK      0xFFFFFFFFFFFFFFFF
00417     #define LDBL_MANL_SHIFT     0
00418 
00419     #define LDBL_MANH_INDEX     1
00420     #define LDBL_MANH_MASK      0x0000FFFFFFFFFFFF
00421     #define LDBL_MANH_SHIFT     0
00422 
00423     #define LDBL_EXP_INDEX      1
00424     #define LDBL_EXP_MASK       0x7FFF000000000000
00425     #define LDBL_EXP_SHIFT      48
00426 
00427     #define LDBL_SIGN_INDEX     1
00428     #define LDBL_SIGN_MASK      0x8000000000000000
00429     #define LDBL_SIGN_SHIFT     63
00430 
00431     #define LDBL_NBIT           0
00432 
00433     typedef npy_uint64 ldouble_man_t;
00434     typedef npy_uint64 ldouble_exp_t;
00435     typedef npy_uint32 ldouble_sign_t;
00436 #endif
00437 
00438 #if !defined(HAVE_LDOUBLE_DOUBLE_DOUBLE_BE) && \
00439     !defined(HAVE_LDOUBLE_DOUBLE_DOUBLE_LE)
00440 /* Get the sign bit of x. x should be of type IEEEl2bitsrep */
00441 #define GET_LDOUBLE_SIGN(x) \
00442     (((x).a[LDBL_SIGN_INDEX] & LDBL_SIGN_MASK) >> LDBL_SIGN_SHIFT)
00443 
00444 /* Set the sign bit of x to v. x should be of type IEEEl2bitsrep */
00445 #define SET_LDOUBLE_SIGN(x, v) \
00446   ((x).a[LDBL_SIGN_INDEX] =                                             \
00447    ((x).a[LDBL_SIGN_INDEX] & ~LDBL_SIGN_MASK) |                         \
00448    (((IEEEl2bitsrep_part)(v) << LDBL_SIGN_SHIFT) & LDBL_SIGN_MASK))
00449 
00450 /* Get the exp bits of x. x should be of type IEEEl2bitsrep */
00451 #define GET_LDOUBLE_EXP(x) \
00452     (((x).a[LDBL_EXP_INDEX] & LDBL_EXP_MASK) >> LDBL_EXP_SHIFT)
00453 
00454 /* Set the exp bit of x to v. x should be of type IEEEl2bitsrep */
00455 #define SET_LDOUBLE_EXP(x, v) \
00456   ((x).a[LDBL_EXP_INDEX] =                                              \
00457    ((x).a[LDBL_EXP_INDEX] & ~LDBL_EXP_MASK) |                           \
00458    (((IEEEl2bitsrep_part)(v) << LDBL_EXP_SHIFT) & LDBL_EXP_MASK))
00459 
00460 /* Get the manl bits of x. x should be of type IEEEl2bitsrep */
00461 #define GET_LDOUBLE_MANL(x) \
00462     (((x).a[LDBL_MANL_INDEX] & LDBL_MANL_MASK) >> LDBL_MANL_SHIFT)
00463 
00464 /* Set the manl bit of x to v. x should be of type IEEEl2bitsrep */
00465 #define SET_LDOUBLE_MANL(x, v) \
00466   ((x).a[LDBL_MANL_INDEX] =                                             \
00467    ((x).a[LDBL_MANL_INDEX] & ~LDBL_MANL_MASK) |                         \
00468    (((IEEEl2bitsrep_part)(v) << LDBL_MANL_SHIFT) & LDBL_MANL_MASK))
00469 
00470 /* Get the manh bits of x. x should be of type IEEEl2bitsrep */
00471 #define GET_LDOUBLE_MANH(x) \
00472     (((x).a[LDBL_MANH_INDEX] & LDBL_MANH_MASK) >> LDBL_MANH_SHIFT)
00473 
00474 /* Set the manh bit of x to v. x should be of type IEEEl2bitsrep */
00475 #define SET_LDOUBLE_MANH(x, v) \
00476     ((x).a[LDBL_MANH_INDEX] = \
00477      ((x).a[LDBL_MANH_INDEX] & ~LDBL_MANH_MASK) |                       \
00478      (((IEEEl2bitsrep_part)(v) << LDBL_MANH_SHIFT) & LDBL_MANH_MASK))
00479 
00480 #endif /* #ifndef HAVE_LDOUBLE_DOUBLE_DOUBLE_BE */
00481 
00482 /*
00483  * Those unions are used to convert a pointer of npy_cdouble to native C99
00484  * complex or our own complex type independently on whether C99 complex
00485  * support is available
00486  */
00487 #ifdef NPY_USE_C99_COMPLEX
00488 
00489 /*
00490  * Microsoft C defines _MSC_VER
00491  * Intel compiler does not use MSVC complex types, but defines _MSC_VER by
00492  * default.
00493  */
00494 #if defined(_MSC_VER) && !defined(__INTEL_COMPILER)
00495 typedef union {
00496         npy_cdouble npy_z;
00497         _Dcomplex c99_z;
00498 } __npy_cdouble_to_c99_cast;
00499 
00500 typedef union {
00501         npy_cfloat npy_z;
00502         _Fcomplex c99_z;
00503 } __npy_cfloat_to_c99_cast;
00504 
00505 typedef union {
00506         npy_clongdouble npy_z;
00507         _Lcomplex c99_z;
00508 } __npy_clongdouble_to_c99_cast;
00509 #else /* !_MSC_VER */
00510 typedef union {
00511         npy_cdouble npy_z;
00512         complex double c99_z;
00513 } __npy_cdouble_to_c99_cast;
00514 
00515 typedef union {
00516         npy_cfloat npy_z;
00517         complex float c99_z;
00518 } __npy_cfloat_to_c99_cast;
00519 
00520 typedef union {
00521         npy_clongdouble npy_z;
00522         complex long double c99_z;
00523 } __npy_clongdouble_to_c99_cast;
00524 #endif /* !_MSC_VER */
00525 
00526 #else /* !NPY_USE_C99_COMPLEX */
00527 typedef union {
00528         npy_cdouble npy_z;
00529         npy_cdouble c99_z;
00530 } __npy_cdouble_to_c99_cast;
00531 
00532 typedef union {
00533         npy_cfloat npy_z;
00534         npy_cfloat c99_z;
00535 } __npy_cfloat_to_c99_cast;
00536 
00537 typedef union {
00538         npy_clongdouble npy_z;
00539         npy_clongdouble c99_z;
00540 } __npy_clongdouble_to_c99_cast;
00541 #endif /* !NPY_USE_C99_COMPLEX */
00542 
00543 
00544 #endif /* !_NPY_MATH_PRIVATE_H_ */