summaryrefslogtreecommitdiff
path: root/libm/powerpc/classic/s_round.c
diff options
context:
space:
mode:
authorJoakim Tjernlund <joakim.tjernlund@transmode.se>2007-03-31 13:28:15 +0000
committerJoakim Tjernlund <joakim.tjernlund@transmode.se>2007-03-31 13:28:15 +0000
commite7bcf43b6440ac9fc61a0eef5591393810daafb5 (patch)
treeb72c3fb15e030b47b2eb02d13169b4548382e855 /libm/powerpc/classic/s_round.c
parent7a40ba19c86e4d2fc7e35f14a0e629ee843b96a9 (diff)
From Steve Papacharalambous:
Add math support for PowerPC e500.
Diffstat (limited to 'libm/powerpc/classic/s_round.c')
-rw-r--r--libm/powerpc/classic/s_round.c115
1 files changed, 115 insertions, 0 deletions
diff --git a/libm/powerpc/classic/s_round.c b/libm/powerpc/classic/s_round.c
new file mode 100644
index 000000000..62d5936d9
--- /dev/null
+++ b/libm/powerpc/classic/s_round.c
@@ -0,0 +1,115 @@
+#include <limits.h>
+#include <math.h>
+#include <endian.h>
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __BIG_ENDIAN)
+ unsigned long int hi;
+ unsigned long int lo;
+#else
+ unsigned long int lo;
+ unsigned long int hi;
+#endif
+ } words;
+ double dbl;
+ } DblInHex;
+
+static const unsigned long int signMask = 0x80000000ul;
+static const double twoTo52 = 4503599627370496.0;
+
+/*******************************************************************************
+* *
+* The function round rounds its double argument to integral value *
+* according to the "add half to the magnitude and truncate" rounding of *
+* Pascal's Round function and FORTRAN's ANINT function and returns the *
+* result in double format. This function signals inexact if an ordered *
+* return value is not equal to the operand. *
+* *
+*******************************************************************************/
+
+libm_hidden_proto(round)
+double round ( double x )
+ {
+ DblInHex argument, OldEnvironment;
+ register double y, z;
+ register unsigned long int xHead;
+ register long int target;
+
+ argument.dbl = x;
+ xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x|
+ target = ( argument.words.hi < signMask ); // flag positive sign
+
+ if ( xHead < 0x43300000ul )
+/*******************************************************************************
+* Is |x| < 2.0^52? *
+*******************************************************************************/
+ {
+ if ( xHead < 0x3ff00000ul )
+/*******************************************************************************
+* Is |x| < 1.0? *
+*******************************************************************************/
+ {
+ asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
+ if ( xHead < 0x3fe00000ul )
+/*******************************************************************************
+* Is |x| < 0.5? *
+*******************************************************************************/
+ {
+ if ( ( xHead | argument.words.lo ) != 0ul )
+ OldEnvironment.words.lo |= 0x02000000ul;
+ asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+ if ( target )
+ return ( 0.0 );
+ else
+ return ( -0.0 );
+ }
+/*******************************************************************************
+* Is 0.5 ² |x| < 1.0? *
+*******************************************************************************/
+ OldEnvironment.words.lo |= 0x02000000ul;
+ asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+ if ( target )
+ return ( 1.0 );
+ else
+ return ( -1.0 );
+ }
+/*******************************************************************************
+* Is 1.0 < |x| < 2.0^52? *
+*******************************************************************************/
+ if ( target )
+ { // positive x
+ y = ( x + twoTo52 ) - twoTo52; // round at binary point
+ if ( y == x ) // exact case
+ return ( x );
+ z = x + 0.5; // inexact case
+ y = ( z + twoTo52 ) - twoTo52; // round at binary point
+ if ( y > z )
+ return ( y - 1.0 );
+ else
+ return ( y );
+ }
+
+/*******************************************************************************
+* Is x < 0? *
+*******************************************************************************/
+ else
+ {
+ y = ( x - twoTo52 ) + twoTo52; // round at binary point
+ if ( y == x )
+ return ( x );
+ z = x - 0.5;
+ y = ( z - twoTo52 ) + twoTo52; // round at binary point
+ if ( y < z )
+ return ( y + 1.0 );
+ else
+ return ( y );
+ }
+ }
+/*******************************************************************************
+* |x| >= 2.0^52 or x is a NaN. *
+*******************************************************************************/
+ return ( x );
+ }
+libm_hidden_def(round)