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