1129203Scognet/* $NetBSD: softfloat.c,v 1.1 2002/05/21 23:51:07 bjh21 Exp $ */ 2129203Scognet 3129203Scognet/* 4129203Scognet * This version hacked for use with gcc -msoft-float by bjh21. 5129203Scognet * (Mostly a case of #ifdefing out things GCC doesn't need or provides 6129203Scognet * itself). 7129203Scognet */ 8129203Scognet 9129203Scognet/* 10129203Scognet * Things you may want to define: 11129203Scognet * 12129203Scognet * SOFTFLOAT_FOR_GCC - build only those functions necessary for GCC (with 13129203Scognet * -msoft-float) to work. Include "softfloat-for-gcc.h" to get them 14129203Scognet * properly renamed. 15129203Scognet */ 16129203Scognet 17129203Scognet/* 18129203Scognet * This differs from the standard bits32/softfloat.c in that float64 19129203Scognet * is defined to be a 64-bit integer rather than a structure. The 20129203Scognet * structure is float64s, with translation between the two going via 21129203Scognet * float64u. 22129203Scognet */ 23129203Scognet 24129203Scognet/* 25129203Scognet=============================================================================== 26129203Scognet 27129203ScognetThis C source file is part of the SoftFloat IEC/IEEE Floating-Point 28129203ScognetArithmetic Package, Release 2a. 29129203Scognet 30129203ScognetWritten by John R. Hauser. This work was made possible in part by the 31129203ScognetInternational Computer Science Institute, located at Suite 600, 1947 Center 32129203ScognetStreet, Berkeley, California 94704. Funding was partially provided by the 33129203ScognetNational Science Foundation under grant MIP-9311980. The original version 34129203Scognetof this code was written as part of a project to build a fixed-point vector 35129203Scognetprocessor in collaboration with the University of California at Berkeley, 36129203Scognetoverseen by Profs. Nelson Morgan and John Wawrzynek. More information 37129203Scognetis available through the Web page `http://HTTP.CS.Berkeley.EDU/~jhauser/ 38129203Scognetarithmetic/SoftFloat.html'. 39129203Scognet 40129203ScognetTHIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort 41129203Scognethas been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT 42129203ScognetTIMES RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO 43129203ScognetPERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY 44129203ScognetAND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE. 45129203Scognet 46129203ScognetDerivative works are acceptable, even for commercial purposes, so long as 47129203Scognet(1) they include prominent notice that the work is derivative, and (2) they 48129203Scognetinclude prominent notice akin to these four paragraphs for those parts of 49129203Scognetthis code that are retained. 50129203Scognet 51129203Scognet=============================================================================== 52129203Scognet*/ 53129203Scognet 54129203Scognet#include <sys/cdefs.h> 55129203Scognet__FBSDID("$FreeBSD$"); 56129203Scognet 57129203Scognet#ifdef SOFTFLOAT_FOR_GCC 58129203Scognet#include "softfloat-for-gcc.h" 59129203Scognet#endif 60129203Scognet 61129203Scognet#include "milieu.h" 62129203Scognet#include "softfloat.h" 63129203Scognet 64129203Scognet/* 65129203Scognet * Conversions between floats as stored in memory and floats as 66129203Scognet * SoftFloat uses them 67129203Scognet */ 68129203Scognet#ifndef FLOAT64_DEMANGLE 69129203Scognet#define FLOAT64_DEMANGLE(a) (a) 70129203Scognet#endif 71129203Scognet#ifndef FLOAT64_MANGLE 72129203Scognet#define FLOAT64_MANGLE(a) (a) 73129203Scognet#endif 74129203Scognet 75129203Scognet/* 76129203Scognet------------------------------------------------------------------------------- 77129203ScognetFloating-point rounding mode and exception flags. 78129203Scognet------------------------------------------------------------------------------- 79129203Scognet*/ 80230189Sdasint float_rounding_mode = float_round_nearest_even; 81230189Sdasint float_exception_flags = 0; 82129203Scognet 83129203Scognet/* 84129203Scognet------------------------------------------------------------------------------- 85129203ScognetPrimitive arithmetic functions, including multi-word arithmetic, and 86129203Scognetdivision and square root approximations. (Can be specialized to target if 87129203Scognetdesired.) 88129203Scognet------------------------------------------------------------------------------- 89129203Scognet*/ 90129203Scognet#include "softfloat-macros" 91129203Scognet 92129203Scognet/* 93129203Scognet------------------------------------------------------------------------------- 94129203ScognetFunctions and definitions to determine: (1) whether tininess for underflow 95129203Scognetis detected before or after rounding by default, (2) what (if anything) 96129203Scognethappens when exceptions are raised, (3) how signaling NaNs are distinguished 97129203Scognetfrom quiet NaNs, (4) the default generated quiet NaNs, and (4) how NaNs 98129203Scognetare propagated from function inputs to output. These details are target- 99129203Scognetspecific. 100129203Scognet------------------------------------------------------------------------------- 101129203Scognet*/ 102129203Scognet#include "softfloat-specialize" 103129203Scognet 104129203Scognet/* 105129203Scognet------------------------------------------------------------------------------- 106129203ScognetReturns the fraction bits of the single-precision floating-point value `a'. 107129203Scognet------------------------------------------------------------------------------- 108129203Scognet*/ 109129203ScognetINLINE bits32 extractFloat32Frac( float32 a ) 110129203Scognet{ 111129203Scognet 112129203Scognet return a & 0x007FFFFF; 113129203Scognet 114129203Scognet} 115129203Scognet 116129203Scognet/* 117129203Scognet------------------------------------------------------------------------------- 118129203ScognetReturns the exponent bits of the single-precision floating-point value `a'. 119129203Scognet------------------------------------------------------------------------------- 120129203Scognet*/ 121129203ScognetINLINE int16 extractFloat32Exp( float32 a ) 122129203Scognet{ 123129203Scognet 124129203Scognet return ( a>>23 ) & 0xFF; 125129203Scognet 126129203Scognet} 127129203Scognet 128129203Scognet/* 129129203Scognet------------------------------------------------------------------------------- 130129203ScognetReturns the sign bit of the single-precision floating-point value `a'. 131129203Scognet------------------------------------------------------------------------------- 132129203Scognet*/ 133129203ScognetINLINE flag extractFloat32Sign( float32 a ) 134129203Scognet{ 135129203Scognet 136129203Scognet return a>>31; 137129203Scognet 138129203Scognet} 139129203Scognet 140129203Scognet/* 141129203Scognet------------------------------------------------------------------------------- 142129203ScognetNormalizes the subnormal single-precision floating-point value represented 143129203Scognetby the denormalized significand `aSig'. The normalized exponent and 144129203Scognetsignificand are stored at the locations pointed to by `zExpPtr' and 145129203Scognet`zSigPtr', respectively. 146129203Scognet------------------------------------------------------------------------------- 147129203Scognet*/ 148129203Scognetstatic void 149129203Scognet normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr ) 150129203Scognet{ 151129203Scognet int8 shiftCount; 152129203Scognet 153129203Scognet shiftCount = countLeadingZeros32( aSig ) - 8; 154129203Scognet *zSigPtr = aSig<<shiftCount; 155129203Scognet *zExpPtr = 1 - shiftCount; 156129203Scognet 157129203Scognet} 158129203Scognet 159129203Scognet/* 160129203Scognet------------------------------------------------------------------------------- 161129203ScognetPacks the sign `zSign', exponent `zExp', and significand `zSig' into a 162129203Scognetsingle-precision floating-point value, returning the result. After being 163129203Scognetshifted into the proper positions, the three fields are simply added 164129203Scognettogether to form the result. This means that any integer portion of `zSig' 165129203Scognetwill be added into the exponent. Since a properly normalized significand 166129203Scognetwill have an integer portion equal to 1, the `zExp' input should be 1 less 167129203Scognetthan the desired result exponent whenever `zSig' is a complete, normalized 168129203Scognetsignificand. 169129203Scognet------------------------------------------------------------------------------- 170129203Scognet*/ 171129203ScognetINLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig ) 172129203Scognet{ 173129203Scognet 174129203Scognet return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig; 175129203Scognet 176129203Scognet} 177129203Scognet 178129203Scognet/* 179129203Scognet------------------------------------------------------------------------------- 180129203ScognetTakes an abstract floating-point value having sign `zSign', exponent `zExp', 181129203Scognetand significand `zSig', and returns the proper single-precision floating- 182129203Scognetpoint value corresponding to the abstract input. Ordinarily, the abstract 183129203Scognetvalue is simply rounded and packed into the single-precision format, with 184129203Scognetthe inexact exception raised if the abstract input cannot be represented 185129203Scognetexactly. However, if the abstract value is too large, the overflow and 186129203Scognetinexact exceptions are raised and an infinity or maximal finite value is 187129203Scognetreturned. If the abstract value is too small, the input value is rounded to 188129203Scogneta subnormal number, and the underflow and inexact exceptions are raised if 189129203Scognetthe abstract input cannot be represented exactly as a subnormal single- 190129203Scognetprecision floating-point number. 191129203Scognet The input significand `zSig' has its binary point between bits 30 192129203Scognetand 29, which is 7 bits to the left of the usual location. This shifted 193129203Scognetsignificand must be normalized or smaller. If `zSig' is not normalized, 194129203Scognet`zExp' must be 0; in that case, the result returned is a subnormal number, 195129203Scognetand it must not require rounding. In the usual case that `zSig' is 196129203Scognetnormalized, `zExp' must be 1 less than the ``true'' floating-point exponent. 197129203ScognetThe handling of underflow and overflow follows the IEC/IEEE Standard for 198129203ScognetBinary Floating-Point Arithmetic. 199129203Scognet------------------------------------------------------------------------------- 200129203Scognet*/ 201129203Scognetstatic float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig ) 202129203Scognet{ 203129203Scognet int8 roundingMode; 204129203Scognet flag roundNearestEven; 205129203Scognet int8 roundIncrement, roundBits; 206129203Scognet flag isTiny; 207129203Scognet 208129203Scognet roundingMode = float_rounding_mode; 209129203Scognet roundNearestEven = roundingMode == float_round_nearest_even; 210129203Scognet roundIncrement = 0x40; 211129203Scognet if ( ! roundNearestEven ) { 212129203Scognet if ( roundingMode == float_round_to_zero ) { 213129203Scognet roundIncrement = 0; 214129203Scognet } 215129203Scognet else { 216129203Scognet roundIncrement = 0x7F; 217129203Scognet if ( zSign ) { 218129203Scognet if ( roundingMode == float_round_up ) roundIncrement = 0; 219129203Scognet } 220129203Scognet else { 221129203Scognet if ( roundingMode == float_round_down ) roundIncrement = 0; 222129203Scognet } 223129203Scognet } 224129203Scognet } 225129203Scognet roundBits = zSig & 0x7F; 226129203Scognet if ( 0xFD <= (bits16) zExp ) { 227129203Scognet if ( ( 0xFD < zExp ) 228129203Scognet || ( ( zExp == 0xFD ) 229129203Scognet && ( (sbits32) ( zSig + roundIncrement ) < 0 ) ) 230129203Scognet ) { 231129203Scognet float_raise( float_flag_overflow | float_flag_inexact ); 232129203Scognet return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 ); 233129203Scognet } 234129203Scognet if ( zExp < 0 ) { 235129203Scognet isTiny = 236129203Scognet ( float_detect_tininess == float_tininess_before_rounding ) 237129203Scognet || ( zExp < -1 ) 238129203Scognet || ( zSig + roundIncrement < 0x80000000 ); 239129203Scognet shift32RightJamming( zSig, - zExp, &zSig ); 240129203Scognet zExp = 0; 241129203Scognet roundBits = zSig & 0x7F; 242129203Scognet if ( isTiny && roundBits ) float_raise( float_flag_underflow ); 243129203Scognet } 244129203Scognet } 245129203Scognet if ( roundBits ) float_exception_flags |= float_flag_inexact; 246129203Scognet zSig = ( zSig + roundIncrement )>>7; 247129203Scognet zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven ); 248129203Scognet if ( zSig == 0 ) zExp = 0; 249129203Scognet return packFloat32( zSign, zExp, zSig ); 250129203Scognet 251129203Scognet} 252129203Scognet 253129203Scognet/* 254129203Scognet------------------------------------------------------------------------------- 255129203ScognetTakes an abstract floating-point value having sign `zSign', exponent `zExp', 256129203Scognetand significand `zSig', and returns the proper single-precision floating- 257129203Scognetpoint value corresponding to the abstract input. This routine is just like 258129203Scognet`roundAndPackFloat32' except that `zSig' does not have to be normalized. 259129203ScognetBit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true'' 260129203Scognetfloating-point exponent. 261129203Scognet------------------------------------------------------------------------------- 262129203Scognet*/ 263129203Scognetstatic float32 264129203Scognet normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig ) 265129203Scognet{ 266129203Scognet int8 shiftCount; 267129203Scognet 268129203Scognet shiftCount = countLeadingZeros32( zSig ) - 1; 269129203Scognet return roundAndPackFloat32( zSign, zExp - shiftCount, zSig<<shiftCount ); 270129203Scognet 271129203Scognet} 272129203Scognet 273129203Scognet/* 274129203Scognet------------------------------------------------------------------------------- 275129203ScognetReturns the least-significant 32 fraction bits of the double-precision 276129203Scognetfloating-point value `a'. 277129203Scognet------------------------------------------------------------------------------- 278129203Scognet*/ 279129203ScognetINLINE bits32 extractFloat64Frac1( float64 a ) 280129203Scognet{ 281129203Scognet 282129203Scognet return FLOAT64_DEMANGLE(a) & LIT64( 0x00000000FFFFFFFF ); 283129203Scognet 284129203Scognet} 285129203Scognet 286129203Scognet/* 287129203Scognet------------------------------------------------------------------------------- 288129203ScognetReturns the most-significant 20 fraction bits of the double-precision 289129203Scognetfloating-point value `a'. 290129203Scognet------------------------------------------------------------------------------- 291129203Scognet*/ 292129203ScognetINLINE bits32 extractFloat64Frac0( float64 a ) 293129203Scognet{ 294129203Scognet 295129203Scognet return ( FLOAT64_DEMANGLE(a)>>32 ) & 0x000FFFFF; 296129203Scognet 297129203Scognet} 298129203Scognet 299129203Scognet/* 300129203Scognet------------------------------------------------------------------------------- 301129203ScognetReturns the exponent bits of the double-precision floating-point value `a'. 302129203Scognet------------------------------------------------------------------------------- 303129203Scognet*/ 304129203ScognetINLINE int16 extractFloat64Exp( float64 a ) 305129203Scognet{ 306129203Scognet 307129203Scognet return ( FLOAT64_DEMANGLE(a)>>52 ) & 0x7FF; 308129203Scognet 309129203Scognet} 310129203Scognet 311129203Scognet/* 312129203Scognet------------------------------------------------------------------------------- 313129203ScognetReturns the sign bit of the double-precision floating-point value `a'. 314129203Scognet------------------------------------------------------------------------------- 315129203Scognet*/ 316129203ScognetINLINE flag extractFloat64Sign( float64 a ) 317129203Scognet{ 318129203Scognet 319129203Scognet return FLOAT64_DEMANGLE(a)>>63; 320129203Scognet 321129203Scognet} 322129203Scognet 323129203Scognet/* 324129203Scognet------------------------------------------------------------------------------- 325129203ScognetNormalizes the subnormal double-precision floating-point value represented 326129203Scognetby the denormalized significand formed by the concatenation of `aSig0' and 327129203Scognet`aSig1'. The normalized exponent is stored at the location pointed to by 328129203Scognet`zExpPtr'. The most significant 21 bits of the normalized significand are 329129203Scognetstored at the location pointed to by `zSig0Ptr', and the least significant 330129203Scognet32 bits of the normalized significand are stored at the location pointed to 331129203Scognetby `zSig1Ptr'. 332129203Scognet------------------------------------------------------------------------------- 333129203Scognet*/ 334129203Scognetstatic void 335129203Scognet normalizeFloat64Subnormal( 336129203Scognet bits32 aSig0, 337129203Scognet bits32 aSig1, 338129203Scognet int16 *zExpPtr, 339129203Scognet bits32 *zSig0Ptr, 340129203Scognet bits32 *zSig1Ptr 341129203Scognet ) 342129203Scognet{ 343129203Scognet int8 shiftCount; 344129203Scognet 345129203Scognet if ( aSig0 == 0 ) { 346129203Scognet shiftCount = countLeadingZeros32( aSig1 ) - 11; 347129203Scognet if ( shiftCount < 0 ) { 348129203Scognet *zSig0Ptr = aSig1>>( - shiftCount ); 349129203Scognet *zSig1Ptr = aSig1<<( shiftCount & 31 ); 350129203Scognet } 351129203Scognet else { 352129203Scognet *zSig0Ptr = aSig1<<shiftCount; 353129203Scognet *zSig1Ptr = 0; 354129203Scognet } 355129203Scognet *zExpPtr = - shiftCount - 31; 356129203Scognet } 357129203Scognet else { 358129203Scognet shiftCount = countLeadingZeros32( aSig0 ) - 11; 359129203Scognet shortShift64Left( aSig0, aSig1, shiftCount, zSig0Ptr, zSig1Ptr ); 360129203Scognet *zExpPtr = 1 - shiftCount; 361129203Scognet } 362129203Scognet 363129203Scognet} 364129203Scognet 365129203Scognet/* 366129203Scognet------------------------------------------------------------------------------- 367129203ScognetPacks the sign `zSign', the exponent `zExp', and the significand formed by 368129203Scognetthe concatenation of `zSig0' and `zSig1' into a double-precision floating- 369129203Scognetpoint value, returning the result. After being shifted into the proper 370129203Scognetpositions, the three fields `zSign', `zExp', and `zSig0' are simply added 371129203Scognettogether to form the most significant 32 bits of the result. This means 372129203Scognetthat any integer portion of `zSig0' will be added into the exponent. Since 373129203Scogneta properly normalized significand will have an integer portion equal to 1, 374129203Scognetthe `zExp' input should be 1 less than the desired result exponent whenever 375129203Scognet`zSig0' and `zSig1' concatenated form a complete, normalized significand. 376129203Scognet------------------------------------------------------------------------------- 377129203Scognet*/ 378129203ScognetINLINE float64 379129203Scognet packFloat64( flag zSign, int16 zExp, bits32 zSig0, bits32 zSig1 ) 380129203Scognet{ 381129203Scognet 382129203Scognet return FLOAT64_MANGLE( ( ( (bits64) zSign )<<63 ) + 383129203Scognet ( ( (bits64) zExp )<<52 ) + 384129203Scognet ( ( (bits64) zSig0 )<<32 ) + zSig1 ); 385129203Scognet 386129203Scognet 387129203Scognet} 388129203Scognet 389129203Scognet/* 390129203Scognet------------------------------------------------------------------------------- 391129203ScognetTakes an abstract floating-point value having sign `zSign', exponent `zExp', 392129203Scognetand extended significand formed by the concatenation of `zSig0', `zSig1', 393129203Scognetand `zSig2', and returns the proper double-precision floating-point value 394129203Scognetcorresponding to the abstract input. Ordinarily, the abstract value is 395129203Scognetsimply rounded and packed into the double-precision format, with the inexact 396129203Scognetexception raised if the abstract input cannot be represented exactly. 397129203ScognetHowever, if the abstract value is too large, the overflow and inexact 398129203Scognetexceptions are raised and an infinity or maximal finite value is returned. 399129203ScognetIf the abstract value is too small, the input value is rounded to a 400129203Scognetsubnormal number, and the underflow and inexact exceptions are raised if the 401129203Scognetabstract input cannot be represented exactly as a subnormal double-precision 402129203Scognetfloating-point number. 403129203Scognet The input significand must be normalized or smaller. If the input 404129203Scognetsignificand is not normalized, `zExp' must be 0; in that case, the result 405129203Scognetreturned is a subnormal number, and it must not require rounding. In the 406129203Scognetusual case that the input significand is normalized, `zExp' must be 1 less 407129203Scognetthan the ``true'' floating-point exponent. The handling of underflow and 408129203Scognetoverflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 409129203Scognet------------------------------------------------------------------------------- 410129203Scognet*/ 411129203Scognetstatic float64 412129203Scognet roundAndPackFloat64( 413129203Scognet flag zSign, int16 zExp, bits32 zSig0, bits32 zSig1, bits32 zSig2 ) 414129203Scognet{ 415129203Scognet int8 roundingMode; 416129203Scognet flag roundNearestEven, increment, isTiny; 417129203Scognet 418129203Scognet roundingMode = float_rounding_mode; 419129203Scognet roundNearestEven = ( roundingMode == float_round_nearest_even ); 420129203Scognet increment = ( (sbits32) zSig2 < 0 ); 421129203Scognet if ( ! roundNearestEven ) { 422129203Scognet if ( roundingMode == float_round_to_zero ) { 423129203Scognet increment = 0; 424129203Scognet } 425129203Scognet else { 426129203Scognet if ( zSign ) { 427129203Scognet increment = ( roundingMode == float_round_down ) && zSig2; 428129203Scognet } 429129203Scognet else { 430129203Scognet increment = ( roundingMode == float_round_up ) && zSig2; 431129203Scognet } 432129203Scognet } 433129203Scognet } 434129203Scognet if ( 0x7FD <= (bits16) zExp ) { 435129203Scognet if ( ( 0x7FD < zExp ) 436129203Scognet || ( ( zExp == 0x7FD ) 437129203Scognet && eq64( 0x001FFFFF, 0xFFFFFFFF, zSig0, zSig1 ) 438129203Scognet && increment 439129203Scognet ) 440129203Scognet ) { 441129203Scognet float_raise( float_flag_overflow | float_flag_inexact ); 442129203Scognet if ( ( roundingMode == float_round_to_zero ) 443129203Scognet || ( zSign && ( roundingMode == float_round_up ) ) 444129203Scognet || ( ! zSign && ( roundingMode == float_round_down ) ) 445129203Scognet ) { 446129203Scognet return packFloat64( zSign, 0x7FE, 0x000FFFFF, 0xFFFFFFFF ); 447129203Scognet } 448129203Scognet return packFloat64( zSign, 0x7FF, 0, 0 ); 449129203Scognet } 450129203Scognet if ( zExp < 0 ) { 451129203Scognet isTiny = 452129203Scognet ( float_detect_tininess == float_tininess_before_rounding ) 453129203Scognet || ( zExp < -1 ) 454129203Scognet || ! increment 455129203Scognet || lt64( zSig0, zSig1, 0x001FFFFF, 0xFFFFFFFF ); 456129203Scognet shift64ExtraRightJamming( 457129203Scognet zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 ); 458129203Scognet zExp = 0; 459129203Scognet if ( isTiny && zSig2 ) float_raise( float_flag_underflow ); 460129203Scognet if ( roundNearestEven ) { 461129203Scognet increment = ( (sbits32) zSig2 < 0 ); 462129203Scognet } 463129203Scognet else { 464129203Scognet if ( zSign ) { 465129203Scognet increment = ( roundingMode == float_round_down ) && zSig2; 466129203Scognet } 467129203Scognet else { 468129203Scognet increment = ( roundingMode == float_round_up ) && zSig2; 469129203Scognet } 470129203Scognet } 471129203Scognet } 472129203Scognet } 473129203Scognet if ( zSig2 ) float_exception_flags |= float_flag_inexact; 474129203Scognet if ( increment ) { 475129203Scognet add64( zSig0, zSig1, 0, 1, &zSig0, &zSig1 ); 476129203Scognet zSig1 &= ~ ( ( zSig2 + zSig2 == 0 ) & roundNearestEven ); 477129203Scognet } 478129203Scognet else { 479129203Scognet if ( ( zSig0 | zSig1 ) == 0 ) zExp = 0; 480129203Scognet } 481129203Scognet return packFloat64( zSign, zExp, zSig0, zSig1 ); 482129203Scognet 483129203Scognet} 484129203Scognet 485129203Scognet/* 486129203Scognet------------------------------------------------------------------------------- 487129203ScognetTakes an abstract floating-point value having sign `zSign', exponent `zExp', 488129203Scognetand significand formed by the concatenation of `zSig0' and `zSig1', and 489129203Scognetreturns the proper double-precision floating-point value corresponding 490129203Scognetto the abstract input. This routine is just like `roundAndPackFloat64' 491129203Scognetexcept that the input significand has fewer bits and does not have to be 492129203Scognetnormalized. In all cases, `zExp' must be 1 less than the ``true'' floating- 493129203Scognetpoint exponent. 494129203Scognet------------------------------------------------------------------------------- 495129203Scognet*/ 496129203Scognetstatic float64 497129203Scognet normalizeRoundAndPackFloat64( 498129203Scognet flag zSign, int16 zExp, bits32 zSig0, bits32 zSig1 ) 499129203Scognet{ 500129203Scognet int8 shiftCount; 501129203Scognet bits32 zSig2; 502129203Scognet 503129203Scognet if ( zSig0 == 0 ) { 504129203Scognet zSig0 = zSig1; 505129203Scognet zSig1 = 0; 506129203Scognet zExp -= 32; 507129203Scognet } 508129203Scognet shiftCount = countLeadingZeros32( zSig0 ) - 11; 509129203Scognet if ( 0 <= shiftCount ) { 510129203Scognet zSig2 = 0; 511129203Scognet shortShift64Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 ); 512129203Scognet } 513129203Scognet else { 514129203Scognet shift64ExtraRightJamming( 515129203Scognet zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 ); 516129203Scognet } 517129203Scognet zExp -= shiftCount; 518129203Scognet return roundAndPackFloat64( zSign, zExp, zSig0, zSig1, zSig2 ); 519129203Scognet 520129203Scognet} 521129203Scognet 522129203Scognet/* 523129203Scognet------------------------------------------------------------------------------- 524129203ScognetReturns the result of converting the 32-bit two's complement integer `a' to 525129203Scognetthe single-precision floating-point format. The conversion is performed 526129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 527129203Scognet------------------------------------------------------------------------------- 528129203Scognet*/ 529129203Scognetfloat32 int32_to_float32( int32 a ) 530129203Scognet{ 531129203Scognet flag zSign; 532129203Scognet 533129203Scognet if ( a == 0 ) return 0; 534129203Scognet if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 ); 535129203Scognet zSign = ( a < 0 ); 536129203Scognet return normalizeRoundAndPackFloat32( zSign, 0x9C, zSign ? - a : a ); 537129203Scognet 538129203Scognet} 539129203Scognet 540129203Scognet/* 541129203Scognet------------------------------------------------------------------------------- 542129203ScognetReturns the result of converting the 32-bit two's complement integer `a' to 543129203Scognetthe double-precision floating-point format. The conversion is performed 544129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 545129203Scognet------------------------------------------------------------------------------- 546129203Scognet*/ 547129203Scognetfloat64 int32_to_float64( int32 a ) 548129203Scognet{ 549129203Scognet flag zSign; 550129203Scognet bits32 absA; 551129203Scognet int8 shiftCount; 552129203Scognet bits32 zSig0, zSig1; 553129203Scognet 554129203Scognet if ( a == 0 ) return packFloat64( 0, 0, 0, 0 ); 555129203Scognet zSign = ( a < 0 ); 556129203Scognet absA = zSign ? - a : a; 557129203Scognet shiftCount = countLeadingZeros32( absA ) - 11; 558129203Scognet if ( 0 <= shiftCount ) { 559129203Scognet zSig0 = absA<<shiftCount; 560129203Scognet zSig1 = 0; 561129203Scognet } 562129203Scognet else { 563129203Scognet shift64Right( absA, 0, - shiftCount, &zSig0, &zSig1 ); 564129203Scognet } 565129203Scognet return packFloat64( zSign, 0x412 - shiftCount, zSig0, zSig1 ); 566129203Scognet 567129203Scognet} 568129203Scognet 569129203Scognet#ifndef SOFTFLOAT_FOR_GCC 570129203Scognet/* 571129203Scognet------------------------------------------------------------------------------- 572129203ScognetReturns the result of converting the single-precision floating-point value 573129203Scognet`a' to the 32-bit two's complement integer format. The conversion is 574129203Scognetperformed according to the IEC/IEEE Standard for Binary Floating-Point 575129203ScognetArithmetic---which means in particular that the conversion is rounded 576129203Scognetaccording to the current rounding mode. If `a' is a NaN, the largest 577129203Scognetpositive integer is returned. Otherwise, if the conversion overflows, the 578129203Scognetlargest integer with the same sign as `a' is returned. 579129203Scognet------------------------------------------------------------------------------- 580129203Scognet*/ 581129203Scognetint32 float32_to_int32( float32 a ) 582129203Scognet{ 583129203Scognet flag aSign; 584129203Scognet int16 aExp, shiftCount; 585129203Scognet bits32 aSig, aSigExtra; 586129203Scognet int32 z; 587129203Scognet int8 roundingMode; 588129203Scognet 589129203Scognet aSig = extractFloat32Frac( a ); 590129203Scognet aExp = extractFloat32Exp( a ); 591129203Scognet aSign = extractFloat32Sign( a ); 592129203Scognet shiftCount = aExp - 0x96; 593129203Scognet if ( 0 <= shiftCount ) { 594129203Scognet if ( 0x9E <= aExp ) { 595129203Scognet if ( a != 0xCF000000 ) { 596129203Scognet float_raise( float_flag_invalid ); 597129203Scognet if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) { 598129203Scognet return 0x7FFFFFFF; 599129203Scognet } 600129203Scognet } 601129203Scognet return (sbits32) 0x80000000; 602129203Scognet } 603129203Scognet z = ( aSig | 0x00800000 )<<shiftCount; 604129203Scognet if ( aSign ) z = - z; 605129203Scognet } 606129203Scognet else { 607129203Scognet if ( aExp < 0x7E ) { 608129203Scognet aSigExtra = aExp | aSig; 609129203Scognet z = 0; 610129203Scognet } 611129203Scognet else { 612129203Scognet aSig |= 0x00800000; 613129203Scognet aSigExtra = aSig<<( shiftCount & 31 ); 614129203Scognet z = aSig>>( - shiftCount ); 615129203Scognet } 616129203Scognet if ( aSigExtra ) float_exception_flags |= float_flag_inexact; 617129203Scognet roundingMode = float_rounding_mode; 618129203Scognet if ( roundingMode == float_round_nearest_even ) { 619129203Scognet if ( (sbits32) aSigExtra < 0 ) { 620129203Scognet ++z; 621129203Scognet if ( (bits32) ( aSigExtra<<1 ) == 0 ) z &= ~1; 622129203Scognet } 623129203Scognet if ( aSign ) z = - z; 624129203Scognet } 625129203Scognet else { 626129203Scognet aSigExtra = ( aSigExtra != 0 ); 627129203Scognet if ( aSign ) { 628129203Scognet z += ( roundingMode == float_round_down ) & aSigExtra; 629129203Scognet z = - z; 630129203Scognet } 631129203Scognet else { 632129203Scognet z += ( roundingMode == float_round_up ) & aSigExtra; 633129203Scognet } 634129203Scognet } 635129203Scognet } 636129203Scognet return z; 637129203Scognet 638129203Scognet} 639129203Scognet#endif 640129203Scognet 641129203Scognet/* 642129203Scognet------------------------------------------------------------------------------- 643129203ScognetReturns the result of converting the single-precision floating-point value 644129203Scognet`a' to the 32-bit two's complement integer format. The conversion is 645129203Scognetperformed according to the IEC/IEEE Standard for Binary Floating-Point 646129203ScognetArithmetic, except that the conversion is always rounded toward zero. 647129203ScognetIf `a' is a NaN, the largest positive integer is returned. Otherwise, if 648129203Scognetthe conversion overflows, the largest integer with the same sign as `a' is 649129203Scognetreturned. 650129203Scognet------------------------------------------------------------------------------- 651129203Scognet*/ 652129203Scognetint32 float32_to_int32_round_to_zero( float32 a ) 653129203Scognet{ 654129203Scognet flag aSign; 655129203Scognet int16 aExp, shiftCount; 656129203Scognet bits32 aSig; 657129203Scognet int32 z; 658129203Scognet 659129203Scognet aSig = extractFloat32Frac( a ); 660129203Scognet aExp = extractFloat32Exp( a ); 661129203Scognet aSign = extractFloat32Sign( a ); 662129203Scognet shiftCount = aExp - 0x9E; 663129203Scognet if ( 0 <= shiftCount ) { 664129203Scognet if ( a != 0xCF000000 ) { 665129203Scognet float_raise( float_flag_invalid ); 666129203Scognet if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF; 667129203Scognet } 668129203Scognet return (sbits32) 0x80000000; 669129203Scognet } 670129203Scognet else if ( aExp <= 0x7E ) { 671129203Scognet if ( aExp | aSig ) float_exception_flags |= float_flag_inexact; 672129203Scognet return 0; 673129203Scognet } 674129203Scognet aSig = ( aSig | 0x00800000 )<<8; 675129203Scognet z = aSig>>( - shiftCount ); 676129203Scognet if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) { 677129203Scognet float_exception_flags |= float_flag_inexact; 678129203Scognet } 679129203Scognet if ( aSign ) z = - z; 680129203Scognet return z; 681129203Scognet 682129203Scognet} 683129203Scognet 684129203Scognet/* 685129203Scognet------------------------------------------------------------------------------- 686129203ScognetReturns the result of converting the single-precision floating-point value 687129203Scognet`a' to the double-precision floating-point format. The conversion is 688129203Scognetperformed according to the IEC/IEEE Standard for Binary Floating-Point 689129203ScognetArithmetic. 690129203Scognet------------------------------------------------------------------------------- 691129203Scognet*/ 692129203Scognetfloat64 float32_to_float64( float32 a ) 693129203Scognet{ 694129203Scognet flag aSign; 695129203Scognet int16 aExp; 696129203Scognet bits32 aSig, zSig0, zSig1; 697129203Scognet 698129203Scognet aSig = extractFloat32Frac( a ); 699129203Scognet aExp = extractFloat32Exp( a ); 700129203Scognet aSign = extractFloat32Sign( a ); 701129203Scognet if ( aExp == 0xFF ) { 702129203Scognet if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a ) ); 703129203Scognet return packFloat64( aSign, 0x7FF, 0, 0 ); 704129203Scognet } 705129203Scognet if ( aExp == 0 ) { 706129203Scognet if ( aSig == 0 ) return packFloat64( aSign, 0, 0, 0 ); 707129203Scognet normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 708129203Scognet --aExp; 709129203Scognet } 710129203Scognet shift64Right( aSig, 0, 3, &zSig0, &zSig1 ); 711129203Scognet return packFloat64( aSign, aExp + 0x380, zSig0, zSig1 ); 712129203Scognet 713129203Scognet} 714129203Scognet 715129203Scognet#ifndef SOFTFLOAT_FOR_GCC 716129203Scognet/* 717129203Scognet------------------------------------------------------------------------------- 718129203ScognetRounds the single-precision floating-point value `a' to an integer, 719129203Scognetand returns the result as a single-precision floating-point value. The 720129203Scognetoperation is performed according to the IEC/IEEE Standard for Binary 721129203ScognetFloating-Point Arithmetic. 722129203Scognet------------------------------------------------------------------------------- 723129203Scognet*/ 724129203Scognetfloat32 float32_round_to_int( float32 a ) 725129203Scognet{ 726129203Scognet flag aSign; 727129203Scognet int16 aExp; 728129203Scognet bits32 lastBitMask, roundBitsMask; 729129203Scognet int8 roundingMode; 730129203Scognet float32 z; 731129203Scognet 732129203Scognet aExp = extractFloat32Exp( a ); 733129203Scognet if ( 0x96 <= aExp ) { 734129203Scognet if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) { 735129203Scognet return propagateFloat32NaN( a, a ); 736129203Scognet } 737129203Scognet return a; 738129203Scognet } 739129203Scognet if ( aExp <= 0x7E ) { 740129203Scognet if ( (bits32) ( a<<1 ) == 0 ) return a; 741129203Scognet float_exception_flags |= float_flag_inexact; 742129203Scognet aSign = extractFloat32Sign( a ); 743129203Scognet switch ( float_rounding_mode ) { 744129203Scognet case float_round_nearest_even: 745129203Scognet if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) { 746129203Scognet return packFloat32( aSign, 0x7F, 0 ); 747129203Scognet } 748129203Scognet break; 749129203Scognet case float_round_to_zero: 750129203Scognet break; 751129203Scognet case float_round_down: 752129203Scognet return aSign ? 0xBF800000 : 0; 753129203Scognet case float_round_up: 754129203Scognet return aSign ? 0x80000000 : 0x3F800000; 755129203Scognet } 756129203Scognet return packFloat32( aSign, 0, 0 ); 757129203Scognet } 758129203Scognet lastBitMask = 1; 759129203Scognet lastBitMask <<= 0x96 - aExp; 760129203Scognet roundBitsMask = lastBitMask - 1; 761129203Scognet z = a; 762129203Scognet roundingMode = float_rounding_mode; 763129203Scognet if ( roundingMode == float_round_nearest_even ) { 764129203Scognet z += lastBitMask>>1; 765129203Scognet if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask; 766129203Scognet } 767129203Scognet else if ( roundingMode != float_round_to_zero ) { 768129203Scognet if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) { 769129203Scognet z += roundBitsMask; 770129203Scognet } 771129203Scognet } 772129203Scognet z &= ~ roundBitsMask; 773129203Scognet if ( z != a ) float_exception_flags |= float_flag_inexact; 774129203Scognet return z; 775129203Scognet 776129203Scognet} 777129203Scognet#endif 778129203Scognet 779129203Scognet/* 780129203Scognet------------------------------------------------------------------------------- 781129203ScognetReturns the result of adding the absolute values of the single-precision 782129203Scognetfloating-point values `a' and `b'. If `zSign' is 1, the sum is negated 783129203Scognetbefore being returned. `zSign' is ignored if the result is a NaN. 784129203ScognetThe addition is performed according to the IEC/IEEE Standard for Binary 785129203ScognetFloating-Point Arithmetic. 786129203Scognet------------------------------------------------------------------------------- 787129203Scognet*/ 788129203Scognetstatic float32 addFloat32Sigs( float32 a, float32 b, flag zSign ) 789129203Scognet{ 790129203Scognet int16 aExp, bExp, zExp; 791129203Scognet bits32 aSig, bSig, zSig; 792129203Scognet int16 expDiff; 793129203Scognet 794129203Scognet aSig = extractFloat32Frac( a ); 795129203Scognet aExp = extractFloat32Exp( a ); 796129203Scognet bSig = extractFloat32Frac( b ); 797129203Scognet bExp = extractFloat32Exp( b ); 798129203Scognet expDiff = aExp - bExp; 799129203Scognet aSig <<= 6; 800129203Scognet bSig <<= 6; 801129203Scognet if ( 0 < expDiff ) { 802129203Scognet if ( aExp == 0xFF ) { 803129203Scognet if ( aSig ) return propagateFloat32NaN( a, b ); 804129203Scognet return a; 805129203Scognet } 806129203Scognet if ( bExp == 0 ) { 807129203Scognet --expDiff; 808129203Scognet } 809129203Scognet else { 810129203Scognet bSig |= 0x20000000; 811129203Scognet } 812129203Scognet shift32RightJamming( bSig, expDiff, &bSig ); 813129203Scognet zExp = aExp; 814129203Scognet } 815129203Scognet else if ( expDiff < 0 ) { 816129203Scognet if ( bExp == 0xFF ) { 817129203Scognet if ( bSig ) return propagateFloat32NaN( a, b ); 818129203Scognet return packFloat32( zSign, 0xFF, 0 ); 819129203Scognet } 820129203Scognet if ( aExp == 0 ) { 821129203Scognet ++expDiff; 822129203Scognet } 823129203Scognet else { 824129203Scognet aSig |= 0x20000000; 825129203Scognet } 826129203Scognet shift32RightJamming( aSig, - expDiff, &aSig ); 827129203Scognet zExp = bExp; 828129203Scognet } 829129203Scognet else { 830129203Scognet if ( aExp == 0xFF ) { 831129203Scognet if ( aSig | bSig ) return propagateFloat32NaN( a, b ); 832129203Scognet return a; 833129203Scognet } 834129203Scognet if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 ); 835129203Scognet zSig = 0x40000000 + aSig + bSig; 836129203Scognet zExp = aExp; 837129203Scognet goto roundAndPack; 838129203Scognet } 839129203Scognet aSig |= 0x20000000; 840129203Scognet zSig = ( aSig + bSig )<<1; 841129203Scognet --zExp; 842129203Scognet if ( (sbits32) zSig < 0 ) { 843129203Scognet zSig = aSig + bSig; 844129203Scognet ++zExp; 845129203Scognet } 846129203Scognet roundAndPack: 847129203Scognet return roundAndPackFloat32( zSign, zExp, zSig ); 848129203Scognet 849129203Scognet} 850129203Scognet 851129203Scognet/* 852129203Scognet------------------------------------------------------------------------------- 853129203ScognetReturns the result of subtracting the absolute values of the single- 854129203Scognetprecision floating-point values `a' and `b'. If `zSign' is 1, the 855129203Scognetdifference is negated before being returned. `zSign' is ignored if the 856129203Scognetresult is a NaN. The subtraction is performed according to the IEC/IEEE 857129203ScognetStandard for Binary Floating-Point Arithmetic. 858129203Scognet------------------------------------------------------------------------------- 859129203Scognet*/ 860129203Scognetstatic float32 subFloat32Sigs( float32 a, float32 b, flag zSign ) 861129203Scognet{ 862129203Scognet int16 aExp, bExp, zExp; 863129203Scognet bits32 aSig, bSig, zSig; 864129203Scognet int16 expDiff; 865129203Scognet 866129203Scognet aSig = extractFloat32Frac( a ); 867129203Scognet aExp = extractFloat32Exp( a ); 868129203Scognet bSig = extractFloat32Frac( b ); 869129203Scognet bExp = extractFloat32Exp( b ); 870129203Scognet expDiff = aExp - bExp; 871129203Scognet aSig <<= 7; 872129203Scognet bSig <<= 7; 873129203Scognet if ( 0 < expDiff ) goto aExpBigger; 874129203Scognet if ( expDiff < 0 ) goto bExpBigger; 875129203Scognet if ( aExp == 0xFF ) { 876129203Scognet if ( aSig | bSig ) return propagateFloat32NaN( a, b ); 877129203Scognet float_raise( float_flag_invalid ); 878129203Scognet return float32_default_nan; 879129203Scognet } 880129203Scognet if ( aExp == 0 ) { 881129203Scognet aExp = 1; 882129203Scognet bExp = 1; 883129203Scognet } 884129203Scognet if ( bSig < aSig ) goto aBigger; 885129203Scognet if ( aSig < bSig ) goto bBigger; 886129203Scognet return packFloat32( float_rounding_mode == float_round_down, 0, 0 ); 887129203Scognet bExpBigger: 888129203Scognet if ( bExp == 0xFF ) { 889129203Scognet if ( bSig ) return propagateFloat32NaN( a, b ); 890129203Scognet return packFloat32( zSign ^ 1, 0xFF, 0 ); 891129203Scognet } 892129203Scognet if ( aExp == 0 ) { 893129203Scognet ++expDiff; 894129203Scognet } 895129203Scognet else { 896129203Scognet aSig |= 0x40000000; 897129203Scognet } 898129203Scognet shift32RightJamming( aSig, - expDiff, &aSig ); 899129203Scognet bSig |= 0x40000000; 900129203Scognet bBigger: 901129203Scognet zSig = bSig - aSig; 902129203Scognet zExp = bExp; 903129203Scognet zSign ^= 1; 904129203Scognet goto normalizeRoundAndPack; 905129203Scognet aExpBigger: 906129203Scognet if ( aExp == 0xFF ) { 907129203Scognet if ( aSig ) return propagateFloat32NaN( a, b ); 908129203Scognet return a; 909129203Scognet } 910129203Scognet if ( bExp == 0 ) { 911129203Scognet --expDiff; 912129203Scognet } 913129203Scognet else { 914129203Scognet bSig |= 0x40000000; 915129203Scognet } 916129203Scognet shift32RightJamming( bSig, expDiff, &bSig ); 917129203Scognet aSig |= 0x40000000; 918129203Scognet aBigger: 919129203Scognet zSig = aSig - bSig; 920129203Scognet zExp = aExp; 921129203Scognet normalizeRoundAndPack: 922129203Scognet --zExp; 923129203Scognet return normalizeRoundAndPackFloat32( zSign, zExp, zSig ); 924129203Scognet 925129203Scognet} 926129203Scognet 927129203Scognet/* 928129203Scognet------------------------------------------------------------------------------- 929129203ScognetReturns the result of adding the single-precision floating-point values `a' 930129203Scognetand `b'. The operation is performed according to the IEC/IEEE Standard for 931129203ScognetBinary Floating-Point Arithmetic. 932129203Scognet------------------------------------------------------------------------------- 933129203Scognet*/ 934129203Scognetfloat32 float32_add( float32 a, float32 b ) 935129203Scognet{ 936129203Scognet flag aSign, bSign; 937129203Scognet 938129203Scognet aSign = extractFloat32Sign( a ); 939129203Scognet bSign = extractFloat32Sign( b ); 940129203Scognet if ( aSign == bSign ) { 941129203Scognet return addFloat32Sigs( a, b, aSign ); 942129203Scognet } 943129203Scognet else { 944129203Scognet return subFloat32Sigs( a, b, aSign ); 945129203Scognet } 946129203Scognet 947129203Scognet} 948129203Scognet 949129203Scognet/* 950129203Scognet------------------------------------------------------------------------------- 951129203ScognetReturns the result of subtracting the single-precision floating-point values 952129203Scognet`a' and `b'. The operation is performed according to the IEC/IEEE Standard 953129203Scognetfor Binary Floating-Point Arithmetic. 954129203Scognet------------------------------------------------------------------------------- 955129203Scognet*/ 956129203Scognetfloat32 float32_sub( float32 a, float32 b ) 957129203Scognet{ 958129203Scognet flag aSign, bSign; 959129203Scognet 960129203Scognet aSign = extractFloat32Sign( a ); 961129203Scognet bSign = extractFloat32Sign( b ); 962129203Scognet if ( aSign == bSign ) { 963129203Scognet return subFloat32Sigs( a, b, aSign ); 964129203Scognet } 965129203Scognet else { 966129203Scognet return addFloat32Sigs( a, b, aSign ); 967129203Scognet } 968129203Scognet 969129203Scognet} 970129203Scognet 971129203Scognet/* 972129203Scognet------------------------------------------------------------------------------- 973129203ScognetReturns the result of multiplying the single-precision floating-point values 974129203Scognet`a' and `b'. The operation is performed according to the IEC/IEEE Standard 975129203Scognetfor Binary Floating-Point Arithmetic. 976129203Scognet------------------------------------------------------------------------------- 977129203Scognet*/ 978129203Scognetfloat32 float32_mul( float32 a, float32 b ) 979129203Scognet{ 980129203Scognet flag aSign, bSign, zSign; 981129203Scognet int16 aExp, bExp, zExp; 982129203Scognet bits32 aSig, bSig, zSig0, zSig1; 983129203Scognet 984129203Scognet aSig = extractFloat32Frac( a ); 985129203Scognet aExp = extractFloat32Exp( a ); 986129203Scognet aSign = extractFloat32Sign( a ); 987129203Scognet bSig = extractFloat32Frac( b ); 988129203Scognet bExp = extractFloat32Exp( b ); 989129203Scognet bSign = extractFloat32Sign( b ); 990129203Scognet zSign = aSign ^ bSign; 991129203Scognet if ( aExp == 0xFF ) { 992129203Scognet if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) { 993129203Scognet return propagateFloat32NaN( a, b ); 994129203Scognet } 995129203Scognet if ( ( bExp | bSig ) == 0 ) { 996129203Scognet float_raise( float_flag_invalid ); 997129203Scognet return float32_default_nan; 998129203Scognet } 999129203Scognet return packFloat32( zSign, 0xFF, 0 ); 1000129203Scognet } 1001129203Scognet if ( bExp == 0xFF ) { 1002129203Scognet if ( bSig ) return propagateFloat32NaN( a, b ); 1003129203Scognet if ( ( aExp | aSig ) == 0 ) { 1004129203Scognet float_raise( float_flag_invalid ); 1005129203Scognet return float32_default_nan; 1006129203Scognet } 1007129203Scognet return packFloat32( zSign, 0xFF, 0 ); 1008129203Scognet } 1009129203Scognet if ( aExp == 0 ) { 1010129203Scognet if ( aSig == 0 ) return packFloat32( zSign, 0, 0 ); 1011129203Scognet normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 1012129203Scognet } 1013129203Scognet if ( bExp == 0 ) { 1014129203Scognet if ( bSig == 0 ) return packFloat32( zSign, 0, 0 ); 1015129203Scognet normalizeFloat32Subnormal( bSig, &bExp, &bSig ); 1016129203Scognet } 1017129203Scognet zExp = aExp + bExp - 0x7F; 1018129203Scognet aSig = ( aSig | 0x00800000 )<<7; 1019129203Scognet bSig = ( bSig | 0x00800000 )<<8; 1020129203Scognet mul32To64( aSig, bSig, &zSig0, &zSig1 ); 1021129203Scognet zSig0 |= ( zSig1 != 0 ); 1022129203Scognet if ( 0 <= (sbits32) ( zSig0<<1 ) ) { 1023129203Scognet zSig0 <<= 1; 1024129203Scognet --zExp; 1025129203Scognet } 1026129203Scognet return roundAndPackFloat32( zSign, zExp, zSig0 ); 1027129203Scognet 1028129203Scognet} 1029129203Scognet 1030129203Scognet/* 1031129203Scognet------------------------------------------------------------------------------- 1032129203ScognetReturns the result of dividing the single-precision floating-point value `a' 1033129203Scognetby the corresponding value `b'. The operation is performed according to the 1034129203ScognetIEC/IEEE Standard for Binary Floating-Point Arithmetic. 1035129203Scognet------------------------------------------------------------------------------- 1036129203Scognet*/ 1037129203Scognetfloat32 float32_div( float32 a, float32 b ) 1038129203Scognet{ 1039129203Scognet flag aSign, bSign, zSign; 1040129203Scognet int16 aExp, bExp, zExp; 1041129203Scognet bits32 aSig, bSig, zSig, rem0, rem1, term0, term1; 1042129203Scognet 1043129203Scognet aSig = extractFloat32Frac( a ); 1044129203Scognet aExp = extractFloat32Exp( a ); 1045129203Scognet aSign = extractFloat32Sign( a ); 1046129203Scognet bSig = extractFloat32Frac( b ); 1047129203Scognet bExp = extractFloat32Exp( b ); 1048129203Scognet bSign = extractFloat32Sign( b ); 1049129203Scognet zSign = aSign ^ bSign; 1050129203Scognet if ( aExp == 0xFF ) { 1051129203Scognet if ( aSig ) return propagateFloat32NaN( a, b ); 1052129203Scognet if ( bExp == 0xFF ) { 1053129203Scognet if ( bSig ) return propagateFloat32NaN( a, b ); 1054129203Scognet float_raise( float_flag_invalid ); 1055129203Scognet return float32_default_nan; 1056129203Scognet } 1057129203Scognet return packFloat32( zSign, 0xFF, 0 ); 1058129203Scognet } 1059129203Scognet if ( bExp == 0xFF ) { 1060129203Scognet if ( bSig ) return propagateFloat32NaN( a, b ); 1061129203Scognet return packFloat32( zSign, 0, 0 ); 1062129203Scognet } 1063129203Scognet if ( bExp == 0 ) { 1064129203Scognet if ( bSig == 0 ) { 1065129203Scognet if ( ( aExp | aSig ) == 0 ) { 1066129203Scognet float_raise( float_flag_invalid ); 1067129203Scognet return float32_default_nan; 1068129203Scognet } 1069129203Scognet float_raise( float_flag_divbyzero ); 1070129203Scognet return packFloat32( zSign, 0xFF, 0 ); 1071129203Scognet } 1072129203Scognet normalizeFloat32Subnormal( bSig, &bExp, &bSig ); 1073129203Scognet } 1074129203Scognet if ( aExp == 0 ) { 1075129203Scognet if ( aSig == 0 ) return packFloat32( zSign, 0, 0 ); 1076129203Scognet normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 1077129203Scognet } 1078129203Scognet zExp = aExp - bExp + 0x7D; 1079129203Scognet aSig = ( aSig | 0x00800000 )<<7; 1080129203Scognet bSig = ( bSig | 0x00800000 )<<8; 1081129203Scognet if ( bSig <= ( aSig + aSig ) ) { 1082129203Scognet aSig >>= 1; 1083129203Scognet ++zExp; 1084129203Scognet } 1085129203Scognet zSig = estimateDiv64To32( aSig, 0, bSig ); 1086129203Scognet if ( ( zSig & 0x3F ) <= 2 ) { 1087129203Scognet mul32To64( bSig, zSig, &term0, &term1 ); 1088129203Scognet sub64( aSig, 0, term0, term1, &rem0, &rem1 ); 1089129203Scognet while ( (sbits32) rem0 < 0 ) { 1090129203Scognet --zSig; 1091129203Scognet add64( rem0, rem1, 0, bSig, &rem0, &rem1 ); 1092129203Scognet } 1093129203Scognet zSig |= ( rem1 != 0 ); 1094129203Scognet } 1095129203Scognet return roundAndPackFloat32( zSign, zExp, zSig ); 1096129203Scognet 1097129203Scognet} 1098129203Scognet 1099129203Scognet#ifndef SOFTFLOAT_FOR_GCC 1100129203Scognet/* 1101129203Scognet------------------------------------------------------------------------------- 1102129203ScognetReturns the remainder of the single-precision floating-point value `a' 1103129203Scognetwith respect to the corresponding value `b'. The operation is performed 1104129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 1105129203Scognet------------------------------------------------------------------------------- 1106129203Scognet*/ 1107129203Scognetfloat32 float32_rem( float32 a, float32 b ) 1108129203Scognet{ 1109129203Scognet flag aSign, bSign, zSign; 1110129203Scognet int16 aExp, bExp, expDiff; 1111129203Scognet bits32 aSig, bSig, q, allZero, alternateASig; 1112129203Scognet sbits32 sigMean; 1113129203Scognet 1114129203Scognet aSig = extractFloat32Frac( a ); 1115129203Scognet aExp = extractFloat32Exp( a ); 1116129203Scognet aSign = extractFloat32Sign( a ); 1117129203Scognet bSig = extractFloat32Frac( b ); 1118129203Scognet bExp = extractFloat32Exp( b ); 1119129203Scognet bSign = extractFloat32Sign( b ); 1120129203Scognet if ( aExp == 0xFF ) { 1121129203Scognet if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) { 1122129203Scognet return propagateFloat32NaN( a, b ); 1123129203Scognet } 1124129203Scognet float_raise( float_flag_invalid ); 1125129203Scognet return float32_default_nan; 1126129203Scognet } 1127129203Scognet if ( bExp == 0xFF ) { 1128129203Scognet if ( bSig ) return propagateFloat32NaN( a, b ); 1129129203Scognet return a; 1130129203Scognet } 1131129203Scognet if ( bExp == 0 ) { 1132129203Scognet if ( bSig == 0 ) { 1133129203Scognet float_raise( float_flag_invalid ); 1134129203Scognet return float32_default_nan; 1135129203Scognet } 1136129203Scognet normalizeFloat32Subnormal( bSig, &bExp, &bSig ); 1137129203Scognet } 1138129203Scognet if ( aExp == 0 ) { 1139129203Scognet if ( aSig == 0 ) return a; 1140129203Scognet normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 1141129203Scognet } 1142129203Scognet expDiff = aExp - bExp; 1143129203Scognet aSig = ( aSig | 0x00800000 )<<8; 1144129203Scognet bSig = ( bSig | 0x00800000 )<<8; 1145129203Scognet if ( expDiff < 0 ) { 1146129203Scognet if ( expDiff < -1 ) return a; 1147129203Scognet aSig >>= 1; 1148129203Scognet } 1149129203Scognet q = ( bSig <= aSig ); 1150129203Scognet if ( q ) aSig -= bSig; 1151129203Scognet expDiff -= 32; 1152129203Scognet while ( 0 < expDiff ) { 1153129203Scognet q = estimateDiv64To32( aSig, 0, bSig ); 1154129203Scognet q = ( 2 < q ) ? q - 2 : 0; 1155129203Scognet aSig = - ( ( bSig>>2 ) * q ); 1156129203Scognet expDiff -= 30; 1157129203Scognet } 1158129203Scognet expDiff += 32; 1159129203Scognet if ( 0 < expDiff ) { 1160129203Scognet q = estimateDiv64To32( aSig, 0, bSig ); 1161129203Scognet q = ( 2 < q ) ? q - 2 : 0; 1162129203Scognet q >>= 32 - expDiff; 1163129203Scognet bSig >>= 2; 1164129203Scognet aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q; 1165129203Scognet } 1166129203Scognet else { 1167129203Scognet aSig >>= 2; 1168129203Scognet bSig >>= 2; 1169129203Scognet } 1170129203Scognet do { 1171129203Scognet alternateASig = aSig; 1172129203Scognet ++q; 1173129203Scognet aSig -= bSig; 1174129203Scognet } while ( 0 <= (sbits32) aSig ); 1175129203Scognet sigMean = aSig + alternateASig; 1176129203Scognet if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) { 1177129203Scognet aSig = alternateASig; 1178129203Scognet } 1179129203Scognet zSign = ( (sbits32) aSig < 0 ); 1180129203Scognet if ( zSign ) aSig = - aSig; 1181129203Scognet return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig ); 1182129203Scognet 1183129203Scognet} 1184129203Scognet#endif 1185129203Scognet 1186129203Scognet#ifndef SOFTFLOAT_FOR_GCC 1187129203Scognet/* 1188129203Scognet------------------------------------------------------------------------------- 1189129203ScognetReturns the square root of the single-precision floating-point value `a'. 1190129203ScognetThe operation is performed according to the IEC/IEEE Standard for Binary 1191129203ScognetFloating-Point Arithmetic. 1192129203Scognet------------------------------------------------------------------------------- 1193129203Scognet*/ 1194129203Scognetfloat32 float32_sqrt( float32 a ) 1195129203Scognet{ 1196129203Scognet flag aSign; 1197129203Scognet int16 aExp, zExp; 1198129203Scognet bits32 aSig, zSig, rem0, rem1, term0, term1; 1199129203Scognet 1200129203Scognet aSig = extractFloat32Frac( a ); 1201129203Scognet aExp = extractFloat32Exp( a ); 1202129203Scognet aSign = extractFloat32Sign( a ); 1203129203Scognet if ( aExp == 0xFF ) { 1204129203Scognet if ( aSig ) return propagateFloat32NaN( a, 0 ); 1205129203Scognet if ( ! aSign ) return a; 1206129203Scognet float_raise( float_flag_invalid ); 1207129203Scognet return float32_default_nan; 1208129203Scognet } 1209129203Scognet if ( aSign ) { 1210129203Scognet if ( ( aExp | aSig ) == 0 ) return a; 1211129203Scognet float_raise( float_flag_invalid ); 1212129203Scognet return float32_default_nan; 1213129203Scognet } 1214129203Scognet if ( aExp == 0 ) { 1215129203Scognet if ( aSig == 0 ) return 0; 1216129203Scognet normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 1217129203Scognet } 1218129203Scognet zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E; 1219129203Scognet aSig = ( aSig | 0x00800000 )<<8; 1220129203Scognet zSig = estimateSqrt32( aExp, aSig ) + 2; 1221129203Scognet if ( ( zSig & 0x7F ) <= 5 ) { 1222129203Scognet if ( zSig < 2 ) { 1223129203Scognet zSig = 0x7FFFFFFF; 1224129203Scognet goto roundAndPack; 1225129203Scognet } 1226129203Scognet else { 1227129203Scognet aSig >>= aExp & 1; 1228129203Scognet mul32To64( zSig, zSig, &term0, &term1 ); 1229129203Scognet sub64( aSig, 0, term0, term1, &rem0, &rem1 ); 1230129203Scognet while ( (sbits32) rem0 < 0 ) { 1231129203Scognet --zSig; 1232129203Scognet shortShift64Left( 0, zSig, 1, &term0, &term1 ); 1233129203Scognet term1 |= 1; 1234129203Scognet add64( rem0, rem1, term0, term1, &rem0, &rem1 ); 1235129203Scognet } 1236129203Scognet zSig |= ( ( rem0 | rem1 ) != 0 ); 1237129203Scognet } 1238129203Scognet } 1239129203Scognet shift32RightJamming( zSig, 1, &zSig ); 1240129203Scognet roundAndPack: 1241129203Scognet return roundAndPackFloat32( 0, zExp, zSig ); 1242129203Scognet 1243129203Scognet} 1244129203Scognet#endif 1245129203Scognet 1246129203Scognet/* 1247129203Scognet------------------------------------------------------------------------------- 1248129203ScognetReturns 1 if the single-precision floating-point value `a' is equal to 1249129203Scognetthe corresponding value `b', and 0 otherwise. The comparison is performed 1250129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 1251129203Scognet------------------------------------------------------------------------------- 1252129203Scognet*/ 1253129203Scognetflag float32_eq( float32 a, float32 b ) 1254129203Scognet{ 1255129203Scognet 1256129203Scognet if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1257129203Scognet || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1258129203Scognet ) { 1259129203Scognet if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) { 1260129203Scognet float_raise( float_flag_invalid ); 1261129203Scognet } 1262129203Scognet return 0; 1263129203Scognet } 1264129203Scognet return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 ); 1265129203Scognet 1266129203Scognet} 1267129203Scognet 1268129203Scognet/* 1269129203Scognet------------------------------------------------------------------------------- 1270129203ScognetReturns 1 if the single-precision floating-point value `a' is less than 1271129203Scognetor equal to the corresponding value `b', and 0 otherwise. The comparison 1272129203Scognetis performed according to the IEC/IEEE Standard for Binary Floating-Point 1273129203ScognetArithmetic. 1274129203Scognet------------------------------------------------------------------------------- 1275129203Scognet*/ 1276129203Scognetflag float32_le( float32 a, float32 b ) 1277129203Scognet{ 1278129203Scognet flag aSign, bSign; 1279129203Scognet 1280129203Scognet if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1281129203Scognet || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1282129203Scognet ) { 1283129203Scognet float_raise( float_flag_invalid ); 1284129203Scognet return 0; 1285129203Scognet } 1286129203Scognet aSign = extractFloat32Sign( a ); 1287129203Scognet bSign = extractFloat32Sign( b ); 1288129203Scognet if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 ); 1289129203Scognet return ( a == b ) || ( aSign ^ ( a < b ) ); 1290129203Scognet 1291129203Scognet} 1292129203Scognet 1293129203Scognet/* 1294129203Scognet------------------------------------------------------------------------------- 1295129203ScognetReturns 1 if the single-precision floating-point value `a' is less than 1296129203Scognetthe corresponding value `b', and 0 otherwise. The comparison is performed 1297129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 1298129203Scognet------------------------------------------------------------------------------- 1299129203Scognet*/ 1300129203Scognetflag float32_lt( float32 a, float32 b ) 1301129203Scognet{ 1302129203Scognet flag aSign, bSign; 1303129203Scognet 1304129203Scognet if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1305129203Scognet || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1306129203Scognet ) { 1307129203Scognet float_raise( float_flag_invalid ); 1308129203Scognet return 0; 1309129203Scognet } 1310129203Scognet aSign = extractFloat32Sign( a ); 1311129203Scognet bSign = extractFloat32Sign( b ); 1312129203Scognet if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 ); 1313129203Scognet return ( a != b ) && ( aSign ^ ( a < b ) ); 1314129203Scognet 1315129203Scognet} 1316129203Scognet 1317129203Scognet#ifndef SOFTFLOAT_FOR_GCC /* Not needed */ 1318129203Scognet/* 1319129203Scognet------------------------------------------------------------------------------- 1320129203ScognetReturns 1 if the single-precision floating-point value `a' is equal to 1321129203Scognetthe corresponding value `b', and 0 otherwise. The invalid exception is 1322129203Scognetraised if either operand is a NaN. Otherwise, the comparison is performed 1323129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 1324129203Scognet------------------------------------------------------------------------------- 1325129203Scognet*/ 1326129203Scognetflag float32_eq_signaling( float32 a, float32 b ) 1327129203Scognet{ 1328129203Scognet 1329129203Scognet if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1330129203Scognet || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1331129203Scognet ) { 1332129203Scognet float_raise( float_flag_invalid ); 1333129203Scognet return 0; 1334129203Scognet } 1335129203Scognet return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 ); 1336129203Scognet 1337129203Scognet} 1338129203Scognet 1339129203Scognet/* 1340129203Scognet------------------------------------------------------------------------------- 1341129203ScognetReturns 1 if the single-precision floating-point value `a' is less than or 1342129203Scognetequal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not 1343129203Scognetcause an exception. Otherwise, the comparison is performed according to the 1344129203ScognetIEC/IEEE Standard for Binary Floating-Point Arithmetic. 1345129203Scognet------------------------------------------------------------------------------- 1346129203Scognet*/ 1347129203Scognetflag float32_le_quiet( float32 a, float32 b ) 1348129203Scognet{ 1349129203Scognet flag aSign, bSign; 1350129203Scognet int16 aExp, bExp; 1351129203Scognet 1352129203Scognet if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1353129203Scognet || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1354129203Scognet ) { 1355129203Scognet if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) { 1356129203Scognet float_raise( float_flag_invalid ); 1357129203Scognet } 1358129203Scognet return 0; 1359129203Scognet } 1360129203Scognet aSign = extractFloat32Sign( a ); 1361129203Scognet bSign = extractFloat32Sign( b ); 1362129203Scognet if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 ); 1363129203Scognet return ( a == b ) || ( aSign ^ ( a < b ) ); 1364129203Scognet 1365129203Scognet} 1366129203Scognet 1367129203Scognet/* 1368129203Scognet------------------------------------------------------------------------------- 1369129203ScognetReturns 1 if the single-precision floating-point value `a' is less than 1370129203Scognetthe corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an 1371129203Scognetexception. Otherwise, the comparison is performed according to the IEC/IEEE 1372129203ScognetStandard for Binary Floating-Point Arithmetic. 1373129203Scognet------------------------------------------------------------------------------- 1374129203Scognet*/ 1375129203Scognetflag float32_lt_quiet( float32 a, float32 b ) 1376129203Scognet{ 1377129203Scognet flag aSign, bSign; 1378129203Scognet 1379129203Scognet if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1380129203Scognet || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1381129203Scognet ) { 1382129203Scognet if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) { 1383129203Scognet float_raise( float_flag_invalid ); 1384129203Scognet } 1385129203Scognet return 0; 1386129203Scognet } 1387129203Scognet aSign = extractFloat32Sign( a ); 1388129203Scognet bSign = extractFloat32Sign( b ); 1389129203Scognet if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 ); 1390129203Scognet return ( a != b ) && ( aSign ^ ( a < b ) ); 1391129203Scognet 1392129203Scognet} 1393129203Scognet#endif /* !SOFTFLOAT_FOR_GCC */ 1394129203Scognet 1395129203Scognet#ifndef SOFTFLOAT_FOR_GCC /* Not needed */ 1396129203Scognet/* 1397129203Scognet------------------------------------------------------------------------------- 1398129203ScognetReturns the result of converting the double-precision floating-point value 1399129203Scognet`a' to the 32-bit two's complement integer format. The conversion is 1400129203Scognetperformed according to the IEC/IEEE Standard for Binary Floating-Point 1401129203ScognetArithmetic---which means in particular that the conversion is rounded 1402129203Scognetaccording to the current rounding mode. If `a' is a NaN, the largest 1403129203Scognetpositive integer is returned. Otherwise, if the conversion overflows, the 1404129203Scognetlargest integer with the same sign as `a' is returned. 1405129203Scognet------------------------------------------------------------------------------- 1406129203Scognet*/ 1407129203Scognetint32 float64_to_int32( float64 a ) 1408129203Scognet{ 1409129203Scognet flag aSign; 1410129203Scognet int16 aExp, shiftCount; 1411129203Scognet bits32 aSig0, aSig1, absZ, aSigExtra; 1412129203Scognet int32 z; 1413129203Scognet int8 roundingMode; 1414129203Scognet 1415129203Scognet aSig1 = extractFloat64Frac1( a ); 1416129203Scognet aSig0 = extractFloat64Frac0( a ); 1417129203Scognet aExp = extractFloat64Exp( a ); 1418129203Scognet aSign = extractFloat64Sign( a ); 1419129203Scognet shiftCount = aExp - 0x413; 1420129203Scognet if ( 0 <= shiftCount ) { 1421129203Scognet if ( 0x41E < aExp ) { 1422129203Scognet if ( ( aExp == 0x7FF ) && ( aSig0 | aSig1 ) ) aSign = 0; 1423129203Scognet goto invalid; 1424129203Scognet } 1425129203Scognet shortShift64Left( 1426129203Scognet aSig0 | 0x00100000, aSig1, shiftCount, &absZ, &aSigExtra ); 1427129203Scognet if ( 0x80000000 < absZ ) goto invalid; 1428129203Scognet } 1429129203Scognet else { 1430129203Scognet aSig1 = ( aSig1 != 0 ); 1431129203Scognet if ( aExp < 0x3FE ) { 1432129203Scognet aSigExtra = aExp | aSig0 | aSig1; 1433129203Scognet absZ = 0; 1434129203Scognet } 1435129203Scognet else { 1436129203Scognet aSig0 |= 0x00100000; 1437129203Scognet aSigExtra = ( aSig0<<( shiftCount & 31 ) ) | aSig1; 1438129203Scognet absZ = aSig0>>( - shiftCount ); 1439129203Scognet } 1440129203Scognet } 1441129203Scognet roundingMode = float_rounding_mode; 1442129203Scognet if ( roundingMode == float_round_nearest_even ) { 1443129203Scognet if ( (sbits32) aSigExtra < 0 ) { 1444129203Scognet ++absZ; 1445129203Scognet if ( (bits32) ( aSigExtra<<1 ) == 0 ) absZ &= ~1; 1446129203Scognet } 1447129203Scognet z = aSign ? - absZ : absZ; 1448129203Scognet } 1449129203Scognet else { 1450129203Scognet aSigExtra = ( aSigExtra != 0 ); 1451129203Scognet if ( aSign ) { 1452129203Scognet z = - ( absZ 1453129203Scognet + ( ( roundingMode == float_round_down ) & aSigExtra ) ); 1454129203Scognet } 1455129203Scognet else { 1456129203Scognet z = absZ + ( ( roundingMode == float_round_up ) & aSigExtra ); 1457129203Scognet } 1458129203Scognet } 1459129203Scognet if ( ( aSign ^ ( z < 0 ) ) && z ) { 1460129203Scognet invalid: 1461129203Scognet float_raise( float_flag_invalid ); 1462129203Scognet return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF; 1463129203Scognet } 1464129203Scognet if ( aSigExtra ) float_exception_flags |= float_flag_inexact; 1465129203Scognet return z; 1466129203Scognet 1467129203Scognet} 1468129203Scognet#endif /* !SOFTFLOAT_FOR_GCC */ 1469129203Scognet 1470129203Scognet/* 1471129203Scognet------------------------------------------------------------------------------- 1472129203ScognetReturns the result of converting the double-precision floating-point value 1473129203Scognet`a' to the 32-bit two's complement integer format. The conversion is 1474129203Scognetperformed according to the IEC/IEEE Standard for Binary Floating-Point 1475129203ScognetArithmetic, except that the conversion is always rounded toward zero. 1476129203ScognetIf `a' is a NaN, the largest positive integer is returned. Otherwise, if 1477129203Scognetthe conversion overflows, the largest integer with the same sign as `a' is 1478129203Scognetreturned. 1479129203Scognet------------------------------------------------------------------------------- 1480129203Scognet*/ 1481129203Scognetint32 float64_to_int32_round_to_zero( float64 a ) 1482129203Scognet{ 1483129203Scognet flag aSign; 1484129203Scognet int16 aExp, shiftCount; 1485129203Scognet bits32 aSig0, aSig1, absZ, aSigExtra; 1486129203Scognet int32 z; 1487129203Scognet 1488129203Scognet aSig1 = extractFloat64Frac1( a ); 1489129203Scognet aSig0 = extractFloat64Frac0( a ); 1490129203Scognet aExp = extractFloat64Exp( a ); 1491129203Scognet aSign = extractFloat64Sign( a ); 1492129203Scognet shiftCount = aExp - 0x413; 1493129203Scognet if ( 0 <= shiftCount ) { 1494129203Scognet if ( 0x41E < aExp ) { 1495129203Scognet if ( ( aExp == 0x7FF ) && ( aSig0 | aSig1 ) ) aSign = 0; 1496129203Scognet goto invalid; 1497129203Scognet } 1498129203Scognet shortShift64Left( 1499129203Scognet aSig0 | 0x00100000, aSig1, shiftCount, &absZ, &aSigExtra ); 1500129203Scognet } 1501129203Scognet else { 1502129203Scognet if ( aExp < 0x3FF ) { 1503129203Scognet if ( aExp | aSig0 | aSig1 ) { 1504129203Scognet float_exception_flags |= float_flag_inexact; 1505129203Scognet } 1506129203Scognet return 0; 1507129203Scognet } 1508129203Scognet aSig0 |= 0x00100000; 1509129203Scognet aSigExtra = ( aSig0<<( shiftCount & 31 ) ) | aSig1; 1510129203Scognet absZ = aSig0>>( - shiftCount ); 1511129203Scognet } 1512129203Scognet z = aSign ? - absZ : absZ; 1513129203Scognet if ( ( aSign ^ ( z < 0 ) ) && z ) { 1514129203Scognet invalid: 1515129203Scognet float_raise( float_flag_invalid ); 1516129203Scognet return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF; 1517129203Scognet } 1518129203Scognet if ( aSigExtra ) float_exception_flags |= float_flag_inexact; 1519129203Scognet return z; 1520129203Scognet 1521129203Scognet} 1522129203Scognet 1523129203Scognet/* 1524129203Scognet------------------------------------------------------------------------------- 1525129203ScognetReturns the result of converting the double-precision floating-point value 1526129203Scognet`a' to the single-precision floating-point format. The conversion is 1527129203Scognetperformed according to the IEC/IEEE Standard for Binary Floating-Point 1528129203ScognetArithmetic. 1529129203Scognet------------------------------------------------------------------------------- 1530129203Scognet*/ 1531129203Scognetfloat32 float64_to_float32( float64 a ) 1532129203Scognet{ 1533129203Scognet flag aSign; 1534129203Scognet int16 aExp; 1535129203Scognet bits32 aSig0, aSig1, zSig; 1536129203Scognet bits32 allZero; 1537129203Scognet 1538129203Scognet aSig1 = extractFloat64Frac1( a ); 1539129203Scognet aSig0 = extractFloat64Frac0( a ); 1540129203Scognet aExp = extractFloat64Exp( a ); 1541129203Scognet aSign = extractFloat64Sign( a ); 1542129203Scognet if ( aExp == 0x7FF ) { 1543129203Scognet if ( aSig0 | aSig1 ) { 1544129203Scognet return commonNaNToFloat32( float64ToCommonNaN( a ) ); 1545129203Scognet } 1546129203Scognet return packFloat32( aSign, 0xFF, 0 ); 1547129203Scognet } 1548129203Scognet shift64RightJamming( aSig0, aSig1, 22, &allZero, &zSig ); 1549129203Scognet if ( aExp ) zSig |= 0x40000000; 1550129203Scognet return roundAndPackFloat32( aSign, aExp - 0x381, zSig ); 1551129203Scognet 1552129203Scognet} 1553129203Scognet 1554129203Scognet#ifndef SOFTFLOAT_FOR_GCC 1555129203Scognet/* 1556129203Scognet------------------------------------------------------------------------------- 1557129203ScognetRounds the double-precision floating-point value `a' to an integer, 1558129203Scognetand returns the result as a double-precision floating-point value. The 1559129203Scognetoperation is performed according to the IEC/IEEE Standard for Binary 1560129203ScognetFloating-Point Arithmetic. 1561129203Scognet------------------------------------------------------------------------------- 1562129203Scognet*/ 1563129203Scognetfloat64 float64_round_to_int( float64 a ) 1564129203Scognet{ 1565129203Scognet flag aSign; 1566129203Scognet int16 aExp; 1567129203Scognet bits32 lastBitMask, roundBitsMask; 1568129203Scognet int8 roundingMode; 1569129203Scognet float64 z; 1570129203Scognet 1571129203Scognet aExp = extractFloat64Exp( a ); 1572129203Scognet if ( 0x413 <= aExp ) { 1573129203Scognet if ( 0x433 <= aExp ) { 1574129203Scognet if ( ( aExp == 0x7FF ) 1575129203Scognet && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) ) { 1576129203Scognet return propagateFloat64NaN( a, a ); 1577129203Scognet } 1578129203Scognet return a; 1579129203Scognet } 1580129203Scognet lastBitMask = 1; 1581129203Scognet lastBitMask = ( lastBitMask<<( 0x432 - aExp ) )<<1; 1582129203Scognet roundBitsMask = lastBitMask - 1; 1583129203Scognet z = a; 1584129203Scognet roundingMode = float_rounding_mode; 1585129203Scognet if ( roundingMode == float_round_nearest_even ) { 1586129203Scognet if ( lastBitMask ) { 1587129203Scognet add64( z.high, z.low, 0, lastBitMask>>1, &z.high, &z.low ); 1588129203Scognet if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask; 1589129203Scognet } 1590129203Scognet else { 1591129203Scognet if ( (sbits32) z.low < 0 ) { 1592129203Scognet ++z.high; 1593129203Scognet if ( (bits32) ( z.low<<1 ) == 0 ) z.high &= ~1; 1594129203Scognet } 1595129203Scognet } 1596129203Scognet } 1597129203Scognet else if ( roundingMode != float_round_to_zero ) { 1598129203Scognet if ( extractFloat64Sign( z ) 1599129203Scognet ^ ( roundingMode == float_round_up ) ) { 1600129203Scognet add64( z.high, z.low, 0, roundBitsMask, &z.high, &z.low ); 1601129203Scognet } 1602129203Scognet } 1603129203Scognet z.low &= ~ roundBitsMask; 1604129203Scognet } 1605129203Scognet else { 1606129203Scognet if ( aExp <= 0x3FE ) { 1607129203Scognet if ( ( ( (bits32) ( a.high<<1 ) ) | a.low ) == 0 ) return a; 1608129203Scognet float_exception_flags |= float_flag_inexact; 1609129203Scognet aSign = extractFloat64Sign( a ); 1610129203Scognet switch ( float_rounding_mode ) { 1611129203Scognet case float_round_nearest_even: 1612129203Scognet if ( ( aExp == 0x3FE ) 1613129203Scognet && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) 1614129203Scognet ) { 1615129203Scognet return packFloat64( aSign, 0x3FF, 0, 0 ); 1616129203Scognet } 1617129203Scognet break; 1618129203Scognet case float_round_down: 1619129203Scognet return 1620129203Scognet aSign ? packFloat64( 1, 0x3FF, 0, 0 ) 1621129203Scognet : packFloat64( 0, 0, 0, 0 ); 1622129203Scognet case float_round_up: 1623129203Scognet return 1624129203Scognet aSign ? packFloat64( 1, 0, 0, 0 ) 1625129203Scognet : packFloat64( 0, 0x3FF, 0, 0 ); 1626129203Scognet } 1627129203Scognet return packFloat64( aSign, 0, 0, 0 ); 1628129203Scognet } 1629129203Scognet lastBitMask = 1; 1630129203Scognet lastBitMask <<= 0x413 - aExp; 1631129203Scognet roundBitsMask = lastBitMask - 1; 1632129203Scognet z.low = 0; 1633129203Scognet z.high = a.high; 1634129203Scognet roundingMode = float_rounding_mode; 1635129203Scognet if ( roundingMode == float_round_nearest_even ) { 1636129203Scognet z.high += lastBitMask>>1; 1637129203Scognet if ( ( ( z.high & roundBitsMask ) | a.low ) == 0 ) { 1638129203Scognet z.high &= ~ lastBitMask; 1639129203Scognet } 1640129203Scognet } 1641129203Scognet else if ( roundingMode != float_round_to_zero ) { 1642129203Scognet if ( extractFloat64Sign( z ) 1643129203Scognet ^ ( roundingMode == float_round_up ) ) { 1644129203Scognet z.high |= ( a.low != 0 ); 1645129203Scognet z.high += roundBitsMask; 1646129203Scognet } 1647129203Scognet } 1648129203Scognet z.high &= ~ roundBitsMask; 1649129203Scognet } 1650129203Scognet if ( ( z.low != a.low ) || ( z.high != a.high ) ) { 1651129203Scognet float_exception_flags |= float_flag_inexact; 1652129203Scognet } 1653129203Scognet return z; 1654129203Scognet 1655129203Scognet} 1656129203Scognet#endif 1657129203Scognet 1658129203Scognet/* 1659129203Scognet------------------------------------------------------------------------------- 1660129203ScognetReturns the result of adding the absolute values of the double-precision 1661129203Scognetfloating-point values `a' and `b'. If `zSign' is 1, the sum is negated 1662129203Scognetbefore being returned. `zSign' is ignored if the result is a NaN. 1663129203ScognetThe addition is performed according to the IEC/IEEE Standard for Binary 1664129203ScognetFloating-Point Arithmetic. 1665129203Scognet------------------------------------------------------------------------------- 1666129203Scognet*/ 1667129203Scognetstatic float64 addFloat64Sigs( float64 a, float64 b, flag zSign ) 1668129203Scognet{ 1669129203Scognet int16 aExp, bExp, zExp; 1670129203Scognet bits32 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2; 1671129203Scognet int16 expDiff; 1672129203Scognet 1673129203Scognet aSig1 = extractFloat64Frac1( a ); 1674129203Scognet aSig0 = extractFloat64Frac0( a ); 1675129203Scognet aExp = extractFloat64Exp( a ); 1676129203Scognet bSig1 = extractFloat64Frac1( b ); 1677129203Scognet bSig0 = extractFloat64Frac0( b ); 1678129203Scognet bExp = extractFloat64Exp( b ); 1679129203Scognet expDiff = aExp - bExp; 1680129203Scognet if ( 0 < expDiff ) { 1681129203Scognet if ( aExp == 0x7FF ) { 1682129203Scognet if ( aSig0 | aSig1 ) return propagateFloat64NaN( a, b ); 1683129203Scognet return a; 1684129203Scognet } 1685129203Scognet if ( bExp == 0 ) { 1686129203Scognet --expDiff; 1687129203Scognet } 1688129203Scognet else { 1689129203Scognet bSig0 |= 0x00100000; 1690129203Scognet } 1691129203Scognet shift64ExtraRightJamming( 1692129203Scognet bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 ); 1693129203Scognet zExp = aExp; 1694129203Scognet } 1695129203Scognet else if ( expDiff < 0 ) { 1696129203Scognet if ( bExp == 0x7FF ) { 1697129203Scognet if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b ); 1698129203Scognet return packFloat64( zSign, 0x7FF, 0, 0 ); 1699129203Scognet } 1700129203Scognet if ( aExp == 0 ) { 1701129203Scognet ++expDiff; 1702129203Scognet } 1703129203Scognet else { 1704129203Scognet aSig0 |= 0x00100000; 1705129203Scognet } 1706129203Scognet shift64ExtraRightJamming( 1707129203Scognet aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 ); 1708129203Scognet zExp = bExp; 1709129203Scognet } 1710129203Scognet else { 1711129203Scognet if ( aExp == 0x7FF ) { 1712129203Scognet if ( aSig0 | aSig1 | bSig0 | bSig1 ) { 1713129203Scognet return propagateFloat64NaN( a, b ); 1714129203Scognet } 1715129203Scognet return a; 1716129203Scognet } 1717129203Scognet add64( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 ); 1718129203Scognet if ( aExp == 0 ) return packFloat64( zSign, 0, zSig0, zSig1 ); 1719129203Scognet zSig2 = 0; 1720129203Scognet zSig0 |= 0x00200000; 1721129203Scognet zExp = aExp; 1722129203Scognet goto shiftRight1; 1723129203Scognet } 1724129203Scognet aSig0 |= 0x00100000; 1725129203Scognet add64( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 ); 1726129203Scognet --zExp; 1727129203Scognet if ( zSig0 < 0x00200000 ) goto roundAndPack; 1728129203Scognet ++zExp; 1729129203Scognet shiftRight1: 1730129203Scognet shift64ExtraRightJamming( zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 ); 1731129203Scognet roundAndPack: 1732129203Scognet return roundAndPackFloat64( zSign, zExp, zSig0, zSig1, zSig2 ); 1733129203Scognet 1734129203Scognet} 1735129203Scognet 1736129203Scognet/* 1737129203Scognet------------------------------------------------------------------------------- 1738129203ScognetReturns the result of subtracting the absolute values of the double- 1739129203Scognetprecision floating-point values `a' and `b'. If `zSign' is 1, the 1740129203Scognetdifference is negated before being returned. `zSign' is ignored if the 1741129203Scognetresult is a NaN. The subtraction is performed according to the IEC/IEEE 1742129203ScognetStandard for Binary Floating-Point Arithmetic. 1743129203Scognet------------------------------------------------------------------------------- 1744129203Scognet*/ 1745129203Scognetstatic float64 subFloat64Sigs( float64 a, float64 b, flag zSign ) 1746129203Scognet{ 1747129203Scognet int16 aExp, bExp, zExp; 1748129203Scognet bits32 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1; 1749129203Scognet int16 expDiff; 1750129203Scognet 1751129203Scognet aSig1 = extractFloat64Frac1( a ); 1752129203Scognet aSig0 = extractFloat64Frac0( a ); 1753129203Scognet aExp = extractFloat64Exp( a ); 1754129203Scognet bSig1 = extractFloat64Frac1( b ); 1755129203Scognet bSig0 = extractFloat64Frac0( b ); 1756129203Scognet bExp = extractFloat64Exp( b ); 1757129203Scognet expDiff = aExp - bExp; 1758129203Scognet shortShift64Left( aSig0, aSig1, 10, &aSig0, &aSig1 ); 1759129203Scognet shortShift64Left( bSig0, bSig1, 10, &bSig0, &bSig1 ); 1760129203Scognet if ( 0 < expDiff ) goto aExpBigger; 1761129203Scognet if ( expDiff < 0 ) goto bExpBigger; 1762129203Scognet if ( aExp == 0x7FF ) { 1763129203Scognet if ( aSig0 | aSig1 | bSig0 | bSig1 ) { 1764129203Scognet return propagateFloat64NaN( a, b ); 1765129203Scognet } 1766129203Scognet float_raise( float_flag_invalid ); 1767129203Scognet return float64_default_nan; 1768129203Scognet } 1769129203Scognet if ( aExp == 0 ) { 1770129203Scognet aExp = 1; 1771129203Scognet bExp = 1; 1772129203Scognet } 1773129203Scognet if ( bSig0 < aSig0 ) goto aBigger; 1774129203Scognet if ( aSig0 < bSig0 ) goto bBigger; 1775129203Scognet if ( bSig1 < aSig1 ) goto aBigger; 1776129203Scognet if ( aSig1 < bSig1 ) goto bBigger; 1777129203Scognet return packFloat64( float_rounding_mode == float_round_down, 0, 0, 0 ); 1778129203Scognet bExpBigger: 1779129203Scognet if ( bExp == 0x7FF ) { 1780129203Scognet if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b ); 1781129203Scognet return packFloat64( zSign ^ 1, 0x7FF, 0, 0 ); 1782129203Scognet } 1783129203Scognet if ( aExp == 0 ) { 1784129203Scognet ++expDiff; 1785129203Scognet } 1786129203Scognet else { 1787129203Scognet aSig0 |= 0x40000000; 1788129203Scognet } 1789129203Scognet shift64RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 ); 1790129203Scognet bSig0 |= 0x40000000; 1791129203Scognet bBigger: 1792129203Scognet sub64( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 ); 1793129203Scognet zExp = bExp; 1794129203Scognet zSign ^= 1; 1795129203Scognet goto normalizeRoundAndPack; 1796129203Scognet aExpBigger: 1797129203Scognet if ( aExp == 0x7FF ) { 1798129203Scognet if ( aSig0 | aSig1 ) return propagateFloat64NaN( a, b ); 1799129203Scognet return a; 1800129203Scognet } 1801129203Scognet if ( bExp == 0 ) { 1802129203Scognet --expDiff; 1803129203Scognet } 1804129203Scognet else { 1805129203Scognet bSig0 |= 0x40000000; 1806129203Scognet } 1807129203Scognet shift64RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 ); 1808129203Scognet aSig0 |= 0x40000000; 1809129203Scognet aBigger: 1810129203Scognet sub64( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 ); 1811129203Scognet zExp = aExp; 1812129203Scognet normalizeRoundAndPack: 1813129203Scognet --zExp; 1814129203Scognet return normalizeRoundAndPackFloat64( zSign, zExp - 10, zSig0, zSig1 ); 1815129203Scognet 1816129203Scognet} 1817129203Scognet 1818129203Scognet/* 1819129203Scognet------------------------------------------------------------------------------- 1820129203ScognetReturns the result of adding the double-precision floating-point values `a' 1821129203Scognetand `b'. The operation is performed according to the IEC/IEEE Standard for 1822129203ScognetBinary Floating-Point Arithmetic. 1823129203Scognet------------------------------------------------------------------------------- 1824129203Scognet*/ 1825129203Scognetfloat64 float64_add( float64 a, float64 b ) 1826129203Scognet{ 1827129203Scognet flag aSign, bSign; 1828129203Scognet 1829129203Scognet aSign = extractFloat64Sign( a ); 1830129203Scognet bSign = extractFloat64Sign( b ); 1831129203Scognet if ( aSign == bSign ) { 1832129203Scognet return addFloat64Sigs( a, b, aSign ); 1833129203Scognet } 1834129203Scognet else { 1835129203Scognet return subFloat64Sigs( a, b, aSign ); 1836129203Scognet } 1837129203Scognet 1838129203Scognet} 1839129203Scognet 1840129203Scognet/* 1841129203Scognet------------------------------------------------------------------------------- 1842129203ScognetReturns the result of subtracting the double-precision floating-point values 1843129203Scognet`a' and `b'. The operation is performed according to the IEC/IEEE Standard 1844129203Scognetfor Binary Floating-Point Arithmetic. 1845129203Scognet------------------------------------------------------------------------------- 1846129203Scognet*/ 1847129203Scognetfloat64 float64_sub( float64 a, float64 b ) 1848129203Scognet{ 1849129203Scognet flag aSign, bSign; 1850129203Scognet 1851129203Scognet aSign = extractFloat64Sign( a ); 1852129203Scognet bSign = extractFloat64Sign( b ); 1853129203Scognet if ( aSign == bSign ) { 1854129203Scognet return subFloat64Sigs( a, b, aSign ); 1855129203Scognet } 1856129203Scognet else { 1857129203Scognet return addFloat64Sigs( a, b, aSign ); 1858129203Scognet } 1859129203Scognet 1860129203Scognet} 1861129203Scognet 1862129203Scognet/* 1863129203Scognet------------------------------------------------------------------------------- 1864129203ScognetReturns the result of multiplying the double-precision floating-point values 1865129203Scognet`a' and `b'. The operation is performed according to the IEC/IEEE Standard 1866129203Scognetfor Binary Floating-Point Arithmetic. 1867129203Scognet------------------------------------------------------------------------------- 1868129203Scognet*/ 1869129203Scognetfloat64 float64_mul( float64 a, float64 b ) 1870129203Scognet{ 1871129203Scognet flag aSign, bSign, zSign; 1872129203Scognet int16 aExp, bExp, zExp; 1873129203Scognet bits32 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3; 1874129203Scognet 1875129203Scognet aSig1 = extractFloat64Frac1( a ); 1876129203Scognet aSig0 = extractFloat64Frac0( a ); 1877129203Scognet aExp = extractFloat64Exp( a ); 1878129203Scognet aSign = extractFloat64Sign( a ); 1879129203Scognet bSig1 = extractFloat64Frac1( b ); 1880129203Scognet bSig0 = extractFloat64Frac0( b ); 1881129203Scognet bExp = extractFloat64Exp( b ); 1882129203Scognet bSign = extractFloat64Sign( b ); 1883129203Scognet zSign = aSign ^ bSign; 1884129203Scognet if ( aExp == 0x7FF ) { 1885129203Scognet if ( ( aSig0 | aSig1 ) 1886129203Scognet || ( ( bExp == 0x7FF ) && ( bSig0 | bSig1 ) ) ) { 1887129203Scognet return propagateFloat64NaN( a, b ); 1888129203Scognet } 1889129203Scognet if ( ( bExp | bSig0 | bSig1 ) == 0 ) goto invalid; 1890129203Scognet return packFloat64( zSign, 0x7FF, 0, 0 ); 1891129203Scognet } 1892129203Scognet if ( bExp == 0x7FF ) { 1893129203Scognet if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b ); 1894129203Scognet if ( ( aExp | aSig0 | aSig1 ) == 0 ) { 1895129203Scognet invalid: 1896129203Scognet float_raise( float_flag_invalid ); 1897129203Scognet return float64_default_nan; 1898129203Scognet } 1899129203Scognet return packFloat64( zSign, 0x7FF, 0, 0 ); 1900129203Scognet } 1901129203Scognet if ( aExp == 0 ) { 1902129203Scognet if ( ( aSig0 | aSig1 ) == 0 ) return packFloat64( zSign, 0, 0, 0 ); 1903129203Scognet normalizeFloat64Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 ); 1904129203Scognet } 1905129203Scognet if ( bExp == 0 ) { 1906129203Scognet if ( ( bSig0 | bSig1 ) == 0 ) return packFloat64( zSign, 0, 0, 0 ); 1907129203Scognet normalizeFloat64Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 ); 1908129203Scognet } 1909129203Scognet zExp = aExp + bExp - 0x400; 1910129203Scognet aSig0 |= 0x00100000; 1911129203Scognet shortShift64Left( bSig0, bSig1, 12, &bSig0, &bSig1 ); 1912129203Scognet mul64To128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3 ); 1913129203Scognet add64( zSig0, zSig1, aSig0, aSig1, &zSig0, &zSig1 ); 1914129203Scognet zSig2 |= ( zSig3 != 0 ); 1915129203Scognet if ( 0x00200000 <= zSig0 ) { 1916129203Scognet shift64ExtraRightJamming( 1917129203Scognet zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 ); 1918129203Scognet ++zExp; 1919129203Scognet } 1920129203Scognet return roundAndPackFloat64( zSign, zExp, zSig0, zSig1, zSig2 ); 1921129203Scognet 1922129203Scognet} 1923129203Scognet 1924129203Scognet/* 1925129203Scognet------------------------------------------------------------------------------- 1926129203ScognetReturns the result of dividing the double-precision floating-point value `a' 1927129203Scognetby the corresponding value `b'. The operation is performed according to the 1928129203ScognetIEC/IEEE Standard for Binary Floating-Point Arithmetic. 1929129203Scognet------------------------------------------------------------------------------- 1930129203Scognet*/ 1931129203Scognetfloat64 float64_div( float64 a, float64 b ) 1932129203Scognet{ 1933129203Scognet flag aSign, bSign, zSign; 1934129203Scognet int16 aExp, bExp, zExp; 1935129203Scognet bits32 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2; 1936129203Scognet bits32 rem0, rem1, rem2, rem3, term0, term1, term2, term3; 1937129203Scognet 1938129203Scognet aSig1 = extractFloat64Frac1( a ); 1939129203Scognet aSig0 = extractFloat64Frac0( a ); 1940129203Scognet aExp = extractFloat64Exp( a ); 1941129203Scognet aSign = extractFloat64Sign( a ); 1942129203Scognet bSig1 = extractFloat64Frac1( b ); 1943129203Scognet bSig0 = extractFloat64Frac0( b ); 1944129203Scognet bExp = extractFloat64Exp( b ); 1945129203Scognet bSign = extractFloat64Sign( b ); 1946129203Scognet zSign = aSign ^ bSign; 1947129203Scognet if ( aExp == 0x7FF ) { 1948129203Scognet if ( aSig0 | aSig1 ) return propagateFloat64NaN( a, b ); 1949129203Scognet if ( bExp == 0x7FF ) { 1950129203Scognet if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b ); 1951129203Scognet goto invalid; 1952129203Scognet } 1953129203Scognet return packFloat64( zSign, 0x7FF, 0, 0 ); 1954129203Scognet } 1955129203Scognet if ( bExp == 0x7FF ) { 1956129203Scognet if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b ); 1957129203Scognet return packFloat64( zSign, 0, 0, 0 ); 1958129203Scognet } 1959129203Scognet if ( bExp == 0 ) { 1960129203Scognet if ( ( bSig0 | bSig1 ) == 0 ) { 1961129203Scognet if ( ( aExp | aSig0 | aSig1 ) == 0 ) { 1962129203Scognet invalid: 1963129203Scognet float_raise( float_flag_invalid ); 1964129203Scognet return float64_default_nan; 1965129203Scognet } 1966129203Scognet float_raise( float_flag_divbyzero ); 1967129203Scognet return packFloat64( zSign, 0x7FF, 0, 0 ); 1968129203Scognet } 1969129203Scognet normalizeFloat64Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 ); 1970129203Scognet } 1971129203Scognet if ( aExp == 0 ) { 1972129203Scognet if ( ( aSig0 | aSig1 ) == 0 ) return packFloat64( zSign, 0, 0, 0 ); 1973129203Scognet normalizeFloat64Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 ); 1974129203Scognet } 1975129203Scognet zExp = aExp - bExp + 0x3FD; 1976129203Scognet shortShift64Left( aSig0 | 0x00100000, aSig1, 11, &aSig0, &aSig1 ); 1977129203Scognet shortShift64Left( bSig0 | 0x00100000, bSig1, 11, &bSig0, &bSig1 ); 1978129203Scognet if ( le64( bSig0, bSig1, aSig0, aSig1 ) ) { 1979129203Scognet shift64Right( aSig0, aSig1, 1, &aSig0, &aSig1 ); 1980129203Scognet ++zExp; 1981129203Scognet } 1982129203Scognet zSig0 = estimateDiv64To32( aSig0, aSig1, bSig0 ); 1983129203Scognet mul64By32To96( bSig0, bSig1, zSig0, &term0, &term1, &term2 ); 1984129203Scognet sub96( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 ); 1985129203Scognet while ( (sbits32) rem0 < 0 ) { 1986129203Scognet --zSig0; 1987129203Scognet add96( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 ); 1988129203Scognet } 1989129203Scognet zSig1 = estimateDiv64To32( rem1, rem2, bSig0 ); 1990129203Scognet if ( ( zSig1 & 0x3FF ) <= 4 ) { 1991129203Scognet mul64By32To96( bSig0, bSig1, zSig1, &term1, &term2, &term3 ); 1992129203Scognet sub96( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 ); 1993129203Scognet while ( (sbits32) rem1 < 0 ) { 1994129203Scognet --zSig1; 1995129203Scognet add96( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 ); 1996129203Scognet } 1997129203Scognet zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 ); 1998129203Scognet } 1999129203Scognet shift64ExtraRightJamming( zSig0, zSig1, 0, 11, &zSig0, &zSig1, &zSig2 ); 2000129203Scognet return roundAndPackFloat64( zSign, zExp, zSig0, zSig1, zSig2 ); 2001129203Scognet 2002129203Scognet} 2003129203Scognet 2004129203Scognet#ifndef SOFTFLOAT_FOR_GCC 2005129203Scognet/* 2006129203Scognet------------------------------------------------------------------------------- 2007129203ScognetReturns the remainder of the double-precision floating-point value `a' 2008129203Scognetwith respect to the corresponding value `b'. The operation is performed 2009129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 2010129203Scognet------------------------------------------------------------------------------- 2011129203Scognet*/ 2012129203Scognetfloat64 float64_rem( float64 a, float64 b ) 2013129203Scognet{ 2014129203Scognet flag aSign, bSign, zSign; 2015129203Scognet int16 aExp, bExp, expDiff; 2016129203Scognet bits32 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2; 2017129203Scognet bits32 allZero, alternateASig0, alternateASig1, sigMean1; 2018129203Scognet sbits32 sigMean0; 2019129203Scognet float64 z; 2020129203Scognet 2021129203Scognet aSig1 = extractFloat64Frac1( a ); 2022129203Scognet aSig0 = extractFloat64Frac0( a ); 2023129203Scognet aExp = extractFloat64Exp( a ); 2024129203Scognet aSign = extractFloat64Sign( a ); 2025129203Scognet bSig1 = extractFloat64Frac1( b ); 2026129203Scognet bSig0 = extractFloat64Frac0( b ); 2027129203Scognet bExp = extractFloat64Exp( b ); 2028129203Scognet bSign = extractFloat64Sign( b ); 2029129203Scognet if ( aExp == 0x7FF ) { 2030129203Scognet if ( ( aSig0 | aSig1 ) 2031129203Scognet || ( ( bExp == 0x7FF ) && ( bSig0 | bSig1 ) ) ) { 2032129203Scognet return propagateFloat64NaN( a, b ); 2033129203Scognet } 2034129203Scognet goto invalid; 2035129203Scognet } 2036129203Scognet if ( bExp == 0x7FF ) { 2037129203Scognet if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b ); 2038129203Scognet return a; 2039129203Scognet } 2040129203Scognet if ( bExp == 0 ) { 2041129203Scognet if ( ( bSig0 | bSig1 ) == 0 ) { 2042129203Scognet invalid: 2043129203Scognet float_raise( float_flag_invalid ); 2044129203Scognet return float64_default_nan; 2045129203Scognet } 2046129203Scognet normalizeFloat64Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 ); 2047129203Scognet } 2048129203Scognet if ( aExp == 0 ) { 2049129203Scognet if ( ( aSig0 | aSig1 ) == 0 ) return a; 2050129203Scognet normalizeFloat64Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 ); 2051129203Scognet } 2052129203Scognet expDiff = aExp - bExp; 2053129203Scognet if ( expDiff < -1 ) return a; 2054129203Scognet shortShift64Left( 2055129203Scognet aSig0 | 0x00100000, aSig1, 11 - ( expDiff < 0 ), &aSig0, &aSig1 ); 2056129203Scognet shortShift64Left( bSig0 | 0x00100000, bSig1, 11, &bSig0, &bSig1 ); 2057129203Scognet q = le64( bSig0, bSig1, aSig0, aSig1 ); 2058129203Scognet if ( q ) sub64( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 ); 2059129203Scognet expDiff -= 32; 2060129203Scognet while ( 0 < expDiff ) { 2061129203Scognet q = estimateDiv64To32( aSig0, aSig1, bSig0 ); 2062129203Scognet q = ( 4 < q ) ? q - 4 : 0; 2063129203Scognet mul64By32To96( bSig0, bSig1, q, &term0, &term1, &term2 ); 2064129203Scognet shortShift96Left( term0, term1, term2, 29, &term1, &term2, &allZero ); 2065129203Scognet shortShift64Left( aSig0, aSig1, 29, &aSig0, &allZero ); 2066129203Scognet sub64( aSig0, 0, term1, term2, &aSig0, &aSig1 ); 2067129203Scognet expDiff -= 29; 2068129203Scognet } 2069129203Scognet if ( -32 < expDiff ) { 2070129203Scognet q = estimateDiv64To32( aSig0, aSig1, bSig0 ); 2071129203Scognet q = ( 4 < q ) ? q - 4 : 0; 2072129203Scognet q >>= - expDiff; 2073129203Scognet shift64Right( bSig0, bSig1, 8, &bSig0, &bSig1 ); 2074129203Scognet expDiff += 24; 2075129203Scognet if ( expDiff < 0 ) { 2076129203Scognet shift64Right( aSig0, aSig1, - expDiff, &aSig0, &aSig1 ); 2077129203Scognet } 2078129203Scognet else { 2079129203Scognet shortShift64Left( aSig0, aSig1, expDiff, &aSig0, &aSig1 ); 2080129203Scognet } 2081129203Scognet mul64By32To96( bSig0, bSig1, q, &term0, &term1, &term2 ); 2082129203Scognet sub64( aSig0, aSig1, term1, term2, &aSig0, &aSig1 ); 2083129203Scognet } 2084129203Scognet else { 2085129203Scognet shift64Right( aSig0, aSig1, 8, &aSig0, &aSig1 ); 2086129203Scognet shift64Right( bSig0, bSig1, 8, &bSig0, &bSig1 ); 2087129203Scognet } 2088129203Scognet do { 2089129203Scognet alternateASig0 = aSig0; 2090129203Scognet alternateASig1 = aSig1; 2091129203Scognet ++q; 2092129203Scognet sub64( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 ); 2093129203Scognet } while ( 0 <= (sbits32) aSig0 ); 2094129203Scognet add64( 2095129203Scognet aSig0, aSig1, alternateASig0, alternateASig1, &sigMean0, &sigMean1 ); 2096129203Scognet if ( ( sigMean0 < 0 ) 2097129203Scognet || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) { 2098129203Scognet aSig0 = alternateASig0; 2099129203Scognet aSig1 = alternateASig1; 2100129203Scognet } 2101129203Scognet zSign = ( (sbits32) aSig0 < 0 ); 2102129203Scognet if ( zSign ) sub64( 0, 0, aSig0, aSig1, &aSig0, &aSig1 ); 2103129203Scognet return 2104129203Scognet normalizeRoundAndPackFloat64( aSign ^ zSign, bExp - 4, aSig0, aSig1 ); 2105129203Scognet 2106129203Scognet} 2107129203Scognet#endif 2108129203Scognet 2109129203Scognet#ifndef SOFTFLOAT_FOR_GCC 2110129203Scognet/* 2111129203Scognet------------------------------------------------------------------------------- 2112129203ScognetReturns the square root of the double-precision floating-point value `a'. 2113129203ScognetThe operation is performed according to the IEC/IEEE Standard for Binary 2114129203ScognetFloating-Point Arithmetic. 2115129203Scognet------------------------------------------------------------------------------- 2116129203Scognet*/ 2117129203Scognetfloat64 float64_sqrt( float64 a ) 2118129203Scognet{ 2119129203Scognet flag aSign; 2120129203Scognet int16 aExp, zExp; 2121129203Scognet bits32 aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0; 2122129203Scognet bits32 rem0, rem1, rem2, rem3, term0, term1, term2, term3; 2123129203Scognet float64 z; 2124129203Scognet 2125129203Scognet aSig1 = extractFloat64Frac1( a ); 2126129203Scognet aSig0 = extractFloat64Frac0( a ); 2127129203Scognet aExp = extractFloat64Exp( a ); 2128129203Scognet aSign = extractFloat64Sign( a ); 2129129203Scognet if ( aExp == 0x7FF ) { 2130129203Scognet if ( aSig0 | aSig1 ) return propagateFloat64NaN( a, a ); 2131129203Scognet if ( ! aSign ) return a; 2132129203Scognet goto invalid; 2133129203Scognet } 2134129203Scognet if ( aSign ) { 2135129203Scognet if ( ( aExp | aSig0 | aSig1 ) == 0 ) return a; 2136129203Scognet invalid: 2137129203Scognet float_raise( float_flag_invalid ); 2138129203Scognet return float64_default_nan; 2139129203Scognet } 2140129203Scognet if ( aExp == 0 ) { 2141129203Scognet if ( ( aSig0 | aSig1 ) == 0 ) return packFloat64( 0, 0, 0, 0 ); 2142129203Scognet normalizeFloat64Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 ); 2143129203Scognet } 2144129203Scognet zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE; 2145129203Scognet aSig0 |= 0x00100000; 2146129203Scognet shortShift64Left( aSig0, aSig1, 11, &term0, &term1 ); 2147129203Scognet zSig0 = ( estimateSqrt32( aExp, term0 )>>1 ) + 1; 2148129203Scognet if ( zSig0 == 0 ) zSig0 = 0x7FFFFFFF; 2149129203Scognet doubleZSig0 = zSig0 + zSig0; 2150129203Scognet shortShift64Left( aSig0, aSig1, 9 - ( aExp & 1 ), &aSig0, &aSig1 ); 2151129203Scognet mul32To64( zSig0, zSig0, &term0, &term1 ); 2152129203Scognet sub64( aSig0, aSig1, term0, term1, &rem0, &rem1 ); 2153129203Scognet while ( (sbits32) rem0 < 0 ) { 2154129203Scognet --zSig0; 2155129203Scognet doubleZSig0 -= 2; 2156129203Scognet add64( rem0, rem1, 0, doubleZSig0 | 1, &rem0, &rem1 ); 2157129203Scognet } 2158129203Scognet zSig1 = estimateDiv64To32( rem1, 0, doubleZSig0 ); 2159129203Scognet if ( ( zSig1 & 0x1FF ) <= 5 ) { 2160129203Scognet if ( zSig1 == 0 ) zSig1 = 1; 2161129203Scognet mul32To64( doubleZSig0, zSig1, &term1, &term2 ); 2162129203Scognet sub64( rem1, 0, term1, term2, &rem1, &rem2 ); 2163129203Scognet mul32To64( zSig1, zSig1, &term2, &term3 ); 2164129203Scognet sub96( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 ); 2165129203Scognet while ( (sbits32) rem1 < 0 ) { 2166129203Scognet --zSig1; 2167129203Scognet shortShift64Left( 0, zSig1, 1, &term2, &term3 ); 2168129203Scognet term3 |= 1; 2169129203Scognet term2 |= doubleZSig0; 2170129203Scognet add96( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 ); 2171129203Scognet } 2172129203Scognet zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 ); 2173129203Scognet } 2174129203Scognet shift64ExtraRightJamming( zSig0, zSig1, 0, 10, &zSig0, &zSig1, &zSig2 ); 2175129203Scognet return roundAndPackFloat64( 0, zExp, zSig0, zSig1, zSig2 ); 2176129203Scognet 2177129203Scognet} 2178129203Scognet#endif 2179129203Scognet 2180129203Scognet/* 2181129203Scognet------------------------------------------------------------------------------- 2182129203ScognetReturns 1 if the double-precision floating-point value `a' is equal to 2183129203Scognetthe corresponding value `b', and 0 otherwise. The comparison is performed 2184129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 2185129203Scognet------------------------------------------------------------------------------- 2186129203Scognet*/ 2187129203Scognetflag float64_eq( float64 a, float64 b ) 2188129203Scognet{ 2189129203Scognet 2190129203Scognet if ( ( ( extractFloat64Exp( a ) == 0x7FF ) 2191129203Scognet && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) ) 2192129203Scognet || ( ( extractFloat64Exp( b ) == 0x7FF ) 2193129203Scognet && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) ) 2194129203Scognet ) { 2195129203Scognet if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) { 2196129203Scognet float_raise( float_flag_invalid ); 2197129203Scognet } 2198129203Scognet return 0; 2199129203Scognet } 2200129203Scognet return ( a == b ) || 2201129203Scognet ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) == 0 ); 2202129203Scognet 2203129203Scognet} 2204129203Scognet 2205129203Scognet/* 2206129203Scognet------------------------------------------------------------------------------- 2207129203ScognetReturns 1 if the double-precision floating-point value `a' is less than 2208129203Scognetor equal to the corresponding value `b', and 0 otherwise. The comparison 2209129203Scognetis performed according to the IEC/IEEE Standard for Binary Floating-Point 2210129203ScognetArithmetic. 2211129203Scognet------------------------------------------------------------------------------- 2212129203Scognet*/ 2213129203Scognetflag float64_le( float64 a, float64 b ) 2214129203Scognet{ 2215129203Scognet flag aSign, bSign; 2216129203Scognet 2217129203Scognet if ( ( ( extractFloat64Exp( a ) == 0x7FF ) 2218129203Scognet && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) ) 2219129203Scognet || ( ( extractFloat64Exp( b ) == 0x7FF ) 2220129203Scognet && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) ) 2221129203Scognet ) { 2222129203Scognet float_raise( float_flag_invalid ); 2223129203Scognet return 0; 2224129203Scognet } 2225129203Scognet aSign = extractFloat64Sign( a ); 2226129203Scognet bSign = extractFloat64Sign( b ); 2227129203Scognet if ( aSign != bSign ) 2228129203Scognet return aSign || 2229129203Scognet ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) == 2230129203Scognet 0 ); 2231129203Scognet return ( a == b ) || 2232129203Scognet ( aSign ^ ( FLOAT64_DEMANGLE(a) < FLOAT64_DEMANGLE(b) ) ); 2233129203Scognet} 2234129203Scognet 2235129203Scognet/* 2236129203Scognet------------------------------------------------------------------------------- 2237129203ScognetReturns 1 if the double-precision floating-point value `a' is less than 2238129203Scognetthe corresponding value `b', and 0 otherwise. The comparison is performed 2239129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 2240129203Scognet------------------------------------------------------------------------------- 2241129203Scognet*/ 2242129203Scognetflag float64_lt( float64 a, float64 b ) 2243129203Scognet{ 2244129203Scognet flag aSign, bSign; 2245129203Scognet 2246129203Scognet if ( ( ( extractFloat64Exp( a ) == 0x7FF ) 2247129203Scognet && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) ) 2248129203Scognet || ( ( extractFloat64Exp( b ) == 0x7FF ) 2249129203Scognet && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) ) 2250129203Scognet ) { 2251129203Scognet float_raise( float_flag_invalid ); 2252129203Scognet return 0; 2253129203Scognet } 2254129203Scognet aSign = extractFloat64Sign( a ); 2255129203Scognet bSign = extractFloat64Sign( b ); 2256129203Scognet if ( aSign != bSign ) 2257129203Scognet return aSign && 2258129203Scognet ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) != 2259129203Scognet 0 ); 2260129203Scognet return ( a != b ) && 2261129203Scognet ( aSign ^ ( FLOAT64_DEMANGLE(a) < FLOAT64_DEMANGLE(b) ) ); 2262129203Scognet 2263129203Scognet} 2264129203Scognet 2265129203Scognet#ifndef SOFTFLOAT_FOR_GCC 2266129203Scognet/* 2267129203Scognet------------------------------------------------------------------------------- 2268129203ScognetReturns 1 if the double-precision floating-point value `a' is equal to 2269129203Scognetthe corresponding value `b', and 0 otherwise. The invalid exception is 2270129203Scognetraised if either operand is a NaN. Otherwise, the comparison is performed 2271129203Scognetaccording to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. 2272129203Scognet------------------------------------------------------------------------------- 2273129203Scognet*/ 2274129203Scognetflag float64_eq_signaling( float64 a, float64 b ) 2275129203Scognet{ 2276129203Scognet 2277129203Scognet if ( ( ( extractFloat64Exp( a ) == 0x7FF ) 2278129203Scognet && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) ) 2279129203Scognet || ( ( extractFloat64Exp( b ) == 0x7FF ) 2280129203Scognet && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) ) 2281129203Scognet ) { 2282129203Scognet float_raise( float_flag_invalid ); 2283129203Scognet return 0; 2284129203Scognet } 2285129203Scognet return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 ); 2286129203Scognet 2287129203Scognet} 2288129203Scognet 2289129203Scognet/* 2290129203Scognet------------------------------------------------------------------------------- 2291129203ScognetReturns 1 if the double-precision floating-point value `a' is less than or 2292129203Scognetequal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not 2293129203Scognetcause an exception. Otherwise, the comparison is performed according to the 2294129203ScognetIEC/IEEE Standard for Binary Floating-Point Arithmetic. 2295129203Scognet------------------------------------------------------------------------------- 2296129203Scognet*/ 2297129203Scognetflag float64_le_quiet( float64 a, float64 b ) 2298129203Scognet{ 2299129203Scognet flag aSign, bSign; 2300129203Scognet 2301129203Scognet if ( ( ( extractFloat64Exp( a ) == 0x7FF ) 2302129203Scognet && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) ) 2303129203Scognet || ( ( extractFloat64Exp( b ) == 0x7FF ) 2304129203Scognet && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) ) 2305129203Scognet ) { 2306129203Scognet if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) { 2307129203Scognet float_raise( float_flag_invalid ); 2308129203Scognet } 2309129203Scognet return 0; 2310129203Scognet } 2311129203Scognet aSign = extractFloat64Sign( a ); 2312129203Scognet bSign = extractFloat64Sign( b ); 2313129203Scognet if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 ); 2314129203Scognet return ( a == b ) || ( aSign ^ ( a < b ) ); 2315129203Scognet 2316129203Scognet} 2317129203Scognet 2318129203Scognet/* 2319129203Scognet------------------------------------------------------------------------------- 2320129203ScognetReturns 1 if the double-precision floating-point value `a' is less than 2321129203Scognetthe corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an 2322129203Scognetexception. Otherwise, the comparison is performed according to the IEC/IEEE 2323129203ScognetStandard for Binary Floating-Point Arithmetic. 2324129203Scognet------------------------------------------------------------------------------- 2325129203Scognet*/ 2326129203Scognetflag float64_lt_quiet( float64 a, float64 b ) 2327129203Scognet{ 2328129203Scognet flag aSign, bSign; 2329129203Scognet 2330129203Scognet if ( ( ( extractFloat64Exp( a ) == 0x7FF ) 2331129203Scognet && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) ) 2332129203Scognet || ( ( extractFloat64Exp( b ) == 0x7FF ) 2333129203Scognet && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) ) 2334129203Scognet ) { 2335129203Scognet if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) { 2336129203Scognet float_raise( float_flag_invalid ); 2337129203Scognet } 2338129203Scognet return 0; 2339129203Scognet } 2340129203Scognet aSign = extractFloat64Sign( a ); 2341129203Scognet bSign = extractFloat64Sign( b ); 2342129203Scognet if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 ); 2343129203Scognet return ( a != b ) && ( aSign ^ ( a < b ) ); 2344129203Scognet 2345129203Scognet} 2346129203Scognet 2347129203Scognet#endif 2348