numpy
2.0.0
|
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_ */