| 1 |
|
|---|
| 2 | /* @(#)fdlibm.h 5.1 93/09/24 */
|
|---|
| 3 | /*
|
|---|
| 4 | * ====================================================
|
|---|
| 5 | * Copyright (C) 1993, 2000 by Sun Microsystems, Inc. All rights reserved.
|
|---|
| 6 | *
|
|---|
| 7 | * Developed at SunPro, a Sun Microsystems, Inc. business.
|
|---|
| 8 | * Permission to use, copy, modify, and distribute this
|
|---|
| 9 | * software is freely granted, provided that this notice
|
|---|
| 10 | * is preserved.
|
|---|
| 11 | * ====================================================
|
|---|
| 12 | */
|
|---|
| 13 |
|
|---|
| 14 |
|
|---|
| 15 | #include <config.h>
|
|---|
| 16 | #include <stdlib.h>
|
|---|
| 17 |
|
|---|
| 18 | /* CYGNUS LOCAL: Include files. */
|
|---|
| 19 | #include "ieeefp.h"
|
|---|
| 20 |
|
|---|
| 21 | #include "mprec.h"
|
|---|
| 22 |
|
|---|
| 23 | /* CYGNUS LOCAL: Default to XOPEN_MODE. */
|
|---|
| 24 | #define _XOPEN_MODE
|
|---|
| 25 |
|
|---|
| 26 | #ifdef __P
|
|---|
| 27 | #undef __P
|
|---|
| 28 | #endif
|
|---|
| 29 |
|
|---|
| 30 | #ifdef __STDC__
|
|---|
| 31 | #define __P(p) p
|
|---|
| 32 | #else
|
|---|
| 33 | #define __P(p) ()
|
|---|
| 34 | #endif
|
|---|
| 35 |
|
|---|
| 36 | #ifndef HUGE
|
|---|
| 37 | #define HUGE ((float)3.40282346638528860e+38)
|
|---|
| 38 | #endif
|
|---|
| 39 |
|
|---|
| 40 | /*
|
|---|
| 41 | * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
|
|---|
| 42 | * (one may replace the following line by "#include <values.h>")
|
|---|
| 43 | */
|
|---|
| 44 |
|
|---|
| 45 | #define X_TLOSS 1.41484755040568800000e+16
|
|---|
| 46 |
|
|---|
| 47 | /* These typedefs are true for the targets running Java. */
|
|---|
| 48 |
|
|---|
| 49 | #define _IEEE_LIBM
|
|---|
| 50 |
|
|---|
| 51 | #ifdef __cplusplus
|
|---|
| 52 | extern "C" {
|
|---|
| 53 | #endif
|
|---|
| 54 |
|
|---|
| 55 | /*
|
|---|
| 56 | * ANSI/POSIX
|
|---|
| 57 | */
|
|---|
| 58 | extern double acos __P((double));
|
|---|
| 59 | extern double asin __P((double));
|
|---|
| 60 | extern double atan __P((double));
|
|---|
| 61 | extern double atan2 __P((double, double));
|
|---|
| 62 | extern double cos __P((double));
|
|---|
| 63 | extern double sin __P((double));
|
|---|
| 64 | extern double tan __P((double));
|
|---|
| 65 |
|
|---|
| 66 | extern double cosh __P((double));
|
|---|
| 67 | extern double sinh __P((double));
|
|---|
| 68 | extern double tanh __P((double));
|
|---|
| 69 |
|
|---|
| 70 | extern double exp __P((double));
|
|---|
| 71 | extern double frexp __P((double, int *));
|
|---|
| 72 | extern double ldexp __P((double, int));
|
|---|
| 73 | extern double log __P((double));
|
|---|
| 74 | extern double log10 __P((double));
|
|---|
| 75 | extern double modf __P((double, double *));
|
|---|
| 76 |
|
|---|
| 77 | extern double pow __P((double, double));
|
|---|
| 78 | extern double sqrt __P((double));
|
|---|
| 79 |
|
|---|
| 80 | extern double ceil __P((double));
|
|---|
| 81 | extern double fabs __P((double));
|
|---|
| 82 | extern double floor __P((double));
|
|---|
| 83 | extern double fmod __P((double, double));
|
|---|
| 84 |
|
|---|
| 85 | extern double erf __P((double));
|
|---|
| 86 | extern double erfc __P((double));
|
|---|
| 87 | extern double gamma __P((double));
|
|---|
| 88 | extern double hypot __P((double, double));
|
|---|
| 89 | extern int isnan __P((double));
|
|---|
| 90 | extern int finite __P((double));
|
|---|
| 91 | extern double j0 __P((double));
|
|---|
| 92 | extern double j1 __P((double));
|
|---|
| 93 | extern double jn __P((int, double));
|
|---|
| 94 | extern double lgamma __P((double));
|
|---|
| 95 | extern double y0 __P((double));
|
|---|
| 96 | extern double y1 __P((double));
|
|---|
| 97 | extern double yn __P((int, double));
|
|---|
| 98 |
|
|---|
| 99 | extern double acosh __P((double));
|
|---|
| 100 | extern double asinh __P((double));
|
|---|
| 101 | extern double atanh __P((double));
|
|---|
| 102 | extern double cbrt __P((double));
|
|---|
| 103 | extern double logb __P((double));
|
|---|
| 104 | extern double nextafter __P((double, double));
|
|---|
| 105 | extern double remainder __P((double, double));
|
|---|
| 106 |
|
|---|
| 107 | /* Functions that are not documented, and are not in <math.h>. */
|
|---|
| 108 |
|
|---|
| 109 | extern double logb __P((double));
|
|---|
| 110 | #ifdef _SCALB_INT
|
|---|
| 111 | extern double scalb __P((double, int));
|
|---|
| 112 | #else
|
|---|
| 113 | extern double scalb __P((double, double));
|
|---|
| 114 | #endif
|
|---|
| 115 | extern double significand __P((double));
|
|---|
| 116 |
|
|---|
| 117 | /* ieee style elementary functions */
|
|---|
| 118 | extern double __ieee754_sqrt __P((double));
|
|---|
| 119 | extern double __ieee754_acos __P((double));
|
|---|
| 120 | extern double __ieee754_acosh __P((double));
|
|---|
| 121 | extern double __ieee754_log __P((double));
|
|---|
| 122 | extern double __ieee754_atanh __P((double));
|
|---|
| 123 | extern double __ieee754_asin __P((double));
|
|---|
| 124 | extern double __ieee754_atan2 __P((double,double));
|
|---|
| 125 | extern double __ieee754_exp __P((double));
|
|---|
| 126 | extern double __ieee754_cosh __P((double));
|
|---|
| 127 | extern double __ieee754_fmod __P((double,double));
|
|---|
| 128 | extern double __ieee754_pow __P((double,double));
|
|---|
| 129 | extern double __ieee754_lgamma_r __P((double,int *));
|
|---|
| 130 | extern double __ieee754_gamma_r __P((double,int *));
|
|---|
| 131 | extern double __ieee754_log10 __P((double));
|
|---|
| 132 | extern double __ieee754_sinh __P((double));
|
|---|
| 133 | extern double __ieee754_hypot __P((double,double));
|
|---|
| 134 | extern double __ieee754_j0 __P((double));
|
|---|
| 135 | extern double __ieee754_j1 __P((double));
|
|---|
| 136 | extern double __ieee754_y0 __P((double));
|
|---|
| 137 | extern double __ieee754_y1 __P((double));
|
|---|
| 138 | extern double __ieee754_jn __P((int,double));
|
|---|
| 139 | extern double __ieee754_yn __P((int,double));
|
|---|
| 140 | extern double __ieee754_remainder __P((double,double));
|
|---|
| 141 | extern int32_t __ieee754_rem_pio2 __P((double,double*));
|
|---|
| 142 | #ifdef _SCALB_INT
|
|---|
| 143 | extern double __ieee754_scalb __P((double,int));
|
|---|
| 144 | #else
|
|---|
| 145 | extern double __ieee754_scalb __P((double,double));
|
|---|
| 146 | #endif
|
|---|
| 147 |
|
|---|
| 148 | /* fdlibm kernel function */
|
|---|
| 149 | extern double __kernel_standard __P((double,double,int));
|
|---|
| 150 | extern double __kernel_sin __P((double,double,int));
|
|---|
| 151 | extern double __kernel_cos __P((double,double));
|
|---|
| 152 | extern double __kernel_tan __P((double,double,int));
|
|---|
| 153 | extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const int32_t*));
|
|---|
| 154 |
|
|---|
| 155 | /* Undocumented float functions. */
|
|---|
| 156 | extern float logbf __P((float));
|
|---|
| 157 | #ifdef _SCALB_INT
|
|---|
| 158 | extern float scalbf __P((float, int));
|
|---|
| 159 | #else
|
|---|
| 160 | extern float scalbf __P((float, float));
|
|---|
| 161 | #endif
|
|---|
| 162 | extern float significandf __P((float));
|
|---|
| 163 |
|
|---|
| 164 | /*
|
|---|
| 165 | * Functions callable from C, intended to support IEEE arithmetic.
|
|---|
| 166 | */
|
|---|
| 167 | extern double copysign __P((double, double));
|
|---|
| 168 | extern int ilogb __P((double));
|
|---|
| 169 | extern double rint __P((double));
|
|---|
| 170 | extern float rintf __P((float));
|
|---|
| 171 | extern double scalbn __P((double, int));
|
|---|
| 172 |
|
|---|
| 173 | /* ieee style elementary float functions */
|
|---|
| 174 | extern float __ieee754_sqrtf __P((float));
|
|---|
| 175 | extern float __ieee754_acosf __P((float));
|
|---|
| 176 | extern float __ieee754_acoshf __P((float));
|
|---|
| 177 | extern float __ieee754_logf __P((float));
|
|---|
| 178 | extern float __ieee754_atanhf __P((float));
|
|---|
| 179 | extern float __ieee754_asinf __P((float));
|
|---|
| 180 | extern float __ieee754_atan2f __P((float,float));
|
|---|
| 181 | extern float __ieee754_expf __P((float));
|
|---|
| 182 | extern float __ieee754_coshf __P((float));
|
|---|
| 183 | extern float __ieee754_fmodf __P((float,float));
|
|---|
| 184 | extern float __ieee754_powf __P((float,float));
|
|---|
| 185 | extern float __ieee754_lgammaf_r __P((float,int *));
|
|---|
| 186 | extern float __ieee754_gammaf_r __P((float,int *));
|
|---|
| 187 | extern float __ieee754_log10f __P((float));
|
|---|
| 188 | extern float __ieee754_sinhf __P((float));
|
|---|
| 189 | extern float __ieee754_hypotf __P((float,float));
|
|---|
| 190 | extern float __ieee754_j0f __P((float));
|
|---|
| 191 | extern float __ieee754_j1f __P((float));
|
|---|
| 192 | extern float __ieee754_y0f __P((float));
|
|---|
| 193 | extern float __ieee754_y1f __P((float));
|
|---|
| 194 | extern float __ieee754_jnf __P((int,float));
|
|---|
| 195 | extern float __ieee754_ynf __P((int,float));
|
|---|
| 196 | extern float __ieee754_remainderf __P((float,float));
|
|---|
| 197 | extern int32_t __ieee754_rem_pio2f __P((float,float*));
|
|---|
| 198 | #ifdef _SCALB_INT
|
|---|
| 199 | extern float __ieee754_scalbf __P((float,int));
|
|---|
| 200 | #else
|
|---|
| 201 | extern float __ieee754_scalbf __P((float,float));
|
|---|
| 202 | #endif
|
|---|
| 203 |
|
|---|
| 204 | /* float versions of fdlibm kernel functions */
|
|---|
| 205 | extern float __kernel_sinf __P((float,float,int));
|
|---|
| 206 | extern float __kernel_cosf __P((float,float));
|
|---|
| 207 | extern float __kernel_tanf __P((float,float,int));
|
|---|
| 208 | extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const int32_t*));
|
|---|
| 209 |
|
|---|
| 210 | /* The original code used statements like
|
|---|
| 211 | n0 = ((*(int*)&one)>>29)^1; * index of high word *
|
|---|
| 212 | ix0 = *(n0+(int*)&x); * high word of x *
|
|---|
| 213 | ix1 = *((1-n0)+(int*)&x); * low word of x *
|
|---|
| 214 | to dig two 32 bit words out of the 64 bit IEEE floating point
|
|---|
| 215 | value. That is non-ANSI, and, moreover, the gcc instruction
|
|---|
| 216 | scheduler gets it wrong. We instead use the following macros.
|
|---|
| 217 | Unlike the original code, we determine the endianness at compile
|
|---|
| 218 | time, not at run time; I don't see much benefit to selecting
|
|---|
| 219 | endianness at run time. */
|
|---|
| 220 |
|
|---|
| 221 | #ifndef __IEEE_BIG_ENDIAN
|
|---|
| 222 | #ifndef __IEEE_LITTLE_ENDIAN
|
|---|
| 223 | #error Must define endianness
|
|---|
| 224 | #endif
|
|---|
| 225 | #endif
|
|---|
| 226 |
|
|---|
| 227 | /* A union which permits us to convert between a double and two 32 bit
|
|---|
| 228 | ints. */
|
|---|
| 229 |
|
|---|
| 230 | #ifdef __IEEE_BIG_ENDIAN
|
|---|
| 231 |
|
|---|
| 232 | typedef union
|
|---|
| 233 | {
|
|---|
| 234 | double value;
|
|---|
| 235 | struct
|
|---|
| 236 | {
|
|---|
| 237 | uint32_t msw;
|
|---|
| 238 | uint32_t lsw;
|
|---|
| 239 | } parts;
|
|---|
| 240 | } ieee_double_shape_type;
|
|---|
| 241 |
|
|---|
| 242 | #endif
|
|---|
| 243 |
|
|---|
| 244 | #ifdef __IEEE_LITTLE_ENDIAN
|
|---|
| 245 |
|
|---|
| 246 | typedef union
|
|---|
| 247 | {
|
|---|
| 248 | double value;
|
|---|
| 249 | struct
|
|---|
| 250 | {
|
|---|
| 251 | uint32_t lsw;
|
|---|
| 252 | uint32_t msw;
|
|---|
| 253 | } parts;
|
|---|
| 254 | } ieee_double_shape_type;
|
|---|
| 255 |
|
|---|
| 256 | #endif
|
|---|
| 257 |
|
|---|
| 258 | /* Get two 32 bit ints from a double. */
|
|---|
| 259 |
|
|---|
| 260 | #define EXTRACT_WORDS(ix0,ix1,d) \
|
|---|
| 261 | do { \
|
|---|
| 262 | ieee_double_shape_type ew_u; \
|
|---|
| 263 | ew_u.value = (d); \
|
|---|
| 264 | (ix0) = ew_u.parts.msw; \
|
|---|
| 265 | (ix1) = ew_u.parts.lsw; \
|
|---|
| 266 | } while (0)
|
|---|
| 267 |
|
|---|
| 268 | /* Get the more significant 32 bit int from a double. */
|
|---|
| 269 |
|
|---|
| 270 | #define GET_HIGH_WORD(i,d) \
|
|---|
| 271 | do { \
|
|---|
| 272 | ieee_double_shape_type gh_u; \
|
|---|
| 273 | gh_u.value = (d); \
|
|---|
| 274 | (i) = gh_u.parts.msw; \
|
|---|
| 275 | } while (0)
|
|---|
| 276 |
|
|---|
| 277 | /* Get the less significant 32 bit int from a double. */
|
|---|
| 278 |
|
|---|
| 279 | #define GET_LOW_WORD(i,d) \
|
|---|
| 280 | do { \
|
|---|
| 281 | ieee_double_shape_type gl_u; \
|
|---|
| 282 | gl_u.value = (d); \
|
|---|
| 283 | (i) = gl_u.parts.lsw; \
|
|---|
| 284 | } while (0)
|
|---|
| 285 |
|
|---|
| 286 | /* Set a double from two 32 bit ints. */
|
|---|
| 287 |
|
|---|
| 288 | #define INSERT_WORDS(d,ix0,ix1) \
|
|---|
| 289 | do { \
|
|---|
| 290 | ieee_double_shape_type iw_u; \
|
|---|
| 291 | iw_u.parts.msw = (ix0); \
|
|---|
| 292 | iw_u.parts.lsw = (ix1); \
|
|---|
| 293 | (d) = iw_u.value; \
|
|---|
| 294 | } while (0)
|
|---|
| 295 |
|
|---|
| 296 | /* Set the more significant 32 bits of a double from an int. */
|
|---|
| 297 |
|
|---|
| 298 | #define SET_HIGH_WORD(d,v) \
|
|---|
| 299 | do { \
|
|---|
| 300 | ieee_double_shape_type sh_u; \
|
|---|
| 301 | sh_u.value = (d); \
|
|---|
| 302 | sh_u.parts.msw = (v); \
|
|---|
| 303 | (d) = sh_u.value; \
|
|---|
| 304 | } while (0)
|
|---|
| 305 |
|
|---|
| 306 | /* Set the less significant 32 bits of a double from an int. */
|
|---|
| 307 |
|
|---|
| 308 | #define SET_LOW_WORD(d,v) \
|
|---|
| 309 | do { \
|
|---|
| 310 | ieee_double_shape_type sl_u; \
|
|---|
| 311 | sl_u.value = (d); \
|
|---|
| 312 | sl_u.parts.lsw = (v); \
|
|---|
| 313 | (d) = sl_u.value; \
|
|---|
| 314 | } while (0)
|
|---|
| 315 |
|
|---|
| 316 | /* A union which permits us to convert between a float and a 32 bit
|
|---|
| 317 | int. */
|
|---|
| 318 |
|
|---|
| 319 | typedef union
|
|---|
| 320 | {
|
|---|
| 321 | float value;
|
|---|
| 322 | uint32_t word;
|
|---|
| 323 | } ieee_float_shape_type;
|
|---|
| 324 |
|
|---|
| 325 | /* Get a 32 bit int from a float. */
|
|---|
| 326 |
|
|---|
| 327 | #define GET_FLOAT_WORD(i,d) \
|
|---|
| 328 | do { \
|
|---|
| 329 | ieee_float_shape_type gf_u; \
|
|---|
| 330 | gf_u.value = (d); \
|
|---|
| 331 | (i) = gf_u.word; \
|
|---|
| 332 | } while (0)
|
|---|
| 333 |
|
|---|
| 334 | /* Set a float from a 32 bit int. */
|
|---|
| 335 |
|
|---|
| 336 | #define SET_FLOAT_WORD(d,i) \
|
|---|
| 337 | do { \
|
|---|
| 338 | ieee_float_shape_type sf_u; \
|
|---|
| 339 | sf_u.word = (i); \
|
|---|
| 340 | (d) = sf_u.value; \
|
|---|
| 341 | } while (0)
|
|---|
| 342 |
|
|---|
| 343 | #ifdef __cplusplus
|
|---|
| 344 | }
|
|---|
| 345 | #endif
|
|---|
| 346 |
|
|---|
