Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | RSS feed

  1.  
  2. /* @(#)fdlibm.h 5.1 93/09/24 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 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. /* REDHAT LOCAL: Include files.  */
  15. #include <math.h>
  16. #include <sys/types.h>
  17. #include <machine/ieeefp.h>
  18.  
  19. /* REDHAT LOCAL: Default to XOPEN_MODE.  */
  20. #define _XOPEN_MODE
  21.  
  22. /* Most routines need to check whether a float is finite, infinite, or not a
  23.    number, and many need to know whether the result of an operation will
  24.    overflow.  These conditions depend on whether the largest exponent is
  25.    used for NaNs & infinities, or whether it's used for finite numbers.  The
  26.    macros below wrap up that kind of information:
  27.  
  28.    FLT_UWORD_IS_FINITE(X)
  29.         True if a positive float with bitmask X is finite.
  30.  
  31.    FLT_UWORD_IS_NAN(X)
  32.         True if a positive float with bitmask X is not a number.
  33.  
  34.    FLT_UWORD_IS_INFINITE(X)
  35.         True if a positive float with bitmask X is +infinity.
  36.  
  37.    FLT_UWORD_MAX
  38.         The bitmask of FLT_MAX.
  39.  
  40.    FLT_UWORD_HALF_MAX
  41.         The bitmask of FLT_MAX/2.
  42.  
  43.    FLT_UWORD_EXP_MAX
  44.         The bitmask of the largest finite exponent (129 if the largest
  45.         exponent is used for finite numbers, 128 otherwise).
  46.  
  47.    FLT_UWORD_LOG_MAX
  48.         The bitmask of log(FLT_MAX), rounded down.  This value is the largest
  49.         input that can be passed to exp() without producing overflow.
  50.  
  51.    FLT_UWORD_LOG_2MAX
  52.         The bitmask of log(2*FLT_MAX), rounded down.  This value is the
  53.         largest input than can be passed to cosh() without producing
  54.         overflow.
  55.  
  56.    FLT_LARGEST_EXP
  57.         The largest biased exponent that can be used for finite numbers
  58.         (255 if the largest exponent is used for finite numbers, 254
  59.         otherwise) */
  60.  
  61. #ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL
  62. #define FLT_UWORD_IS_FINITE(x) 1
  63. #define FLT_UWORD_IS_NAN(x) 0
  64. #define FLT_UWORD_IS_INFINITE(x) 0
  65. #define FLT_UWORD_MAX 0x7fffffff
  66. #define FLT_UWORD_EXP_MAX 0x43010000
  67. #define FLT_UWORD_LOG_MAX 0x42b2d4fc
  68. #define FLT_UWORD_LOG_2MAX 0x42b437e0
  69. #define HUGE ((float)0X1.FFFFFEP128)
  70. #else
  71. #define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)
  72. #define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)
  73. #define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)
  74. #define FLT_UWORD_MAX 0x7f7fffffL
  75. #define FLT_UWORD_EXP_MAX 0x43000000
  76. #define FLT_UWORD_LOG_MAX 0x42b17217
  77. #define FLT_UWORD_LOG_2MAX 0x42b2d4fc
  78. #define HUGE ((float)3.40282346638528860e+38)
  79. #endif
  80. #define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23))
  81. #define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)
  82.  
  83. /* Many routines check for zero and subnormal numbers.  Such things depend
  84.    on whether the target supports denormals or not:
  85.  
  86.    FLT_UWORD_IS_ZERO(X)
  87.         True if a positive float with bitmask X is +0.  Without denormals,
  88.         any float with a zero exponent is a +0 representation.  With
  89.         denormals, the only +0 representation is a 0 bitmask.
  90.  
  91.    FLT_UWORD_IS_SUBNORMAL(X)
  92.         True if a non-zero positive float with bitmask X is subnormal.
  93.         (Routines should check for zeros first.)
  94.  
  95.    FLT_UWORD_MIN
  96.         The bitmask of the smallest float above +0.  Call this number
  97.         REAL_FLT_MIN...
  98.  
  99.    FLT_UWORD_EXP_MIN
  100.         The bitmask of the float representation of REAL_FLT_MIN's exponent.
  101.  
  102.    FLT_UWORD_LOG_MIN
  103.         The bitmask of |log(REAL_FLT_MIN)|, rounding down.
  104.  
  105.    FLT_SMALLEST_EXP
  106.         REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
  107.         -22 if they are).
  108. */
  109.  
  110. #ifdef _FLT_NO_DENORMALS
  111. #define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)
  112. #define FLT_UWORD_IS_SUBNORMAL(x) 0
  113. #define FLT_UWORD_MIN 0x00800000
  114. #define FLT_UWORD_EXP_MIN 0x42fc0000
  115. #define FLT_UWORD_LOG_MIN 0x42aeac50
  116. #define FLT_SMALLEST_EXP 1
  117. #else
  118. #define FLT_UWORD_IS_ZERO(x) ((x)==0)
  119. #define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)
  120. #define FLT_UWORD_MIN 0x00000001
  121. #define FLT_UWORD_EXP_MIN 0x43160000
  122. #define FLT_UWORD_LOG_MIN 0x42cff1b5
  123. #define FLT_SMALLEST_EXP -22
  124. #endif
  125.  
  126. #ifdef __STDC__
  127. #undef __P
  128. #define __P(p)  p
  129. #else
  130. #define __P(p)  ()
  131. #endif
  132.  
  133. /*
  134.  * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
  135.  * (one may replace the following line by "#include <values.h>")
  136.  */
  137.  
  138. #define X_TLOSS         1.41484755040568800000e+16
  139.  
  140. /* Functions that are not documented, and are not in <math.h>.  */
  141.  
  142. #ifdef _SCALB_INT
  143. extern double scalb __P((double, int));
  144. #else
  145. extern double scalb __P((double, double));
  146. #endif
  147. extern double significand __P((double));
  148.  
  149. /* ieee style elementary functions */
  150. extern double __ieee754_sqrt __P((double));                    
  151. extern double __ieee754_acos __P((double));                    
  152. extern double __ieee754_acosh __P((double));                   
  153. extern double __ieee754_log __P((double));                     
  154. extern double __ieee754_atanh __P((double));                   
  155. extern double __ieee754_asin __P((double));                    
  156. extern double __ieee754_atan2 __P((double,double));                    
  157. extern double __ieee754_exp __P((double));
  158. extern double __ieee754_cosh __P((double));
  159. extern double __ieee754_fmod __P((double,double));
  160. extern double __ieee754_pow __P((double,double));
  161. extern double __ieee754_lgamma_r __P((double,int *));
  162. extern double __ieee754_gamma_r __P((double,int *));
  163. extern double __ieee754_log10 __P((double));
  164. extern double __ieee754_sinh __P((double));
  165. extern double __ieee754_hypot __P((double,double));
  166. extern double __ieee754_j0 __P((double));
  167. extern double __ieee754_j1 __P((double));
  168. extern double __ieee754_y0 __P((double));
  169. extern double __ieee754_y1 __P((double));
  170. extern double __ieee754_jn __P((int,double));
  171. extern double __ieee754_yn __P((int,double));
  172. extern double __ieee754_remainder __P((double,double));
  173. extern __int32_t __ieee754_rem_pio2 __P((double,double*));
  174. #ifdef _SCALB_INT
  175. extern double __ieee754_scalb __P((double,int));
  176. #else
  177. extern double __ieee754_scalb __P((double,double));
  178. #endif
  179.  
  180. /* fdlibm kernel function */
  181. extern double __kernel_standard __P((double,double,int));
  182. extern double __kernel_sin __P((double,double,int));
  183. extern double __kernel_cos __P((double,double));
  184. extern double __kernel_tan __P((double,double,int));
  185. extern int    __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));
  186.  
  187. /* Undocumented float functions.  */
  188. #ifdef _SCALB_INT
  189. extern float scalbf __P((float, int));
  190. #else
  191. extern float scalbf __P((float, float));
  192. #endif
  193. extern float significandf __P((float));
  194.  
  195. /* ieee style elementary float functions */
  196. extern float __ieee754_sqrtf __P((float));                     
  197. extern float __ieee754_acosf __P((float));                     
  198. extern float __ieee754_acoshf __P((float));                    
  199. extern float __ieee754_logf __P((float));                      
  200. extern float __ieee754_atanhf __P((float));                    
  201. extern float __ieee754_asinf __P((float));                     
  202. extern float __ieee754_atan2f __P((float,float));                      
  203. extern float __ieee754_expf __P((float));
  204. extern float __ieee754_coshf __P((float));
  205. extern float __ieee754_fmodf __P((float,float));
  206. extern float __ieee754_powf __P((float,float));
  207. extern float __ieee754_lgammaf_r __P((float,int *));
  208. extern float __ieee754_gammaf_r __P((float,int *));
  209. extern float __ieee754_log10f __P((float));
  210. extern float __ieee754_sinhf __P((float));
  211. extern float __ieee754_hypotf __P((float,float));
  212. extern float __ieee754_j0f __P((float));
  213. extern float __ieee754_j1f __P((float));
  214. extern float __ieee754_y0f __P((float));
  215. extern float __ieee754_y1f __P((float));
  216. extern float __ieee754_jnf __P((int,float));
  217. extern float __ieee754_ynf __P((int,float));
  218. extern float __ieee754_remainderf __P((float,float));
  219. extern __int32_t __ieee754_rem_pio2f __P((float,float*));
  220. #ifdef _SCALB_INT
  221. extern float __ieee754_scalbf __P((float,int));
  222. #else
  223. extern float __ieee754_scalbf __P((float,float));
  224. #endif
  225.  
  226. /* float versions of fdlibm kernel functions */
  227. extern float __kernel_sinf __P((float,float,int));
  228. extern float __kernel_cosf __P((float,float));
  229. extern float __kernel_tanf __P((float,float,int));
  230. extern int   __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));
  231.  
  232. /* The original code used statements like
  233.         n0 = ((*(int*)&one)>>29)^1;             * index of high word *
  234.         ix0 = *(n0+(int*)&x);                   * high word of x *
  235.         ix1 = *((1-n0)+(int*)&x);               * low word of x *
  236.    to dig two 32 bit words out of the 64 bit IEEE floating point
  237.    value.  That is non-ANSI, and, moreover, the gcc instruction
  238.    scheduler gets it wrong.  We instead use the following macros.
  239.    Unlike the original code, we determine the endianness at compile
  240.    time, not at run time; I don't see much benefit to selecting
  241.    endianness at run time.  */
  242.  
  243. #ifndef __IEEE_BIG_ENDIAN
  244. #ifndef __IEEE_LITTLE_ENDIAN
  245.  #error Must define endianness
  246. #endif
  247. #endif
  248.  
  249. /* A union which permits us to convert between a double and two 32 bit
  250.    ints.  */
  251.  
  252. #ifdef __IEEE_BIG_ENDIAN
  253.  
  254. typedef union
  255. {
  256.   double value;
  257.   struct
  258.   {
  259.     __uint32_t msw;
  260.     __uint32_t lsw;
  261.   } parts;
  262. } ieee_double_shape_type;
  263.  
  264. #endif
  265.  
  266. #ifdef __IEEE_LITTLE_ENDIAN
  267.  
  268. typedef union
  269. {
  270.   double value;
  271.   struct
  272.   {
  273.     __uint32_t lsw;
  274.     __uint32_t msw;
  275.   } parts;
  276. } ieee_double_shape_type;
  277.  
  278. #endif
  279.  
  280. /* Get two 32 bit ints from a double.  */
  281.  
  282. #define EXTRACT_WORDS(ix0,ix1,d)                                \
  283. do {                                                            \
  284.   ieee_double_shape_type ew_u;                                  \
  285.   ew_u.value = (d);                                             \
  286.   (ix0) = ew_u.parts.msw;                                       \
  287.   (ix1) = ew_u.parts.lsw;                                       \
  288. } while (0)
  289.  
  290. /* Get the more significant 32 bit int from a double.  */
  291.  
  292. #define GET_HIGH_WORD(i,d)                                      \
  293. do {                                                            \
  294.   ieee_double_shape_type gh_u;                                  \
  295.   gh_u.value = (d);                                             \
  296.   (i) = gh_u.parts.msw;                                         \
  297. } while (0)
  298.  
  299. /* Get the less significant 32 bit int from a double.  */
  300.  
  301. #define GET_LOW_WORD(i,d)                                       \
  302. do {                                                            \
  303.   ieee_double_shape_type gl_u;                                  \
  304.   gl_u.value = (d);                                             \
  305.   (i) = gl_u.parts.lsw;                                         \
  306. } while (0)
  307.  
  308. /* Set a double from two 32 bit ints.  */
  309.  
  310. #define INSERT_WORDS(d,ix0,ix1)                                 \
  311. do {                                                            \
  312.   ieee_double_shape_type iw_u;                                  \
  313.   iw_u.parts.msw = (ix0);                                       \
  314.   iw_u.parts.lsw = (ix1);                                       \
  315.   (d) = iw_u.value;                                             \
  316. } while (0)
  317.  
  318. /* Set the more significant 32 bits of a double from an int.  */
  319.  
  320. #define SET_HIGH_WORD(d,v)                                      \
  321. do {                                                            \
  322.   ieee_double_shape_type sh_u;                                  \
  323.   sh_u.value = (d);                                             \
  324.   sh_u.parts.msw = (v);                                         \
  325.   (d) = sh_u.value;                                             \
  326. } while (0)
  327.  
  328. /* Set the less significant 32 bits of a double from an int.  */
  329.  
  330. #define SET_LOW_WORD(d,v)                                       \
  331. do {                                                            \
  332.   ieee_double_shape_type sl_u;                                  \
  333.   sl_u.value = (d);                                             \
  334.   sl_u.parts.lsw = (v);                                         \
  335.   (d) = sl_u.value;                                             \
  336. } while (0)
  337.  
  338. /* A union which permits us to convert between a float and a 32 bit
  339.    int.  */
  340.  
  341. typedef union
  342. {
  343.   float value;
  344.   __uint32_t word;
  345. } ieee_float_shape_type;
  346.  
  347. /* Get a 32 bit int from a float.  */
  348.  
  349. #define GET_FLOAT_WORD(i,d)                                     \
  350. do {                                                            \
  351.   ieee_float_shape_type gf_u;                                   \
  352.   gf_u.value = (d);                                             \
  353.   (i) = gf_u.word;                                              \
  354. } while (0)
  355.  
  356. /* Set a float from a 32 bit int.  */
  357.  
  358. #define SET_FLOAT_WORD(d,i)                                     \
  359. do {                                                            \
  360.   ieee_float_shape_type sf_u;                                   \
  361.   sf_u.word = (i);                                              \
  362.   (d) = sf_u.value;                                             \
  363. } while (0)
  364.  
  365. /* Macros to avoid undefined behaviour that can arise if the amount
  366.    of a shift is exactly equal to the size of the shifted operand.  */
  367.  
  368. #define SAFE_LEFT_SHIFT(op,amt)                                 \
  369.   (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0)
  370.  
  371. #define SAFE_RIGHT_SHIFT(op,amt)                                \
  372.   (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0)
  373.  
  374. #ifdef  _COMPLEX_H
  375.  
  376. /*
  377.  * Quoting from ISO/IEC 9899:TC2:
  378.  *
  379.  * 6.2.5.13 Types
  380.  * Each complex type has the same representation and alignment requirements as
  381.  * an array type containing exactly two elements of the corresponding real type;
  382.  * the first element is equal to the real part, and the second element to the
  383.  * imaginary part, of the complex number.
  384.  */
  385. typedef union {
  386.         float complex z;
  387.         float parts[2];
  388. } float_complex;
  389.  
  390. typedef union {
  391.         double complex z;
  392.         double parts[2];
  393. } double_complex;
  394.  
  395. typedef union {
  396.         long double complex z;
  397.         long double parts[2];
  398. } long_double_complex;
  399.  
  400. #define REAL_PART(z)    ((z).parts[0])
  401. #define IMAG_PART(z)    ((z).parts[1])
  402.  
  403. #endif  /* _COMPLEX_H */
  404.  
  405.