libm/powerpc/rndint.c

   1 /*******************************************************************************
   2 **      File:   rndint.c
   3 **
   4 **      Contains: C source code for implementations of floating-point
   5 **                functions which round to integral value or format, as
   6 **                defined in header <fp.h>.  In particular, this file
   7 **                contains implementations of functions rint, nearbyint,
   8 **                rinttol, round, roundtol, trunc, modf and modfl.  This file
   9 **                targets PowerPC or Power platforms.
  10 **
  11 **      Written by: A. Sazegari, Apple AltiVec Group
  12 **          Created originally by Jon Okada, Apple Numerics Group
  13 **
  14 **      Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
  15 **
  16 **      Change History (most recent first):
  17 **
  18 **      13 Jul 01  ram  replaced --setflm calls with inline assembly
  19 **      03 Mar 01  ali  first port to os x using gcc, added the crucial __setflm
  20 **                      definition.
  21 **                              1. removed double_t, put in double for now.
  22 **                              2. removed iclass from nearbyint.
  23 **                              3. removed wrong comments intrunc.
  24 **                              4.
  25 **      13 May 97  ali  made performance improvements in rint, rinttol, roundtol
  26 **                      and trunc by folding some of the taligent ideas into this
  27 **                      implementation.  nearbyint is faster than the one in taligent,
  28 **                      rint is more elegant, but slower by %30 than the taligent one.
  29 **      09 Apr 97  ali  deleted modfl and deferred to AuxiliaryDD.c
  30 **      15 Sep 94  ali  Major overhaul and performance improvements of all functions.
  31 **      20 Jul 94  PAF  New faster version
  32 **      16 Jul 93  ali  Added the modfl function.
  33 **      18 Feb 93  ali  Changed the return value of fenv functions
  34 **                      feclearexcept and feraiseexcept to their new
  35 **                      NCEG X3J11.1/93-001 definitions.
  36 **      16 Dec 92  JPO  Removed __itrunc implementation to a
  37 **                      separate file.
  38 **      15 Dec 92  JPO  Added __itrunc implementation and modified
  39 **                      rinttol to include conversion from double
  40 **                      to long int format.  Modified roundtol to
  41 **                      call __itrunc.
  42 **      10 Dec 92  JPO  Added modf (double) implementation.
  43 **      04 Dec 92  JPO  First created.
  44 **
  45 *******************************************************************************/
  46
  47 #include <limits.h>
  48 #include <math.h>
  49
  50 #define      SET_INVALID      0x01000000UL
  51
  52 typedef union
  53       {
  54       struct {
  55 #if defined(__BIG_ENDIAN__)
  56         unsigned long int hi;
  57         unsigned long int lo;
  58 #else
  59         unsigned long int lo;
  60         unsigned long int hi;
  61 #endif
  62       } words;
  63       double dbl;
  64       } DblInHex;
  65
  66 static const unsigned long int signMask = 0x80000000ul;
  67 static const double twoTo52      = 4503599627370496.0;
  68 static const double doubleToLong = 4503603922337792.0;              // 2^52
  69 static const DblInHex Huge       = {{ 0x7FF00000, 0x00000000 }};
  70 static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
  71
  72 /*******************************************************************************
  73 *                                                                              *
  74 *     The function rint rounds its double argument to integral value           *
  75 *     according to the current rounding direction and returns the result in    *
  76 *     double format.  This function signals inexact if an ordered return       *
  77 *     value is not equal to the operand.                                       *
  78 *                                                                              *
  79 ********************************************************************************
  80 *                                                                              *
  81 *     This function calls:  fabs.                                                  *
  82 *                                                                              *
  83 *******************************************************************************/
  84
  85 /*******************************************************************************
  86 *     First, an elegant implementation.                                        *
  87 ********************************************************************************
  88 *
  89 *double rint ( double x )
  90 *      {
  91 *      double y;
  92 *
  93 *      y = twoTo52.fval;
  94 *
  95 *      if ( fabs ( x ) >= y )                          // huge case is exact
  96 *            return x;
  97 *      if ( x < 0 ) y = -y;                            // negative case
  98 *      y = ( x + y ) - y;                              // force rounding
  99 *      if ( y == 0.0 )                                 // zero results mirror sign of x
 100 *            y = copysign ( y, x );
 101 *      return ( y );
 102 *      }
 103 ********************************************************************************
 104 *     Now a bit twidling version that is about %30 faster.                     *
 105 *******************************************************************************/
 106
 107 double rint ( double x )
 108       {
 109       DblInHex argument;
 110       register double y;
 111       unsigned long int xHead;
 112       register long int target;
 113
 114       argument.dbl = x;
 115       xHead = argument.words.hi & 0x7fffffffUL;          // xHead <- high half of |x|
 116       target = ( argument.words.hi < signMask );         // flags positive sign
 117
 118       if ( xHead < 0x43300000ul )
 119 /*******************************************************************************
 120 *     Is |x| < 2.0^52?                                                         *
 121 *******************************************************************************/
 122             {
 123             if ( xHead < 0x3ff00000ul )
 124 /*******************************************************************************
 125 *     Is |x| < 1.0?                                                            *
 126 *******************************************************************************/
 127                   {
 128                   if ( target )
 129                         y = ( x + twoTo52 ) - twoTo52;  // round at binary point
 130                   else
 131                         y = ( x - twoTo52 ) + twoTo52;  // round at binary point
 132                   if ( y == 0.0 )
 133                         {                               // fix sign of zero result
 134                         if ( target )
 135                               return ( 0.0 );
 136                         else
 137                               return ( -0.0 );
 138                         }
 139                   return y;
 140                   }
 141
 142 /*******************************************************************************
 143 *     Is 1.0 < |x| < 2.0^52?                                                   *
 144 *******************************************************************************/
 145
 146             if ( target )
 147                   return ( ( x + twoTo52 ) - twoTo52 ); //   round at binary pt.
 148             else
 149                   return ( ( x - twoTo52 ) + twoTo52 );
 150             }
 151
 152 /*******************************************************************************
 153 *     |x| >= 2.0^52 or x is a NaN.                                             *
 154 *******************************************************************************/
 155       return ( x );
 156       }
 157
 158 /*******************************************************************************
 159 *                                                                              *
 160 *     The function nearbyint rounds its double argument to integral value      *
 161 *     according to the current rounding direction and returns the result in    *
 162 *     double format.  This function does not signal inexact.                   *
 163 *                                                                              *
 164 ********************************************************************************
 165 *                                                                              *
 166 *     This function calls fabs and copysign.                                     *
 167 *                                                                              *
 168 *******************************************************************************/
 169
 170 double nearbyint ( double x )
 171       {
 172         double y;
 173         double OldEnvironment;
 174
 175         y = twoTo52;
 176
 177         asm ("mffs %0" : "=f" (OldEnvironment));        /* get the environement */
 178
 179       if ( fabs ( x ) >= y )                          /* huge case is exact */
 180             return x;
 181       if ( x < 0 ) y = -y;                                   /* negative case */
 182       y = ( x + y ) - y;                                    /* force rounding */
 183       if ( y == 0.0 )                        /* zero results mirror sign of x */
 184             y = copysign ( y, x );
 185 //      restore old flags
 186         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
 187       return ( y );
 188         }
 189
 190 /*******************************************************************************
 191 *                                                                              *
 192 *     The function rinttol converts its double argument to integral value      *
 193 *     according to the current rounding direction and returns the result in    *
 194 *     long int format.  This conversion signals invalid if the argument is a   *
 195 *     NaN or the rounded intermediate result is out of range of the            *
 196 *     destination long int format, and it delivers an unspecified result in    *
 197 *     this case.  This function signals inexact if the rounded result is       *
 198 *     within range of the long int format but unequal to the operand.          *
 199 *                                                                              *
 200 *******************************************************************************/
 201
 202 long int rinttol ( double x )
 203       {
 204       register double y;
 205       DblInHex argument, OldEnvironment;
 206       unsigned long int xHead;
 207       register long int target;
 208
 209       argument.dbl = x;
 210       target = ( argument.words.hi < signMask );        // flag positive sign
 211       xHead = argument.words.hi & 0x7ffffffful;         // high 32 bits of x
 212
 213       if ( target )
 214 /*******************************************************************************
 215 *    Sign of x is positive.                                                    *
 216 *******************************************************************************/
 217             {
 218             if ( xHead < 0x41dffffful )
 219                   {                                    // x is safely in long range
 220                   y = ( x + twoTo52 ) - twoTo52;       // round at binary point
 221                   argument.dbl = y + doubleToLong;     // force result into argument.words.lo
 222                   return ( ( long ) argument.words.lo );
 223                   }
 224
 225                 asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
 226
 227             if ( xHead > 0x41dffffful )
 228                   {                                    // x is safely out of long range
 229                   OldEnvironment.words.lo |= SET_INVALID;
 230                         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 231                   return ( LONG_MAX );
 232                   }
 233
 234 /*******************************************************************************
 235 *     x > 0.0 and may or may not be out of range of long.                      *
 236 *******************************************************************************/
 237
 238             y = ( x + twoTo52 ) - twoTo52;             // do rounding
 239             if ( y > ( double ) LONG_MAX )
 240                   {                                    // out of range of long
 241                   OldEnvironment.words.lo |= SET_INVALID;
 242                         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 243                   return ( LONG_MAX );
 244                   }
 245             argument.dbl = y + doubleToLong;           // in range
 246             return ( ( long ) argument.words.lo );      // return result & flags
 247             }
 248
 249 /*******************************************************************************
 250 *    Sign of x is negative.                                                    *
 251 *******************************************************************************/
 252       if ( xHead < 0x41e00000ul )
 253             {                                          // x is safely in long range
 254             y = ( x - twoTo52 ) + twoTo52;
 255             argument.dbl = y + doubleToLong;
 256             return ( ( long ) argument.words.lo );
 257             }
 258
 259         asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
 260
 261       if ( xHead > 0x41e00000ul )
 262             {                                          // x is safely out of long range
 263             OldEnvironment.words.lo |= SET_INVALID;
 264                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 265             return ( LONG_MIN );
 266             }
 267
 268 /*******************************************************************************
 269 *    x < 0.0 and may or may not be out of range of long.                       *
 270 *******************************************************************************/
 271
 272       y = ( x - twoTo52 ) + twoTo52;                   // do rounding
 273       if ( y < ( double ) LONG_MIN )
 274             {                                          // out of range of long
 275             OldEnvironment.words.lo |= SET_INVALID;
 276                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 277             return ( LONG_MIN );
 278             }
 279       argument.dbl = y + doubleToLong;                       // in range
 280       return ( ( long ) argument.words.lo );           // return result & flags
 281       }
 282
 283 /*******************************************************************************
 284 *                                                                              *
 285 *     The function round rounds its double argument to integral value          *
 286 *     according to the "add half to the magnitude and truncate" rounding of    *
 287 *     Pascal's Round function and FORTRAN's ANINT function and returns the     *
 288 *     result in double format.  This function signals inexact if an ordered    *
 289 *     return value is not equal to the operand.                                *
 290 *                                                                              *
 291 *******************************************************************************/
 292
 293 double round ( double x )
 294       {
 295       DblInHex argument, OldEnvironment;
 296       register double y, z;
 297       register unsigned long int xHead;
 298       register long int target;
 299
 300       argument.dbl = x;
 301       xHead = argument.words.hi & 0x7fffffffUL;      // xHead <- high half of |x|
 302       target = ( argument.words.hi < signMask );     // flag positive sign
 303
 304       if ( xHead < 0x43300000ul )
 305 /*******************************************************************************
 306 *     Is |x| < 2.0^52?                                                        *
 307 *******************************************************************************/
 308             {
 309             if ( xHead < 0x3ff00000ul )
 310 /*******************************************************************************
 311 *     Is |x| < 1.0?                                                           *
 312 *******************************************************************************/
 313                   {
 314                         asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
 315                   if ( xHead < 0x3fe00000ul )
 316 /*******************************************************************************
 317 *     Is |x| < 0.5?                                                           *
 318 *******************************************************************************/
 319                         {
 320                         if ( ( xHead | argument.words.lo ) != 0ul )
 321                               OldEnvironment.words.lo |= 0x02000000ul;
 322                                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 323                         if ( target )
 324                               return ( 0.0 );
 325                         else
 326                               return ( -0.0 );
 327                         }
 328 /*******************************************************************************
 329 *     Is 0.5 ² |x| < 1.0?                                                      *
 330 *******************************************************************************/
 331                   OldEnvironment.words.lo |= 0x02000000ul;
 332                         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 333                   if ( target )
 334                         return ( 1.0 );
 335                   else
 336                         return ( -1.0 );
 337                   }
 338 /*******************************************************************************
 339 *     Is 1.0 < |x| < 2.0^52?                                                   *
 340 *******************************************************************************/
 341             if ( target )
 342                   {                                     // positive x
 343                   y = ( x + twoTo52 ) - twoTo52;        // round at binary point
 344                   if ( y == x )                         // exact case
 345                         return ( x );
 346                   z = x + 0.5;                          // inexact case
 347                   y = ( z + twoTo52 ) - twoTo52;        // round at binary point
 348                   if ( y > z )
 349                         return ( y - 1.0 );
 350                   else
 351                         return ( y );
 352                   }
 353
 354 /*******************************************************************************
 355 *     Is x < 0?                                                                *
 356 *******************************************************************************/
 357             else
 358                   {
 359                   y = ( x - twoTo52 ) + twoTo52;        // round at binary point
 360                   if ( y == x )
 361                         return ( x );
 362                   z = x - 0.5;
 363                   y = ( z - twoTo52 ) + twoTo52;        // round at binary point
 364                   if ( y < z )
 365                         return ( y + 1.0 );
 366                   else
 367                   return ( y );
 368                   }
 369             }
 370 /*******************************************************************************
 371 *      |x| >= 2.0^52 or x is a NaN.                                            *
 372 *******************************************************************************/
 373       return ( x );
 374       }
 375
 376 /*******************************************************************************
 377 *                                                                              *
 378 *     The function roundtol converts its double argument to integral format    *
 379 *     according to the "add half to the magnitude and chop" rounding mode of   *
 380 *     Pascal's Round function and FORTRAN's NINT function.  This conversion    *
 381 *     signals invalid if the argument is a NaN or the rounded intermediate     *
 382 *     result is out of range of the destination long int format, and it        *
 383 *     delivers an unspecified result in this case.  This function signals      *
 384 *     inexact if the rounded result is within range of the long int format but *
 385 *     unequal to the operand.                                                  *
 386 *                                                                              *
 387 *******************************************************************************/
 388
 389 long int roundtol ( double x )
 390         {
 391         register double y, z;
 392         DblInHex argument, OldEnvironment;
 393         register unsigned long int xhi;
 394         register long int target;
 395         const DblInHex kTZ = {{ 0x0, 0x1 }};
 396         const DblInHex kUP = {{ 0x0, 0x2 }};
 397
 398         argument.dbl = x;
 399         xhi = argument.words.hi & 0x7ffffffful;                 // high 32 bits of x
 400         target = ( argument.words.hi < signMask );              // flag positive sign
 401
 402         if ( xhi > 0x41e00000ul )
 403 /*******************************************************************************
 404 *     Is x is out of long range or NaN?                                        *
 405 *******************************************************************************/
 406                 {
 407                 asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // get environment
 408                 OldEnvironment.words.lo |= SET_INVALID;
 409                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 410                 if ( target )                                   // pin result
 411                         return ( LONG_MAX );
 412                 else
 413                         return ( LONG_MIN );
 414                 }
 415
 416         if ( target )
 417 /*******************************************************************************
 418 *     Is sign of x is "+"?                                                     *
 419 *******************************************************************************/
 420                 {
 421                 if ( x < 2147483647.5 )
 422 /*******************************************************************************
 423 *     x is in the range of a long.                                             *
 424 *******************************************************************************/
 425                         {
 426                         y = ( x + doubleToLong ) - doubleToLong;        // round at binary point
 427                         if ( y != x )
 428                                 {                                       // inexact case
 429                                 asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // save environment
 430                                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding
 431                                 z = x + 0.5;                            // truncate x + 0.5
 432                                 argument.dbl = z + doubleToLong;
 433                                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 434                                 return ( ( long ) argument.words.lo );
 435                                 }
 436
 437                         argument.dbl = y + doubleToLong;                // force result into argument.words.lo
 438                         return ( ( long ) argument.words.lo );  // return long result
 439                         }
 440 /*******************************************************************************
 441 *     Rounded positive x is out of the range of a long.                        *
 442 *******************************************************************************/
 443                 asm ("mffs %0" : "=f" (OldEnvironment.dbl));
 444                 OldEnvironment.words.lo |= SET_INVALID;
 445                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 446                 return ( LONG_MAX );                            // return pinned result
 447                 }
 448 /*******************************************************************************
 449 *     x < 0.0 and may or may not be out of the range of a long.                *
 450 *******************************************************************************/
 451         if ( x > -2147483648.5 )
 452 /*******************************************************************************
 453 *     x is in the range of a long.                                             *
 454 *******************************************************************************/
 455                 {
 456                 y = ( x + doubleToLong ) - doubleToLong;                // round at binary point
 457                 if ( y != x )
 458                         {                                               // inexact case
 459                         asm ("mffs %0" : "=f" (OldEnvironment.dbl));    // save environment
 460                         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up
 461                         z = x - 0.5;                            // truncate x - 0.5
 462                         argument.dbl = z + doubleToLong;
 463                         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 464                         return ( ( long ) argument.words.lo );
 465                         }
 466
 467                 argument.dbl = y + doubleToLong;
 468                 return ( ( long ) argument.words.lo );          //  return long result
 469                 }
 470 /*******************************************************************************
 471 *     Rounded negative x is out of the range of a long.                        *
 472 *******************************************************************************/
 473         asm ("mffs %0" : "=f" (OldEnvironment.dbl));
 474         OldEnvironment.words.lo |= SET_INVALID;
 475         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 476         return ( LONG_MIN );                                    // return pinned result
 477         }
 478
 479 /*******************************************************************************
 480 *                                                                              *
 481 *     The function trunc truncates its double argument to integral value       *
 482 *     and returns the result in double format.  This function signals          *
 483 *     inexact if an ordered return value is not equal to the operand.          *
 484 *                                                                              *
 485 *******************************************************************************/
 486
 487 double trunc ( double x )
 488       {
 489         DblInHex argument,OldEnvironment;
 490         register double y;
 491         register unsigned long int xhi;
 492         register long int target;
 493
 494         argument.dbl = x;
 495         xhi = argument.words.hi & 0x7fffffffUL;         // xhi <- high half of |x|
 496         target = ( argument.words.hi < signMask );              // flag positive sign
 497
 498         if ( xhi < 0x43300000ul )
 499 /*******************************************************************************
 500 *     Is |x| < 2.0^53?                                                         *
 501 *******************************************************************************/
 502                 {
 503                 if ( xhi < 0x3ff00000ul )
 504 /*******************************************************************************
 505 *     Is |x| < 1.0?                                                            *
 506 *******************************************************************************/
 507                         {
 508                         if ( ( xhi | argument.words.lo ) != 0ul )
 509                                 {                               // raise deserved INEXACT
 510                                 asm ("mffs %0" : "=f" (OldEnvironment.dbl));
 511                                 OldEnvironment.words.lo |= 0x02000000ul;
 512                                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 513                                 }
 514                         if ( target )                           // return properly signed zero
 515                                 return ( 0.0 );
 516                         else
 517                                 return ( -0.0 );
 518                         }
 519 /*******************************************************************************
 520 *     Is 1.0 < |x| < 2.0^52?                                                   *
 521 *******************************************************************************/
 522                 if ( target )
 523                         {
 524                         y = ( x + twoTo52 ) - twoTo52;          // round at binary point
 525                         if ( y > x )
 526                                 return ( y - 1.0 );
 527                         else
 528                                 return ( y );
 529                         }
 530
 531                 else
 532                         {
 533                         y = ( x - twoTo52 ) + twoTo52;          // round at binary point.
 534                         if ( y < x )
 535                                 return ( y + 1.0 );
 536                         else
 537                                 return ( y );
 538                         }
 539                 }
 540 /*******************************************************************************
 541 *      Is |x| >= 2.0^52 or x is a NaN.                                         *
 542 *******************************************************************************/
 543         return ( x );
 544         }
 545
 546 /*******************************************************************************
 547 *     The modf family of functions separate a floating-point number into its   *
 548 *     fractional and integral parts, returning the fractional part and writing *
 549 *     the integral part in floating-point format to the object pointed to by a *
 550 *     pointer argument.  If the input argument is integral or infinite in      *
 551 *     value, the return value is a zero with the sign of the input argument.   *
 552 *     The modf family of functions raises no floating-point exceptions. older  *
 553 *     implemenation set the INVALID flag due to signaling NaN input.           *
 554 *                                                                              *
 555 *******************************************************************************/
 556
 557 /*******************************************************************************
 558 *     modf is the double implementation.                                       *
 559 *******************************************************************************/
 560
 561 double modf ( double x, double *iptr )
 562       {
 563       register double OldEnvironment, xtrunc;
 564       register unsigned long int xHead, signBit;
 565       DblInHex argument;
 566
 567       argument.dbl = x;
 568       xHead = argument.words.hi & 0x7ffffffful;            // |x| high bit pattern
 569       signBit = ( argument.words.hi & 0x80000000ul );      // isolate sign bit
 570           if (xHead == 0x7ff81fe0)
 571                         signBit = signBit | 0;
 572
 573       if ( xHead < 0x43300000ul )
 574 /*******************************************************************************
 575 *     Is |x| < 2.0^53?                                                         *
 576 *******************************************************************************/
 577             {
 578             if ( xHead < 0x3ff00000ul )
 579 /*******************************************************************************
 580 *     Is |x| < 1.0?                                                            *
 581 *******************************************************************************/
 582                   {
 583                   argument.words.hi = signBit;             // truncate to zero
 584                   argument.words.lo = 0ul;
 585                   *iptr = argument.dbl;
 586                   return ( x );
 587                   }
 588 /*******************************************************************************
 589 *     Is 1.0 < |x| < 2.0^52?                                                   *
 590 *******************************************************************************/
 591                         asm ("mffs %0" : "=f" (OldEnvironment));        // save environment
 592                         // round toward zero
 593                         asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl ));
 594             if ( signBit == 0ul )                         // truncate to integer
 595                   xtrunc = ( x + twoTo52 ) - twoTo52;
 596             else
 597                   xtrunc = ( x - twoTo52 ) + twoTo52;
 598                 // restore caller's env
 599                 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
 600             *iptr = xtrunc;                               // store integral part
 601             if ( x != xtrunc )                            // nonzero fraction
 602                   return ( x - xtrunc );
 603             else
 604                   {                                       // zero with x's sign
 605                   argument.words.hi = signBit;
 606                   argument.words.lo = 0ul;
 607                   return ( argument.dbl );
 608                   }
 609             }
 610
 611       *iptr = x;                                          // x is integral or NaN
 612       if ( x != x )                                       // NaN is returned
 613             return x;
 614       else
 615             {                                             // zero with x's sign
 616             argument.words.hi = signBit;
 617             argument.words.lo = 0ul;
 618             return ( argument.dbl );
 619             }
 620       }