diff options
Diffstat (limited to 'sysdeps/powerpc/fpu')
| -rw-r--r-- | sysdeps/powerpc/fpu/fegetround.c | 4 | ||||
| -rw-r--r-- | sysdeps/powerpc/fpu/fenv_libc.h | 37 | ||||
| -rw-r--r-- | sysdeps/powerpc/fpu/fesetround.c | 18 | ||||
| -rw-r--r-- | sysdeps/powerpc/fpu/math_ldbl.h | 189 |
4 files changed, 230 insertions, 18 deletions
diff --git a/sysdeps/powerpc/fpu/fegetround.c b/sysdeps/powerpc/fpu/fegetround.c index ae521fdf10..93663ad546 100644 --- a/sysdeps/powerpc/fpu/fegetround.c +++ b/sysdeps/powerpc/fpu/fegetround.c @@ -23,7 +23,5 @@ int fegetround (void) { - int result; - asm ("mcrfs 7,7 ; mfcr %0" : "=r"(result) : : "cr7"); \ - return result & 3; + return __fegetround(); } diff --git a/sysdeps/powerpc/fpu/fenv_libc.h b/sysdeps/powerpc/fpu/fenv_libc.h index 7ae12a7d2b..fd5fc0c767 100644 --- a/sysdeps/powerpc/fpu/fenv_libc.h +++ b/sysdeps/powerpc/fpu/fenv_libc.h @@ -1,5 +1,5 @@ /* Internal libc stuff for floating point environment routines. - Copyright (C) 1997 Free Software Foundation, Inc. + Copyright (C) 1997, 2006 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -54,6 +54,41 @@ typedef union unsigned int l[2]; } fenv_union_t; + +static inline int +__fegetround (void) +{ + int result; + asm volatile ("mcrfs 7,7\n\t" + "mfcr %0" : "=r"(result) : : "cr7"); + return result & 3; +} +#define fegetround() __fegetround() + +static inline int +__fesetround (int round) +{ + if ((unsigned int) round < 2) + { + asm volatile ("mtfsb0 30"); + if ((unsigned int) round == 0) + asm volatile ("mtfsb0 31"); + else + asm volatile ("mtfsb1 31"); + } + else + { + asm volatile ("mtfsb1 30"); + if ((unsigned int) round == 2) + asm volatile ("mtfsb0 31"); + else + asm volatile ("mtfsb1 31"); + } + + return 0; +} +#define fesetround(mode) __fesetround(mode) + /* Definitions of all the FPSCR bit numbers */ enum { FPSCR_FX = 0, /* exception summary */ diff --git a/sysdeps/powerpc/fpu/fesetround.c b/sysdeps/powerpc/fpu/fesetround.c index 67518d0df4..a7efa3bbb0 100644 --- a/sysdeps/powerpc/fpu/fesetround.c +++ b/sysdeps/powerpc/fpu/fesetround.c @@ -1,5 +1,5 @@ /* Set current rounding direction. - Copyright (C) 1997, 2005 Free Software Foundation, Inc. + Copyright (C) 1997, 2005, 2006 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,23 +20,13 @@ #include <fenv_libc.h> +#undef fesetround int fesetround (int round) { - fenv_union_t u; - if ((unsigned int) round > 3) return 1; - - /* Get the current state. */ - u.fenv = fegetenv_register (); - - /* Set the relevant bits. */ - u.l[1] = (u.l[1] & ~3) | (round & 3); - - /* Put the new state in effect. */ - fesetenv_register (u.fenv); - - return 0; + else + return __fesetround(round); } libm_hidden_def (fesetround) diff --git a/sysdeps/powerpc/fpu/math_ldbl.h b/sysdeps/powerpc/fpu/math_ldbl.h new file mode 100644 index 0000000000..6cd6d0bdfe --- /dev/null +++ b/sysdeps/powerpc/fpu/math_ldbl.h @@ -0,0 +1,189 @@ +#ifndef _MATH_PRIVATE_H_ +#error "Never use <math_ldbl.h> directly; include <math_private.h> instead." +#endif + +#include <sysdeps/ieee754/ldbl-128/math_ldbl.h> +#include <ieee754.h> + +static inline void +ldbl_extract_mantissa (int64_t *hi64, u_int64_t *lo64, int *exp, long double x) +{ + /* We have 105 bits of mantissa plus one implicit digit. Since + 106 bits are representable we use the first implicit digit for + the number before the decimal point and the second implicit bit + as bit 53 of the mantissa. */ + unsigned long long hi, lo; + int ediff; + union ibm_extended_long_double eldbl; + eldbl.d = x; + *exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS; + + lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3; + hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1; + /* If the lower double is not a denomal or zero then set the hidden + 53rd bit. */ + if (eldbl.ieee.exponent2 > 0x001) + { + lo |= (1ULL << 52); + lo = lo << 7; /* pre-shift lo to match ieee854. */ + /* The lower double is normalized separately from the upper. We + may need to adjust the lower manitissa to reflect this. */ + ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2; + if (ediff > 53) + lo = lo >> (ediff-53); + } + hi |= (1ULL << 52); + + if ((eldbl.ieee.negative != eldbl.ieee.negative2) + && ((eldbl.ieee.exponent2 != 0) && (lo != 0LL))) + { + hi--; + lo = (1ULL << 60) - lo; + if (hi < (1ULL << 52)) + { + /* we have a borrow from the hidden bit, so shift left 1. */ + hi = (hi << 1) | (lo >> 59); + lo = 0xfffffffffffffffLL & (lo << 1); + *exp = *exp - 1; + } + } + *lo64 = (hi << 60) | lo; + *hi64 = hi >> 4; +} + +static inline long double +ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64) +{ + union ibm_extended_long_double u; + unsigned long hidden2, lzcount; + unsigned long long hi, lo; + + u.ieee.negative = sign; + u.ieee.negative2 = sign; + u.ieee.exponent = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS; + u.ieee.exponent2 = exp-53 + IBM_EXTENDED_LONG_DOUBLE_BIAS; + /* Expect 113 bits (112 bits + hidden) right justified in two longs. + The low order 53 bits (52 + hidden) go into the lower double */ + lo = (lo64 >> 7)& ((1ULL << 53) - 1); + hidden2 = (lo64 >> 59) & 1ULL; + /* The high order 53 bits (52 + hidden) go into the upper double */ + hi = (lo64 >> 60) & ((1ULL << 11) - 1); + hi |= (hi64 << 4); + + if (lo != 0LL) + { + /* hidden2 bit of low double controls rounding of the high double. + If hidden2 is '1' then round up hi and adjust lo (2nd mantissa) + plus change the sign of the low double to compensate. */ + if (hidden2) + { + hi++; + u.ieee.negative2 = !sign; + lo = (1ULL << 53) - lo; + } + /* The hidden bit of the lo mantissa is zero so we need to + normalize the it for the low double. Shift it left until the + hidden bit is '1' then adjust the 2nd exponent accordingly. */ + + if (sizeof (lo) == sizeof (long)) + lzcount = __builtin_clzl (lo); + else if ((lo >> 32) != 0) + lzcount = __builtin_clzl ((long) (lo >> 32)); + else + lzcount = __builtin_clzl ((long) lo) + 32; + lzcount = lzcount - 11; + if (lzcount > 0) + { + int expnt2 = u.ieee.exponent2 - lzcount; + if (expnt2 >= 1) + { + /* Not denormal. Normalize and set low exponent. */ + lo = lo << lzcount; + u.ieee.exponent2 = expnt2; + } + else + { + /* Is denormal. */ + lo = lo << (lzcount + expnt2); + u.ieee.exponent2 = 0; + } + } + } + else + { + u.ieee.negative2 = 0; + u.ieee.exponent2 = 0; + } + + u.ieee.mantissa3 = lo & ((1ULL << 32) - 1); + u.ieee.mantissa2 = (lo >> 32) & ((1ULL << 20) - 1); + u.ieee.mantissa1 = hi & ((1ULL << 32) - 1); + u.ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1); + return u.d; +} + +/* gcc generates disgusting code to pack and unpack long doubles. + This tells gcc that pack/unpack is really a nop. We use fr1/fr2 + because those are the regs used to pass/return a single + long double arg. */ +static inline long double +ldbl_pack (double a, double aa) +{ + register long double x __asm__ ("fr1"); + register double xh __asm__ ("fr1"); + register double xl __asm__ ("fr2"); + xh = a; + xl = aa; + __asm__ ("" : "=f" (x) : "f" (xh), "f" (xl)); + return x; +} + +static inline void +ldbl_unpack (long double l, double *a, double *aa) +{ + register long double x __asm__ ("fr1"); + register double xh __asm__ ("fr1"); + register double xl __asm__ ("fr2"); + x = l; + __asm__ ("" : "=f" (xh), "=f" (xl) : "f" (x)); + *a = xh; + *aa = xl; +} + + +/* Convert a finite long double to canonical form. + Does not handle +/-Inf properly. */ +static inline void +ldbl_canonicalize (double *a, double *aa) +{ + double xh, xl; + + xh = *a + *aa; + xl = (*a - xh) + *aa; + *a = xh; + *aa = xl; +} + +/* Simple inline nearbyint (double) function . + Only works in the default rounding mode + but is useful in long double rounding functions. */ +static inline double +ldbl_nearbyint (double a) +{ + double two52 = 0x10000000000000LL; + + if (__builtin_expect ((__builtin_fabs (a) < two52), 1)) + { + if (__builtin_expect ((a > 0.0), 1)) + { + a += two52; + a -= two52; + } + else if (__builtin_expect ((a < 0.0), 1)) + { + a = two52 - a; + a = -(a - two52); + } + } + return a; +} |