base/math/rndint.c

Go to the documentation of this file.
00001 /*******************************************************************************
00002 **      File:   rndint.c
00003 **      
00004 **      Contains: C source code for implementations of floating-point
00005 **                functions which round to integral value or format, as
00006 **                defined in header <fp.h>.  In particular, this file
00007 **                contains implementations of functions rint, nearbyint,
00008 **                rinttol, round, roundtol, trunc, modf and modfl.  This file
00009 **                targets PowerPC or Power platforms.
00010 **                        
00011 **      Written by: A. Sazegari, Apple AltiVec Group
00012 **      Created originally by Jon Okada, Apple Numerics Group
00013 **      
00014 **      Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
00015 **      
00016 **      Change History (most recent first):
00017 **
00018 **      13 Jul 01  ram  replaced --setflm calls with inline assembly
00019 **      03 Mar 01  ali  first port to os x using gcc, added the crucial __setflm
00020 **                      definition.
00021 **              1. removed double_t, put in double for now.
00022 **              2. removed iclass from nearbyint.
00023 **              3. removed wrong comments intrunc.
00024 **              4. 
00025 **      13 May 97  ali  made performance improvements in rint, rinttol, roundtol
00026 **                      and trunc by folding some of the taligent ideas into this
00027 **                      implementation.  nearbyint is faster than the one in taligent,
00028 **                      rint is more elegant, but slower by %30 than the taligent one.
00029 **      09 Apr 97  ali  deleted modfl and deferred to AuxiliaryDD.c
00030 **      15 Sep 94  ali  Major overhaul and performance improvements of all functions.
00031 **      20 Jul 94  PAF  New faster version
00032 **      16 Jul 93  ali  Added the modfl function.
00033 **      18 Feb 93  ali  Changed the return value of fenv functions
00034 **                      feclearexcept and feraiseexcept to their new
00035 **                      NCEG X3J11.1/93-001 definitions.
00036 **      16 Dec 92  JPO  Removed __itrunc implementation to a 
00037 **                      separate file.
00038 **      15 Dec 92  JPO  Added __itrunc implementation and modified
00039 **                      rinttol to include conversion from double
00040 **                      to long int format.  Modified roundtol to
00041 **                      call __itrunc.
00042 **      10 Dec 92  JPO  Added modf (double) implementation.
00043 **      04 Dec 92  JPO  First created.
00044 **                        
00045 *******************************************************************************/
00046 
00047 #include <limits.h>
00048 #include <math.h>
00049 
00050 #if !defined(__ppc__)
00051 #define asm(x)
00052 #endif
00053 
00054 #define      SET_INVALID      0x01000000UL
00055 
00056 typedef union
00057       {
00058       struct {
00059 #if defined(__BIG_ENDIAN__)
00060         unsigned long int hi;
00061         unsigned long int lo;
00062 #else
00063         unsigned long int lo;
00064         unsigned long int hi;
00065 #endif
00066       } words;
00067       double dbl;
00068       } DblInHex;
00069 
00070 static const unsigned long int signMask = 0x80000000ul;
00071 static const double twoTo52      = 4503599627370496.0;
00072 static const double doubleToLong = 4503603922337792.0;              // 2^52
00073 static const DblInHex Huge       = {{ 0x7FF00000, 0x00000000 }};
00074 static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
00075 
00076 /*******************************************************************************
00077 *                                                                              *
00078 *     The function rint rounds its double argument to integral value           *
00079 *     according to the current rounding direction and returns the result in    *
00080 *     double format.  This function signals inexact if an ordered return       * 
00081 *     value is not equal to the operand.                                       *
00082 *                                                                              *
00083 ********************************************************************************
00084 *                                                                              *
00085 *     This function calls:  fabs.                                              *
00086 *                                                                              *
00087 *******************************************************************************/
00088 
00089 /*******************************************************************************
00090 *     First, an elegant implementation.                                        *
00091 ********************************************************************************
00092 *
00093 *double rint ( double x )
00094 *      {
00095 *      double y;
00096 *      
00097 *      y = twoTo52.fval;
00098 *      
00099 *      if ( fabs ( x ) >= y )                          // huge case is exact 
00100 *            return x;
00101 *      if ( x < 0 ) y = -y;                            // negative case 
00102 *      y = ( x + y ) - y;                              // force rounding 
00103 *      if ( y == 0.0 )                                 // zero results mirror sign of x 
00104 *            y = copysign ( y, x );
00105 *      return ( y );      
00106 *      }
00107 ********************************************************************************
00108 *     Now a bit twidling version that is about %30 faster.                     *
00109 *******************************************************************************/
00110 
00111 #if defined(__ppc__)
00112 double rint ( double x )
00113       {
00114       DblInHex argument;
00115       register double y;
00116       unsigned long int xHead;
00117       register long int target;
00118       
00119       argument.dbl = x;
00120       xHead = argument.words.hi & 0x7fffffffUL;          // xHead <- high half of |x|
00121       target = ( argument.words.hi < signMask );         // flags positive sign
00122       
00123       if ( xHead < 0x43300000ul ) 
00124 /*******************************************************************************
00125 *     Is |x| < 2.0^52?                                                         *
00126 *******************************************************************************/
00127             {
00128             if ( xHead < 0x3ff00000ul ) 
00129 /*******************************************************************************
00130 *     Is |x| < 1.0?                                                            *
00131 *******************************************************************************/
00132                   {
00133                   if ( target )
00134                         y = ( x + twoTo52 ) - twoTo52;  // round at binary point
00135                   else
00136                         y = ( x - twoTo52 ) + twoTo52;  // round at binary point
00137                   if ( y == 0.0 ) 
00138                         {                               // fix sign of zero result
00139                         if ( target )
00140                               return ( 0.0 );
00141                         else
00142                               return ( -0.0 );
00143                         }
00144                   return y;
00145                   }
00146             
00147 /*******************************************************************************
00148 *     Is 1.0 < |x| < 2.0^52?                                                   *
00149 *******************************************************************************/
00150 
00151             if ( target )
00152                   return ( ( x + twoTo52 ) - twoTo52 ); //   round at binary pt.
00153             else
00154                   return ( ( x - twoTo52 ) + twoTo52 );
00155             }
00156       
00157 /*******************************************************************************
00158 *     |x| >= 2.0^52 or x is a NaN.                                             *
00159 *******************************************************************************/
00160       return ( x );
00161       }
00162 #endif /* __ppc__ */
00163 
00164 /*******************************************************************************
00165 *                                                                              *
00166 *     The function nearbyint rounds its double argument to integral value      *
00167 *     according to the current rounding direction and returns the result in    *
00168 *     double format.  This function does not signal inexact.                   *
00169 *                                                                              *
00170 ********************************************************************************
00171 *                                                                              *
00172 *     This function calls fabs and copysign.                                 *
00173 *                                                                              *
00174 *******************************************************************************/
00175    
00176 double nearbyint ( double x )
00177       {
00178     double y;
00179 #if defined(__ppc__)
00180     double OldEnvironment;
00181 #endif /* __ppc__ */
00182       
00183     y = twoTo52;
00184     
00185     asm ("mffs %0" : "=f" (OldEnvironment));    /* get the environement */
00186 
00187       if ( fabs ( x ) >= y )                          /* huge case is exact */
00188             return x;
00189       if ( x < 0 ) y = -y;                                   /* negative case */
00190       y = ( x + y ) - y;                                    /* force rounding */
00191       if ( y == 0.0 )                        /* zero results mirror sign of x */
00192             y = copysign ( y, x );
00193 //  restore old flags
00194     asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
00195       return ( y );      
00196     }
00197       
00198 /*******************************************************************************
00199 *                                                                              *
00200 *     The function rinttol converts its double argument to integral value      *
00201 *     according to the current rounding direction and returns the result in    *
00202 *     long int format.  This conversion signals invalid if the argument is a   *
00203 *     NaN or the rounded intermediate result is out of range of the            *
00204 *     destination long int format, and it delivers an unspecified result in    *
00205 *     this case.  This function signals inexact if the rounded result is       *
00206 *     within range of the long int format but unequal to the operand.          *
00207 *                                                                              *
00208 *******************************************************************************/
00209 
00210 long int rinttol ( double x )
00211       {
00212       register double y;
00213       DblInHex argument, OldEnvironment;
00214       unsigned long int xHead;
00215       register long int target;
00216       
00217       argument.dbl = x;
00218       target = ( argument.words.hi < signMask );        // flag positive sign
00219       xHead = argument.words.hi & 0x7ffffffful;         // high 32 bits of x
00220       
00221       if ( target ) 
00222 /*******************************************************************************
00223 *    Sign of x is positive.                                                    *
00224 *******************************************************************************/
00225             {
00226             if ( xHead < 0x41dffffful ) 
00227                   {                                    // x is safely in long range
00228                   y = ( x + twoTo52 ) - twoTo52;       // round at binary point
00229                   argument.dbl = y + doubleToLong;     // force result into argument.words.lo
00230                   return ( ( long ) argument.words.lo );
00231                   }
00232             
00233         asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
00234             
00235             if ( xHead > 0x41dffffful ) 
00236                   {                                    // x is safely out of long range
00237                   OldEnvironment.words.lo |= SET_INVALID;
00238             asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00239                   return ( LONG_MAX );
00240                   }
00241             
00242 /*******************************************************************************
00243 *     x > 0.0 and may or may not be out of range of long.                      *
00244 *******************************************************************************/
00245 
00246             y = ( x + twoTo52 ) - twoTo52;             // do rounding
00247             if ( y > ( double ) LONG_MAX ) 
00248                   {                                    // out of range of long
00249                   OldEnvironment.words.lo |= SET_INVALID;
00250             asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00251                   return ( LONG_MAX );
00252                   }
00253             argument.dbl = y + doubleToLong;           // in range
00254             return ( ( long ) argument.words.lo );      // return result & flags
00255             }
00256       
00257 /*******************************************************************************
00258 *    Sign of x is negative.                                                    *
00259 *******************************************************************************/
00260       if ( xHead < 0x41e00000ul ) 
00261             {                                          // x is safely in long range
00262             y = ( x - twoTo52 ) + twoTo52;
00263             argument.dbl = y + doubleToLong;
00264             return ( ( long ) argument.words.lo );
00265             }
00266       
00267     asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
00268       
00269       if ( xHead > 0x41e00000ul ) 
00270             {                                          // x is safely out of long range
00271             OldEnvironment.words.lo |= SET_INVALID;
00272         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00273             return ( LONG_MIN );
00274             }
00275       
00276 /*******************************************************************************
00277 *    x < 0.0 and may or may not be out of range of long.                       *
00278 *******************************************************************************/
00279 
00280       y = ( x - twoTo52 ) + twoTo52;                   // do rounding
00281       if ( y < ( double ) LONG_MIN ) 
00282             {                                          // out of range of long
00283             OldEnvironment.words.lo |= SET_INVALID;
00284         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00285             return ( LONG_MIN );
00286             }
00287       argument.dbl = y + doubleToLong;                       // in range
00288       return ( ( long ) argument.words.lo );           // return result & flags
00289       }
00290 
00291 /*******************************************************************************
00292 *                                                                              *
00293 *     The function round rounds its double argument to integral value          *
00294 *     according to the "add half to the magnitude and truncate" rounding of    *
00295 *     Pascal's Round function and FORTRAN's ANINT function and returns the     *
00296 *     result in double format.  This function signals inexact if an ordered    *
00297 *     return value is not equal to the operand.                                *
00298 *                                                                              *
00299 *******************************************************************************/
00300    
00301 double round ( double x )
00302       {      
00303       DblInHex argument, OldEnvironment;
00304       register double y, z;
00305       register unsigned long int xHead;
00306       register long int target;
00307       
00308       argument.dbl = x;
00309       xHead = argument.words.hi & 0x7fffffffUL;      // xHead <- high half of |x|
00310       target = ( argument.words.hi < signMask );     // flag positive sign
00311       
00312       if ( xHead < 0x43300000ul ) 
00313 /*******************************************************************************
00314 *     Is |x| < 2.0^52?                                                        *
00315 *******************************************************************************/
00316             {
00317             if ( xHead < 0x3ff00000ul ) 
00318 /*******************************************************************************
00319 *     Is |x| < 1.0?                                                           *
00320 *******************************************************************************/
00321                   {
00322             asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
00323                   if ( xHead < 0x3fe00000ul ) 
00324 /*******************************************************************************
00325 *     Is |x| < 0.5?                                                           *
00326 *******************************************************************************/
00327                         {
00328                         if ( ( xHead | argument.words.lo ) != 0ul )
00329                               OldEnvironment.words.lo |= 0x02000000ul;
00330                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00331                         if ( target ) 
00332                               return ( 0.0 );
00333                         else
00334                               return ( -0.0 );
00335                         }
00336 /*******************************************************************************
00337 *     Is 0.5 ² |x| < 1.0?                                                      *
00338 *******************************************************************************/
00339                   OldEnvironment.words.lo |= 0x02000000ul;
00340             asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00341                   if ( target )
00342                         return ( 1.0 );
00343                   else
00344                         return ( -1.0 );
00345                   }
00346 /*******************************************************************************
00347 *     Is 1.0 < |x| < 2.0^52?                                                   *
00348 *******************************************************************************/
00349             if ( target ) 
00350                   {                                     // positive x
00351                   y = ( x + twoTo52 ) - twoTo52;        // round at binary point
00352                   if ( y == x )                         // exact case
00353                         return ( x );
00354                   z = x + 0.5;                          // inexact case
00355                   y = ( z + twoTo52 ) - twoTo52;        // round at binary point
00356                   if ( y > z )
00357                         return ( y - 1.0 );
00358                   else
00359                         return ( y );
00360                   }
00361             
00362 /*******************************************************************************
00363 *     Is x < 0?                                                                *
00364 *******************************************************************************/
00365             else 
00366                   {
00367                   y = ( x - twoTo52 ) + twoTo52;        // round at binary point
00368                   if ( y == x )
00369                         return ( x );
00370                   z = x - 0.5;
00371                   y = ( z - twoTo52 ) + twoTo52;        // round at binary point
00372                   if ( y < z )
00373                         return ( y + 1.0 );
00374                   else
00375                   return ( y );
00376                   }
00377             }
00378 /*******************************************************************************
00379 *      |x| >= 2.0^52 or x is a NaN.                                            *
00380 *******************************************************************************/
00381       return ( x );
00382       }
00383 
00384 /*******************************************************************************
00385 *                                                                              *
00386 *     The function roundtol converts its double argument to integral format    *
00387 *     according to the "add half to the magnitude and chop" rounding mode of   *
00388 *     Pascal's Round function and FORTRAN's NINT function.  This conversion    *
00389 *     signals invalid if the argument is a NaN or the rounded intermediate     *
00390 *     result is out of range of the destination long int format, and it        *
00391 *     delivers an unspecified result in this case.  This function signals      *
00392 *     inexact if the rounded result is within range of the long int format but *
00393 *     unequal to the operand.                                                  *
00394 *                                                                              *
00395 *******************************************************************************/
00396 
00397 long int roundtol ( double x )
00398     {   
00399     register double y, z;
00400     DblInHex argument, OldEnvironment;
00401     register unsigned long int xhi;
00402     register long int target;
00403 #if defined(__ppc__)
00404     const DblInHex kTZ = {{ 0x0, 0x1 }};
00405     const DblInHex kUP = {{ 0x0, 0x2 }};
00406 #endif /* __ppc__ */
00407     
00408     argument.dbl = x;
00409     xhi = argument.words.hi & 0x7ffffffful;             // high 32 bits of x
00410     target = ( argument.words.hi < signMask );          // flag positive sign
00411     
00412     if ( xhi > 0x41e00000ul ) 
00413 /*******************************************************************************
00414 *     Is x is out of long range or NaN?                                        *
00415 *******************************************************************************/
00416         {
00417         asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
00418         OldEnvironment.words.lo |= SET_INVALID;
00419         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00420         if ( target )                           // pin result
00421             return ( LONG_MAX );
00422         else
00423             return ( LONG_MIN );
00424         }
00425     
00426     if ( target ) 
00427 /*******************************************************************************
00428 *     Is sign of x is "+"?                                                     *
00429 *******************************************************************************/
00430         {
00431         if ( x < 2147483647.5 ) 
00432 /*******************************************************************************
00433 *     x is in the range of a long.                                             *
00434 *******************************************************************************/
00435             {
00436             y = ( x + doubleToLong ) - doubleToLong;    // round at binary point
00437             if ( y != x )   
00438                 {                               // inexact case
00439                 asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // save environment
00440                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding
00441                 z = x + 0.5;                    // truncate x + 0.5
00442                 argument.dbl = z + doubleToLong;
00443                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00444                 return ( ( long ) argument.words.lo );
00445                 }
00446             
00447             argument.dbl = y + doubleToLong;        // force result into argument.words.lo
00448             return ( ( long ) argument.words.lo );  // return long result
00449             }
00450 /*******************************************************************************
00451 *     Rounded positive x is out of the range of a long.                        *
00452 *******************************************************************************/
00453         asm ("mffs %0" : "=f" (OldEnvironment.dbl));
00454         OldEnvironment.words.lo |= SET_INVALID;
00455         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00456         return ( LONG_MAX );                        // return pinned result
00457         }
00458 /*******************************************************************************
00459 *     x < 0.0 and may or may not be out of the range of a long.                *
00460 *******************************************************************************/
00461     if ( x > -2147483648.5 ) 
00462 /*******************************************************************************
00463 *     x is in the range of a long.                                             *
00464 *******************************************************************************/
00465         {
00466         y = ( x + doubleToLong ) - doubleToLong;        // round at binary point
00467         if ( y != x ) 
00468             {                                   // inexact case
00469             asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // save environment
00470             asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up
00471             z = x - 0.5;                        // truncate x - 0.5
00472             argument.dbl = z + doubleToLong;
00473             asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00474             return ( ( long ) argument.words.lo );
00475             }
00476         
00477         argument.dbl = y + doubleToLong;
00478         return ( ( long ) argument.words.lo );      //  return long result
00479         }
00480 /*******************************************************************************
00481 *     Rounded negative x is out of the range of a long.                        *
00482 *******************************************************************************/
00483     asm ("mffs %0" : "=f" (OldEnvironment.dbl));
00484     OldEnvironment.words.lo |= SET_INVALID;
00485     asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00486     return ( LONG_MIN );                            // return pinned result
00487     }
00488 
00489 /*******************************************************************************
00490 *                                                                              *
00491 *     The function trunc truncates its double argument to integral value       *
00492 *     and returns the result in double format.  This function signals          *
00493 *     inexact if an ordered return value is not equal to the operand.          *
00494 *                                                                              *
00495 *******************************************************************************/
00496    
00497 double trunc ( double x )
00498       { 
00499     DblInHex argument,OldEnvironment;
00500     register double y;
00501     register unsigned long int xhi;
00502     register long int target;
00503     
00504     argument.dbl = x;
00505     xhi = argument.words.hi & 0x7fffffffUL;         // xhi <- high half of |x|
00506     target = ( argument.words.hi < signMask );          // flag positive sign
00507     
00508     if ( xhi < 0x43300000ul ) 
00509 /*******************************************************************************
00510 *     Is |x| < 2.0^53?                                                         *
00511 *******************************************************************************/
00512         {
00513         if ( xhi < 0x3ff00000ul ) 
00514 /*******************************************************************************
00515 *     Is |x| < 1.0?                                                            *
00516 *******************************************************************************/
00517             {
00518             if ( ( xhi | argument.words.lo ) != 0ul ) 
00519                 {                               // raise deserved INEXACT
00520                 asm ("mffs %0" : "=f" (OldEnvironment.dbl));
00521                 OldEnvironment.words.lo |= 0x02000000ul;
00522                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
00523                 }
00524             if ( target )                       // return properly signed zero
00525                 return ( 0.0 );
00526             else
00527                 return ( -0.0 );
00528             }
00529 /*******************************************************************************
00530 *     Is 1.0 < |x| < 2.0^52?                                                   *
00531 *******************************************************************************/
00532         if ( target ) 
00533             {
00534             y = ( x + twoTo52 ) - twoTo52;          // round at binary point
00535             if ( y > x )
00536                 return ( y - 1.0 );
00537             else
00538                 return ( y );
00539             }
00540         
00541         else 
00542             {
00543             y = ( x - twoTo52 ) + twoTo52;          // round at binary point.
00544             if ( y < x )
00545                 return ( y + 1.0 );
00546             else
00547                 return ( y );
00548             }
00549         }
00550 /*******************************************************************************
00551 *      Is |x| >= 2.0^52 or x is a NaN.                                         *
00552 *******************************************************************************/
00553     return ( x );
00554     }
00555 
00556 /*******************************************************************************
00557 *     The modf family of functions separate a floating-point number into its   *
00558 *     fractional and integral parts, returning the fractional part and writing *
00559 *     the integral part in floating-point format to the object pointed to by a *
00560 *     pointer argument.  If the input argument is integral or infinite in      *
00561 *     value, the return value is a zero with the sign of the input argument.   *
00562 *     The modf family of functions raises no floating-point exceptions. older  *
00563 *     implemenation set the INVALID flag due to signaling NaN input.           *
00564 *                                                                              *
00565 *******************************************************************************/
00566 
00567 /*******************************************************************************
00568 *     modf is the double implementation.                                       *                             
00569 *******************************************************************************/
00570 
00571 #if defined(__ppc__)
00572 double modf ( double x, double *iptr )
00573       {
00574       register double OldEnvironment, xtrunc;
00575       register unsigned long int xHead, signBit;
00576       DblInHex argument;
00577       
00578       argument.dbl = x;
00579       xHead = argument.words.hi & 0x7ffffffful;            // |x| high bit pattern
00580       signBit = ( argument.words.hi & 0x80000000ul );      // isolate sign bit
00581       if (xHead == 0x7ff81fe0)
00582             signBit = signBit | 0;
00583       
00584       if ( xHead < 0x43300000ul ) 
00585 /*******************************************************************************
00586 *     Is |x| < 2.0^53?                                                         *
00587 *******************************************************************************/
00588             {
00589             if ( xHead < 0x3ff00000ul )      
00590 /*******************************************************************************
00591 *     Is |x| < 1.0?                                                            *
00592 *******************************************************************************/
00593                   {
00594                   argument.words.hi = signBit;             // truncate to zero
00595                   argument.words.lo = 0ul;
00596                   *iptr = argument.dbl;
00597                   return ( x );
00598                   }
00599 /*******************************************************************************
00600 *     Is 1.0 < |x| < 2.0^52?                                                   *
00601 *******************************************************************************/
00602             asm ("mffs %0" : "=f" (OldEnvironment));    // save environment
00603             // round toward zero
00604             asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl ));
00605             if ( signBit == 0ul )                         // truncate to integer
00606                   xtrunc = ( x + twoTo52 ) - twoTo52;
00607             else
00608                   xtrunc = ( x - twoTo52 ) + twoTo52;
00609         // restore caller's env
00610         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
00611             *iptr = xtrunc;                               // store integral part
00612             if ( x != xtrunc )                            // nonzero fraction
00613                   return ( x - xtrunc );
00614             else 
00615                   {                                       // zero with x's sign
00616                   argument.words.hi = signBit;
00617                   argument.words.lo = 0ul;
00618                   return ( argument.dbl );
00619                   }
00620             }
00621       
00622       *iptr = x;                                          // x is integral or NaN
00623       if ( x != x )                                       // NaN is returned
00624             return x;
00625       else 
00626             {                                             // zero with x's sign
00627             argument.words.hi = signBit;
00628             argument.words.lo = 0ul;
00629             return ( argument.dbl );
00630             }
00631       }
00632 #endif /* __ppc__ */

Generated on Tue Feb 2 17:46:05 2010 for RTAI API by  doxygen 1.4.7