git.sesse.net Git - pistorm/blob - softfloat/softfloat.c

   1
   2 #define SOFTFLOAT_68K
   3
   4 #include <stdint.h>
   5 #include <stdlib.h>
   6 #include "softfloat/softfloat.h"
   7
   8
   9 /*
  10  * QEMU float support
  11  *
  12  * The code in this source file is derived from release 2a of the SoftFloat
  13  * IEC/IEEE Floating-point Arithmetic Package. Those parts of the code (and
  14  * some later contributions) are provided under that license, as detailed below.
  15  * It has subsequently been modified by contributors to the QEMU Project,
  16  * so some portions are provided under:
  17  *  the SoftFloat-2a license
  18  *  the BSD license
  19  *  GPL-v2-or-later
  20  *
  21  * Any future contributions to this file after December 1st 2014 will be
  22  * taken to be licensed under the Softfloat-2a license unless specifically
  23  * indicated otherwise.
  24  */
  25
  26 /*
  27 ===============================================================================
  28 This C source file is part of the SoftFloat IEC/IEEE Floating-point
  29 Arithmetic Package, Release 2a.
  30
  31 Written by John R. Hauser.  This work was made possible in part by the
  32 International Computer Science Institute, located at Suite 600, 1947 Center
  33 Street, Berkeley, California 94704.  Funding was partially provided by the
  34 National Science Foundation under grant MIP-9311980.  The original version
  35 of this code was written as part of a project to build a fixed-point vector
  36 processor in collaboration with the University of California at Berkeley,
  37 overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
  38 is available through the Web page `http://HTTP.CS.Berkeley.EDU/~jhauser/
  39 arithmetic/SoftFloat.html'.
  40
  41 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
  42 has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
  43 TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
  44 PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
  45 AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
  46
  47 Derivative works are acceptable, even for commercial purposes, so long as
  48 (1) they include prominent notice that the work is derivative, and (2) they
  49 include prominent notice akin to these four paragraphs for those parts of
  50 this code that are retained.
  51
  52 ===============================================================================
  53 */
  54
  55 /* BSD licensing:
  56  * Copyright (c) 2006, Fabrice Bellard
  57  * All rights reserved.
  58  *
  59  * Redistribution and use in source and binary forms, with or without
  60  * modification, are permitted provided that the following conditions are met:
  61  *
  62  * 1. Redistributions of source code must retain the above copyright notice,
  63  * this list of conditions and the following disclaimer.
  64  *
  65  * 2. Redistributions in binary form must reproduce the above copyright notice,
  66  * this list of conditions and the following disclaimer in the documentation
  67  * and/or other materials provided with the distribution.
  68  *
  69  * 3. Neither the name of the copyright holder nor the names of its contributors
  70  * may be used to endorse or promote products derived from this software without
  71  * specific prior written permission.
  72  *
  73  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  74  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  75  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  76  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  77  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  78  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  79  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  80  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  81  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  82  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
  83  * THE POSSIBILITY OF SUCH DAMAGE.
  84  */
  85
  86 /* Portions of this work are licensed under the terms of the GNU GPL,
  87  * version 2 or later. See the COPYING file in the top-level directory.
  88  */
  89
  90 /* We only need stdlib for abort() */
  91
  92 /*----------------------------------------------------------------------------
  93 | Primitive arithmetic functions, including multi-word arithmetic, and
  94 | division and square root approximations.  (Can be specialized to target if
  95 | desired.)
  96 *----------------------------------------------------------------------------*/
  97 #include "softfloat-macros.h"
  98
  99 /*----------------------------------------------------------------------------
 100  | Variables for storing sign, exponent and significand of internal extended
 101  | double-precision floating-point value for external use.
 102  *----------------------------------------------------------------------------*/
 103 flag floatx80_internal_sign = 0;
 104 int32_t floatx80_internal_exp = 0;
 105 uint64_t floatx80_internal_sig = 0;
 106 int32_t floatx80_internal_exp0 = 0;
 107 uint64_t floatx80_internal_sig0 = 0;
 108 uint64_t floatx80_internal_sig1 = 0;
 109 int8_t floatx80_internal_precision = 80;
 110 int8_t floatx80_internal_mode = float_round_nearest_even;
 111
 112 /*----------------------------------------------------------------------------
 113  | Functions for storing sign, exponent and significand of extended
 114  | double-precision floating-point intermediate result for external use.
 115  *----------------------------------------------------------------------------*/
 116 floatx80 roundSaveFloatx80Internal( int8_t roundingPrecision, flag zSign, int32_t zExp, uint64_t zSig0, uint64_t zSig1, float_status *status )
 117 {
 118     uint64_t roundIncrement, roundMask, roundBits;
 119     flag increment;
 120
 121     if ( roundingPrecision == 80 ) {
 122         goto precision80;
 123     } else if ( roundingPrecision == 64 ) {
 124         roundIncrement = LIT64( 0x0000000000000400 );
 125         roundMask = LIT64( 0x00000000000007FF );
 126     } else if ( roundingPrecision == 32 ) {
 127         roundIncrement = LIT64( 0x0000008000000000 );
 128         roundMask = LIT64( 0x000000FFFFFFFFFF );
 129     } else {
 130         goto precision80;
 131     }
 132
 133     zSig0 |= ( zSig1 != 0 );
 134     if ( status->float_rounding_mode != float_round_nearest_even ) {
 135         if ( status->float_rounding_mode == float_round_to_zero ) {
 136             roundIncrement = 0;
 137         } else {
 138             roundIncrement = roundMask;
 139             if ( zSign ) {
 140                 if ( status->float_rounding_mode == float_round_up ) roundIncrement = 0;
 141             } else {
 142                 if ( status->float_rounding_mode == float_round_down ) roundIncrement = 0;
 143             }
 144         }
 145     }
 146
 147     roundBits = zSig0 & roundMask;
 148
 149     zSig0 += roundIncrement;
 150     if ( zSig0 < roundIncrement ) {
 151         ++zExp;
 152         zSig0 = LIT64( 0x8000000000000000 );
 153     }
 154     roundIncrement = roundMask + 1;
 155     if ( status->float_rounding_mode == float_round_nearest_even && ( roundBits<<1 == roundIncrement ) ) {
 156         roundMask |= roundIncrement;
 157     }
 158     zSig0 &= ~ roundMask;
 159     if ( zSig0 == 0 ) zExp = 0;
 160     return packFloatx80( zSign, zExp, zSig0 );
 161
 162 precision80:
 163     increment = ( (int64_t) zSig1 < 0 );
 164     if ( status->float_rounding_mode != float_round_nearest_even ) {
 165         if ( status->float_rounding_mode == float_round_to_zero ) {
 166             increment = 0;
 167         } else {
 168             if ( zSign ) {
 169                 increment = ( status->float_rounding_mode == float_round_down ) && zSig1;
 170             } else {
 171                 increment = ( status->float_rounding_mode == float_round_up ) && zSig1;
 172             }
 173         }
 174     }
 175     if ( increment ) {
 176         ++zSig0;
 177         if ( zSig0 == 0 ) {
 178             ++zExp;
 179             zSig0 = LIT64( 0x8000000000000000 );
 180         } else {
 181             if ((zSig1 << 1) == 0 && status->float_rounding_mode == float_round_nearest_even)
 182                 zSig0 &= ~1;
 183         }
 184     } else {
 185         if ( zSig0 == 0 ) zExp = 0;
 186     }
 187     return packFloatx80( zSign, zExp, zSig0 );
 188 }
 189
 190 static void saveFloatx80Internal( int8_t prec, flag zSign, int32_t zExp, uint64_t zSig0, uint64_t zSig1, float_status *status )
 191 {
 192     floatx80_internal_sign = zSign;
 193     floatx80_internal_exp = zExp;
 194     floatx80_internal_sig0 = zSig0;
 195     floatx80_internal_sig1 = zSig1;
 196     floatx80_internal_precision = prec;
 197     floatx80_internal_mode = status->float_rounding_mode;
 198 }
 199
 200 static void saveFloat64Internal( flag zSign, int16_t zExp, uint64_t zSig, float_status *status )
 201 {
 202     floatx80_internal_sign = zSign;
 203     floatx80_internal_exp = zExp + 0x3C01;
 204     floatx80_internal_sig0 = zSig<<1;
 205     floatx80_internal_sig1 = 0;
 206     floatx80_internal_precision = 64;
 207     floatx80_internal_mode = status->float_rounding_mode;
 208 }
 209
 210 static void saveFloat32Internal( flag zSign, int16_t zExp, uint32_t zSig, float_status *status )
 211 {
 212         floatx80 z = roundSaveFloatx80Internal( 32, zSign, zExp + 0x3F81, ( (uint64_t) zSig )<<33, 0, status );
 213
 214     floatx80_internal_sign = zSign;
 215     floatx80_internal_exp = extractFloatx80Exp( z );
 216     floatx80_internal_sig = extractFloatx80Frac( z );
 217     floatx80_internal_exp0 = zExp + 0x3F81;
 218     floatx80_internal_sig0 = ( (uint64_t) zSig )<<33;
 219     floatx80_internal_sig1 = 0;
 220 }
 221
 222 /*----------------------------------------------------------------------------
 223  | Functions for returning sign, exponent and significand of extended
 224  | double-precision floating-point intermediate result for external use.
 225  *----------------------------------------------------------------------------*/
 226
 227 void getRoundedFloatInternal( int8_t roundingPrecision, flag *pzSign, int32_t *pzExp, uint64_t *pzSig )
 228 {
 229     uint64_t roundIncrement, roundMask, roundBits;
 230     flag increment;
 231
 232     flag zSign = floatx80_internal_sign;
 233     int32_t zExp = floatx80_internal_exp;
 234     uint64_t zSig0 = floatx80_internal_sig0;
 235     uint64_t zSig1 = floatx80_internal_sig1;
 236
 237     if ( roundingPrecision == 80 ) {
 238         goto precision80;
 239     } else if ( roundingPrecision == 64 ) {
 240         roundIncrement = LIT64( 0x0000000000000400 );
 241         roundMask = LIT64( 0x00000000000007FF );
 242     } else if ( roundingPrecision == 32 ) {
 243         roundIncrement = LIT64( 0x0000008000000000 );
 244         roundMask = LIT64( 0x000000FFFFFFFFFF );
 245     } else {
 246         goto precision80;
 247     }
 248
 249     zSig0 |= ( zSig1 != 0 );
 250     if ( floatx80_internal_mode != float_round_nearest_even ) {
 251         if ( floatx80_internal_mode == float_round_to_zero ) {
 252             roundIncrement = 0;
 253         } else {
 254             roundIncrement = roundMask;
 255             if ( zSign ) {
 256                 if ( floatx80_internal_mode == float_round_up ) roundIncrement = 0;
 257             } else {
 258                 if ( floatx80_internal_mode == float_round_down ) roundIncrement = 0;
 259             }
 260         }
 261     }
 262
 263     roundBits = zSig0 & roundMask;
 264
 265     zSig0 += roundIncrement;
 266     if ( zSig0 < roundIncrement ) {
 267         ++zExp;
 268         zSig0 = LIT64( 0x8000000000000000 );
 269     }
 270     roundIncrement = roundMask + 1;
 271     if ( floatx80_internal_mode == float_round_nearest_even && ( roundBits<<1 == roundIncrement ) ) {
 272         roundMask |= roundIncrement;
 273     }
 274     zSig0 &= ~ roundMask;
 275     if ( zSig0 == 0 ) zExp = 0;
 276
 277     *pzSign = zSign;
 278     *pzExp = zExp;
 279     *pzSig = zSig0;
 280     return;
 281
 282 precision80:
 283     increment = ( (int64_t) zSig1 < 0 );
 284     if ( floatx80_internal_mode != float_round_nearest_even ) {
 285         if ( floatx80_internal_mode == float_round_to_zero ) {
 286             increment = 0;
 287         } else {
 288             if ( zSign ) {
 289                 increment = ( floatx80_internal_mode == float_round_down ) && zSig1;
 290             } else {
 291                 increment = ( floatx80_internal_mode == float_round_up ) && zSig1;
 292             }
 293         }
 294     }
 295     if ( increment ) {
 296         ++zSig0;
 297         if ( zSig0 == 0 ) {
 298             ++zExp;
 299             zSig0 = LIT64( 0x8000000000000000 );
 300         } else {
 301             if ((zSig1 << 1) == 0 && floatx80_internal_mode == float_round_nearest_even)
 302                 zSig0 &= ~1;
 303         }
 304     } else {
 305         if ( zSig0 == 0 ) zExp = 0;
 306     }
 307
 308     *pzSign = zSign;
 309     *pzExp = zExp;
 310     *pzSig = zSig0;
 311 }
 312
 313 floatx80 getFloatInternalOverflow( void )
 314 {
 315     flag zSign;
 316     int32_t zExp;
 317     uint64_t zSig;
 318
 319     getRoundedFloatInternal( floatx80_internal_precision, &zSign, &zExp, &zSig );
 320
 321     if (zExp > (0x7fff + 0x6000)) { // catastrophic
 322         zExp = 0;
 323     } else {
 324         zExp -= 0x6000;
 325     }
 326
 327     return packFloatx80( zSign, zExp, zSig );
 328
 329 }
 330
 331 floatx80 getFloatInternalUnderflow( void )
 332 {
 333     flag zSign;
 334     int32_t zExp;
 335     uint64_t zSig;
 336
 337     getRoundedFloatInternal( floatx80_internal_precision, &zSign, &zExp, &zSig );
 338
 339     if (zExp < (0x0000 - 0x6000)) { // catastrophic
 340         zExp = 0;
 341     } else {
 342         zExp += 0x6000;
 343     }
 344
 345     return packFloatx80( zSign, zExp, zSig );
 346
 347 }
 348
 349 floatx80 getFloatInternalRoundedAll( void )
 350 {
 351     flag zSign;
 352     int32_t zExp;
 353     uint64_t zSig, zSig32, zSig64, zSig80;
 354
 355     if (floatx80_internal_precision == 80) {
 356         getRoundedFloatInternal( 80, &zSign, &zExp, &zSig80 );
 357         zSig = zSig80;
 358     } else if (floatx80_internal_precision == 64) {
 359         getRoundedFloatInternal( 80, &zSign, &zExp, &zSig80 );
 360         getRoundedFloatInternal( 64, &zSign, &zExp, &zSig64 );
 361         zSig = zSig64;
 362         zSig |= zSig80 & LIT64( 0x00000000000007FF );
 363     } else {
 364         getRoundedFloatInternal( 80, &zSign, &zExp, &zSig80 );
 365         getRoundedFloatInternal( 64, &zSign, &zExp, &zSig64 );
 366         getRoundedFloatInternal( 32, &zSign, &zExp, &zSig32 );
 367         zSig = zSig32;
 368         zSig |= zSig64 & LIT64( 0x000000FFFFFFFFFF );
 369         zSig |= zSig80 & LIT64( 0x00000000000007FF );
 370     }
 371
 372     return packFloatx80( zSign, zExp & 0x7FFF, zSig );
 373
 374 }
 375
 376 floatx80 getFloatInternalRoundedSome( void )
 377 {
 378     flag zSign;
 379     int32_t zExp;
 380     uint64_t zSig, zSig32, zSig64, zSig80;
 381
 382     if (floatx80_internal_precision == 80) {
 383         getRoundedFloatInternal( 80, &zSign, &zExp, &zSig80 );
 384         zSig = zSig80;
 385     } else if (floatx80_internal_precision == 64) {
 386         getRoundedFloatInternal( 64, &zSign, &zExp, &zSig64 );
 387         zSig80 = floatx80_internal_sig0;
 388         if (zSig64 != (zSig80 & LIT64( 0xFFFFFFFFFFFFF800 ))) {
 389             zSig80++;
 390         }
 391         zSig = zSig64;
 392         zSig |= zSig80 & LIT64( 0x00000000000007FF );
 393     } else {
 394         getRoundedFloatInternal( 32, &zSign, &zExp, &zSig32 );
 395         zSig80 = floatx80_internal_sig0;
 396         if (zSig32 != (zSig80 & LIT64( 0xFFFFFF0000000000 ))) {
 397            zSig80++;
 398         }
 399         zSig = zSig32;
 400         zSig |= zSig80 & LIT64( 0x000000FFFFFFFFFF );
 401     }
 402
 403     return packFloatx80( zSign, zExp & 0x7FFF, zSig );
 404
 405 }
 406
 407 floatx80 getFloatInternalFloatx80( void )
 408 {
 409     flag zSign;
 410     int32_t zExp;
 411     uint64_t zSig;
 412
 413     getRoundedFloatInternal( 80, &zSign, &zExp, &zSig );
 414
 415     return packFloatx80( zSign, zExp & 0x7FFF, zSig );
 416
 417 }
 418
 419 floatx80 getFloatInternalUnrounded( void )
 420 {
 421     flag zSign = floatx80_internal_sign;
 422     int32_t zExp = floatx80_internal_exp;
 423     uint64_t zSig = floatx80_internal_sig0;
 424
 425     return packFloatx80( zSign, zExp & 0x7FFF, zSig );
 426
 427 }
 428
 429 uint64_t getFloatInternalGRS( void )
 430 {
 431 #if 1
 432     if (floatx80_internal_sig1)
 433         return 5;
 434
 435     if (floatx80_internal_precision == 64 &&
 436         floatx80_internal_sig0 & LIT64( 0x00000000000007FF )) {
 437         return 1;
 438     }
 439     if (floatx80_internal_precision == 32 &&
 440         floatx80_internal_sig0 & LIT64( 0x000000FFFFFFFFFF )) {
 441         return 1;
 442     }
 443
 444     return 0;
 445 #else
 446     uint64_t roundbits;
 447     shift64RightJamming(floatx80_internal_sig1, 61, &roundbits);
 448
 449     return roundbits;
 450 #endif
 451 }
 452
 453 /*----------------------------------------------------------------------------
 454 | Functions and definitions to determine:  (1) whether tininess for underflow
 455 | is detected before or after rounding by default, (2) what (if anything)
 456 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
 457 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
 458 | are propagated from function inputs to output.  These details are target-
 459 | specific.
 460 *----------------------------------------------------------------------------*/
 461 #include "softfloat-specialize.h"
 462
 463 /*----------------------------------------------------------------------------
 464 | Raises the exceptions specified by `flags'.  Floating-point traps can be
 465 | defined here if desired.  It is currently not possible for such a trap
 466 | to substitute a result value.  If traps are not implemented, this routine
 467 | should be simply `float_exception_flags |= flags;'.
 468 *----------------------------------------------------------------------------*/
 469
 470 void float_raise(uint8_t flags, float_status *status)
 471 {
 472     status->float_exception_flags |= flags;
 473 }
 474
 475
 476 /*----------------------------------------------------------------------------
 477 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
 478 | and 7, and returns the properly rounded 32-bit integer corresponding to the
 479 | input.  If `zSign' is 1, the input is negated before being converted to an
 480 | integer.  Bit 63 of `absZ' must be zero.  Ordinarily, the fixed-point input
 481 | is simply rounded to an integer, with the inexact exception raised if the
 482 | input cannot be represented exactly as an integer.  However, if the fixed-
 483 | point input is too large, the invalid exception is raised and the largest
 484 | positive or negative integer is returned.
 485 *----------------------------------------------------------------------------*/
 486
 487 static int32_t roundAndPackInt32(flag zSign, uint64_t absZ, float_status *status)
 488 {
 489     int8_t roundingMode;
 490     flag roundNearestEven;
 491     int8_t roundIncrement, roundBits;
 492     int32_t z;
 493
 494     roundingMode = status->float_rounding_mode;
 495     roundNearestEven = ( roundingMode == float_round_nearest_even );
 496     switch (roundingMode) {
 497     case float_round_nearest_even:
 498     case float_round_ties_away:
 499         roundIncrement = 0x40;
 500         break;
 501     case float_round_to_zero:
 502         roundIncrement = 0;
 503         break;
 504     case float_round_up:
 505         roundIncrement = zSign ? 0 : 0x7f;
 506         break;
 507     case float_round_down:
 508         roundIncrement = zSign ? 0x7f : 0;
 509         break;
 510     default:
 511         abort();
 512     }
 513     roundBits = absZ & 0x7F;
 514     absZ = ( absZ + roundIncrement )>>7;
 515     absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
 516     z = absZ;
 517     if ( zSign ) z = - z;
 518     if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
 519         float_raise(float_flag_invalid, status);
 520         return zSign ? (int32_t) 0x80000000 : 0x7FFFFFFF;
 521     }
 522     if (roundBits) {
 523         status->float_exception_flags |= float_flag_inexact;
 524     }
 525     return z;
 526
 527 }
 528
 529
 530 #ifdef SOFTFLOAT_68K // 30-01-2017: Added for Previous
 531 static int16_t roundAndPackInt16( flag zSign, uint64_t absZ, float_status *status )
 532 {
 533     int8_t roundingMode;
 534     flag roundNearestEven;
 535     int8_t roundIncrement, roundBits;
 536     int16_t z;
 537
 538     roundingMode = status->float_rounding_mode;
 539     roundNearestEven = ( roundingMode == float_round_nearest_even );
 540     roundIncrement = 0x40;
 541     if ( ! roundNearestEven ) {
 542         if ( roundingMode == float_round_to_zero ) {
 543             roundIncrement = 0;
 544         }
 545         else {
 546             roundIncrement = 0x7F;
 547             if ( zSign ) {
 548                 if ( roundingMode == float_round_up ) roundIncrement = 0;
 549             }
 550             else {
 551                 if ( roundingMode == float_round_down ) roundIncrement = 0;
 552             }
 553         }
 554     }
 555     roundBits = absZ & 0x7F;
 556     absZ = ( absZ + roundIncrement )>>7;
 557     absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
 558     z = absZ;
 559     if ( zSign ) z = - z;
 560     z = (int16_t) z;
 561     if ( ( absZ>>16 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
 562         float_raise( float_flag_invalid, status );
 563         return zSign ? (int16_t) 0x8000 : 0x7FFF;
 564     }
 565     if ( roundBits ) status->float_exception_flags |= float_flag_inexact;
 566     return z;
 567
 568 }
 569
 570 static int8_t roundAndPackInt8( flag zSign, uint64_t absZ, float_status *status )
 571 {
 572     int8_t roundingMode;
 573     flag roundNearestEven;
 574     int8_t roundIncrement, roundBits;
 575     int8_t z;
 576
 577     roundingMode = status->float_rounding_mode;
 578     roundNearestEven = ( roundingMode == float_round_nearest_even );
 579     roundIncrement = 0x40;
 580     if ( ! roundNearestEven ) {
 581         if ( roundingMode == float_round_to_zero ) {
 582             roundIncrement = 0;
 583         }
 584         else {
 585             roundIncrement = 0x7F;
 586             if ( zSign ) {
 587                 if ( roundingMode == float_round_up ) roundIncrement = 0;
 588             }
 589             else {
 590                 if ( roundingMode == float_round_down ) roundIncrement = 0;
 591             }
 592         }
 593     }
 594     roundBits = absZ & 0x7F;
 595     absZ = ( absZ + roundIncrement )>>7;
 596     absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
 597     z = absZ;
 598     if ( zSign ) z = - z;
 599     z = (int8_t) z;
 600     if ( ( absZ>>8 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
 601         float_raise( float_flag_invalid, status );
 602         return zSign ? (int8_t) 0x80 : 0x7F;
 603     }
 604     if ( roundBits ) status->float_exception_flags |= float_flag_inexact;
 605     return z;
 606
 607 }
 608 #endif // End of addition for Previous
 609
 610 /*----------------------------------------------------------------------------
 611 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
 612 | `absZ1', with binary point between bits 63 and 64 (between the input words),
 613 | and returns the properly rounded 64-bit integer corresponding to the input.
 614 | If `zSign' is 1, the input is negated before being converted to an integer.
 615 | Ordinarily, the fixed-point input is simply rounded to an integer, with
 616 | the inexact exception raised if the input cannot be represented exactly as
 617 | an integer.  However, if the fixed-point input is too large, the invalid
 618 | exception is raised and the largest positive or negative integer is
 619 | returned.
 620 *----------------------------------------------------------------------------*/
 621
 622 static int64_t roundAndPackInt64(flag zSign, uint64_t absZ0, uint64_t absZ1,
 623                                float_status *status)
 624 {
 625     int8_t roundingMode;
 626     flag roundNearestEven, increment;
 627     int64_t z;
 628
 629     roundingMode = status->float_rounding_mode;
 630     roundNearestEven = ( roundingMode == float_round_nearest_even );
 631     switch (roundingMode) {
 632     case float_round_nearest_even:
 633     case float_round_ties_away:
 634         increment = ((int64_t) absZ1 < 0);
 635         break;
 636     case float_round_to_zero:
 637         increment = 0;
 638         break;
 639     case float_round_up:
 640         increment = !zSign && absZ1;
 641         break;
 642     case float_round_down:
 643         increment = zSign && absZ1;
 644         break;
 645     default:
 646         abort();
 647     }
 648     if ( increment ) {
 649         ++absZ0;
 650         if ( absZ0 == 0 ) goto overflow;
 651         absZ0 &= ~ ( ( (uint64_t) ( absZ1<<1 ) == 0 ) & roundNearestEven );
 652     }
 653     z = absZ0;
 654     if ( zSign ) z = - z;
 655     if ( z && ( ( z < 0 ) ^ zSign ) ) {
 656  overflow:
 657         float_raise(float_flag_invalid, status);
 658         return
 659               zSign ? (uint64_t) LIT64( 0x8000000000000000 )
 660             : LIT64( 0x7FFFFFFFFFFFFFFF );
 661     }
 662     if (absZ1) {
 663         status->float_exception_flags |= float_flag_inexact;
 664     }
 665     return z;
 666
 667 }
 668
 669 /*----------------------------------------------------------------------------
 670 | Returns the fraction bits of the single-precision floating-point value `a'.
 671 *----------------------------------------------------------------------------*/
 672
 673 static inline uint32_t extractFloat32Frac( float32 a )
 674 {
 675
 676     return float32_val(a) & 0x007FFFFF;
 677
 678 }
 679
 680 /*----------------------------------------------------------------------------
 681 | Returns the exponent bits of the single-precision floating-point value `a'.
 682 *----------------------------------------------------------------------------*/
 683
 684 static inline int extractFloat32Exp(float32 a)
 685 {
 686
 687     return ( float32_val(a)>>23 ) & 0xFF;
 688
 689 }
 690
 691 /*----------------------------------------------------------------------------
 692 | Returns the sign bit of the single-precision floating-point value `a'.
 693 *----------------------------------------------------------------------------*/
 694
 695 static inline flag extractFloat32Sign( float32 a )
 696 {
 697
 698     return float32_val(a)>>31;
 699
 700 }
 701
 702 /*----------------------------------------------------------------------------
 703 | Normalizes the subnormal single-precision floating-point value represented
 704 | by the denormalized significand `aSig'.  The normalized exponent and
 705 | significand are stored at the locations pointed to by `zExpPtr' and
 706 | `zSigPtr', respectively.
 707 *----------------------------------------------------------------------------*/
 708
 709 static void
 710  normalizeFloat32Subnormal(uint32_t aSig, int *zExpPtr, uint32_t *zSigPtr)
 711 {
 712     int8_t shiftCount;
 713
 714     shiftCount = countLeadingZeros32( aSig ) - 8;
 715     *zSigPtr = aSig<<shiftCount;
 716     *zExpPtr = 1 - shiftCount;
 717
 718 }
 719
 720 /*----------------------------------------------------------------------------
 721 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
 722 | single-precision floating-point value, returning the result.  After being
 723 | shifted into the proper positions, the three fields are simply added
 724 | together to form the result.  This means that any integer portion of `zSig'
 725 | will be added into the exponent.  Since a properly normalized significand
 726 | will have an integer portion equal to 1, the `zExp' input should be 1 less
 727 | than the desired result exponent whenever `zSig' is a complete, normalized
 728 | significand.
 729 *----------------------------------------------------------------------------*/
 730
 731 static inline float32 packFloat32(flag zSign, int zExp, uint32_t zSig)
 732 {
 733
 734     return make_float32(
 735           ( ( (uint32_t) zSign )<<31 ) + ( ( (uint32_t) zExp )<<23 ) + zSig);
 736
 737 }
 738
 739 /*----------------------------------------------------------------------------
 740 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 741 | and significand `zSig', and returns the proper single-precision floating-
 742 | point value corresponding to the abstract input.  Ordinarily, the abstract
 743 | value is simply rounded and packed into the single-precision format, with
 744 | the inexact exception raised if the abstract input cannot be represented
 745 | exactly.  However, if the abstract value is too large, the overflow and
 746 | inexact exceptions are raised and an infinity or maximal finite value is
 747 | returned.  If the abstract value is too small, the input value is rounded to
 748 | a subnormal number, and the underflow and inexact exceptions are raised if
 749 | the abstract input cannot be represented exactly as a subnormal single-
 750 | precision floating-point number.
 751 |     The input significand `zSig' has its binary point between bits 30
 752 | and 29, which is 7 bits to the left of the usual location.  This shifted
 753 | significand must be normalized or smaller.  If `zSig' is not normalized,
 754 | `zExp' must be 0; in that case, the result returned is a subnormal number,
 755 | and it must not require rounding.  In the usual case that `zSig' is
 756 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
 757 | The handling of underflow and overflow follows the IEC/IEEE Standard for
 758 | Binary Floating-Point Arithmetic.
 759 *----------------------------------------------------------------------------*/
 760
 761 static float32 roundAndPackFloat32(flag zSign, int zExp, uint32_t zSig,
 762                                    float_status *status)
 763 {
 764     int8_t roundingMode;
 765     flag roundNearestEven;
 766     int8_t roundIncrement, roundBits;
 767     flag isTiny;
 768
 769     roundingMode = status->float_rounding_mode;
 770     roundNearestEven = ( roundingMode == float_round_nearest_even );
 771     switch (roundingMode) {
 772     case float_round_nearest_even:
 773     case float_round_ties_away:
 774         roundIncrement = 0x40;
 775         break;
 776     case float_round_to_zero:
 777         roundIncrement = 0;
 778         break;
 779     case float_round_up:
 780         roundIncrement = zSign ? 0 : 0x7f;
 781         break;
 782     case float_round_down:
 783         roundIncrement = zSign ? 0x7f : 0;
 784         break;
 785     default:
 786         abort();
 787         break;
 788     }
 789     roundBits = zSig & 0x7F;
 790     if ( 0xFD <= (uint16_t) zExp ) {
 791         if (    ( 0xFD < zExp )
 792              || (    ( zExp == 0xFD )
 793                   && ( (int32_t) ( zSig + roundIncrement ) < 0 ) )
 794            ) {
 795 #ifdef SOFTFLOAT_68K
 796             float_raise( float_flag_overflow, status );
 797                         saveFloat32Internal( zSign, zExp, zSig, status );
 798             if ( roundBits ) float_raise( float_flag_inexact, status );
 799 #else
 800             float_raise(float_flag_overflow | float_flag_inexact, status);
 801 #endif
 802             return packFloat32( zSign, 0xFF, - ( roundIncrement == 0 ));
 803         }
 804         if ( zExp < 0 ) {
 805             if (status->flush_to_zero) {
 806                 //float_raise(float_flag_output_denormal, status);
 807                 return packFloat32(zSign, 0, 0);
 808             }
 809             isTiny =
 810                 (status->float_detect_tininess
 811                  == float_tininess_before_rounding)
 812                 || ( zExp < -1 )
 813                 || ( zSig + roundIncrement < 0x80000000 );
 814 #ifdef SOFTFLOAT_68K
 815             if ( isTiny ) {
 816                 float_raise( float_flag_underflow, status );
 817                 saveFloat32Internal( zSign, zExp, zSig, status );
 818             }
 819 #endif
 820             shift32RightJamming( zSig, - zExp, &zSig );
 821             zExp = 0;
 822             roundBits = zSig & 0x7F;
 823 #ifndef SOFTFLOAT_68K
 824             if (isTiny && roundBits)
 825                 float_raise(float_flag_underflow, status);
 826 #endif
 827         }
 828     }
 829     if (roundBits) {
 830         status->float_exception_flags |= float_flag_inexact;
 831     }
 832     zSig = ( zSig + roundIncrement )>>7;
 833     zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
 834     if ( zSig == 0 ) zExp = 0;
 835     return packFloat32( zSign, zExp, zSig );
 836
 837 }
 838
 839 /*----------------------------------------------------------------------------
 840 | Returns the fraction bits of the double-precision floating-point value `a'.
 841 *----------------------------------------------------------------------------*/
 842
 843 static inline uint64_t extractFloat64Frac( float64 a )
 844 {
 845
 846     return float64_val(a) & LIT64( 0x000FFFFFFFFFFFFF );
 847
 848 }
 849
 850 /*----------------------------------------------------------------------------
 851 | Returns the exponent bits of the double-precision floating-point value `a'.
 852 *----------------------------------------------------------------------------*/
 853
 854 static inline int extractFloat64Exp(float64 a)
 855 {
 856
 857     return ( float64_val(a)>>52 ) & 0x7FF;
 858
 859 }
 860
 861 /*----------------------------------------------------------------------------
 862 | Returns the sign bit of the double-precision floating-point value `a'.
 863 *----------------------------------------------------------------------------*/
 864
 865 static inline flag extractFloat64Sign( float64 a )
 866 {
 867
 868     return float64_val(a)>>63;
 869
 870 }
 871
 872 /*----------------------------------------------------------------------------
 873 | If `a' is denormal and we are in flush-to-zero mode then set the
 874 | input-denormal exception and return zero. Otherwise just return the value.
 875 *----------------------------------------------------------------------------*/
 876 float64 float64_squash_input_denormal(float64 a, float_status *status)
 877 {
 878     if (status->flush_inputs_to_zero) {
 879         if (extractFloat64Exp(a) == 0 && extractFloat64Frac(a) != 0) {
 880             //float_raise(float_flag_input_denormal, status);
 881             return make_float64(float64_val(a) & (1ULL << 63));
 882         }
 883     }
 884     return a;
 885 }
 886
 887 /*----------------------------------------------------------------------------
 888 | Normalizes the subnormal double-precision floating-point value represented
 889 | by the denormalized significand `aSig'.  The normalized exponent and
 890 | significand are stored at the locations pointed to by `zExpPtr' and
 891 | `zSigPtr', respectively.
 892 *----------------------------------------------------------------------------*/
 893
 894 static void
 895  normalizeFloat64Subnormal(uint64_t aSig, int *zExpPtr, uint64_t *zSigPtr)
 896 {
 897     int8_t shiftCount;
 898
 899     shiftCount = countLeadingZeros64( aSig ) - 11;
 900     *zSigPtr = aSig<<shiftCount;
 901     *zExpPtr = 1 - shiftCount;
 902
 903 }
 904
 905 /*----------------------------------------------------------------------------
 906 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
 907 | double-precision floating-point value, returning the result.  After being
 908 | shifted into the proper positions, the three fields are simply added
 909 | together to form the result.  This means that any integer portion of `zSig'
 910 | will be added into the exponent.  Since a properly normalized significand
 911 | will have an integer portion equal to 1, the `zExp' input should be 1 less
 912 | than the desired result exponent whenever `zSig' is a complete, normalized
 913 | significand.
 914 *----------------------------------------------------------------------------*/
 915
 916 static inline float64 packFloat64(flag zSign, int zExp, uint64_t zSig)
 917 {
 918
 919     return make_float64(
 920         ( ( (uint64_t) zSign )<<63 ) + ( ( (uint64_t) zExp )<<52 ) + zSig);
 921
 922 }
 923
 924 /*----------------------------------------------------------------------------
 925 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 926 | and significand `zSig', and returns the proper double-precision floating-
 927 | point value corresponding to the abstract input.  Ordinarily, the abstract
 928 | value is simply rounded and packed into the double-precision format, with
 929 | the inexact exception raised if the abstract input cannot be represented
 930 | exactly.  However, if the abstract value is too large, the overflow and
 931 | inexact exceptions are raised and an infinity or maximal finite value is
 932 | returned.  If the abstract value is too small, the input value is rounded to
 933 | a subnormal number, and the underflow and inexact exceptions are raised if
 934 | the abstract input cannot be represented exactly as a subnormal double-
 935 | precision floating-point number.
 936 |     The input significand `zSig' has its binary point between bits 62
 937 | and 61, which is 10 bits to the left of the usual location.  This shifted
 938 | significand must be normalized or smaller.  If `zSig' is not normalized,
 939 | `zExp' must be 0; in that case, the result returned is a subnormal number,
 940 | and it must not require rounding.  In the usual case that `zSig' is
 941 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
 942 | The handling of underflow and overflow follows the IEC/IEEE Standard for
 943 | Binary Floating-Point Arithmetic.
 944 *----------------------------------------------------------------------------*/
 945
 946 static float64 roundAndPackFloat64(flag zSign, int zExp, uint64_t zSig,
 947                                    float_status *status)
 948 {
 949     int8_t roundingMode;
 950     flag roundNearestEven;
 951     int roundIncrement, roundBits;
 952     flag isTiny;
 953
 954     roundingMode = status->float_rounding_mode;
 955     roundNearestEven = ( roundingMode == float_round_nearest_even );
 956     switch (roundingMode) {
 957     case float_round_nearest_even:
 958     case float_round_ties_away:
 959         roundIncrement = 0x200;
 960         break;
 961     case float_round_to_zero:
 962         roundIncrement = 0;
 963         break;
 964     case float_round_up:
 965         roundIncrement = zSign ? 0 : 0x3ff;
 966         break;
 967     case float_round_down:
 968         roundIncrement = zSign ? 0x3ff : 0;
 969         break;
 970     default:
 971         abort();
 972     }
 973     roundBits = zSig & 0x3FF;
 974     if ( 0x7FD <= (uint16_t) zExp ) {
 975         if (    ( 0x7FD < zExp )
 976              || (    ( zExp == 0x7FD )
 977                   && ( (int64_t) ( zSig + roundIncrement ) < 0 ) )
 978            ) {
 979 #ifdef SOFTFLOAT_68K
 980                         float_raise( float_flag_overflow, status );
 981                         saveFloat64Internal( zSign, zExp, zSig, status );
 982             if ( roundBits ) float_raise( float_flag_inexact, status );
 983 #else
 984             float_raise(float_flag_overflow | float_flag_inexact, status);
 985 #endif
 986             return packFloat64( zSign, 0x7FF, - ( roundIncrement == 0 ));
 987         }
 988         if ( zExp < 0 ) {
 989             if (status->flush_to_zero) {
 990                 //float_raise(float_flag_output_denormal, status);
 991                 return packFloat64(zSign, 0, 0);
 992             }
 993             isTiny =
 994                    (status->float_detect_tininess
 995                     == float_tininess_before_rounding)
 996                 || ( zExp < -1 )
 997                 || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );
 998 #ifdef SOFTFLOAT_68K
 999             if ( isTiny ) {
1000                 float_raise( float_flag_underflow, status );
1001                 saveFloat64Internal( zSign, zExp, zSig, status );
1002             }
1003 #endif
1004             shift64RightJamming( zSig, - zExp, &zSig );
1005             zExp = 0;
1006             roundBits = zSig & 0x3FF;
1007 #ifndef SOFTFLOAT_68K
1008             if (isTiny && roundBits)
1009                 float_raise(float_flag_underflow, status);
1010 #endif
1011         }
1012     }
1013     if (roundBits) {
1014         status->float_exception_flags |= float_flag_inexact;
1015     }
1016     zSig = ( zSig + roundIncrement )>>10;
1017     zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );
1018     if ( zSig == 0 ) zExp = 0;
1019     return packFloat64( zSign, zExp, zSig );
1020
1021 }
1022
1023 /*----------------------------------------------------------------------------
1024 | Returns the fraction bits of the extended double-precision floating-point
1025 | value `a'.
1026 *----------------------------------------------------------------------------*/
1027
1028 uint64_t extractFloatx80Frac( floatx80 a )
1029 {
1030
1031     return a.low;
1032
1033 }
1034
1035 /*----------------------------------------------------------------------------
1036 | Returns the exponent bits of the extended double-precision floating-point
1037 | value `a'.
1038 *----------------------------------------------------------------------------*/
1039
1040 int32_t extractFloatx80Exp( floatx80 a )
1041 {
1042
1043     return a.high & 0x7FFF;
1044
1045 }
1046
1047 /*----------------------------------------------------------------------------
1048 | Returns the sign bit of the extended double-precision floating-point value
1049 | `a'.
1050 *----------------------------------------------------------------------------*/
1051
1052 flag extractFloatx80Sign( floatx80 a )
1053 {
1054
1055     return a.high>>15;
1056
1057 }
1058
1059 /*----------------------------------------------------------------------------
1060 | Normalizes the subnormal extended double-precision floating-point value
1061 | represented by the denormalized significand `aSig'.  The normalized exponent
1062 | and significand are stored at the locations pointed to by `zExpPtr' and
1063 | `zSigPtr', respectively.
1064 *----------------------------------------------------------------------------*/
1065
1066 void normalizeFloatx80Subnormal( uint64_t aSig, int32_t *zExpPtr, uint64_t *zSigPtr )
1067 {
1068     int8_t shiftCount;
1069
1070     shiftCount = countLeadingZeros64( aSig );
1071     *zSigPtr = aSig<<shiftCount;
1072 #ifdef SOFTFLOAT_68K
1073         *zExpPtr = -shiftCount;
1074 #else
1075         *zExpPtr = 1 - shiftCount;
1076 #endif
1077 }
1078
1079 /*----------------------------------------------------------------------------
1080 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
1081 | extended double-precision floating-point value, returning the result.
1082 *----------------------------------------------------------------------------*/
1083
1084 floatx80 packFloatx80( flag zSign, int32_t zExp, uint64_t zSig )
1085 {
1086     floatx80 z;
1087
1088     z.low = zSig;
1089     z.high = ( ( (uint16_t) zSign )<<15 ) + zExp;
1090     return z;
1091
1092 }
1093
1094 /*----------------------------------------------------------------------------
1095 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1096 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
1097 | and returns the proper extended double-precision floating-point value
1098 | corresponding to the abstract input.  Ordinarily, the abstract value is
1099 | rounded and packed into the extended double-precision format, with the
1100 | inexact exception raised if the abstract input cannot be represented
1101 | exactly.  However, if the abstract value is too large, the overflow and
1102 | inexact exceptions are raised and an infinity or maximal finite value is
1103 | returned.  If the abstract value is too small, the input value is rounded to
1104 | a subnormal number, and the underflow and inexact exceptions are raised if
1105 | the abstract input cannot be represented exactly as a subnormal extended
1106 | double-precision floating-point number.
1107 |     If `roundingPrecision' is 32 or 64, the result is rounded to the same
1108 | number of bits as single or double precision, respectively.  Otherwise, the
1109 | result is rounded to the full precision of the extended double-precision
1110 | format.
1111 |     The input significand must be normalized or smaller.  If the input
1112 | significand is not normalized, `zExp' must be 0; in that case, the result
1113 | returned is a subnormal number, and it must not require rounding.  The
1114 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
1115 | Floating-Point Arithmetic.
1116 *----------------------------------------------------------------------------*/
1117
1118 #ifndef SOFTFLOAT_68K
1119 floatx80 roundAndPackFloatx80(int8_t roundingPrecision, flag zSign,
1120                                      int32_t zExp, uint64_t zSig0, uint64_t zSig1,
1121                                      float_status *status)
1122 {
1123     int8_t roundingMode;
1124     flag roundNearestEven, increment, isTiny;
1125     int64_t roundIncrement, roundMask, roundBits;
1126
1127     roundingMode = status->float_rounding_mode;
1128     roundNearestEven = ( roundingMode == float_round_nearest_even );
1129     if ( roundingPrecision == 80 ) goto precision80;
1130     if ( roundingPrecision == 64 ) {
1131         roundIncrement = LIT64( 0x0000000000000400 );
1132         roundMask = LIT64( 0x00000000000007FF );
1133     }
1134     else if ( roundingPrecision == 32 ) {
1135         roundIncrement = LIT64( 0x0000008000000000 );
1136         roundMask = LIT64( 0x000000FFFFFFFFFF );
1137     }
1138     else {
1139         goto precision80;
1140     }
1141     zSig0 |= ( zSig1 != 0 );
1142     switch (roundingMode) {
1143     case float_round_nearest_even:
1144     case float_round_ties_away:
1145         break;
1146     case float_round_to_zero:
1147         roundIncrement = 0;
1148         break;
1149     case float_round_up:
1150         roundIncrement = zSign ? 0 : roundMask;
1151         break;
1152     case float_round_down:
1153         roundIncrement = zSign ? roundMask : 0;
1154         break;
1155     default:
1156         abort();
1157     }
1158     roundBits = zSig0 & roundMask;
1159 #ifdef SOFTFLOAT_68K
1160         if ( 0x7FFE <= (uint32_t) zExp ) {
1161 #else
1162         if ( 0x7FFD <= (uint32_t) ( zExp - 1 ) ) {
1163 #endif
1164         if (    ( 0x7FFE < zExp )
1165              || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
1166            ) {
1167             goto overflow;
1168         }
1169 #ifdef SOFTFLOAT_68K
1170         if ( zExp < 0 ) {
1171 #else
1172                 if ( zExp <= 0 ) {
1173 #endif
1174             if (status->flush_to_zero) {
1175                 //float_raise(float_flag_output_denormal, status);
1176                 return packFloatx80(zSign, 0, 0);
1177             }
1178             isTiny =
1179                    (status->float_detect_tininess
1180                     == float_tininess_before_rounding)
1181 #ifdef SOFTFLOAT_68K
1182                                 || ( zExp < -1 )
1183 #else
1184                                 || ( zExp < 0 )
1185 #endif
1186                                 || ( zSig0 <= zSig0 + roundIncrement );
1187 #ifdef SOFTFLOAT_68K
1188                         if ( isTiny ) {
1189                                 float_raise( float_flag_underflow, status );
1190                                 saveFloatx80Internal( zSign, zExp, zSig0, zSig1, status );
1191                         }
1192                         shift64RightJamming( zSig0, -zExp, &zSig0 );
1193 #else
1194                         shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
1195 #endif
1196             zExp = 0;
1197             roundBits = zSig0 & roundMask;
1198 #ifdef SOFTFLOAT_68K
1199             if ( isTiny ) float_raise( float_flag_underflow, status );
1200 #else
1201             if (isTiny && roundBits) {
1202                 float_raise(float_flag_underflow, status);
1203             }
1204 #endif
1205 if (roundBits) {
1206                 status->float_exception_flags |= float_flag_inexact;
1207             }
1208             zSig0 += roundIncrement;
1209 #ifndef SOFTFLOAT_68K
1210                         if ( (int64_t) zSig0 < 0 ) zExp = 1;
1211 #endif
1212             roundIncrement = roundMask + 1;
1213             if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
1214                 roundMask |= roundIncrement;
1215             }
1216             zSig0 &= ~ roundMask;
1217             return packFloatx80( zSign, zExp, zSig0 );
1218         }
1219     }
1220     if (roundBits) {
1221         status->float_exception_flags |= float_flag_inexact;
1222     }
1223     zSig0 += roundIncrement;
1224     if ( zSig0 < roundIncrement ) {
1225         ++zExp;
1226         zSig0 = LIT64( 0x8000000000000000 );
1227     }
1228     roundIncrement = roundMask + 1;
1229     if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
1230         roundMask |= roundIncrement;
1231     }
1232     zSig0 &= ~ roundMask;
1233     if ( zSig0 == 0 ) zExp = 0;
1234     return packFloatx80( zSign, zExp, zSig0 );
1235  precision80:
1236     switch (roundingMode) {
1237     case float_round_nearest_even:
1238     case float_round_ties_away:
1239         increment = ((int64_t)zSig1 < 0);
1240         break;
1241     case float_round_to_zero:
1242         increment = 0;
1243         break;
1244     case float_round_up:
1245         increment = !zSign && zSig1;
1246         break;
1247     case float_round_down:
1248         increment = zSign && zSig1;
1249         break;
1250     default:
1251         abort();
1252     }
1253 #ifdef SOFTFLOAT_68K
1254         if ( 0x7FFE <= (uint32_t) zExp ) {
1255 #else
1256         if ( 0x7FFD <= (uint32_t) ( zExp - 1 ) ) {
1257 #endif
1258         if (    ( 0x7FFE < zExp )
1259              || (    ( zExp == 0x7FFE )
1260                   && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
1261                   && increment
1262                 )
1263            ) {
1264             roundMask = 0;
1265  overflow:
1266 #ifndef SOFTFLOAT_68K
1267             float_raise(float_flag_overflow | float_flag_inexact, status);
1268 #else
1269             float_raise( float_flag_overflow, status );
1270                         saveFloatx80Internal( zSign, zExp, zSig0, zSig1, status );
1271             if ( ( zSig0 & roundMask ) || zSig1 ) float_raise( float_flag_inexact, status );
1272 #endif
1273                         if (    ( roundingMode == float_round_to_zero )
1274                  || ( zSign && ( roundingMode == float_round_up ) )
1275                  || ( ! zSign && ( roundingMode == float_round_down ) )
1276                ) {
1277                 return packFloatx80( zSign, 0x7FFE, ~ roundMask );
1278             }
1279             return packFloatx80( zSign, 0x7FFF, floatx80_default_infinity_low );
1280         }
1281 #ifdef SOFTFLOAT_68K
1282                 if ( zExp < 0 ) {
1283 #else
1284                 if ( zExp <= 0 ) {
1285 #endif
1286             isTiny =
1287                    (status->float_detect_tininess
1288                     == float_tininess_before_rounding)
1289 #ifdef SOFTFLOAT_68K
1290                                 || ( zExp < -1 )
1291 #else
1292                                 || ( zExp < 0 )
1293 #endif
1294                 || ! increment
1295                 || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );
1296 #ifdef SOFTFLOAT_68K
1297                         if ( isTiny ) {
1298                                 float_raise( float_flag_underflow, status );
1299                                 saveFloatx80Internal( zSign, zExp, zSig0, zSig1, status );
1300                         }
1301                         shift64ExtraRightJamming( zSig0, zSig1, -zExp, &zSig0, &zSig1 );
1302 #else
1303                         shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
1304 #endif
1305             zExp = 0;
1306 #ifndef SOFTFLOAT_68K
1307                         if ( isTiny && zSig1 ) float_raise( float_flag_underflow, status );
1308 #endif
1309             if (zSig1) float_raise(float_flag_inexact, status);
1310             switch (roundingMode) {
1311             case float_round_nearest_even:
1312             case float_round_ties_away:
1313                 increment = ((int64_t)zSig1 < 0);
1314                 break;
1315             case float_round_to_zero:
1316                 increment = 0;
1317                 break;
1318             case float_round_up:
1319                 increment = !zSign && zSig1;
1320                 break;
1321             case float_round_down:
1322                 increment = zSign && zSig1;
1323                 break;
1324             default:
1325                 abort();
1326             }
1327             if ( increment ) {
1328                 ++zSig0;
1329                 zSig0 &=
1330                     ~ ( ( (uint64_t) ( zSig1<<1 ) == 0 ) & roundNearestEven );
1331 #ifndef SOFTFLOAT_68K
1332                                 if ( (int64_t) zSig0 < 0 ) zExp = 1;
1333 #endif
1334             }
1335             return packFloatx80( zSign, zExp, zSig0 );
1336         }
1337     }
1338     if (zSig1) {
1339         status->float_exception_flags |= float_flag_inexact;
1340     }
1341     if ( increment ) {
1342         ++zSig0;
1343         if ( zSig0 == 0 ) {
1344             ++zExp;
1345             zSig0 = LIT64( 0x8000000000000000 );
1346         }
1347         else {
1348             zSig0 &= ~ ( ( (uint64_t) ( zSig1<<1 ) == 0 ) & roundNearestEven );
1349         }
1350     }
1351     else {
1352         if ( zSig0 == 0 ) zExp = 0;
1353     }
1354     return packFloatx80( zSign, zExp, zSig0 );
1355
1356 }
1357
1358 #else // SOFTFLOAT_68K
1359
1360 floatx80 roundAndPackFloatx80( int8_t roundingPrecision, flag zSign, int32_t zExp, uint64_t zSig0, uint64_t zSig1, float_status *status )
1361 {
1362     int8_t roundingMode;
1363     flag roundNearestEven, increment;
1364     uint64_t roundIncrement, roundMask, roundBits;
1365     int32_t expOffset;
1366
1367     roundingMode = status->float_rounding_mode;
1368     roundNearestEven = ( roundingMode == float_round_nearest_even );
1369     if ( roundingPrecision == 80 ) goto precision80;
1370     if ( roundingPrecision == 64 ) {
1371         roundIncrement = LIT64( 0x0000000000000400 );
1372         roundMask = LIT64( 0x00000000000007FF );
1373         expOffset = 0x3C00;
1374     } else if ( roundingPrecision == 32 ) {
1375         roundIncrement = LIT64( 0x0000008000000000 );
1376         roundMask = LIT64( 0x000000FFFFFFFFFF );
1377         expOffset = 0x3F80;
1378     } else {
1379         goto precision80;
1380     }
1381     zSig0 |= ( zSig1 != 0 );
1382     if ( ! roundNearestEven ) {
1383         if ( roundingMode == float_round_to_zero ) {
1384             roundIncrement = 0;
1385         } else {
1386             roundIncrement = roundMask;
1387             if ( zSign ) {
1388                 if ( roundingMode == float_round_up ) roundIncrement = 0;
1389             } else {
1390                 if ( roundingMode == float_round_down ) roundIncrement = 0;
1391             }
1392         }
1393     }
1394     roundBits = zSig0 & roundMask;
1395     if ( ( ( 0x7FFE - expOffset ) < zExp ) ||
1396         ( ( zExp == ( 0x7FFE - expOffset ) ) && ( zSig0 + roundIncrement < zSig0 ) ) ) {
1397         float_raise( float_flag_overflow, status );
1398         saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status );
1399         if ( zSig0 & roundMask ) float_raise( float_flag_inexact, status );
1400         if (    ( roundingMode == float_round_to_zero )
1401             || ( zSign && ( roundingMode == float_round_up ) )
1402             || ( ! zSign && ( roundingMode == float_round_down ) )
1403             ) {
1404             return packFloatx80( zSign, 0x7FFE - expOffset, ~ roundMask );
1405         }
1406         return packFloatx80( zSign, 0x7FFF, floatx80_default_infinity_low );
1407     }
1408     if ( zExp < ( expOffset + 1 ) ) {
1409         float_raise( float_flag_underflow, status );
1410         saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status );
1411         shift64RightJamming( zSig0, -( zExp - ( expOffset + 1 ) ), &zSig0 );
1412         zExp = expOffset + 1;
1413         roundBits = zSig0 & roundMask;
1414         if ( roundBits ) float_raise( float_flag_inexact, status );
1415         zSig0 += roundIncrement;
1416         roundIncrement = roundMask + 1;
1417         if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
1418             roundMask |= roundIncrement;
1419         }
1420         zSig0 &= ~ roundMask;
1421         return packFloatx80( zSign, zExp, zSig0 );
1422     }
1423     if ( roundBits ) {
1424         float_raise( float_flag_inexact, status );
1425         saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status);
1426     }
1427     zSig0 += roundIncrement;
1428     if ( zSig0 < roundIncrement ) {
1429         ++zExp;
1430         zSig0 = LIT64( 0x8000000000000000 );
1431     }
1432     roundIncrement = roundMask + 1;
1433     if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
1434         roundMask |= roundIncrement;
1435     }
1436     zSig0 &= ~ roundMask;
1437     if ( zSig0 == 0 ) zExp = 0;
1438     return packFloatx80( zSign, zExp, zSig0 );
1439 precision80:
1440     increment = ( (int64_t) zSig1 < 0 );
1441     if ( ! roundNearestEven ) {
1442         if ( roundingMode == float_round_to_zero ) {
1443             increment = 0;
1444         } else {
1445             if ( zSign ) {
1446                 increment = ( roundingMode == float_round_down ) && zSig1;
1447             } else {
1448                 increment = ( roundingMode == float_round_up ) && zSig1;
1449             }
1450         }
1451     }
1452     if ( 0x7FFE <= (uint32_t) zExp ) {
1453         if ( ( 0x7FFE < zExp ) ||
1454             ( ( zExp == 0x7FFE ) && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) ) && increment )
1455             ) {
1456             roundMask = 0;
1457             float_raise( float_flag_overflow, status );
1458             saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status );
1459             if ( ( zSig0 & roundMask ) || zSig1 ) float_raise( float_flag_inexact, status );
1460             if (    ( roundingMode == float_round_to_zero )
1461                 || ( zSign && ( roundingMode == float_round_up ) )
1462                 || ( ! zSign && ( roundingMode == float_round_down ) )
1463                 ) {
1464                 return packFloatx80( zSign, 0x7FFE, ~ roundMask );
1465             }
1466             return packFloatx80( zSign, 0x7FFF, floatx80_default_infinity_low );
1467         }
1468         if ( zExp < 0 ) {
1469             float_raise( float_flag_underflow, status );
1470             saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status);
1471             shift64ExtraRightJamming( zSig0, zSig1, -zExp, &zSig0, &zSig1 );
1472             zExp = 0;
1473             if ( zSig1 ) float_raise( float_flag_inexact, status );
1474             if ( roundNearestEven ) {
1475                 increment = ( (int64_t) zSig1 < 0 );
1476             } else {
1477                 if ( zSign ) {
1478                     increment = ( roundingMode == float_round_down ) && zSig1;
1479                 } else {
1480                     increment = ( roundingMode == float_round_up ) && zSig1;
1481                 }
1482             }
1483             if ( increment ) {
1484                 ++zSig0;
1485                 zSig0 &=
1486                 ~ ( ( (uint64_t) ( zSig1<<1 ) == 0 ) & roundNearestEven );
1487             }
1488             return packFloatx80( zSign, zExp, zSig0 );
1489         }
1490     }
1491     if ( zSig1 ) {
1492         float_raise( float_flag_inexact, status );
1493         saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status );
1494     }
1495     if ( increment ) {
1496         ++zSig0;
1497         if ( zSig0 == 0 ) {
1498             ++zExp;
1499             zSig0 = LIT64( 0x8000000000000000 );
1500         } else {
1501             zSig0 &= ~ ( ( (uint64_t) ( zSig1<<1 ) == 0 ) & roundNearestEven );
1502         }
1503     } else {
1504         if ( zSig0 == 0 ) zExp = 0;
1505     }
1506     return packFloatx80( zSign, zExp, zSig0 );
1507
1508 }
1509
1510 #endif
1511
1512 #ifdef SOFTFLOAT_68K // 21-01-2017: Added for Previous
1513 floatx80 roundSigAndPackFloatx80( int8_t roundingPrecision, flag zSign, int32_t zExp, uint64_t zSig0, uint64_t zSig1, float_status *status )
1514 {
1515     int8_t roundingMode;
1516     flag roundNearestEven, isTiny;
1517     uint64_t roundIncrement, roundMask, roundBits;
1518
1519     roundingMode = status->float_rounding_mode;
1520     roundNearestEven = ( roundingMode == float_round_nearest_even );
1521     if ( roundingPrecision == 32 ) {
1522         roundIncrement = LIT64( 0x0000008000000000 );
1523         roundMask = LIT64( 0x000000FFFFFFFFFF );
1524     } else if ( roundingPrecision == 64 ) {
1525         roundIncrement = LIT64( 0x0000000000000400 );
1526         roundMask = LIT64( 0x00000000000007FF );
1527     } else {
1528         return roundAndPackFloatx80( 80, zSign, zExp, zSig0, zSig1, status );
1529     }
1530     zSig0 |= ( zSig1 != 0 );
1531     if ( ! roundNearestEven ) {
1532         if ( roundingMode == float_round_to_zero ) {
1533             roundIncrement = 0;
1534         }
1535         else {
1536             roundIncrement = roundMask;
1537             if ( zSign ) {
1538                 if ( roundingMode == float_round_up ) roundIncrement = 0;
1539             }
1540             else {
1541                 if ( roundingMode == float_round_down ) roundIncrement = 0;
1542             }
1543         }
1544     }
1545     roundBits = zSig0 & roundMask;
1546
1547     if ( 0x7FFE <= (uint32_t) zExp ) {
1548         if (    ( 0x7FFE < zExp )
1549             || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
1550             ) {
1551             float_raise( float_flag_overflow, status );
1552                         saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status);
1553                         if ( zSig0 & roundMask ) float_raise( float_flag_inexact, status );
1554             if (    ( roundingMode == float_round_to_zero )
1555                 || ( zSign && ( roundingMode == float_round_up ) )
1556                 || ( ! zSign && ( roundingMode == float_round_down ) )
1557                 ) {
1558                 return packFloatx80( zSign, 0x7FFE, LIT64( 0xFFFFFFFFFFFFFFFF ) );
1559             }
1560             return packFloatx80( zSign, 0x7FFF, floatx80_default_infinity_low );
1561         }
1562
1563         if ( zExp < 0 ) {
1564             isTiny =
1565             ( status->float_detect_tininess == float_tininess_before_rounding )
1566             || ( zExp < -1 )
1567             || ( zSig0 <= zSig0 + roundIncrement );
1568                         if ( isTiny ) {
1569                                 float_raise( float_flag_underflow, status );
1570                                 saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status );
1571                         }
1572             shift64RightJamming( zSig0, -zExp, &zSig0 );
1573             zExp = 0;
1574             roundBits = zSig0 & roundMask;
1575             if ( roundBits ) float_raise ( float_flag_inexact, status );
1576             zSig0 += roundIncrement;
1577             if ( roundNearestEven && ( roundBits == roundIncrement ) ) {
1578                 roundMask |= roundIncrement<<1;
1579             }
1580             zSig0 &= ~roundMask;
1581             return packFloatx80( zSign, zExp, zSig0 );
1582         }
1583     }
1584     if ( roundBits ) {
1585         float_raise( float_flag_inexact, status );
1586         saveFloatx80Internal( roundingPrecision, zSign, zExp, zSig0, zSig1, status );
1587     }
1588     zSig0 += roundIncrement;
1589     if ( zSig0 < roundIncrement ) {
1590         ++zExp;
1591         zSig0 = LIT64( 0x8000000000000000 );
1592     }
1593     roundIncrement = roundMask + 1;
1594     if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
1595         roundMask |= roundIncrement;
1596     }
1597     zSig0 &= ~ roundMask;
1598     if ( zSig0 == 0 ) zExp = 0;
1599     return packFloatx80( zSign, zExp, zSig0 );
1600
1601 }
1602 #endif // End of Addition for Previous
1603
1604
1605 /*----------------------------------------------------------------------------
1606 | Takes an abstract floating-point value having sign `zSign', exponent
1607 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
1608 | and returns the proper extended double-precision floating-point value
1609 | corresponding to the abstract input.  This routine is just like
1610 | `roundAndPackFloatx80' except that the input significand does not have to be
1611 | normalized.
1612 *----------------------------------------------------------------------------*/
1613
1614 static floatx80 normalizeRoundAndPackFloatx80(int8_t roundingPrecision,
1615                                               flag zSign, int32_t zExp,
1616                                               uint64_t zSig0, uint64_t zSig1,
1617                                               float_status *status)
1618 {
1619     int8_t shiftCount;
1620
1621     if ( zSig0 == 0 ) {
1622         zSig0 = zSig1;
1623         zSig1 = 0;
1624         zExp -= 64;
1625     }
1626     shiftCount = countLeadingZeros64( zSig0 );
1627     shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
1628     zExp -= shiftCount;
1629     return roundAndPackFloatx80(roundingPrecision, zSign, zExp,
1630                                 zSig0, zSig1, status);
1631
1632 }
1633
1634 /*----------------------------------------------------------------------------
1635 | Returns the result of converting the 32-bit two's complement integer `a'
1636 | to the extended double-precision floating-point format.  The conversion
1637 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1638 | Arithmetic.
1639 *----------------------------------------------------------------------------*/
1640
1641 floatx80 int32_to_floatx80(int32_t a)
1642 {
1643     flag zSign;
1644     uint32_t absA;
1645     int8_t shiftCount;
1646     uint64_t zSig;
1647
1648     if ( a == 0 ) return packFloatx80( 0, 0, 0 );
1649     zSign = ( a < 0 );
1650     absA = zSign ? - a : a;
1651     shiftCount = countLeadingZeros32( absA ) + 32;
1652     zSig = absA;
1653     return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );
1654
1655 }
1656
1657 /*----------------------------------------------------------------------------
1658 | Returns the result of converting the single-precision floating-point value
1659 | `a' to the extended double-precision floating-point format.  The conversion
1660 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1661 | Arithmetic.
1662 *----------------------------------------------------------------------------*/
1663
1664 floatx80 float32_to_floatx80(float32 a, float_status *status)
1665 {
1666     flag aSign;
1667     int aExp;
1668     uint32_t aSig;
1669
1670     aSig = extractFloat32Frac( a );
1671     aExp = extractFloat32Exp( a );
1672     aSign = extractFloat32Sign( a );
1673     if ( aExp == 0xFF ) {
1674                 if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a, status ), status );
1675                 return packFloatx80( aSign, 0x7FFF, floatx80_default_infinity_low );
1676     }
1677     if ( aExp == 0 ) {
1678         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
1679         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1680     }
1681     aSig |= 0x00800000;
1682     return packFloatx80( aSign, aExp + 0x3F80, ( (uint64_t) aSig )<<40 );
1683
1684 }
1685
1686 #ifdef SOFTFLOAT_68K // 31-12-2016: Added for Previous
1687 floatx80 float32_to_floatx80_allowunnormal(float32 a , float_status *status)
1688 {
1689         (void)status;
1690     flag aSign;
1691     int16_t aExp;
1692     uint32_t aSig;
1693
1694     aSig = extractFloat32Frac(a);
1695     aExp = extractFloat32Exp(a);
1696     aSign = extractFloat32Sign(a);
1697     if (aExp == 0xFF) {
1698         return packFloatx80( aSign, 0x7FFF, ( (uint64_t) aSig )<<40 );
1699     }
1700     if (aExp == 0) {
1701         if (aSig == 0) return packFloatx80(aSign, 0, 0);
1702         return packFloatx80(aSign, 0x3F81, ((uint64_t) aSig) << 40);
1703     }
1704     aSig |= 0x00800000;
1705     return packFloatx80(aSign, aExp + 0x3F80, ((uint64_t)aSig) << 40);
1706
1707 }
1708 #endif // end of addition for Previous
1709
1710 /*----------------------------------------------------------------------------
1711 | Returns the result of converting the double-precision floating-point value
1712 | `a' to the extended double-precision floating-point format.  The conversion
1713 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1714 | Arithmetic.
1715 *----------------------------------------------------------------------------*/
1716
1717 floatx80 float64_to_floatx80(float64 a, float_status *status)
1718 {
1719     flag aSign;
1720     int aExp;
1721     uint64_t aSig;
1722
1723     aSig = extractFloat64Frac( a );
1724     aExp = extractFloat64Exp( a );
1725     aSign = extractFloat64Sign( a );
1726     if ( aExp == 0x7FF ) {
1727         if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a, status ), status );
1728         return packFloatx80( aSign, 0x7FFF, floatx80_default_infinity_low );
1729     }
1730     if ( aExp == 0 ) {
1731         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
1732         normalizeFloat64Subnormal( aSig, &aExp, &aSig );
1733     }
1734     return
1735         packFloatx80(
1736             aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );
1737
1738 }
1739
1740 #ifdef SOFTFLOAT_68K // 31-12-2016: Added for Previous
1741 floatx80 float64_to_floatx80_allowunnormal( float64 a, float_status *status )
1742 {
1743         (void)status;
1744     flag aSign;
1745     int16_t aExp;
1746     uint64_t aSig;
1747
1748     aSig = extractFloat64Frac( a );
1749     aExp = extractFloat64Exp( a );
1750     aSign = extractFloat64Sign( a );
1751     if ( aExp == 0x7FF ) {
1752         return packFloatx80( aSign, 0x7FFF, aSig<<11 );
1753     }
1754     if ( aExp == 0 ) {
1755         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
1756         return packFloatx80( aSign, 0x3C01, aSig<<11 );
1757     }
1758     return
1759     packFloatx80(
1760                  aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );
1761
1762 }
1763 #endif // end of addition for Previous
1764
1765 /*----------------------------------------------------------------------------
1766 | Returns the result of converting the extended double-precision floating-
1767 | point value `a' to the 32-bit two's complement integer format.  The
1768 | conversion is performed according to the IEC/IEEE Standard for Binary
1769 | Floating-Point Arithmetic---which means in particular that the conversion
1770 | is rounded according to the current rounding mode.  If `a' is a NaN, the
1771 | largest positive integer is returned.  Otherwise, if the conversion
1772 | overflows, the largest integer with the same sign as `a' is returned.
1773 *----------------------------------------------------------------------------*/
1774
1775 int32_t floatx80_to_int32(floatx80 a, float_status *status)
1776 {
1777     flag aSign;
1778     int32_t aExp, shiftCount;
1779     uint64_t aSig;
1780
1781     aSig = extractFloatx80Frac( a );
1782     aExp = extractFloatx80Exp( a );
1783     aSign = extractFloatx80Sign( a );
1784 #ifdef SOFTFLOAT_68K
1785     if ( aExp == 0x7FFF ) {
1786                 if ( (uint64_t) ( aSig<<1 ) ) {
1787             a = propagateFloatx80NaNOneArg( a, status );
1788             if ( a.low == aSig ) float_raise( float_flag_invalid, status );
1789             return (int32_t)(a.low>>32);
1790         }
1791                 float_raise( float_flag_invalid, status );
1792         return aSign ? (int32_t) 0x80000000 : 0x7FFFFFFF;
1793     }
1794 #else
1795         if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
1796 #endif
1797     shiftCount = 0x4037 - aExp;
1798     if ( shiftCount <= 0 ) shiftCount = 1;
1799     shift64RightJamming( aSig, shiftCount, &aSig );
1800     return roundAndPackInt32(aSign, aSig, status);
1801
1802 }
1803
1804 #ifdef SOFTFLOAT_68K // 30-01-2017: Addition for Previous
1805 int16_t floatx80_to_int16( floatx80 a, float_status *status)
1806 {
1807     flag aSign;
1808     int32_t aExp, shiftCount;
1809     uint64_t aSig;
1810
1811     aSig = extractFloatx80Frac( a );
1812     aExp = extractFloatx80Exp( a );
1813     aSign = extractFloatx80Sign( a );
1814     if ( aExp == 0x7FFF ) {
1815         float_raise( float_flag_invalid, status );
1816         if ( (uint64_t) ( aSig<<1 ) ) {
1817             a = propagateFloatx80NaNOneArg( a, status );
1818             if ( a.low == aSig ) float_raise( float_flag_invalid, status );
1819             return (int16_t)(a.low>>48);
1820         }
1821         return aSign ? (int16_t) 0x8000 : 0x7FFF;
1822     }
1823     shiftCount = 0x4037 - aExp;
1824     if ( shiftCount <= 0 ) shiftCount = 1;
1825     shift64RightJamming( aSig, shiftCount, &aSig );
1826     return roundAndPackInt16( aSign, aSig, status );
1827
1828 }
1829 int8_t floatx80_to_int8( floatx80 a, float_status *status)
1830 {
1831     flag aSign;
1832     int32_t aExp, shiftCount;
1833     uint64_t aSig;
1834
1835     aSig = extractFloatx80Frac( a );
1836     aExp = extractFloatx80Exp( a );
1837     aSign = extractFloatx80Sign( a );
1838     if ( aExp == 0x7FFF ) {
1839          if ( (uint64_t) ( aSig<<1 ) ) {
1840             a = propagateFloatx80NaNOneArg( a, status );
1841             if ( a.low == aSig ) float_raise( float_flag_invalid, status );
1842             return (int8_t)(a.low>>56);
1843         }
1844         float_raise( float_flag_invalid, status );
1845         return aSign ? (int8_t) 0x80 : 0x7F;
1846     }
1847     shiftCount = 0x4037 - aExp;
1848     if ( shiftCount <= 0 ) shiftCount = 1;
1849     shift64RightJamming( aSig, shiftCount, &aSig );
1850     return roundAndPackInt8( aSign, aSig, status );
1851
1852 }
1853 #endif // End of addition for Previous
1854
1855
1856 /*----------------------------------------------------------------------------
1857 | Returns the result of converting the extended double-precision floating-
1858 | point value `a' to the 32-bit two's complement integer format.  The
1859 | conversion is performed according to the IEC/IEEE Standard for Binary
1860 | Floating-Point Arithmetic, except that the conversion is always rounded
1861 | toward zero.  If `a' is a NaN, the largest positive integer is returned.
1862 | Otherwise, if the conversion overflows, the largest integer with the same
1863 | sign as `a' is returned.
1864 *----------------------------------------------------------------------------*/
1865
1866 int32_t floatx80_to_int32_round_to_zero(floatx80 a, float_status *status)
1867 {
1868     flag aSign;
1869     int32_t aExp, shiftCount;
1870     uint64_t aSig, savedASig;
1871     int32_t z;
1872
1873     if (floatx80_invalid_encoding(a)) {
1874         float_raise(float_flag_invalid, status);
1875         return 1 << 31;
1876     }
1877     aSig = extractFloatx80Frac( a );
1878     aExp = extractFloatx80Exp( a );
1879     aSign = extractFloatx80Sign( a );
1880     if ( 0x401E < aExp ) {
1881         if ( ( aExp == 0x7FFF ) && (uint64_t) ( aSig<<1 ) ) aSign = 0;
1882         goto invalid;
1883     }
1884     else if ( aExp < 0x3FFF ) {
1885         if (aExp || aSig) {
1886             status->float_exception_flags |= float_flag_inexact;
1887         }
1888         return 0;
1889     }
1890     shiftCount = 0x403E - aExp;
1891     savedASig = aSig;
1892     aSig >>= shiftCount;
1893     z = aSig;
1894     if ( aSign ) z = - z;
1895     if ( ( z < 0 ) ^ aSign ) {
1896  invalid:
1897         float_raise(float_flag_invalid, status);
1898         return aSign ? (int32_t) 0x80000000 : 0x7FFFFFFF;
1899     }
1900     if ( ( aSig<<shiftCount ) != savedASig ) {
1901         status->float_exception_flags |= float_flag_inexact;
1902     }
1903     return z;
1904
1905 }
1906
1907 /*----------------------------------------------------------------------------
1908 | Returns the result of converting the extended double-precision floating-
1909 | point value `a' to the 64-bit two's complement integer format.  The
1910 | conversion is performed according to the IEC/IEEE Standard for Binary
1911 | Floating-Point Arithmetic---which means in particular that the conversion
1912 | is rounded according to the current rounding mode.  If `a' is a NaN,
1913 | the largest positive integer is returned.  Otherwise, if the conversion
1914 | overflows, the largest integer with the same sign as `a' is returned.
1915 *----------------------------------------------------------------------------*/
1916
1917 int64_t floatx80_to_int64(floatx80 a, float_status *status)
1918 {
1919     flag aSign;
1920     int32_t aExp, shiftCount;
1921     uint64_t aSig, aSigExtra;
1922
1923     if (floatx80_invalid_encoding(a)) {
1924         float_raise(float_flag_invalid, status);
1925         return 1ULL << 63;
1926     }
1927     aSig = extractFloatx80Frac( a );
1928     aExp = extractFloatx80Exp( a );
1929     aSign = extractFloatx80Sign( a );
1930     shiftCount = 0x403E - aExp;
1931     if ( shiftCount <= 0 ) {
1932         if ( shiftCount ) {
1933             float_raise(float_flag_invalid, status);
1934             if (    ! aSign
1935                  || (    ( aExp == 0x7FFF )
1936                       && ( aSig != LIT64( 0x8000000000000000 ) ) )
1937                ) {
1938                 return LIT64( 0x7FFFFFFFFFFFFFFF );
1939             }
1940             return (int64_t) LIT64( 0x8000000000000000 );
1941         }
1942         aSigExtra = 0;
1943     }
1944     else {
1945         shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
1946     }
1947     return roundAndPackInt64(aSign, aSig, aSigExtra, status);
1948
1949 }
1950
1951 /*----------------------------------------------------------------------------
1952 | Returns the result of converting the extended double-precision floating-
1953 | point value `a' to the single-precision floating-point format.  The
1954 | conversion is performed according to the IEC/IEEE Standard for Binary
1955 | Floating-Point Arithmetic.
1956 *----------------------------------------------------------------------------*/
1957
1958 float32 floatx80_to_float32(floatx80 a, float_status *status)
1959 {
1960     flag aSign;
1961     int32_t aExp;
1962     uint64_t aSig;
1963
1964     aSig = extractFloatx80Frac( a );
1965     aExp = extractFloatx80Exp( a );
1966     aSign = extractFloatx80Sign( a );
1967     if ( aExp == 0x7FFF ) {
1968         if ( (uint64_t) ( aSig<<1 ) ) {
1969             return commonNaNToFloat32(floatx80ToCommonNaN(a, status));
1970         }
1971         return packFloat32( aSign, 0xFF, 0 );
1972     }
1973 #ifdef SOFTFLOAT_68K
1974     if ( aExp == 0 ) {
1975         if ( aSig == 0) return packFloat32( aSign, 0, 0 );
1976         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
1977     }
1978     shift64RightJamming( aSig, 33, &aSig );
1979     aExp -= 0x3F81;
1980 #else
1981     shift64RightJamming( aSig, 33, &aSig );
1982     if ( aExp || aSig ) aExp -= 0x3F81;
1983 #endif
1984     return roundAndPackFloat32(aSign, aExp, aSig, status);
1985
1986 }
1987
1988 /*----------------------------------------------------------------------------
1989 | Returns the result of converting the extended double-precision floating-
1990 | point value `a' to the double-precision floating-point format.  The
1991 | conversion is performed according to the IEC/IEEE Standard for Binary
1992 | Floating-Point Arithmetic.
1993 *----------------------------------------------------------------------------*/
1994
1995 float64 floatx80_to_float64(floatx80 a, float_status *status)
1996 {
1997     flag aSign;
1998     int32_t aExp;
1999     uint64_t aSig, zSig;
2000
2001     aSig = extractFloatx80Frac( a );
2002     aExp = extractFloatx80Exp( a );
2003     aSign = extractFloatx80Sign( a );
2004     if ( aExp == 0x7FFF ) {
2005         if ( (uint64_t) ( aSig<<1 ) ) {
2006             return commonNaNToFloat64(floatx80ToCommonNaN(a, status), status);
2007         }
2008         return packFloat64( aSign, 0x7FF, 0 );
2009     }
2010 #ifdef SOFTFLOAT_68K
2011     if ( aExp == 0 ) {
2012         if ( aSig == 0) return packFloat64( aSign, 0, 0 );
2013         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2014     }
2015     shift64RightJamming( aSig, 1, &zSig );
2016     aExp -= 0x3C01;
2017 #else
2018     shift64RightJamming( aSig, 1, &zSig );
2019     if ( aExp || aSig ) aExp -= 0x3C01;
2020 #endif
2021     return roundAndPackFloat64(aSign, aExp, zSig, status);
2022
2023 }
2024
2025 #ifdef SOFTFLOAT_68K // 31-01-2017
2026 /*----------------------------------------------------------------------------
2027  | Returns the result of converting the extended double-precision floating-
2028  | point value `a' to the extended double-precision floating-point format.
2029  | The conversion is performed according to the IEC/IEEE Standard for Binary
2030  | Floating-Point Arithmetic.
2031  *----------------------------------------------------------------------------*/
2032
2033 floatx80 floatx80_to_floatx80( floatx80 a, float_status *status )
2034 {
2035     flag aSign;
2036     int32_t aExp;
2037     uint64_t aSig;
2038
2039     aSig = extractFloatx80Frac( a );
2040     aExp = extractFloatx80Exp( a );
2041     aSign = extractFloatx80Sign( a );
2042
2043     if ( aExp == 0x7FFF && (uint64_t) ( aSig<<1 ) ) {
2044         return propagateFloatx80NaNOneArg( a, status );
2045     }
2046     if ( aExp == 0 && aSig != 0 ) {
2047         return normalizeRoundAndPackFloatx80( status->floatx80_rounding_precision, aSign, aExp, aSig, 0, status );
2048     }
2049     return a;
2050
2051 }
2052 #endif
2053
2054 #ifdef SOFTFLOAT_68K // 30-01-2016: Added for Previous
2055 floatx80 floatx80_round32( floatx80 a, float_status *status )
2056 {
2057     flag aSign;
2058     int32_t aExp;
2059     uint64_t aSig;
2060
2061     aSig = extractFloatx80Frac( a );
2062     aExp = extractFloatx80Exp( a );
2063     aSign = extractFloatx80Sign( a );
2064
2065     if ( aExp == 0x7FFF || aSig == 0 ) {
2066         return a;
2067     }
2068     if ( aExp == 0 ) {
2069         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2070     }
2071
2072     return roundSigAndPackFloatx80( 32, aSign, aExp, aSig, 0, status );
2073 }
2074
2075 floatx80 floatx80_round64( floatx80 a, float_status *status )
2076 {
2077     flag aSign;
2078     int32_t aExp;
2079     uint64_t aSig;
2080
2081     aSig = extractFloatx80Frac( a );
2082     aExp = extractFloatx80Exp( a );
2083     aSign = extractFloatx80Sign( a );
2084
2085     if ( aExp == 0x7FFF || aSig == 0 ) {
2086         return a;
2087     }
2088     if ( aExp == 0 ) {
2089         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2090     }
2091
2092     return roundSigAndPackFloatx80( 64, aSign, aExp, aSig, 0, status );
2093 }
2094
2095 floatx80 floatx80_round_to_float32( floatx80 a, float_status *status )
2096 {
2097         flag aSign;
2098     int32_t aExp;
2099     uint64_t aSig;
2100
2101         aSign = extractFloatx80Sign( a );
2102     aSig = extractFloatx80Frac( a );
2103     aExp = extractFloatx80Exp( a );
2104
2105     if ( aExp == 0x7FFF ) {
2106         if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaNOneArg( a, status );
2107         return a;
2108     }
2109     if ( aExp == 0 ) {
2110         if ( aSig == 0 ) return a;
2111         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2112     }
2113
2114     return roundAndPackFloatx80( 32, aSign, aExp, aSig, 0, status );
2115 }
2116
2117 floatx80 floatx80_round_to_float64( floatx80 a, float_status *status )
2118 {
2119     flag aSign;
2120     int32_t aExp;
2121     uint64_t aSig;
2122
2123         aSign = extractFloatx80Sign( a );
2124     aSig = extractFloatx80Frac( a );
2125     aExp = extractFloatx80Exp( a );
2126
2127     if ( aExp == 0x7FFF ) {
2128         if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaNOneArg( a, status );
2129         return a;
2130     }
2131     if ( aExp == 0 ) {
2132         if ( aSig == 0 ) return a;
2133         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2134     }
2135
2136     return roundAndPackFloatx80( 64, aSign, aExp, aSig, 0, status );
2137 }
2138
2139
2140 floatx80 floatx80_normalize( floatx80 a )
2141 {
2142     flag aSign;
2143     int16_t aExp;
2144     uint64_t aSig;
2145     int8_t shiftCount;
2146
2147     aSig = extractFloatx80Frac( a );
2148     aExp = extractFloatx80Exp( a );
2149     aSign = extractFloatx80Sign( a );
2150
2151     if ( aExp == 0x7FFF || aExp == 0 ) return a;
2152     if ( aSig == 0 ) return packFloatx80(aSign, 0, 0);
2153
2154     shiftCount = countLeadingZeros64( aSig );
2155
2156     if ( shiftCount > aExp ) shiftCount = aExp;
2157
2158     aExp -= shiftCount;
2159     aSig <<= shiftCount;
2160
2161     return packFloatx80( aSign, aExp, aSig );
2162 }
2163 #endif // end of addition for Previous
2164
2165 /*----------------------------------------------------------------------------
2166 | Rounds the extended double-precision floating-point value `a' to an integer,
2167 | and returns the result as an extended quadruple-precision floating-point
2168 | value.  The operation is performed according to the IEC/IEEE Standard for
2169 | Binary Floating-Point Arithmetic.
2170 *----------------------------------------------------------------------------*/
2171
2172 floatx80 floatx80_round_to_int(floatx80 a, float_status *status)
2173 {
2174     flag aSign;
2175     int32_t aExp;
2176     uint64_t lastBitMask, roundBitsMask;
2177 //      int8_t roundingMode;
2178         floatx80 z;
2179
2180 //      roundingMode = status->float_rounding_mode;
2181         aSign = extractFloatx80Sign(a);
2182     aExp = extractFloatx80Exp( a );
2183     if ( 0x403E <= aExp ) {
2184         if ( aExp == 0x7FFF ) {
2185                         if ((uint64_t) ( extractFloatx80Frac( a )<<1 ) )
2186                                 return propagateFloatx80NaNOneArg(a, status);
2187                         return inf_clear_intbit(status) ? packFloatx80(aSign, aExp, 0) : a;
2188         }
2189         return a;
2190     }
2191     if ( aExp < 0x3FFF ) {
2192         if (    ( aExp == 0 )
2193  #ifdef SOFTFLOAT_68K
2194                                 && ( (uint64_t) extractFloatx80Frac( a ) == 0 ) ) {
2195 #else
2196                                 && ( (uint64_t) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
2197 #endif
2198            return a;
2199         }
2200         status->float_exception_flags |= float_flag_inexact;
2201         switch (status->float_rounding_mode) {
2202          case float_round_nearest_even:
2203             if ( ( aExp == 0x3FFE ) && (uint64_t) ( extractFloatx80Frac( a )<<1 )
2204                ) {
2205                 return
2206                     packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
2207             }
2208             break;
2209         case float_round_ties_away:
2210             if (aExp == 0x3FFE) {
2211                 return packFloatx80(aSign, 0x3FFF, LIT64(0x8000000000000000));
2212             }
2213             break;
2214          case float_round_down:
2215             return
2216                   aSign ?
2217                       packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
2218                 : packFloatx80( 0, 0, 0 );
2219          case float_round_up:
2220             return
2221                   aSign ? packFloatx80( 1, 0, 0 )
2222                 : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
2223         }
2224         return packFloatx80( aSign, 0, 0 );
2225     }
2226     lastBitMask = 1;
2227     lastBitMask <<= 0x403E - aExp;
2228     roundBitsMask = lastBitMask - 1;
2229     z = a;
2230     switch (status->float_rounding_mode) {
2231     case float_round_nearest_even:
2232         z.low += lastBitMask>>1;
2233         if ((z.low & roundBitsMask) == 0) {
2234             z.low &= ~lastBitMask;
2235         }
2236         break;
2237     case float_round_ties_away:
2238         z.low += lastBitMask >> 1;
2239         break;
2240     case float_round_to_zero:
2241         break;
2242     case float_round_up:
2243         if (!extractFloatx80Sign(z)) {
2244             z.low += roundBitsMask;
2245         }
2246         break;
2247     case float_round_down:
2248         if (extractFloatx80Sign(z)) {
2249             z.low += roundBitsMask;
2250         }
2251         break;
2252     default:
2253         abort();
2254     }
2255     z.low &= ~ roundBitsMask;
2256     if ( z.low == 0 ) {
2257         ++z.high;
2258         z.low = LIT64( 0x8000000000000000 );
2259     }
2260     if (z.low != a.low) {
2261         status->float_exception_flags |= float_flag_inexact;
2262     }
2263     return z;
2264
2265 }
2266
2267 #ifdef SOFTFLOAT_68K // 09-01-2017: Added for Previous
2268 floatx80 floatx80_round_to_int_toward_zero( floatx80 a, float_status *status)
2269 {
2270     flag aSign;
2271     int32_t aExp;
2272     uint64_t lastBitMask, roundBitsMask;
2273     floatx80 z;
2274
2275         aSign = extractFloatx80Sign(a);
2276         aExp = extractFloatx80Exp( a );
2277     if ( 0x403E <= aExp ) {
2278         if ( aExp == 0x7FFF ) {
2279                         if ( (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
2280                     return propagateFloatx80NaNOneArg( a, status );
2281                         return inf_clear_intbit(status) ? packFloatx80(aSign, aExp, 0) : a;
2282         }
2283         return a;
2284     }
2285     if ( aExp < 0x3FFF ) {
2286         if (    ( aExp == 0 )
2287 #ifdef SOFTFLOAT_68K
2288             && ( (uint64_t) extractFloatx80Frac( a ) == 0 ) ) {
2289 #else
2290             && ( (uint64_t) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
2291 #endif
2292             return a;
2293         }
2294         status->float_exception_flags |= float_flag_inexact;
2295         return packFloatx80( aSign, 0, 0 );
2296     }
2297     lastBitMask = 1;
2298     lastBitMask <<= 0x403E - aExp;
2299     roundBitsMask = lastBitMask - 1;
2300     z = a;
2301     z.low &= ~ roundBitsMask;
2302     if ( z.low == 0 ) {
2303         ++z.high;
2304         z.low = LIT64( 0x8000000000000000 );
2305     }
2306     if ( z.low != a.low ) status->float_exception_flags |= float_flag_inexact;
2307     return z;
2308
2309 }
2310 #endif // End of addition for Previous
2311
2312 /*----------------------------------------------------------------------------
2313 | Returns the result of adding the absolute values of the extended double-
2314 | precision floating-point values `a' and `b'.  If `zSign' is 1, the sum is
2315 | negated before being returned.  `zSign' is ignored if the result is a NaN.
2316 | The addition is performed according to the IEC/IEEE Standard for Binary
2317 | Floating-Point Arithmetic.
2318 *----------------------------------------------------------------------------*/
2319
2320 static floatx80 addFloatx80Sigs(floatx80 a, floatx80 b, flag zSign,
2321                                 float_status *status)
2322 {
2323     int32_t aExp, bExp, zExp;
2324     uint64_t aSig, bSig, zSig0, zSig1;
2325     int32_t expDiff;
2326
2327     aSig = extractFloatx80Frac( a );
2328     aExp = extractFloatx80Exp( a );
2329     bSig = extractFloatx80Frac( b );
2330     bExp = extractFloatx80Exp( b );
2331 #ifdef SOFTFLOAT_68K
2332         if ( aExp == 0 ) {
2333                 normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2334         }
2335         if ( bExp == 0 ) {
2336                 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
2337         }
2338 #endif
2339     expDiff = aExp - bExp;
2340     if ( 0 < expDiff ) {
2341         if ( aExp == 0x7FFF ) {
2342             if ((uint64_t)(aSig << 1))
2343                 return propagateFloatx80NaN(a, b, status);
2344                         return inf_clear_intbit(status) ? packFloatx80(extractFloatx80Sign(a), aExp, 0) : a;
2345         }
2346 #ifndef SOFTFLOAT_68K
2347                 if ( bExp == 0 ) --expDiff;
2348 #endif
2349         shift64ExtraRightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
2350         zExp = aExp;
2351     }
2352     else if ( expDiff < 0 ) {
2353         if ( bExp == 0x7FFF ) {
2354             if ((uint64_t)(bSig << 1))
2355                 return propagateFloatx80NaN(a, b, status);
2356                         if (inf_clear_intbit(status)) bSig = 0;
2357             return packFloatx80( zSign, bExp, bSig );
2358         }
2359 #ifndef SOFTFLOAT_68K
2360                 if ( aExp == 0 ) ++expDiff;
2361 #endif
2362         shift64ExtraRightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
2363         zExp = bExp;
2364     }
2365     else {
2366         if ( aExp == 0x7FFF ) {
2367             if ( (uint64_t) ( ( aSig | bSig )<<1 ) ) {
2368                 return propagateFloatx80NaN(a, b, status);
2369             }
2370                         if (inf_clear_intbit(status)) return packFloatx80(extractFloatx80Sign(a), aExp, 0);
2371                         return faddsub_swap_inf(status) ? b : a;
2372         }
2373         zSig1 = 0;
2374         zSig0 = aSig + bSig;
2375  #ifndef SOFTFLOAT_68K
2376                 if ( aExp == 0 ) {
2377                         normalizeFloatx80Subnormal( zSig0, &zExp, &zSig0 );
2378                         goto roundAndPack;
2379                 }
2380 #endif
2381         zExp = aExp;
2382 #ifdef SOFTFLOAT_68K
2383                 if ( aSig == 0 && bSig == 0 ) return packFloatx80( zSign, 0, 0 );
2384         if ( aSig == 0 || bSig == 0 ) goto roundAndPack;
2385 #endif
2386                 goto shiftRight1;
2387     }
2388     zSig0 = aSig + bSig;
2389     if ( (int64_t) zSig0 < 0 ) goto roundAndPack;
2390  shiftRight1:
2391     shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 );
2392     zSig0 |= LIT64( 0x8000000000000000 );
2393     ++zExp;
2394  roundAndPack:
2395     return roundAndPackFloatx80(status->floatx80_rounding_precision,
2396                                 zSign, zExp, zSig0, zSig1, status);
2397 }
2398
2399 /*----------------------------------------------------------------------------
2400 | Returns the result of subtracting the absolute values of the extended
2401 | double-precision floating-point values `a' and `b'.  If `zSign' is 1, the
2402 | difference is negated before being returned.  `zSign' is ignored if the
2403 | result is a NaN.  The subtraction is performed according to the IEC/IEEE
2404 | Standard for Binary Floating-Point Arithmetic.
2405 *----------------------------------------------------------------------------*/
2406
2407 static floatx80 subFloatx80Sigs(floatx80 a, floatx80 b, flag zSign,
2408                                 float_status *status)
2409 {
2410     int32_t aExp, bExp, zExp;
2411     uint64_t aSig, bSig, zSig0, zSig1;
2412     int32_t expDiff;
2413
2414     aSig = extractFloatx80Frac( a );
2415     aExp = extractFloatx80Exp( a );
2416     bSig = extractFloatx80Frac( b );
2417     bExp = extractFloatx80Exp( b );
2418     expDiff = aExp - bExp;
2419     if ( 0 < expDiff ) goto aExpBigger;
2420     if ( expDiff < 0 ) goto bExpBigger;
2421     if ( aExp == 0x7FFF ) {
2422         if ( (uint64_t) ( ( aSig | bSig )<<1 ) ) {
2423             return propagateFloatx80NaN(a, b, status);
2424         }
2425         float_raise(float_flag_invalid, status);
2426         return floatx80_default_nan(status);
2427     }
2428  #ifndef SOFTFLOAT_68K
2429         if ( aExp == 0 ) {
2430                 aExp = 1;
2431                 bExp = 1;
2432         }
2433 #endif
2434     zSig1 = 0;
2435     if ( bSig < aSig ) goto aBigger;
2436     if ( aSig < bSig ) goto bBigger;
2437     return packFloatx80(status->float_rounding_mode == float_round_down, 0, 0);
2438  bExpBigger:
2439     if ( bExp == 0x7FFF ) {
2440                 if ((uint64_t)(bSig << 1)) return propagateFloatx80NaN(a, b, status);
2441                 if (inf_clear_intbit(status)) bSig = 0;
2442                 return packFloatx80(zSign ^ 1, bExp, bSig);
2443         }
2444 #ifndef SOFTFLOAT_68K
2445         if ( aExp == 0 ) ++expDiff;
2446 #endif
2447     shift128RightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
2448  bBigger:
2449     sub128( bSig, 0, aSig, zSig1, &zSig0, &zSig1 );
2450     zExp = bExp;
2451     zSign ^= 1;
2452     goto normalizeRoundAndPack;
2453  aExpBigger:
2454     if ( aExp == 0x7FFF ) {
2455                 if ((uint64_t)(aSig << 1)) return propagateFloatx80NaN(a, b, status);
2456                 return inf_clear_intbit(status) ? packFloatx80(extractFloatx80Sign(a), aExp, 0) : a;
2457         }
2458 #ifndef SOFTFLOAT_68K
2459         if ( bExp == 0 ) --expDiff;
2460 #endif
2461     shift128RightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
2462  aBigger:
2463     sub128( aSig, 0, bSig, zSig1, &zSig0, &zSig1 );
2464     zExp = aExp;
2465  normalizeRoundAndPack:
2466     return normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision,
2467                                          zSign, zExp, zSig0, zSig1, status);
2468 }
2469
2470 /*----------------------------------------------------------------------------
2471 | Returns the result of adding the extended double-precision floating-point
2472 | values `a' and `b'.  The operation is performed according to the IEC/IEEE
2473 | Standard for Binary Floating-Point Arithmetic.
2474 *----------------------------------------------------------------------------*/
2475
2476 floatx80 floatx80_add(floatx80 a, floatx80 b, float_status *status)
2477 {
2478     flag aSign, bSign;
2479
2480     if (floatx80_invalid_encoding(a) || floatx80_invalid_encoding(b)) {
2481         float_raise(float_flag_invalid, status);
2482         return floatx80_default_nan(status);
2483     }
2484     aSign = extractFloatx80Sign( a );
2485     bSign = extractFloatx80Sign( b );
2486     if ( aSign == bSign ) {
2487         return addFloatx80Sigs(a, b, aSign, status);
2488     }
2489     else {
2490         return subFloatx80Sigs(a, b, aSign, status);
2491     }
2492
2493 }
2494
2495 /*----------------------------------------------------------------------------
2496 | Returns the result of subtracting the extended double-precision floating-
2497 | point values `a' and `b'.  The operation is performed according to the
2498 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2499 *----------------------------------------------------------------------------*/
2500
2501 floatx80 floatx80_sub(floatx80 a, floatx80 b, float_status *status)
2502 {
2503     flag aSign, bSign;
2504
2505     if (floatx80_invalid_encoding(a) || floatx80_invalid_encoding(b)) {
2506         float_raise(float_flag_invalid, status);
2507         return floatx80_default_nan(status);
2508     }
2509     aSign = extractFloatx80Sign( a );
2510     bSign = extractFloatx80Sign( b );
2511     if ( aSign == bSign ) {
2512         return subFloatx80Sigs(a, b, aSign, status);
2513     }
2514     else {
2515         return addFloatx80Sigs(a, b, aSign, status);
2516     }
2517
2518 }
2519
2520 /*----------------------------------------------------------------------------
2521 | Returns the result of multiplying the extended double-precision floating-
2522 | point values `a' and `b'.  The operation is performed according to the
2523 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2524 *----------------------------------------------------------------------------*/
2525
2526 floatx80 floatx80_mul(floatx80 a, floatx80 b, float_status *status)
2527 {
2528     flag aSign, bSign, zSign;
2529     int32_t aExp, bExp, zExp;
2530     uint64_t aSig, bSig, zSig0, zSig1;
2531
2532     if (floatx80_invalid_encoding(a) || floatx80_invalid_encoding(b)) {
2533         float_raise(float_flag_invalid, status);
2534         return floatx80_default_nan(status);
2535     }
2536     aSig = extractFloatx80Frac( a );
2537     aExp = extractFloatx80Exp( a );
2538     aSign = extractFloatx80Sign( a );
2539     bSig = extractFloatx80Frac( b );
2540     bExp = extractFloatx80Exp( b );
2541     bSign = extractFloatx80Sign( b );
2542     zSign = aSign ^ bSign;
2543     if ( aExp == 0x7FFF ) {
2544         if (    (uint64_t) ( aSig<<1 )
2545              || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) {
2546             return propagateFloatx80NaN(a, b, status);
2547         }
2548         if ( ( bExp | bSig ) == 0 ) goto invalid;
2549                 if (inf_clear_intbit(status)) aSig = 0;
2550                 return packFloatx80(zSign, aExp, aSig);
2551         }
2552     if ( bExp == 0x7FFF ) {
2553         if ((uint64_t)(bSig << 1)) {
2554             return propagateFloatx80NaN(a, b, status);
2555         }
2556         if ( ( aExp | aSig ) == 0 ) {
2557  invalid:
2558             float_raise(float_flag_invalid, status);
2559             return floatx80_default_nan(status);
2560         }
2561                 if (inf_clear_intbit(status)) bSig = 0;
2562                 return packFloatx80(zSign, bExp, bSig);
2563         }
2564     if ( aExp == 0 ) {
2565         if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
2566         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2567     }
2568     if ( bExp == 0 ) {
2569         if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 );
2570         normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
2571     }
2572     zExp = aExp + bExp - 0x3FFE;
2573     mul64To128( aSig, bSig, &zSig0, &zSig1 );
2574     if ( 0 < (int64_t) zSig0 ) {
2575         shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
2576         --zExp;
2577     }
2578     return roundAndPackFloatx80(status->floatx80_rounding_precision,
2579                                 zSign, zExp, zSig0, zSig1, status);
2580 }
2581
2582 #ifdef SOFTFLOAT_68K // 21-01-2017: Added for Previous
2583 floatx80 floatx80_sglmul( floatx80 a, floatx80 b, float_status *status )
2584 {
2585         flag aSign, bSign, zSign;
2586         int32_t aExp, bExp, zExp;
2587         uint64_t aSig, bSig, zSig0, zSig1;
2588
2589         aSig = extractFloatx80Frac( a );
2590         aExp = extractFloatx80Exp( a );
2591         aSign = extractFloatx80Sign( a );
2592         bSig = extractFloatx80Frac( b );
2593         bExp = extractFloatx80Exp( b );
2594         bSign = extractFloatx80Sign( b );
2595         zSign = aSign ^ bSign;
2596         if ( aExp == 0x7FFF ) {
2597                 if (    (uint64_t) ( aSig<<1 )
2598                         || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) {
2599                         return propagateFloatx80NaN( a, b, status );
2600                 }
2601                 if ( ( bExp | bSig ) == 0 ) goto invalid;
2602                 if (inf_clear_intbit(status)) aSig = 0;
2603                 return packFloatx80(zSign, aExp, aSig);
2604         }
2605         if ( bExp == 0x7FFF ) {
2606                 if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b, status );
2607                 if ( ( aExp | aSig ) == 0 ) {
2608                 invalid:
2609                         float_raise( float_flag_invalid, status );
2610                         return floatx80_default_nan(status);
2611                 }
2612                 if (inf_clear_intbit(status)) bSig = 0;
2613                 return packFloatx80(zSign, bExp, bSig);
2614         }
2615         if ( aExp == 0 ) {
2616                 if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
2617                 normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2618         }
2619         if ( bExp == 0 ) {
2620                 if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 );
2621                 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
2622         }
2623         aSig &= LIT64( 0xFFFFFF0000000000 );
2624         bSig &= LIT64( 0xFFFFFF0000000000 );
2625         zExp = aExp + bExp - 0x3FFE;
2626         mul64To128( aSig, bSig, &zSig0, &zSig1 );
2627         if ( 0 < (int64_t) zSig0 ) {
2628                 shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
2629                 --zExp;
2630         }
2631         return roundSigAndPackFloatx80( 32, zSign, zExp, zSig0, zSig1, status);
2632
2633 }
2634 #endif // End of addition for Previous
2635
2636
2637 /*----------------------------------------------------------------------------
2638 | Returns the result of dividing the extended double-precision floating-point
2639 | value `a' by the corresponding value `b'.  The operation is performed
2640 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2641 *----------------------------------------------------------------------------*/
2642
2643 floatx80 floatx80_div(floatx80 a, floatx80 b, float_status *status)
2644 {
2645     flag aSign, bSign, zSign;
2646     int32_t aExp, bExp, zExp;
2647     uint64_t aSig, bSig, zSig0, zSig1;
2648     uint64_t rem0, rem1, rem2, term0, term1, term2;
2649
2650     if (floatx80_invalid_encoding(a) || floatx80_invalid_encoding(b)) {
2651         float_raise(float_flag_invalid, status);
2652         return floatx80_default_nan(status);
2653     }
2654     aSig = extractFloatx80Frac( a );
2655     aExp = extractFloatx80Exp( a );
2656     aSign = extractFloatx80Sign( a );
2657     bSig = extractFloatx80Frac( b );
2658     bExp = extractFloatx80Exp( b );
2659     bSign = extractFloatx80Sign( b );
2660     zSign = aSign ^ bSign;
2661     if ( aExp == 0x7FFF ) {
2662         if ((uint64_t)(aSig << 1)) {
2663             return propagateFloatx80NaN(a, b, status);
2664         }
2665         if ( bExp == 0x7FFF ) {
2666             if ((uint64_t)(bSig << 1)) {
2667                 return propagateFloatx80NaN(a, b, status);
2668             }
2669             goto invalid;
2670         }
2671                 if (inf_clear_intbit(status)) aSig = 0;
2672                 return packFloatx80(zSign, aExp, aSig);
2673         }
2674     if ( bExp == 0x7FFF ) {
2675         if ((uint64_t)(bSig << 1)) {
2676             return propagateFloatx80NaN(a, b, status);
2677         }
2678         return packFloatx80( zSign, 0, 0 );
2679     }
2680     if ( bExp == 0 ) {
2681         if ( bSig == 0 ) {
2682             if ( ( aExp | aSig ) == 0 ) {
2683  invalid:
2684                 float_raise(float_flag_invalid, status);
2685                 return floatx80_default_nan(status);
2686             }
2687             float_raise(float_flag_divbyzero, status);
2688             return packFloatx80( zSign, 0x7FFF, floatx80_default_infinity_low );
2689         }
2690         normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
2691     }
2692     if ( aExp == 0 ) {
2693         if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
2694         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2695     }
2696     zExp = aExp - bExp + 0x3FFE;
2697     rem1 = 0;
2698     if ( bSig <= aSig ) {
2699         shift128Right( aSig, 0, 1, &aSig, &rem1 );
2700         ++zExp;
2701     }
2702     zSig0 = estimateDiv128To64( aSig, rem1, bSig );
2703     mul64To128( bSig, zSig0, &term0, &term1 );
2704     sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
2705     while ( (int64_t) rem0 < 0 ) {
2706         --zSig0;
2707         add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
2708     }
2709     zSig1 = estimateDiv128To64( rem1, 0, bSig );
2710     if ( (uint64_t) ( zSig1<<1 ) <= 8 ) {
2711         mul64To128( bSig, zSig1, &term1, &term2 );
2712         sub128( rem1, 0, term1, term2, &rem1, &rem2 );
2713         while ( (int64_t) rem1 < 0 ) {
2714             --zSig1;
2715             add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
2716         }
2717         zSig1 |= ( ( rem1 | rem2 ) != 0 );
2718     }
2719     return roundAndPackFloatx80(status->floatx80_rounding_precision,
2720                                 zSign, zExp, zSig0, zSig1, status);
2721 }
2722
2723 #ifdef SOFTFLOAT_68K // 21-01-2017: Addition for Previous
2724 floatx80 floatx80_sgldiv( floatx80 a, floatx80 b, float_status *status )
2725 {
2726         flag aSign, bSign, zSign;
2727         int32_t aExp, bExp, zExp;
2728         uint64_t aSig, bSig, zSig0, zSig1;
2729         uint64_t rem0, rem1, rem2, term0, term1, term2;
2730
2731         aSig = extractFloatx80Frac( a );
2732         aExp = extractFloatx80Exp( a );
2733         aSign = extractFloatx80Sign( a );
2734         bSig = extractFloatx80Frac( b );
2735         bExp = extractFloatx80Exp( b );
2736         bSign = extractFloatx80Sign( b );
2737         zSign = aSign ^ bSign;
2738         if ( aExp == 0x7FFF ) {
2739                 if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b, status );
2740                 if ( bExp == 0x7FFF ) {
2741                         if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b, status );
2742                         goto invalid;
2743                 }
2744                 if (inf_clear_intbit(status)) aSig = 0;
2745                 return packFloatx80(zSign, aExp, aSig);
2746         }
2747         if ( bExp == 0x7FFF ) {
2748                 if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b, status );
2749                 return packFloatx80( zSign, 0, 0 );
2750         }
2751         if ( bExp == 0 ) {
2752                 if ( bSig == 0 ) {
2753                         if ( ( aExp | aSig ) == 0 ) {
2754                         invalid:
2755                                 float_raise( float_flag_invalid, status );
2756                                 return floatx80_default_nan(status);
2757                         }
2758                         float_raise( float_flag_divbyzero, status );
2759                         return packFloatx80( zSign, 0x7FFF, floatx80_default_infinity_low );
2760                 }
2761                 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
2762         }
2763         if ( aExp == 0 ) {
2764                 if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
2765                 normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
2766         }
2767
2768         zExp = aExp - bExp + 0x3FFE;
2769         rem1 = 0;
2770         if ( bSig <= aSig ) {
2771                 shift128Right( aSig, 0, 1, &aSig, &rem1 );
2772                 ++zExp;
2773         }
2774         zSig0 = estimateDiv128To64( aSig, rem1, bSig );
2775         mul64To128( bSig, zSig0, &term0, &term1 );
2776         sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
2777         while ( (int64_t) rem0 < 0 ) {
2778                 --zSig0;
2779                 add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
2780         }
2781         zSig1 = estimateDiv128To64( rem1, 0, bSig );
2782         if ( (uint64_t) ( zSig1<<1 ) <= 8 ) {
2783                 mul64To128( bSig, zSig1, &term1, &term2 );
2784                 sub128( rem1, 0, term1, term2, &rem1, &rem2 );
2785                 while ( (int64_t) rem1 < 0 ) {
2786                         --zSig1;
2787                         add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
2788                 }
2789                 zSig1 |= ( ( rem1 | rem2 ) != 0 );
2790         }
2791         return roundSigAndPackFloatx80( 32, zSign, zExp, zSig0, zSig1, status);
2792
2793 }
2794 #endif // End of addition for Previous
2795
2796
2797 /*----------------------------------------------------------------------------
2798 | Returns the remainder of the extended double-precision floating-point value
2799 | `a' with respect to the corresponding value `b'.  The operation is performed
2800 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2801 *----------------------------------------------------------------------------*/
2802
2803 #ifndef SOFTFLOAT_68K
2804 floatx80 floatx80_rem(floatx80 a, floatx80 b, float_status *status)
2805 {
2806     flag aSign, zSign;
2807     int32_t aExp, bExp, expDiff;
2808     uint64_t aSig0, aSig1, bSig;
2809     uint64_t q, term0, term1, alternateASig0, alternateASig1;
2810
2811     if (floatx80_invalid_encoding(a) || floatx80_invalid_encoding(b)) {
2812         float_raise(float_flag_invalid, status);
2813         return floatx80_default_nan(status);
2814     }
2815     aSig0 = extractFloatx80Frac( a );
2816     aExp = extractFloatx80Exp( a );
2817     aSign = extractFloatx80Sign( a );
2818     bSig = extractFloatx80Frac( b );
2819     bExp = extractFloatx80Exp( b );
2820     if ( aExp == 0x7FFF ) {
2821         if (    (uint64_t) ( aSig0<<1 )
2822              || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) {
2823             return propagateFloatx80NaN(a, b, status);
2824         }
2825         goto invalid;
2826     }
2827     if ( bExp == 0x7FFF ) {
2828         if ((uint64_t)(bSig << 1)) {
2829             return propagateFloatx80NaN(a, b, status);
2830         }
2831         return a;
2832     }
2833     if ( bExp == 0 ) {
2834         if ( bSig == 0 ) {
2835  invalid:
2836             float_raise(float_flag_invalid, status);
2837             return floatx80_default_nan(status);
2838         }
2839         normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
2840     }
2841     if ( aExp == 0 ) {
2842         if ( (uint64_t) ( aSig0<<1 ) == 0 ) return a;
2843         normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
2844     }
2845     bSig |= LIT64( 0x8000000000000000 );
2846     zSign = aSign;
2847     expDiff = aExp - bExp;
2848     aSig1 = 0;
2849     if ( expDiff < 0 ) {
2850         if ( expDiff < -1 ) return a;
2851         shift128Right( aSig0, 0, 1, &aSig0, &aSig1 );
2852         expDiff = 0;
2853     }
2854     q = ( bSig <= aSig0 );
2855     if ( q ) aSig0 -= bSig;
2856     expDiff -= 64;
2857     while ( 0 < expDiff ) {
2858         q = estimateDiv128To64( aSig0, aSig1, bSig );
2859         q = ( 2 < q ) ? q - 2 : 0;
2860         mul64To128( bSig, q, &term0, &term1 );
2861         sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
2862         shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
2863         expDiff -= 62;
2864     }
2865     expDiff += 64;
2866     if ( 0 < expDiff ) {
2867         q = estimateDiv128To64( aSig0, aSig1, bSig );
2868         q = ( 2 < q ) ? q - 2 : 0;
2869         q >>= 64 - expDiff;
2870         mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 );
2871         sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
2872         shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
2873         while ( le128( term0, term1, aSig0, aSig1 ) ) {
2874             ++q;
2875             sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
2876         }
2877     }
2878     else {
2879         term1 = 0;
2880         term0 = bSig;
2881     }
2882     sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 );
2883     if (    lt128( alternateASig0, alternateASig1, aSig0, aSig1 )
2884          || (    eq128( alternateASig0, alternateASig1, aSig0, aSig1 )
2885               && ( q & 1 ) )
2886        ) {
2887         aSig0 = alternateASig0;
2888         aSig1 = alternateASig1;
2889         zSign = ! zSign;
2890     }
2891     return
2892         normalizeRoundAndPackFloatx80(
2893             80, zSign, bExp + expDiff, aSig0, aSig1, status);
2894
2895 }
2896 #else // 09-01-2017: Modified version for Previous
2897 floatx80 floatx80_rem( floatx80 a, floatx80 b, uint64_t *q, flag *s, float_status *status )
2898 {
2899     flag aSign, bSign, zSign;
2900     int32_t aExp, bExp, expDiff;
2901     uint64_t aSig0, aSig1, bSig;
2902     uint64_t qTemp, term0, term1, alternateASig0, alternateASig1;
2903
2904     aSig0 = extractFloatx80Frac( a );
2905     aExp = extractFloatx80Exp( a );
2906     aSign = extractFloatx80Sign( a );
2907     bSig = extractFloatx80Frac( b );
2908     bExp = extractFloatx80Exp( b );
2909     bSign = extractFloatx80Sign( b );
2910
2911         if ( aExp == 0x7FFF ) {
2912         if (    (uint64_t) ( aSig0<<1 )
2913             || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) {
2914             return propagateFloatx80NaN( a, b, status );
2915         }
2916         goto invalid;
2917     }
2918     if ( bExp == 0x7FFF ) {
2919         if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b, status );
2920         *s = (aSign != bSign);
2921         *q = 0;
2922         return a;
2923     }
2924     if ( bExp == 0 ) {
2925         if ( bSig == 0 ) {
2926         invalid:
2927             float_raise( float_flag_invalid, status );
2928                         return floatx80_default_nan(status);
2929         }
2930         normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
2931     }
2932     if ( aExp == 0 ) {
2933 #ifdef SOFTFLOAT_68K
2934         if ( aSig0 == 0 ) {
2935             *s = (aSign != bSign);
2936             *q = 0;
2937             return a;
2938         }
2939 #else
2940         if ( (uint64_t) ( aSig0<<1 ) == 0 ) return a;
2941 #endif
2942         normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
2943     }
2944     bSig |= LIT64( 0x8000000000000000 );
2945     zSign = aSign;
2946     expDiff = aExp - bExp;
2947     *s = (aSign != bSign);
2948     aSig1 = 0;
2949     if ( expDiff < 0 ) {
2950         if ( expDiff < -1 ) return a;
2951         shift128Right( aSig0, 0, 1, &aSig0, &aSig1 );
2952         expDiff = 0;
2953     }
2954     qTemp = ( bSig <= aSig0 );
2955     if ( qTemp ) aSig0 -= bSig;
2956     *q = ( expDiff > 63 ) ? 0 : ( qTemp<<expDiff );
2957     expDiff -= 64;
2958     while ( 0 < expDiff ) {
2959         qTemp = estimateDiv128To64( aSig0, aSig1, bSig );
2960         qTemp = ( 2 < qTemp ) ? qTemp - 2 : 0;
2961         mul64To128( bSig, qTemp, &term0, &term1 );
2962         sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
2963         shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
2964         *q = ( expDiff > 63 ) ? 0 : ( qTemp<<expDiff );
2965         expDiff -= 62;
2966     }
2967     expDiff += 64;
2968     if ( 0 < expDiff ) {
2969         qTemp = estimateDiv128To64( aSig0, aSig1, bSig );
2970         qTemp = ( 2 < qTemp ) ? qTemp - 2 : 0;
2971         qTemp >>= 64 - expDiff;
2972         mul64To128( bSig, qTemp<<( 64 - expDiff ), &term0, &term1 );
2973         sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
2974         shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
2975         while ( le128( term0, term1, aSig0, aSig1 ) ) {
2976             ++qTemp;
2977             sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
2978         }
2979         *q += qTemp;
2980     }
2981     else {
2982         term1 = 0;
2983         term0 = bSig;
2984     }
2985     sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 );
2986     if (    lt128( alternateASig0, alternateASig1, aSig0, aSig1 )
2987         || (    eq128( alternateASig0, alternateASig1, aSig0, aSig1 )
2988             && ( qTemp & 1 ) )
2989         ) {
2990         aSig0 = alternateASig0;
2991         aSig1 = alternateASig1;
2992         zSign = ! zSign;
2993         ++*q;
2994     }
2995     return
2996     normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision,
2997                                   zSign, bExp + expDiff, aSig0, aSig1, status );
2998
2999 }
3000 #endif // End of modification
3001
3002
3003 #ifdef SOFTFLOAT_68K // 08-01-2017: Added for Previous
3004 /*----------------------------------------------------------------------------
3005  | Returns the modulo remainder of the extended double-precision floating-point
3006  | value `a' with respect to the corresponding value `b'.
3007  *----------------------------------------------------------------------------*/
3008
3009 floatx80 floatx80_mod( floatx80 a, floatx80 b, uint64_t *q, flag *s, float_status *status )
3010 {
3011     flag aSign, bSign, zSign;
3012     int32_t aExp, bExp, expDiff;
3013     uint64_t aSig0, aSig1, bSig;
3014     uint64_t qTemp, term0, term1;
3015
3016     aSig0 = extractFloatx80Frac( a );
3017     aExp = extractFloatx80Exp( a );
3018     aSign = extractFloatx80Sign( a );
3019     bSig = extractFloatx80Frac( b );
3020     bExp = extractFloatx80Exp( b );
3021     bSign = extractFloatx80Sign( b );
3022
3023         if ( aExp == 0x7FFF ) {
3024         if (    (uint64_t) ( aSig0<<1 )
3025             || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) {
3026             return propagateFloatx80NaN( a, b, status );
3027         }
3028         goto invalid;
3029     }
3030     if ( bExp == 0x7FFF ) {
3031         if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b, status );
3032         *s = (aSign != bSign);
3033         *q = 0;
3034         return a;
3035     }
3036     if ( bExp == 0 ) {
3037         if ( bSig == 0 ) {
3038         invalid:
3039             float_raise( float_flag_invalid, status );
3040                         return floatx80_default_nan(status);
3041                 }
3042         normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3043     }
3044     if ( aExp == 0 ) {
3045 #ifdef SOFTFLOAT_68K
3046         if ( aSig0 == 0 ) {
3047             *s = (aSign != bSign);
3048             *q = 0;
3049             return a;
3050         }
3051 #else
3052         if ( (uint64_t) ( aSig0<<1 ) == 0 ) return a;
3053 #endif
3054         normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
3055     }
3056     bSig |= LIT64( 0x8000000000000000 );
3057     zSign = aSign;
3058     expDiff = aExp - bExp;
3059     *s = (aSign != bSign);
3060     aSig1 = 0;
3061     if ( expDiff < 0 ) return a;
3062     qTemp = ( bSig <= aSig0 );
3063     if ( qTemp ) aSig0 -= bSig;
3064     *q = ( expDiff > 63 ) ? 0 : ( qTemp<<expDiff );
3065     expDiff -= 64;
3066     while ( 0 < expDiff ) {
3067         qTemp = estimateDiv128To64( aSig0, aSig1, bSig );
3068         qTemp = ( 2 < qTemp ) ? qTemp - 2 : 0;
3069         mul64To128( bSig, qTemp, &term0, &term1 );
3070         sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3071         shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
3072         *q = ( expDiff > 63 ) ? 0 : ( qTemp<<expDiff );
3073                 expDiff -= 62;
3074     }
3075     expDiff += 64;
3076     if ( 0 < expDiff ) {
3077         qTemp = estimateDiv128To64( aSig0, aSig1, bSig );
3078         qTemp = ( 2 < qTemp ) ? qTemp - 2 : 0;
3079         qTemp >>= 64 - expDiff;
3080         mul64To128( bSig, qTemp<<( 64 - expDiff ), &term0, &term1 );
3081         sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3082         shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
3083         while ( le128( term0, term1, aSig0, aSig1 ) ) {
3084             ++qTemp;
3085             sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3086         }
3087         *q += qTemp;
3088     }
3089     return
3090         normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision,
3091             zSign, bExp + expDiff, aSig0, aSig1, status );
3092
3093 }
3094 #endif // end of addition for Previous
3095
3096
3097 /*----------------------------------------------------------------------------
3098 | Returns the square root of the extended double-precision floating-point
3099 | value `a'.  The operation is performed according to the IEC/IEEE Standard
3100 | for Binary Floating-Point Arithmetic.
3101 *----------------------------------------------------------------------------*/
3102
3103 floatx80 floatx80_sqrt(floatx80 a, float_status *status)
3104 {
3105     flag aSign;
3106     int32_t aExp, zExp;
3107     uint64_t aSig0, aSig1, zSig0, zSig1, doubleZSig0;
3108     uint64_t rem0, rem1, rem2, rem3, term0, term1, term2, term3;
3109
3110     if (floatx80_invalid_encoding(a)) {
3111         float_raise(float_flag_invalid, status);
3112         return floatx80_default_nan(status);
3113     }
3114     aSig0 = extractFloatx80Frac( a );
3115     aExp = extractFloatx80Exp( a );
3116     aSign = extractFloatx80Sign( a );
3117     if ( aExp == 0x7FFF ) {
3118         if ((uint64_t)(aSig0 << 1))
3119             return propagateFloatx80NaNOneArg(a, status);
3120                 if (!aSign) return inf_clear_intbit(status) ? packFloatx80(aSign, aExp, 0) : a;
3121         goto invalid;
3122     }
3123     if ( aSign ) {
3124         if ( ( aExp | aSig0 ) == 0 ) return a;
3125  invalid:
3126         float_raise(float_flag_invalid, status);
3127         return floatx80_default_nan(status);
3128     }
3129     if ( aExp == 0 ) {
3130         if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 );
3131         normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
3132     }
3133     zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF;
3134     zSig0 = estimateSqrt32( aExp, aSig0>>32 );
3135     shift128Right( aSig0, 0, 2 + ( aExp & 1 ), &aSig0, &aSig1 );
3136     zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
3137     doubleZSig0 = zSig0<<1;
3138     mul64To128( zSig0, zSig0, &term0, &term1 );
3139     sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
3140     while ( (int64_t) rem0 < 0 ) {
3141         --zSig0;
3142         doubleZSig0 -= 2;
3143         add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
3144     }
3145     zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
3146     if ( ( zSig1 & LIT64( 0x3FFFFFFFFFFFFFFF ) ) <= 5 ) {
3147         if ( zSig1 == 0 ) zSig1 = 1;
3148         mul64To128( doubleZSig0, zSig1, &term1, &term2 );
3149         sub128( rem1, 0, term1, term2, &rem1, &rem2 );
3150         mul64To128( zSig1, zSig1, &term2, &term3 );
3151         sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
3152         while ( (int64_t) rem1 < 0 ) {
3153             --zSig1;
3154             shortShift128Left( 0, zSig1, 1, &term2, &term3 );
3155             term3 |= 1;
3156             term2 |= doubleZSig0;
3157             add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
3158         }
3159         zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
3160     }
3161     shortShift128Left( 0, zSig1, 1, &zSig0, &zSig1 );
3162     zSig0 |= doubleZSig0;
3163     return roundAndPackFloatx80(status->floatx80_rounding_precision,
3164                                 0, zExp, zSig0, zSig1, status);
3165 }
3166
3167
3168 #ifdef SOFTFLOAT_68K // 07-01-2017: Added for Previous
3169 /*----------------------------------------------------------------------------
3170  | Returns the mantissa of the extended double-precision floating-point
3171  | value `a'.
3172  *----------------------------------------------------------------------------*/
3173
3174 floatx80 floatx80_getman( floatx80 a, float_status *status)
3175 {
3176     flag aSign;
3177     int32_t aExp;
3178     uint64_t aSig;
3179
3180     aSig = extractFloatx80Frac( a );
3181     aExp = extractFloatx80Exp( a );
3182     aSign = extractFloatx80Sign( a );
3183
3184     if ( aExp == 0x7FFF ) {
3185         if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaNOneArg( a, status );
3186         float_raise( float_flag_invalid, status );
3187                 return floatx80_default_nan(status);
3188     }
3189
3190     if ( aExp == 0 ) {
3191         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
3192         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3193     }
3194
3195     return packFloatx80(aSign, 0x3fff, aSig);
3196 }
3197
3198 /*----------------------------------------------------------------------------
3199  | Returns the exponent of the extended double-precision floating-point
3200  | value `a' as an extended double-precision value.
3201  *----------------------------------------------------------------------------*/
3202
3203 floatx80 floatx80_getexp( floatx80 a, float_status *status)
3204 {
3205     flag aSign;
3206     int32_t aExp;
3207     uint64_t aSig;
3208
3209     aSig = extractFloatx80Frac( a );
3210     aExp = extractFloatx80Exp( a );
3211     aSign = extractFloatx80Sign( a );
3212
3213     if ( aExp == 0x7FFF ) {
3214         if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaNOneArg( a, status );
3215         float_raise( float_flag_invalid, status );
3216                 return floatx80_default_nan(status);
3217     }
3218
3219     if ( aExp == 0 ) {
3220         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
3221         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3222     }
3223
3224     return int32_to_floatx80(aExp - 0x3FFF);
3225 }
3226
3227 /*----------------------------------------------------------------------------
3228  | Scales extended double-precision floating-point value in operand `a' by
3229  | value `b'. The function truncates the value in the second operand 'b' to
3230  | an integral value and adds that value to the exponent of the operand 'a'.
3231  | The operation performed according to the IEC/IEEE Standard for Binary
3232  | Floating-Point Arithmetic.
3233  *----------------------------------------------------------------------------*/
3234
3235 floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status)
3236 {
3237     flag aSign, bSign;
3238     int32_t aExp, bExp, shiftCount;
3239     uint64_t aSig, bSig;
3240
3241     aSig = extractFloatx80Frac(a);
3242     aExp = extractFloatx80Exp(a);
3243     aSign = extractFloatx80Sign(a);
3244     bSig = extractFloatx80Frac(b);
3245     bExp = extractFloatx80Exp(b);
3246     bSign = extractFloatx80Sign(b);
3247
3248     if ( bExp == 0x7FFF ) {
3249         if ( (uint64_t) ( bSig<<1 ) ||
3250             ( ( aExp == 0x7FFF ) && (uint64_t) ( aSig<<1 ) ) ) {
3251             return propagateFloatx80NaN( a, b, status );
3252         }
3253         float_raise( float_flag_invalid, status );
3254                 return floatx80_default_nan(status);
3255     }
3256     if ( aExp == 0x7FFF ) {
3257         if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b, status );
3258         return a;
3259     }
3260     if ( aExp == 0 ) {
3261         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0);
3262         if ( bExp < 0x3FFF ) return a;
3263         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3264     }
3265
3266     if (bExp < 0x3FFF) {
3267         return roundAndPackFloatx80(
3268             status->floatx80_rounding_precision, aSign, aExp, aSig, 0, status);
3269     }
3270
3271     if ( 0x400F < bExp ) {
3272         aExp = bSign ? -0x6001 : 0xE000;
3273         return roundAndPackFloatx80(
3274                     status->floatx80_rounding_precision, aSign, aExp, aSig, 0, status );
3275     }
3276
3277     shiftCount = 0x403E - bExp;
3278     bSig >>= shiftCount;
3279     aExp = bSign ? ( aExp - bSig ) : ( aExp + bSig );
3280
3281     return roundAndPackFloatx80(
3282                 status->floatx80_rounding_precision, aSign, aExp, aSig, 0, status);
3283
3284 }
3285
3286 /*-----------------------------------------------------------------------------
3287  | Calculates the absolute value of the extended double-precision floating-point
3288  | value `a'.  The operation is performed according to the IEC/IEEE Standard
3289  | for Binary Floating-Point Arithmetic.
3290  *----------------------------------------------------------------------------*/
3291
3292 floatx80 floatx80_abs(floatx80 a, float_status *status)
3293 {
3294     int32_t aExp;
3295     uint64_t aSig;
3296
3297     aSig = extractFloatx80Frac(a);
3298     aExp = extractFloatx80Exp(a);
3299
3300     if ( aExp == 0x7FFF ) {
3301         if ( (uint64_t) ( aSig<<1 ) )
3302             return propagateFloatx80NaNOneArg( a, status );
3303                 if (inf_clear_intbit(status)) aSig = 0;
3304                 return packFloatx80(0, aExp, aSig);
3305         }
3306
3307     if ( aExp == 0 ) {
3308         if ( aSig == 0 ) return packFloatx80( 0, 0, 0 );
3309         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3310     }
3311
3312     return roundAndPackFloatx80(
3313                 status->floatx80_rounding_precision, 0, aExp, aSig, 0, status );
3314
3315 }
3316
3317 /*-----------------------------------------------------------------------------
3318  | Changes the sign of the extended double-precision floating-point value 'a'.
3319  | The operation is performed according to the IEC/IEEE Standard for Binary
3320  | Floating-Point Arithmetic.
3321  *----------------------------------------------------------------------------*/
3322
3323 floatx80 floatx80_neg(floatx80 a, float_status *status)
3324 {
3325     flag aSign;
3326     int32_t aExp;
3327     uint64_t aSig;
3328
3329     aSig = extractFloatx80Frac(a);
3330     aExp = extractFloatx80Exp(a);
3331     aSign = extractFloatx80Sign(a);
3332
3333     if ( aExp == 0x7FFF ) {
3334         if ( (uint64_t) ( aSig<<1 ) )
3335             return propagateFloatx80NaNOneArg( a, status );
3336                 if (inf_clear_intbit(status)) aSig = 0;
3337                 return packFloatx80(!aSign, aExp, aSig);
3338         }
3339
3340         aSign = !aSign;
3341
3342         if ( aExp == 0 ) {
3343         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
3344         normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3345     }
3346
3347     return roundAndPackFloatx80(
3348                 status->floatx80_rounding_precision, aSign, aExp, aSig, 0, status );
3349
3350 }
3351
3352 /*----------------------------------------------------------------------------
3353  | Returns the result of comparing the extended double-precision floating-
3354  | point values `a' and `b'.  The result is abstracted for matching the
3355  | corresponding condition codes.
3356  *----------------------------------------------------------------------------*/
3357
3358 floatx80 floatx80_cmp( floatx80 a, floatx80 b, float_status *status )
3359 {
3360     flag aSign, bSign;
3361     int32_t aExp, bExp;
3362     uint64_t aSig, bSig;
3363
3364     aSig = extractFloatx80Frac( a );
3365     aExp = extractFloatx80Exp( a );
3366     aSign = extractFloatx80Sign( a );
3367     bSig = extractFloatx80Frac( b );
3368     bExp = extractFloatx80Exp( b );
3369     bSign = extractFloatx80Sign( b );
3370
3371     if ( ( aExp == 0x7FFF && (uint64_t) ( aSig<<1 ) ) ||
3372          ( bExp == 0x7FFF && (uint64_t) ( bSig<<1 ) ) ) {
3373                 // 68040 FCMP -NaN return N flag set
3374                 if (fcmp_signed_nan(status))
3375                 return propagateFloatx80NaN(a, b, status );
3376                 return propagateFloatx80NaN(packFloatx80(0, aExp, aSig),
3377                         packFloatx80(0, bExp, bSig), status);
3378         }
3379
3380     if ( bExp < aExp ) return packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
3381     if ( aExp < bExp ) return packFloatx80( bSign ^ 1, 0x3FFF, LIT64( 0x8000000000000000 ) );
3382
3383     if ( aExp == 0x7FFF ) {
3384         if ( aSign == bSign ) return packFloatx80( aSign, 0, 0 );
3385         return packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
3386     }
3387
3388     if ( bSig < aSig ) return packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
3389     if ( aSig < bSig ) return packFloatx80( bSign ^ 1, 0x3FFF, LIT64( 0x8000000000000000 ) );
3390
3391     if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
3392
3393     if ( aSign == bSign ) return packFloatx80( 0, 0, 0 );
3394
3395     return packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
3396
3397 }
3398
3399 floatx80 floatx80_tst( floatx80 a, float_status *status )
3400 {
3401     int32_t aExp;
3402     uint64_t aSig;
3403
3404     aSig = extractFloatx80Frac( a );
3405     aExp = extractFloatx80Exp( a );
3406
3407     if ( aExp == 0x7FFF && (uint64_t) ( aSig<<1 ) )
3408         return propagateFloatx80NaNOneArg( a, status );
3409     return a;
3410 }
3411
3412 floatx80 floatx80_move( floatx80 a, float_status *status )
3413 {
3414     flag aSign;
3415     int32_t aExp;
3416     uint64_t aSig;
3417
3418     aSig = extractFloatx80Frac( a );
3419     aExp = extractFloatx80Exp( a );
3420     aSign = extractFloatx80Sign( a );
3421
3422     if ( aExp == 0x7FFF ) {
3423                 if ((uint64_t)(aSig << 1)) return propagateFloatx80NaNOneArg(a, status);
3424                 return inf_clear_intbit(status) ? packFloatx80(aSign, aExp, 0) : a;
3425     }
3426     if ( aExp == 0 ) {
3427         if ( aSig == 0 ) return a;
3428         return normalizeRoundAndPackFloatx80( status->floatx80_rounding_precision, aSign, aExp, aSig, 0, status );
3429     }
3430     return roundAndPackFloatx80( status->floatx80_rounding_precision, aSign, aExp, aSig, 0, status );
3431 }
3432
3433 floatx80 floatx80_denormalize( floatx80 a, flag eSign)
3434 {
3435         flag aSign;
3436         int32_t aExp;
3437         uint64_t aSig;
3438         int32_t shiftCount;
3439
3440         aSig = extractFloatx80Frac( a );
3441         aExp = extractFloatx80Exp( a );
3442         aSign = extractFloatx80Sign( a );
3443
3444         if ( eSign ) {
3445                 shiftCount = 0x8000 - aExp;
3446                 aExp = 0;
3447                 if (shiftCount > 63) {
3448                         aSig = 0;
3449                 } else {
3450                         aSig >>= shiftCount;
3451                 }
3452         }
3453         return packFloatx80(aSign, aExp, aSig);
3454 }
3455
3456 #endif // End of addition for Previous
3457
3458 /*----------------------------------------------------------------------------
3459 | Returns 1 if the extended double-precision floating-point value `a' is
3460 | equal to the corresponding value `b', and 0 otherwise.  The comparison is
3461 | performed according to the IEC/IEEE Standard for Binary Floating-Point
3462 | Arithmetic.
3463 *----------------------------------------------------------------------------*/
3464
3465 flag floatx80_eq( floatx80 a, floatx80 b, float_status *status )
3466 {
3467         if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
3468                                 && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
3469                         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
3470                                 && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
3471                 ) {
3472                 if (    floatx80_is_signaling_nan( a )
3473                                 || floatx80_is_signaling_nan( b ) ) {
3474                         float_raise( float_flag_invalid, status );
3475                 }
3476                 return 0;
3477         }
3478         return
3479                         ( a.low == b.low )
3480                 && (    ( a.high == b.high )
3481                                 || (    ( a.low == 0 )
3482                                         && ( (uint16_t) ( ( a.high | b.high )<<1 ) == 0 ) )
3483                         );
3484
3485 }
3486
3487 /*----------------------------------------------------------------------------
3488 | Returns 1 if the extended double-precision floating-point value `a' is
3489 | less than or equal to the corresponding value `b', and 0 otherwise.  The
3490 | comparison is performed according to the IEC/IEEE Standard for Binary
3491 | Floating-Point Arithmetic.
3492 *----------------------------------------------------------------------------*/
3493
3494 flag floatx80_le( floatx80 a, floatx80 b, float_status *status )
3495 {
3496         flag aSign, bSign;
3497
3498         if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
3499                                 && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
3500                         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
3501                                 && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
3502                 ) {
3503                 float_raise( float_flag_invalid, status );
3504                 return 0;
3505         }
3506         aSign = extractFloatx80Sign( a );
3507         bSign = extractFloatx80Sign( b );
3508         if ( aSign != bSign ) {
3509                 return
3510                                 aSign
3511                         || (    ( ( (uint16_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
3512                                         == 0 );
3513         }
3514         return
3515                         aSign ? le128( b.high, b.low, a.high, a.low )
3516                 : le128( a.high, a.low, b.high, b.low );
3517 }
3518
3519 /*----------------------------------------------------------------------------
3520 | Returns 1 if the extended double-precision floating-point value `a' is
3521 | less than the corresponding value `b', and 0 otherwise.  The comparison
3522 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
3523 | Arithmetic.
3524 *----------------------------------------------------------------------------*/
3525
3526 flag floatx80_lt( floatx80 a, floatx80 b, float_status *status )
3527 {
3528         flag aSign, bSign;
3529
3530         if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
3531                                 && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
3532                         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
3533                                 && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
3534                 ) {
3535                 float_raise( float_flag_invalid, status );
3536                 return 0;
3537         }
3538         aSign = extractFloatx80Sign( a );
3539         bSign = extractFloatx80Sign( b );
3540         if ( aSign != bSign ) {
3541                 return
3542                                 aSign
3543                         && (    ( ( (uint16_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
3544                                         != 0 );
3545         }
3546         return
3547                         aSign ? lt128( b.high, b.low, a.high, a.low )
3548                 : lt128( a.high, a.low, b.high, b.low );
3549
3550 }
3551
3552
3553 /*----------------------------------------------------------------------------
3554 | Returns the result of converting the 64-bit two's complement integer `a'
3555 | to the extended double-precision floating-point format.  The conversion
3556 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
3557 | Arithmetic.
3558 *----------------------------------------------------------------------------*/
3559
3560 floatx80 int64_to_floatx80( int64_t a )
3561 {
3562         flag zSign;
3563         uint64_t absA;
3564         int8_t shiftCount;
3565
3566         if ( a == 0 ) return packFloatx80( 0, 0, 0 );
3567         zSign = ( a < 0 );
3568         absA = zSign ? - a : a;
3569         shiftCount = countLeadingZeros64( absA );
3570         return packFloatx80( zSign, 0x403E - shiftCount, absA<<shiftCount );
3571 }