1 /* Native implementation of soft float functions */

     2 #include <math.h>

     3 

     4 #if (defined(_BSD) && !defined(__APPLE__)) || defined(HOST_SOLARIS)

     5 #include <ieeefp.h>

     6 #define fabsf(f) ((float)fabs(f))

     7 #else

     8 #include <fenv.h>

     9 #endif

    10 

    11 #if defined(__OpenBSD__) || defined(__NetBSD__)

    12 #include <sys/param.h>

    13 #endif

    14 

    15 /*

    16  * Define some C99-7.12.3 classification macros and

    17  *        some C99-.12.4 for Solaris systems OS less than 10,

    18  *        or Solaris 10 systems running GCC 3.x or less.

    19  *   Solaris 10 with GCC4 does not need these macros as they

    20  *   are defined in <iso/math_c99.h> with a compiler directive

    21  */

    22 #if defined(HOST_SOLARIS) && (( HOST_SOLARIS <= 9 ) || ((HOST_SOLARIS >= 10) \

    23                                                         && (__GNUC__ <= 4))) \

    24     || (defined(__OpenBSD__) && (OpenBSD < 200811))

    25 /*

    26  * C99 7.12.3 classification macros

    27  * and

    28  * C99 7.12.14 comparison macros

    29  *

    30  * ... do not work on Solaris 10 using GNU CC 3.4.x.

    31  * Try to workaround the missing / broken C99 math macros.

    32  */

    33 #if defined(__OpenBSD__)

    34 #define unordered(x, y) (isnan(x) || isnan(y))

    35 #endif

    36 

    37 #ifdef __NetBSD__

    38 #ifndef isgreater

    39 #define isgreater(x, y)		__builtin_isgreater(x, y)

    40 #endif

    41 #ifndef isgreaterequal

    42 #define isgreaterequal(x, y)	__builtin_isgreaterequal(x, y)

    43 #endif

    44 #ifndef isless

    45 #define isless(x, y)		__builtin_isless(x, y)

    46 #endif

    47 #ifndef islessequal

    48 #define islessequal(x, y)	__builtin_islessequal(x, y)

    49 #endif

    50 #ifndef isunordered

    51 #define isunordered(x, y)	__builtin_isunordered(x, y)

    52 #endif

    53 #endif

    54 

    55 

    56 #define isnormal(x)             (fpclass(x) >= FP_NZERO)

    57 #define isgreater(x, y)         ((!unordered(x, y)) && ((x) > (y)))

    58 #define isgreaterequal(x, y)    ((!unordered(x, y)) && ((x) >= (y)))

    59 #define isless(x, y)            ((!unordered(x, y)) && ((x) < (y)))

    60 #define islessequal(x, y)       ((!unordered(x, y)) && ((x) <= (y)))

    61 #define isunordered(x,y)        unordered(x, y)

    62 #endif

    63 

    64 #if defined(__sun__) && !defined(NEED_LIBSUNMATH)

    65 

    66 #ifndef isnan

    67 # define isnan(x) \

    68     (sizeof (x) == sizeof (long double) ? isnan_ld (x) \

    69      : sizeof (x) == sizeof (double) ? isnan_d (x) \

    70      : isnan_f (x))

    71 static inline int isnan_f  (float       x) { return x != x; }

    72 static inline int isnan_d  (double      x) { return x != x; }

    73 static inline int isnan_ld (long double x) { return x != x; }

    74 #endif

    75 

    76 #ifndef isinf

    77 # define isinf(x) \

    78     (sizeof (x) == sizeof (long double) ? isinf_ld (x) \

    79      : sizeof (x) == sizeof (double) ? isinf_d (x) \

    80      : isinf_f (x))

    81 static inline int isinf_f  (float       x) { return isnan (x - x); }

    82 static inline int isinf_d  (double      x) { return isnan (x - x); }

    83 static inline int isinf_ld (long double x) { return isnan (x - x); }

    84 #endif

    85 #endif

    86 

    87 typedef float float32;

    88 typedef double float64;

    89 #ifdef FLOATX80

    90 typedef long double floatx80;

    91 #endif

    92 

    93 typedef union {

    94     float32 f;

    95     uint32_t i;

    96 } float32u;

    97 typedef union {

    98     float64 f;

    99     uint64_t i;

   100 } float64u;

   101 #ifdef FLOATX80

   102 typedef union {

   103     floatx80 f;

   104     struct {

   105         uint64_t low;

   106         uint16_t high;

   107     } i;

   108 } floatx80u;

   109 #endif

   110 

   111 /*----------------------------------------------------------------------------

   112 | Software IEC/IEEE floating-point rounding mode.

   113 *----------------------------------------------------------------------------*/

   114 #if (defined(_BSD) && !defined(__APPLE__)) || defined(HOST_SOLARIS)

   115 #if defined(__OpenBSD__)

   116 #define FE_RM FP_RM

   117 #define FE_RP FP_RP

   118 #define FE_RZ FP_RZ

   119 #endif

   120 enum {

   121     float_round_nearest_even = FP_RN,

   122     float_round_down         = FP_RM,

   123     float_round_up           = FP_RP,

   124     float_round_to_zero      = FP_RZ

   125 };

   126 #elif defined(__arm__)

   127 enum {

   128     float_round_nearest_even = 0,

   129     float_round_down         = 1,

   130     float_round_up           = 2,

   131     float_round_to_zero      = 3

   132 };

   133 #else

   134 enum {

   135     float_round_nearest_even = FE_TONEAREST,

   136     float_round_down         = FE_DOWNWARD,

   137     float_round_up           = FE_UPWARD,

   138     float_round_to_zero      = FE_TOWARDZERO

   139 };

   140 #endif

   141 

   142 typedef struct float_status {

   143     int float_rounding_mode;

   144 #ifdef FLOATX80

   145     int floatx80_rounding_precision;

   146 #endif

   147 } float_status;

   148 

   149 void set_float_rounding_mode(int val STATUS_PARAM);

   150 #ifdef FLOATX80

   151 void set_floatx80_rounding_precision(int val STATUS_PARAM);

   152 #endif

   153 

   154 /*----------------------------------------------------------------------------

   155 | Software IEC/IEEE integer-to-floating-point conversion routines.

   156 *----------------------------------------------------------------------------*/

   157 float32 int32_to_float32( int STATUS_PARAM);

   158 float32 uint32_to_float32( unsigned int STATUS_PARAM);

   159 float64 int32_to_float64( int STATUS_PARAM);

   160 float64 uint32_to_float64( unsigned int STATUS_PARAM);

   161 #ifdef FLOATX80

   162 floatx80 int32_to_floatx80( int STATUS_PARAM);

   163 #endif

   164 #ifdef FLOAT128

   165 float128 int32_to_float128( int STATUS_PARAM);

   166 #endif

   167 float32 int64_to_float32( int64_t STATUS_PARAM);

   168 float32 uint64_to_float32( uint64_t STATUS_PARAM);

   169 float64 int64_to_float64( int64_t STATUS_PARAM);

   170 float64 uint64_to_float64( uint64_t v STATUS_PARAM);

   171 #ifdef FLOATX80

   172 floatx80 int64_to_floatx80( int64_t STATUS_PARAM);

   173 #endif

   174 #ifdef FLOAT128

   175 float128 int64_to_float128( int64_t STATUS_PARAM);

   176 #endif

   177 

   178 /*----------------------------------------------------------------------------

   179 | Software IEC/IEEE single-precision conversion routines.

   180 *----------------------------------------------------------------------------*/

   181 int float32_to_int32( float32  STATUS_PARAM);

   182 int float32_to_int32_round_to_zero( float32  STATUS_PARAM);

   183 unsigned int float32_to_uint32( float32 a STATUS_PARAM);

   184 unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM);

   185 int64_t float32_to_int64( float32  STATUS_PARAM);

   186 int64_t float32_to_int64_round_to_zero( float32  STATUS_PARAM);

   187 float64 float32_to_float64( float32  STATUS_PARAM);

   188 #ifdef FLOATX80

   189 floatx80 float32_to_floatx80( float32  STATUS_PARAM);

   190 #endif

   191 #ifdef FLOAT128

   192 float128 float32_to_float128( float32  STATUS_PARAM);

   193 #endif

   194 

   195 /*----------------------------------------------------------------------------

   196 | Software IEC/IEEE single-precision operations.

   197 *----------------------------------------------------------------------------*/

   198 float32 float32_round_to_int( float32  STATUS_PARAM);

   199 INLINE float32 float32_add( float32 a, float32 b STATUS_PARAM)

   200 {

   201     return a + b;

   202 }

   203 INLINE float32 float32_sub( float32 a, float32 b STATUS_PARAM)

   204 {

   205     return a - b;

   206 }

   207 INLINE float32 float32_mul( float32 a, float32 b STATUS_PARAM)

   208 {

   209     return a * b;

   210 }

   211 INLINE float32 float32_div( float32 a, float32 b STATUS_PARAM)

   212 {

   213     return a / b;

   214 }

   215 float32 float32_rem( float32, float32  STATUS_PARAM);

   216 float32 float32_sqrt( float32  STATUS_PARAM);

   217 INLINE int float32_eq( float32 a, float32 b STATUS_PARAM)

   218 {

   219     return a == b;

   220 }

   221 INLINE int float32_le( float32 a, float32 b STATUS_PARAM)

   222 {

   223     return a <= b;

   224 }

   225 INLINE int float32_lt( float32 a, float32 b STATUS_PARAM)

   226 {

   227     return a < b;

   228 }

   229 INLINE int float32_eq_signaling( float32 a, float32 b STATUS_PARAM)

   230 {

   231     return a <= b && a >= b;

   232 }

   233 INLINE int float32_le_quiet( float32 a, float32 b STATUS_PARAM)

   234 {

   235     return islessequal(a, b);

   236 }

   237 INLINE int float32_lt_quiet( float32 a, float32 b STATUS_PARAM)

   238 {

   239     return isless(a, b);

   240 }

   241 INLINE int float32_unordered( float32 a, float32 b STATUS_PARAM)

   242 {

   243     return isunordered(a, b);

   244 

   245 }

   246 int float32_compare( float32, float32 STATUS_PARAM );

   247 int float32_compare_quiet( float32, float32 STATUS_PARAM );

   248 int float32_is_signaling_nan( float32 );

   249 int float32_is_nan( float32 );

   250 

   251 INLINE float32 float32_abs(float32 a)

   252 {

   253     return fabsf(a);

   254 }

   255 

   256 INLINE float32 float32_chs(float32 a)

   257 {

   258     return -a;

   259 }

   260 

   261 INLINE float32 float32_is_infinity(float32 a)

   262 {

   263     return fpclassify(a) == FP_INFINITE;

   264 }

   265 

   266 INLINE float32 float32_is_neg(float32 a)

   267 {

   268     float32u u;

   269     u.f = a;

   270     return u.i >> 31;

   271 }

   272 

   273 INLINE float32 float32_is_zero(float32 a)

   274 {

   275     return fpclassify(a) == FP_ZERO;

   276 }

   277 

   278 INLINE float32 float32_scalbn(float32 a, int n)

   279 {

   280     return scalbnf(a, n);

   281 }

   282 

   283 /*----------------------------------------------------------------------------

   284 | Software IEC/IEEE double-precision conversion routines.

   285 *----------------------------------------------------------------------------*/

   286 int float64_to_int32( float64 STATUS_PARAM );

   287 int float64_to_int32_round_to_zero( float64 STATUS_PARAM );

   288 unsigned int float64_to_uint32( float64 STATUS_PARAM );

   289 unsigned int float64_to_uint32_round_to_zero( float64 STATUS_PARAM );

   290 int64_t float64_to_int64( float64 STATUS_PARAM );

   291 int64_t float64_to_int64_round_to_zero( float64 STATUS_PARAM );

   292 uint64_t float64_to_uint64( float64 STATUS_PARAM );

   293 uint64_t float64_to_uint64_round_to_zero( float64 STATUS_PARAM );

   294 float32 float64_to_float32( float64 STATUS_PARAM );

   295 #ifdef FLOATX80

   296 floatx80 float64_to_floatx80( float64 STATUS_PARAM );

   297 #endif

   298 #ifdef FLOAT128

   299 float128 float64_to_float128( float64 STATUS_PARAM );

   300 #endif

   301 

   302 /*----------------------------------------------------------------------------

   303 | Software IEC/IEEE double-precision operations.

   304 *----------------------------------------------------------------------------*/

   305 float64 float64_round_to_int( float64 STATUS_PARAM );

   306 float64 float64_trunc_to_int( float64 STATUS_PARAM );

   307 INLINE float64 float64_add( float64 a, float64 b STATUS_PARAM)

   308 {

   309     return a + b;

   310 }

   311 INLINE float64 float64_sub( float64 a, float64 b STATUS_PARAM)

   312 {

   313     return a - b;

   314 }

   315 INLINE float64 float64_mul( float64 a, float64 b STATUS_PARAM)

   316 {

   317     return a * b;

   318 }

   319 INLINE float64 float64_div( float64 a, float64 b STATUS_PARAM)

   320 {

   321     return a / b;

   322 }

   323 float64 float64_rem( float64, float64 STATUS_PARAM );

   324 float64 float64_sqrt( float64 STATUS_PARAM );

   325 INLINE int float64_eq( float64 a, float64 b STATUS_PARAM)

   326 {

   327     return a == b;

   328 }

   329 INLINE int float64_le( float64 a, float64 b STATUS_PARAM)

   330 {

   331     return a <= b;

   332 }

   333 INLINE int float64_lt( float64 a, float64 b STATUS_PARAM)

   334 {

   335     return a < b;

   336 }

   337 INLINE int float64_eq_signaling( float64 a, float64 b STATUS_PARAM)

   338 {

   339     return a <= b && a >= b;

   340 }

   341 INLINE int float64_le_quiet( float64 a, float64 b STATUS_PARAM)

   342 {

   343     return islessequal(a, b);

   344 }

   345 INLINE int float64_lt_quiet( float64 a, float64 b STATUS_PARAM)

   346 {

   347     return isless(a, b);

   348 

   349 }

   350 INLINE int float64_unordered( float64 a, float64 b STATUS_PARAM)

   351 {

   352     return isunordered(a, b);

   353 

   354 }

   355 int float64_compare( float64, float64 STATUS_PARAM );

   356 int float64_compare_quiet( float64, float64 STATUS_PARAM );

   357 int float64_is_signaling_nan( float64 );

   358 int float64_is_nan( float64 );

   359 

   360 INLINE float64 float64_abs(float64 a)

   361 {

   362     return fabs(a);

   363 }

   364 

   365 INLINE float64 float64_chs(float64 a)

   366 {

   367     return -a;

   368 }

   369 

   370 INLINE float64 float64_is_infinity(float64 a)

   371 {

   372     return fpclassify(a) == FP_INFINITE;

   373 }

   374 

   375 INLINE float64 float64_is_neg(float64 a)

   376 {

   377     float64u u;

   378     u.f = a;

   379     return u.i >> 63;

   380 }

   381 

   382 INLINE float64 float64_is_zero(float64 a)

   383 {

   384     return fpclassify(a) == FP_ZERO;

   385 }

   386 

   387 INLINE float64 float64_scalbn(float64 a, int n)

   388 {

   389     return scalbn(a, n);

   390 }

   391 

   392 #ifdef FLOATX80

   393 

   394 /*----------------------------------------------------------------------------

   395 | Software IEC/IEEE extended double-precision conversion routines.

   396 *----------------------------------------------------------------------------*/

   397 int floatx80_to_int32( floatx80 STATUS_PARAM );

   398 int floatx80_to_int32_round_to_zero( floatx80 STATUS_PARAM );

   399 int64_t floatx80_to_int64( floatx80 STATUS_PARAM);

   400 int64_t floatx80_to_int64_round_to_zero( floatx80 STATUS_PARAM);

   401 float32 floatx80_to_float32( floatx80 STATUS_PARAM );

   402 float64 floatx80_to_float64( floatx80 STATUS_PARAM );

   403 #ifdef FLOAT128

   404 float128 floatx80_to_float128( floatx80 STATUS_PARAM );

   405 #endif

   406 

   407 /*----------------------------------------------------------------------------

   408 | Software IEC/IEEE extended double-precision operations.

   409 *----------------------------------------------------------------------------*/

   410 floatx80 floatx80_round_to_int( floatx80 STATUS_PARAM );

   411 INLINE floatx80 floatx80_add( floatx80 a, floatx80 b STATUS_PARAM)

   412 {

   413     return a + b;

   414 }

   415 INLINE floatx80 floatx80_sub( floatx80 a, floatx80 b STATUS_PARAM)

   416 {

   417     return a - b;

   418 }

   419 INLINE floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM)

   420 {

   421     return a * b;

   422 }

   423 INLINE floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM)

   424 {

   425     return a / b;

   426 }

   427 floatx80 floatx80_rem( floatx80, floatx80 STATUS_PARAM );

   428 floatx80 floatx80_sqrt( floatx80 STATUS_PARAM );

   429 INLINE int floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM)

   430 {

   431     return a == b;

   432 }

   433 INLINE int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM)

   434 {

   435     return a <= b;

   436 }

   437 INLINE int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM)

   438 {

   439     return a < b;

   440 }

   441 INLINE int floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM)

   442 {

   443     return a <= b && a >= b;

   444 }

   445 INLINE int floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM)

   446 {

   447     return islessequal(a, b);

   448 }

   449 INLINE int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM)

   450 {

   451     return isless(a, b);

   452 

   453 }

   454 INLINE int floatx80_unordered( floatx80 a, floatx80 b STATUS_PARAM)

   455 {

   456     return isunordered(a, b);

   457 

   458 }

   459 int floatx80_compare( floatx80, floatx80 STATUS_PARAM );

   460 int floatx80_compare_quiet( floatx80, floatx80 STATUS_PARAM );

   461 int floatx80_is_signaling_nan( floatx80 );

   462 int floatx80_is_nan( floatx80 );

   463 

   464 INLINE floatx80 floatx80_abs(floatx80 a)

   465 {

   466     return fabsl(a);

   467 }

   468 

   469 INLINE floatx80 floatx80_chs(floatx80 a)

   470 {

   471     return -a;

   472 }

   473 

   474 INLINE floatx80 floatx80_is_infinity(floatx80 a)

   475 {

   476     return fpclassify(a) == FP_INFINITE;

   477 }

   478 

   479 INLINE floatx80 floatx80_is_neg(floatx80 a)

   480 {

   481     floatx80u u;

   482     u.f = a;

   483     return u.i.high >> 15;

   484 }

   485 

   486 INLINE floatx80 floatx80_is_zero(floatx80 a)

   487 {

   488     return fpclassify(a) == FP_ZERO;

   489 }

   490 

   491 INLINE floatx80 floatx80_scalbn(floatx80 a, int n)

   492 {

   493     return scalbnl(a, n);

   494 }

   495 

   496 #endif