Home | History | Annotate | Download | only in fpu
      1 /* Native implementation of soft float functions */
      2 #include <math.h>
      3 
      4 #if (defined(CONFIG_BSD) && !defined(__APPLE__) && !defined(__GLIBC__)) \
      5     || defined(CONFIG_SOLARIS)
      6 #include <ieeefp.h>
      7 #define fabsf(f) ((float)fabs(f))
      8 #else
      9 #include <fenv.h>
     10 #endif
     11 
     12 #if defined(__OpenBSD__) || defined(__NetBSD__)
     13 #include <sys/param.h>
     14 #endif
     15 
     16 /*
     17  * Define some C99-7.12.3 classification macros and
     18  *        some C99-.12.4 for Solaris systems OS less than 10,
     19  *        or Solaris 10 systems running GCC 3.x or less.
     20  *   Solaris 10 with GCC4 does not need these macros as they
     21  *   are defined in <iso/math_c99.h> with a compiler directive
     22  */
     23 #if defined(CONFIG_SOLARIS) && \
     24            ((CONFIG_SOLARIS_VERSION <= 9 ) || \
     25            ((CONFIG_SOLARIS_VERSION == 10) && (__GNUC__ < 4))) \
     26     || (defined(__OpenBSD__) && (OpenBSD < 200811))
     27 /*
     28  * C99 7.12.3 classification macros
     29  * and
     30  * C99 7.12.14 comparison macros
     31  *
     32  * ... do not work on Solaris 10 using GNU CC 3.4.x.
     33  * Try to workaround the missing / broken C99 math macros.
     34  */
     35 #if defined(__OpenBSD__)
     36 #define unordered(x, y) (isnan(x) || isnan(y))
     37 #endif
     38 
     39 #ifdef __NetBSD__
     40 #ifndef isgreater
     41 #define isgreater(x, y)		__builtin_isgreater(x, y)
     42 #endif
     43 #ifndef isgreaterequal
     44 #define isgreaterequal(x, y)	__builtin_isgreaterequal(x, y)
     45 #endif
     46 #ifndef isless
     47 #define isless(x, y)		__builtin_isless(x, y)
     48 #endif
     49 #ifndef islessequal
     50 #define islessequal(x, y)	__builtin_islessequal(x, y)
     51 #endif
     52 #ifndef isunordered
     53 #define isunordered(x, y)	__builtin_isunordered(x, y)
     54 #endif
     55 #endif
     56 
     57 
     58 #define isnormal(x)             (fpclass(x) >= FP_NZERO)
     59 #define isgreater(x, y)         ((!unordered(x, y)) && ((x) > (y)))
     60 #define isgreaterequal(x, y)    ((!unordered(x, y)) && ((x) >= (y)))
     61 #define isless(x, y)            ((!unordered(x, y)) && ((x) < (y)))
     62 #define islessequal(x, y)       ((!unordered(x, y)) && ((x) <= (y)))
     63 #define isunordered(x,y)        unordered(x, y)
     64 #endif
     65 
     66 #if defined(__sun__) && !defined(CONFIG_NEEDS_LIBSUNMATH)
     67 
     68 #ifndef isnan
     69 # define isnan(x) \
     70     (sizeof (x) == sizeof (long double) ? isnan_ld (x) \
     71      : sizeof (x) == sizeof (double) ? isnan_d (x) \
     72      : isnan_f (x))
     73 static inline int isnan_f  (float       x) { return x != x; }
     74 static inline int isnan_d  (double      x) { return x != x; }
     75 static inline int isnan_ld (long double x) { return x != x; }
     76 #endif
     77 
     78 #ifndef isinf
     79 # define isinf(x) \
     80     (sizeof (x) == sizeof (long double) ? isinf_ld (x) \
     81      : sizeof (x) == sizeof (double) ? isinf_d (x) \
     82      : isinf_f (x))
     83 static inline int isinf_f  (float       x) { return isnan (x - x); }
     84 static inline int isinf_d  (double      x) { return isnan (x - x); }
     85 static inline int isinf_ld (long double x) { return isnan (x - x); }
     86 #endif
     87 #endif
     88 
     89 typedef float float32;
     90 typedef double float64;
     91 #ifdef FLOATX80
     92 typedef long double floatx80;
     93 #endif
     94 
     95 typedef union {
     96     float32 f;
     97     uint32_t i;
     98 } float32u;
     99 typedef union {
    100     float64 f;
    101     uint64_t i;
    102 } float64u;
    103 #ifdef FLOATX80
    104 typedef union {
    105     floatx80 f;
    106     struct {
    107         uint64_t low;
    108         uint16_t high;
    109     } i;
    110 } floatx80u;
    111 #endif
    112 
    113 /*----------------------------------------------------------------------------
    114 | Software IEC/IEEE floating-point rounding mode.
    115 *----------------------------------------------------------------------------*/
    116 #if (defined(CONFIG_BSD) && !defined(__APPLE__) && !defined(__GLIBC__)) \
    117     || defined(CONFIG_SOLARIS)
    118 #if defined(__OpenBSD__)
    119 #define FE_RM FP_RM
    120 #define FE_RP FP_RP
    121 #define FE_RZ FP_RZ
    122 #endif
    123 enum {
    124     float_round_nearest_even = FP_RN,
    125     float_round_down         = FP_RM,
    126     float_round_up           = FP_RP,
    127     float_round_to_zero      = FP_RZ
    128 };
    129 #else
    130 enum {
    131     float_round_nearest_even = FE_TONEAREST,
    132     float_round_down         = FE_DOWNWARD,
    133     float_round_up           = FE_UPWARD,
    134     float_round_to_zero      = FE_TOWARDZERO
    135 };
    136 #endif
    137 
    138 typedef struct float_status {
    139     int float_rounding_mode;
    140 #ifdef FLOATX80
    141     int floatx80_rounding_precision;
    142 #endif
    143 } float_status;
    144 
    145 void set_float_rounding_mode(int val STATUS_PARAM);
    146 #ifdef FLOATX80
    147 void set_floatx80_rounding_precision(int val STATUS_PARAM);
    148 #endif
    149 
    150 /*----------------------------------------------------------------------------
    151 | Software IEC/IEEE integer-to-floating-point conversion routines.
    152 *----------------------------------------------------------------------------*/
    153 float32 int32_to_float32( int STATUS_PARAM);
    154 float32 uint32_to_float32( unsigned int STATUS_PARAM);
    155 float64 int32_to_float64( int STATUS_PARAM);
    156 float64 uint32_to_float64( unsigned int STATUS_PARAM);
    157 #ifdef FLOATX80
    158 floatx80 int32_to_floatx80( int STATUS_PARAM);
    159 #endif
    160 #ifdef FLOAT128
    161 float128 int32_to_float128( int STATUS_PARAM);
    162 #endif
    163 float32 int64_to_float32( int64_t STATUS_PARAM);
    164 float32 uint64_to_float32( uint64_t STATUS_PARAM);
    165 float64 int64_to_float64( int64_t STATUS_PARAM);
    166 float64 uint64_to_float64( uint64_t v STATUS_PARAM);
    167 #ifdef FLOATX80
    168 floatx80 int64_to_floatx80( int64_t STATUS_PARAM);
    169 #endif
    170 #ifdef FLOAT128
    171 float128 int64_to_float128( int64_t STATUS_PARAM);
    172 #endif
    173 
    174 /*----------------------------------------------------------------------------
    175 | Software IEC/IEEE single-precision conversion routines.
    176 *----------------------------------------------------------------------------*/
    177 int float32_to_int32( float32  STATUS_PARAM);
    178 int float32_to_int32_round_to_zero( float32  STATUS_PARAM);
    179 unsigned int float32_to_uint32( float32 a STATUS_PARAM);
    180 unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM);
    181 int64_t float32_to_int64( float32  STATUS_PARAM);
    182 int64_t float32_to_int64_round_to_zero( float32  STATUS_PARAM);
    183 float64 float32_to_float64( float32  STATUS_PARAM);
    184 #ifdef FLOATX80
    185 floatx80 float32_to_floatx80( float32  STATUS_PARAM);
    186 #endif
    187 #ifdef FLOAT128
    188 float128 float32_to_float128( float32  STATUS_PARAM);
    189 #endif
    190 
    191 /*----------------------------------------------------------------------------
    192 | Software IEC/IEEE single-precision operations.
    193 *----------------------------------------------------------------------------*/
    194 float32 float32_round_to_int( float32  STATUS_PARAM);
    195 INLINE float32 float32_add( float32 a, float32 b STATUS_PARAM)
    196 {
    197     return a + b;
    198 }
    199 INLINE float32 float32_sub( float32 a, float32 b STATUS_PARAM)
    200 {
    201     return a - b;
    202 }
    203 INLINE float32 float32_mul( float32 a, float32 b STATUS_PARAM)
    204 {
    205     return a * b;
    206 }
    207 INLINE float32 float32_div( float32 a, float32 b STATUS_PARAM)
    208 {
    209     return a / b;
    210 }
    211 float32 float32_rem( float32, float32  STATUS_PARAM);
    212 float32 float32_sqrt( float32  STATUS_PARAM);
    213 INLINE int float32_eq( float32 a, float32 b STATUS_PARAM)
    214 {
    215     return a == b;
    216 }
    217 INLINE int float32_le( float32 a, float32 b STATUS_PARAM)
    218 {
    219     return a <= b;
    220 }
    221 INLINE int float32_lt( float32 a, float32 b STATUS_PARAM)
    222 {
    223     return a < b;
    224 }
    225 INLINE int float32_eq_signaling( float32 a, float32 b STATUS_PARAM)
    226 {
    227     return a <= b && a >= b;
    228 }
    229 INLINE int float32_le_quiet( float32 a, float32 b STATUS_PARAM)
    230 {
    231     return islessequal(a, b);
    232 }
    233 INLINE int float32_lt_quiet( float32 a, float32 b STATUS_PARAM)
    234 {
    235     return isless(a, b);
    236 }
    237 INLINE int float32_unordered( float32 a, float32 b STATUS_PARAM)
    238 {
    239     return isunordered(a, b);
    240 
    241 }
    242 int float32_compare( float32, float32 STATUS_PARAM );
    243 int float32_compare_quiet( float32, float32 STATUS_PARAM );
    244 int float32_is_signaling_nan( float32 );
    245 int float32_is_nan( float32 );
    246 
    247 INLINE float32 float32_abs(float32 a)
    248 {
    249     return fabsf(a);
    250 }
    251 
    252 INLINE float32 float32_chs(float32 a)
    253 {
    254     return -a;
    255 }
    256 
    257 INLINE float32 float32_is_infinity(float32 a)
    258 {
    259     return fpclassify(a) == FP_INFINITE;
    260 }
    261 
    262 INLINE float32 float32_is_neg(float32 a)
    263 {
    264     float32u u;
    265     u.f = a;
    266     return u.i >> 31;
    267 }
    268 
    269 INLINE float32 float32_is_zero(float32 a)
    270 {
    271     return fpclassify(a) == FP_ZERO;
    272 }
    273 
    274 INLINE float32 float32_scalbn(float32 a, int n)
    275 {
    276     return scalbnf(a, n);
    277 }
    278 
    279 /*----------------------------------------------------------------------------
    280 | Software IEC/IEEE double-precision conversion routines.
    281 *----------------------------------------------------------------------------*/
    282 int float64_to_int32( float64 STATUS_PARAM );
    283 int float64_to_int32_round_to_zero( float64 STATUS_PARAM );
    284 unsigned int float64_to_uint32( float64 STATUS_PARAM );
    285 unsigned int float64_to_uint32_round_to_zero( float64 STATUS_PARAM );
    286 int64_t float64_to_int64( float64 STATUS_PARAM );
    287 int64_t float64_to_int64_round_to_zero( float64 STATUS_PARAM );
    288 uint64_t float64_to_uint64( float64 STATUS_PARAM );
    289 uint64_t float64_to_uint64_round_to_zero( float64 STATUS_PARAM );
    290 float32 float64_to_float32( float64 STATUS_PARAM );
    291 #ifdef FLOATX80
    292 floatx80 float64_to_floatx80( float64 STATUS_PARAM );
    293 #endif
    294 #ifdef FLOAT128
    295 float128 float64_to_float128( float64 STATUS_PARAM );
    296 #endif
    297 
    298 /*----------------------------------------------------------------------------
    299 | Software IEC/IEEE double-precision operations.
    300 *----------------------------------------------------------------------------*/
    301 float64 float64_round_to_int( float64 STATUS_PARAM );
    302 float64 float64_trunc_to_int( float64 STATUS_PARAM );
    303 INLINE float64 float64_add( float64 a, float64 b STATUS_PARAM)
    304 {
    305     return a + b;
    306 }
    307 INLINE float64 float64_sub( float64 a, float64 b STATUS_PARAM)
    308 {
    309     return a - b;
    310 }
    311 INLINE float64 float64_mul( float64 a, float64 b STATUS_PARAM)
    312 {
    313     return a * b;
    314 }
    315 INLINE float64 float64_div( float64 a, float64 b STATUS_PARAM)
    316 {
    317     return a / b;
    318 }
    319 float64 float64_rem( float64, float64 STATUS_PARAM );
    320 float64 float64_sqrt( float64 STATUS_PARAM );
    321 INLINE int float64_eq( float64 a, float64 b STATUS_PARAM)
    322 {
    323     return a == b;
    324 }
    325 INLINE int float64_le( float64 a, float64 b STATUS_PARAM)
    326 {
    327     return a <= b;
    328 }
    329 INLINE int float64_lt( float64 a, float64 b STATUS_PARAM)
    330 {
    331     return a < b;
    332 }
    333 INLINE int float64_eq_signaling( float64 a, float64 b STATUS_PARAM)
    334 {
    335     return a <= b && a >= b;
    336 }
    337 INLINE int float64_le_quiet( float64 a, float64 b STATUS_PARAM)
    338 {
    339     return islessequal(a, b);
    340 }
    341 INLINE int float64_lt_quiet( float64 a, float64 b STATUS_PARAM)
    342 {
    343     return isless(a, b);
    344 
    345 }
    346 INLINE int float64_unordered( float64 a, float64 b STATUS_PARAM)
    347 {
    348     return isunordered(a, b);
    349 
    350 }
    351 int float64_compare( float64, float64 STATUS_PARAM );
    352 int float64_compare_quiet( float64, float64 STATUS_PARAM );
    353 int float64_is_signaling_nan( float64 );
    354 int float64_is_nan( float64 );
    355 
    356 INLINE float64 float64_abs(float64 a)
    357 {
    358     return fabs(a);
    359 }
    360 
    361 INLINE float64 float64_chs(float64 a)
    362 {
    363     return -a;
    364 }
    365 
    366 INLINE float64 float64_is_infinity(float64 a)
    367 {
    368     return fpclassify(a) == FP_INFINITE;
    369 }
    370 
    371 INLINE float64 float64_is_neg(float64 a)
    372 {
    373     float64u u;
    374     u.f = a;
    375     return u.i >> 63;
    376 }
    377 
    378 INLINE float64 float64_is_zero(float64 a)
    379 {
    380     return fpclassify(a) == FP_ZERO;
    381 }
    382 
    383 INLINE float64 float64_scalbn(float64 a, int n)
    384 {
    385     return scalbn(a, n);
    386 }
    387 
    388 #ifdef FLOATX80
    389 
    390 /*----------------------------------------------------------------------------
    391 | Software IEC/IEEE extended double-precision conversion routines.
    392 *----------------------------------------------------------------------------*/
    393 int floatx80_to_int32( floatx80 STATUS_PARAM );
    394 int floatx80_to_int32_round_to_zero( floatx80 STATUS_PARAM );
    395 int64_t floatx80_to_int64( floatx80 STATUS_PARAM);
    396 int64_t floatx80_to_int64_round_to_zero( floatx80 STATUS_PARAM);
    397 float32 floatx80_to_float32( floatx80 STATUS_PARAM );
    398 float64 floatx80_to_float64( floatx80 STATUS_PARAM );
    399 #ifdef FLOAT128
    400 float128 floatx80_to_float128( floatx80 STATUS_PARAM );
    401 #endif
    402 
    403 /*----------------------------------------------------------------------------
    404 | Software IEC/IEEE extended double-precision operations.
    405 *----------------------------------------------------------------------------*/
    406 floatx80 floatx80_round_to_int( floatx80 STATUS_PARAM );
    407 INLINE floatx80 floatx80_add( floatx80 a, floatx80 b STATUS_PARAM)
    408 {
    409     return a + b;
    410 }
    411 INLINE floatx80 floatx80_sub( floatx80 a, floatx80 b STATUS_PARAM)
    412 {
    413     return a - b;
    414 }
    415 INLINE floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM)
    416 {
    417     return a * b;
    418 }
    419 INLINE floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM)
    420 {
    421     return a / b;
    422 }
    423 floatx80 floatx80_rem( floatx80, floatx80 STATUS_PARAM );
    424 floatx80 floatx80_sqrt( floatx80 STATUS_PARAM );
    425 INLINE int floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM)
    426 {
    427     return a == b;
    428 }
    429 INLINE int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM)
    430 {
    431     return a <= b;
    432 }
    433 INLINE int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM)
    434 {
    435     return a < b;
    436 }
    437 INLINE int floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM)
    438 {
    439     return a <= b && a >= b;
    440 }
    441 INLINE int floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM)
    442 {
    443     return islessequal(a, b);
    444 }
    445 INLINE int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM)
    446 {
    447     return isless(a, b);
    448 
    449 }
    450 INLINE int floatx80_unordered( floatx80 a, floatx80 b STATUS_PARAM)
    451 {
    452     return isunordered(a, b);
    453 
    454 }
    455 int floatx80_compare( floatx80, floatx80 STATUS_PARAM );
    456 int floatx80_compare_quiet( floatx80, floatx80 STATUS_PARAM );
    457 int floatx80_is_signaling_nan( floatx80 );
    458 int floatx80_is_nan( floatx80 );
    459 
    460 INLINE floatx80 floatx80_abs(floatx80 a)
    461 {
    462     return fabsl(a);
    463 }
    464 
    465 INLINE floatx80 floatx80_chs(floatx80 a)
    466 {
    467     return -a;
    468 }
    469 
    470 INLINE floatx80 floatx80_is_infinity(floatx80 a)
    471 {
    472     return fpclassify(a) == FP_INFINITE;
    473 }
    474 
    475 INLINE floatx80 floatx80_is_neg(floatx80 a)
    476 {
    477     floatx80u u;
    478     u.f = a;
    479     return u.i.high >> 15;
    480 }
    481 
    482 INLINE floatx80 floatx80_is_zero(floatx80 a)
    483 {
    484     return fpclassify(a) == FP_ZERO;
    485 }
    486 
    487 INLINE floatx80 floatx80_scalbn(floatx80 a, int n)
    488 {
    489     return scalbnl(a, n);
    490 }
    491 
    492 #endif
    493