libm/powerpc/classic/s_round.c

   1 #include <limits.h>
   2 #include <math.h>
   3 #include <endian.h>
   4
   5 typedef union
   6       {
   7       struct {
   8 #if (__BYTE_ORDER == __BIG_ENDIAN)
   9         unsigned long int hi;
  10         unsigned long int lo;
  11 #else
  12         unsigned long int lo;
  13         unsigned long int hi;
  14 #endif
  15       } words;
  16       double dbl;
  17       } DblInHex;
  18
  19 static const unsigned long int signMask = 0x80000000ul;
  20 static const double twoTo52      = 4503599627370496.0;
  21
  22 /*******************************************************************************
  23 *                                                                              *
  24 *     The function round rounds its double argument to integral value          *
  25 *     according to the "add half to the magnitude and truncate" rounding of    *
  26 *     Pascal's Round function and FORTRAN's ANINT function and returns the     *
  27 *     result in double format.  This function signals inexact if an ordered    *
  28 *     return value is not equal to the operand.                                *
  29 *                                                                              *
  30 *******************************************************************************/
  31
  32 libm_hidden_proto(round)
  33 double round ( double x )
  34       {
  35       DblInHex argument, OldEnvironment;
  36       register double y, z;
  37       register unsigned long int xHead;
  38       register long int target;
  39
  40       argument.dbl = x;
  41       xHead = argument.words.hi & 0x7fffffffUL;      // xHead <- high half of |x|
  42       target = ( argument.words.hi < signMask );     // flag positive sign
  43
  44       if ( xHead < 0x43300000ul )
  45 /*******************************************************************************
  46 *     Is |x| < 2.0^52?                                                        *
  47 *******************************************************************************/
  48             {
  49             if ( xHead < 0x3ff00000ul )
  50 /*******************************************************************************
  51 *     Is |x| < 1.0?                                                           *
  52 *******************************************************************************/
  53                   {
  54                         asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
  55                   if ( xHead < 0x3fe00000ul )
  56 /*******************************************************************************
  57 *     Is |x| < 0.5?                                                           *
  58 *******************************************************************************/
  59                         {
  60                         if ( ( xHead | argument.words.lo ) != 0ul )
  61                               OldEnvironment.words.lo |= 0x02000000ul;
  62                                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  63                         if ( target )
  64                               return ( 0.0 );
  65                         else
  66                               return ( -0.0 );
  67                         }
  68 /*******************************************************************************
  69 *     Is 0.5 ² |x| < 1.0?                                                      *
  70 *******************************************************************************/
  71                   OldEnvironment.words.lo |= 0x02000000ul;
  72                         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  73                   if ( target )
  74                         return ( 1.0 );
  75                   else
  76                         return ( -1.0 );
  77                   }
  78 /*******************************************************************************
  79 *     Is 1.0 < |x| < 2.0^52?                                                   *
  80 *******************************************************************************/
  81             if ( target )
  82                   {                                     // positive x
  83                   y = ( x + twoTo52 ) - twoTo52;        // round at binary point
  84                   if ( y == x )                         // exact case
  85                         return ( x );
  86                   z = x + 0.5;                          // inexact case
  87                   y = ( z + twoTo52 ) - twoTo52;        // round at binary point
  88                   if ( y > z )
  89                         return ( y - 1.0 );
  90                   else
  91                         return ( y );
  92                   }
  93
  94 /*******************************************************************************
  95 *     Is x < 0?                                                                *
  96 *******************************************************************************/
  97             else
  98                   {
  99                   y = ( x - twoTo52 ) + twoTo52;        // round at binary point
 100                   if ( y == x )
 101                         return ( x );
 102                   z = x - 0.5;
 103                   y = ( z - twoTo52 ) + twoTo52;        // round at binary point
 104                   if ( y < z )
 105                         return ( y + 1.0 );
 106                   else
 107                   return ( y );
 108                   }
 109             }
 110 /*******************************************************************************
 111 *      |x| >= 2.0^52 or x is a NaN.                                            *
 112 *******************************************************************************/
 113       return ( x );
 114       }
 115 libm_hidden_def(round)