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)) 73static inline int isnan_f (float x) { return x != x; } 74static inline int isnan_d (double x) { return x != x; } 75static 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)) 83static inline int isinf_f (float x) { return isnan (x - x); } 84static inline int isinf_d (double x) { return isnan (x - x); } 85static inline int isinf_ld (long double x) { return isnan (x - x); } 86#endif 87#endif 88 89typedef float float32; 90typedef double float64; 91#ifdef FLOATX80 92typedef long double floatx80; 93#endif 94 95typedef union { 96 float32 f; 97 uint32_t i; 98} float32u; 99typedef union { 100 float64 f; 101 uint64_t i; 102} float64u; 103#ifdef FLOATX80 104typedef 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 123enum { 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 130enum { 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 138typedef struct float_status { 139 int float_rounding_mode; 140#ifdef FLOATX80 141 int floatx80_rounding_precision; 142#endif 143} float_status; 144 145void set_float_rounding_mode(int val STATUS_PARAM); 146#ifdef FLOATX80 147void set_floatx80_rounding_precision(int val STATUS_PARAM); 148#endif 149 150/*---------------------------------------------------------------------------- 151| Software IEC/IEEE integer-to-floating-point conversion routines. 152*----------------------------------------------------------------------------*/ 153float32 int32_to_float32( int STATUS_PARAM); 154float32 uint32_to_float32( unsigned int STATUS_PARAM); 155float64 int32_to_float64( int STATUS_PARAM); 156float64 uint32_to_float64( unsigned int STATUS_PARAM); 157#ifdef FLOATX80 158floatx80 int32_to_floatx80( int STATUS_PARAM); 159#endif 160#ifdef FLOAT128 161float128 int32_to_float128( int STATUS_PARAM); 162#endif 163float32 int64_to_float32( int64_t STATUS_PARAM); 164float32 uint64_to_float32( uint64_t STATUS_PARAM); 165float64 int64_to_float64( int64_t STATUS_PARAM); 166float64 uint64_to_float64( uint64_t v STATUS_PARAM); 167#ifdef FLOATX80 168floatx80 int64_to_floatx80( int64_t STATUS_PARAM); 169#endif 170#ifdef FLOAT128 171float128 int64_to_float128( int64_t STATUS_PARAM); 172#endif 173 174/*---------------------------------------------------------------------------- 175| Software IEC/IEEE single-precision conversion constants. 176*----------------------------------------------------------------------------*/ 177#define float32_zero (0.0) 178#define float32_one (1.0) 179#define float32_ln2 (0.6931471) 180#define float32_pi (3.1415926) 181#define float32_half (0.5) 182 183/*---------------------------------------------------------------------------- 184| Software IEC/IEEE single-precision conversion routines. 185*----------------------------------------------------------------------------*/ 186int float32_to_int32( float32 STATUS_PARAM); 187int float32_to_int32_round_to_zero( float32 STATUS_PARAM); 188unsigned int float32_to_uint32( float32 a STATUS_PARAM); 189unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM); 190int64_t float32_to_int64( float32 STATUS_PARAM); 191int64_t float32_to_int64_round_to_zero( float32 STATUS_PARAM); 192float64 float32_to_float64( float32 STATUS_PARAM); 193#ifdef FLOATX80 194floatx80 float32_to_floatx80( float32 STATUS_PARAM); 195#endif 196#ifdef FLOAT128 197float128 float32_to_float128( float32 STATUS_PARAM); 198#endif 199 200/*---------------------------------------------------------------------------- 201| Software IEC/IEEE single-precision operations. 202*----------------------------------------------------------------------------*/ 203float32 float32_round_to_int( float32 STATUS_PARAM); 204INLINE float32 float32_add( float32 a, float32 b STATUS_PARAM) 205{ 206 return a + b; 207} 208INLINE float32 float32_sub( float32 a, float32 b STATUS_PARAM) 209{ 210 return a - b; 211} 212INLINE float32 float32_mul( float32 a, float32 b STATUS_PARAM) 213{ 214 return a * b; 215} 216INLINE float32 float32_div( float32 a, float32 b STATUS_PARAM) 217{ 218 return a / b; 219} 220float32 float32_rem( float32, float32 STATUS_PARAM); 221float32 float32_sqrt( float32 STATUS_PARAM); 222INLINE int float32_eq_quiet( float32 a, float32 b STATUS_PARAM) 223{ 224 return a == b; 225} 226INLINE int float32_le( float32 a, float32 b STATUS_PARAM) 227{ 228 return a <= b; 229} 230INLINE int float32_lt( float32 a, float32 b STATUS_PARAM) 231{ 232 return a < b; 233} 234INLINE int float32_eq( float32 a, float32 b STATUS_PARAM) 235{ 236 return a <= b && a >= b; 237} 238INLINE int float32_le_quiet( float32 a, float32 b STATUS_PARAM) 239{ 240 return islessequal(a, b); 241} 242INLINE int float32_lt_quiet( float32 a, float32 b STATUS_PARAM) 243{ 244 return isless(a, b); 245} 246INLINE int float32_unordered( float32 a, float32 b STATUS_PARAM) 247{ 248 return isunordered(a, b); 249} 250INLINE int float32_unordered_quiet( float32 a, float32 b STATUS_PARAM) 251{ 252 return isunordered(a, b); 253} 254int float32_compare( float32, float32 STATUS_PARAM ); 255int float32_compare_quiet( float32, float32 STATUS_PARAM ); 256int float32_is_signaling_nan( float32 ); 257int float32_is_quiet_nan( float32 ); 258int float32_is_any_nan( float32 ); 259 260INLINE float32 float32_abs(float32 a) 261{ 262 return fabsf(a); 263} 264 265INLINE float32 float32_chs(float32 a) 266{ 267 return -a; 268} 269 270INLINE float32 float32_is_infinity(float32 a) 271{ 272 return fpclassify(a) == FP_INFINITE; 273} 274 275INLINE float32 float32_is_neg(float32 a) 276{ 277 float32u u; 278 u.f = a; 279 return u.i >> 31; 280} 281 282INLINE float32 float32_is_zero(float32 a) 283{ 284 return fpclassify(a) == FP_ZERO; 285} 286 287INLINE float32 float32_scalbn(float32 a, int n STATUS_PARAM) 288{ 289 return scalbnf(a, n); 290} 291 292/*---------------------------------------------------------------------------- 293| Software IEC/IEEE double-precision conversion constants. 294*----------------------------------------------------------------------------*/ 295#define float64_zero (0.0) 296#define float64_one (1.0) 297#define float64_ln2 (0.693147180559945) 298#define float64_pi (3.141592653589793) 299#define float64_half (0.5) 300 301/*---------------------------------------------------------------------------- 302| Software IEC/IEEE double-precision conversion routines. 303*----------------------------------------------------------------------------*/ 304int float64_to_int32( float64 STATUS_PARAM ); 305int float64_to_int32_round_to_zero( float64 STATUS_PARAM ); 306unsigned int float64_to_uint32( float64 STATUS_PARAM ); 307unsigned int float64_to_uint32_round_to_zero( float64 STATUS_PARAM ); 308int64_t float64_to_int64( float64 STATUS_PARAM ); 309int64_t float64_to_int64_round_to_zero( float64 STATUS_PARAM ); 310uint64_t float64_to_uint64( float64 STATUS_PARAM ); 311uint64_t float64_to_uint64_round_to_zero( float64 STATUS_PARAM ); 312float32 float64_to_float32( float64 STATUS_PARAM ); 313#ifdef FLOATX80 314floatx80 float64_to_floatx80( float64 STATUS_PARAM ); 315#endif 316#ifdef FLOAT128 317float128 float64_to_float128( float64 STATUS_PARAM ); 318#endif 319 320/*---------------------------------------------------------------------------- 321| Software IEC/IEEE double-precision operations. 322*----------------------------------------------------------------------------*/ 323float64 float64_round_to_int( float64 STATUS_PARAM ); 324float64 float64_trunc_to_int( float64 STATUS_PARAM ); 325INLINE float64 float64_add( float64 a, float64 b STATUS_PARAM) 326{ 327 return a + b; 328} 329INLINE float64 float64_sub( float64 a, float64 b STATUS_PARAM) 330{ 331 return a - b; 332} 333INLINE float64 float64_mul( float64 a, float64 b STATUS_PARAM) 334{ 335 return a * b; 336} 337INLINE float64 float64_div( float64 a, float64 b STATUS_PARAM) 338{ 339 return a / b; 340} 341float64 float64_rem( float64, float64 STATUS_PARAM ); 342float64 float64_sqrt( float64 STATUS_PARAM ); 343INLINE int float64_eq_quiet( float64 a, float64 b STATUS_PARAM) 344{ 345 return a == b; 346} 347INLINE int float64_le( float64 a, float64 b STATUS_PARAM) 348{ 349 return a <= b; 350} 351INLINE int float64_lt( float64 a, float64 b STATUS_PARAM) 352{ 353 return a < b; 354} 355INLINE int float64_eq( float64 a, float64 b STATUS_PARAM) 356{ 357 return a <= b && a >= b; 358} 359INLINE int float64_le_quiet( float64 a, float64 b STATUS_PARAM) 360{ 361 return islessequal(a, b); 362} 363INLINE int float64_lt_quiet( float64 a, float64 b STATUS_PARAM) 364{ 365 return isless(a, b); 366 367} 368INLINE int float64_unordered( float64 a, float64 b STATUS_PARAM) 369{ 370 return isunordered(a, b); 371} 372INLINE int float64_unordered_quiet( float64 a, float64 b STATUS_PARAM) 373{ 374 return isunordered(a, b); 375} 376int float64_compare( float64, float64 STATUS_PARAM ); 377int float64_compare_quiet( float64, float64 STATUS_PARAM ); 378int float64_is_signaling_nan( float64 ); 379int float64_is_any_nan( float64 ); 380int float64_is_quiet_nan( float64 ); 381 382INLINE float64 float64_abs(float64 a) 383{ 384 return fabs(a); 385} 386 387INLINE float64 float64_chs(float64 a) 388{ 389 return -a; 390} 391 392INLINE float64 float64_is_infinity(float64 a) 393{ 394 return fpclassify(a) == FP_INFINITE; 395} 396 397INLINE float64 float64_is_neg(float64 a) 398{ 399 float64u u; 400 u.f = a; 401 return u.i >> 63; 402} 403 404INLINE float64 float64_is_zero(float64 a) 405{ 406 return fpclassify(a) == FP_ZERO; 407} 408 409INLINE float64 float64_scalbn(float64 a, int n STATUS_PARAM) 410{ 411 return scalbn(a, n); 412} 413 414#ifdef FLOATX80 415 416/*---------------------------------------------------------------------------- 417| Software IEC/IEEE extended double-precision conversion constants. 418*----------------------------------------------------------------------------*/ 419#define floatx80_zero (0.0L) 420#define floatx80_one (1.0L) 421#define floatx80_ln2 (0.69314718055994530943L) 422#define floatx80_pi (3.14159265358979323851L) 423#define floatx80_half (0.5L) 424 425/*---------------------------------------------------------------------------- 426| Software IEC/IEEE extended double-precision conversion routines. 427*----------------------------------------------------------------------------*/ 428int floatx80_to_int32( floatx80 STATUS_PARAM ); 429int floatx80_to_int32_round_to_zero( floatx80 STATUS_PARAM ); 430int64_t floatx80_to_int64( floatx80 STATUS_PARAM); 431int64_t floatx80_to_int64_round_to_zero( floatx80 STATUS_PARAM); 432float32 floatx80_to_float32( floatx80 STATUS_PARAM ); 433float64 floatx80_to_float64( floatx80 STATUS_PARAM ); 434#ifdef FLOAT128 435float128 floatx80_to_float128( floatx80 STATUS_PARAM ); 436#endif 437 438/*---------------------------------------------------------------------------- 439| Software IEC/IEEE extended double-precision operations. 440*----------------------------------------------------------------------------*/ 441floatx80 floatx80_round_to_int( floatx80 STATUS_PARAM ); 442INLINE floatx80 floatx80_add( floatx80 a, floatx80 b STATUS_PARAM) 443{ 444 return a + b; 445} 446INLINE floatx80 floatx80_sub( floatx80 a, floatx80 b STATUS_PARAM) 447{ 448 return a - b; 449} 450INLINE floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM) 451{ 452 return a * b; 453} 454INLINE floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM) 455{ 456 return a / b; 457} 458floatx80 floatx80_rem( floatx80, floatx80 STATUS_PARAM ); 459floatx80 floatx80_sqrt( floatx80 STATUS_PARAM ); 460INLINE int floatx80_eq_quiet( floatx80 a, floatx80 b STATUS_PARAM) 461{ 462 return a == b; 463} 464INLINE int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM) 465{ 466 return a <= b; 467} 468INLINE int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM) 469{ 470 return a < b; 471} 472INLINE int floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM) 473{ 474 return a <= b && a >= b; 475} 476INLINE int floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM) 477{ 478 return islessequal(a, b); 479} 480INLINE int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM) 481{ 482 return isless(a, b); 483 484} 485INLINE int floatx80_unordered( floatx80 a, floatx80 b STATUS_PARAM) 486{ 487 return isunordered(a, b); 488} 489INLINE int floatx80_unordered_quiet( floatx80 a, floatx80 b STATUS_PARAM) 490{ 491 return isunordered(a, b); 492} 493int floatx80_compare( floatx80, floatx80 STATUS_PARAM ); 494int floatx80_compare_quiet( floatx80, floatx80 STATUS_PARAM ); 495int floatx80_is_signaling_nan( floatx80 ); 496int floatx80_is_quiet_nan( floatx80 ); 497int floatx80_is_any_nan( floatx80 ); 498 499INLINE floatx80 floatx80_abs(floatx80 a) 500{ 501 return fabsl(a); 502} 503 504INLINE floatx80 floatx80_chs(floatx80 a) 505{ 506 return -a; 507} 508 509INLINE floatx80 floatx80_is_infinity(floatx80 a) 510{ 511 return fpclassify(a) == FP_INFINITE; 512} 513 514INLINE floatx80 floatx80_is_neg(floatx80 a) 515{ 516 floatx80u u; 517 u.f = a; 518 return u.i.high >> 15; 519} 520 521INLINE floatx80 floatx80_is_zero(floatx80 a) 522{ 523 return fpclassify(a) == FP_ZERO; 524} 525 526INLINE floatx80 floatx80_scalbn(floatx80 a, int n STATUS_PARAM) 527{ 528 return scalbnl(a, n); 529} 530 531#endif 532