git.oblomov.eu Git - linux-2.6/blob - arch/sh/kernel/cpu/sh4/softfloat.c

   1 /*
   2  * Floating point emulation support for subnormalised numbers on SH4
   3  * architecture This file is derived from the SoftFloat IEC/IEEE
   4  * Floating-point Arithmetic Package, Release 2 the original license of
   5  * which is reproduced below.
   6  *
   7  * ========================================================================
   8  *
   9  * This C source file is part of the SoftFloat IEC/IEEE Floating-point
  10  * Arithmetic Package, Release 2.
  11  *
  12  * Written by John R. Hauser.  This work was made possible in part by the
  13  * International Computer Science Institute, located at Suite 600, 1947 Center
  14  * Street, Berkeley, California 94704.  Funding was partially provided by the
  15  * National Science Foundation under grant MIP-9311980.  The original version
  16  * of this code was written as part of a project to build a fixed-point vector
  17  * processor in collaboration with the University of California at Berkeley,
  18  * overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
  19  * is available through the web page `http://HTTP.CS.Berkeley.EDU/~jhauser/
  20  * arithmetic/softfloat.html'.
  21  *
  22  * THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
  23  * has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
  24  * TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
  25  * PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
  26  * AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
  27  *
  28  * Derivative works are acceptable, even for commercial purposes, so long as
  29  * (1) they include prominent notice that the work is derivative, and (2) they
  30  * include prominent notice akin to these three paragraphs for those parts of
  31  * this code that are retained.
  32  *
  33  * ========================================================================
  34  *
  35  * SH4 modifications by Ismail Dhaoui <ismail.dhaoui@st.com>
  36  * and Kamel Khelifi <kamel.khelifi@st.com>
  37  */
  38 #include <linux/kernel.h>
  39 #include <cpu/fpu.h>
  40
  41 #define LIT64( a ) a##LL
  42
  43 typedef char flag;
  44 typedef unsigned char uint8;
  45 typedef signed char int8;
  46 typedef int uint16;
  47 typedef int int16;
  48 typedef unsigned int uint32;
  49 typedef signed int int32;
  50
  51 typedef unsigned long long int bits64;
  52 typedef signed long long int sbits64;
  53
  54 typedef unsigned char bits8;
  55 typedef signed char sbits8;
  56 typedef unsigned short int bits16;
  57 typedef signed short int sbits16;
  58 typedef unsigned int bits32;
  59 typedef signed int sbits32;
  60
  61 typedef unsigned long long int uint64;
  62 typedef signed long long int int64;
  63
  64 typedef unsigned long int float32;
  65 typedef unsigned long long float64;
  66
  67 extern void float_raise(unsigned int flags);    /* in fpu.c */
  68 extern int float_rounding_mode(void);   /* in fpu.c */
  69
  70 inline bits64 extractFloat64Frac(float64 a);
  71 inline flag extractFloat64Sign(float64 a);
  72 inline int16 extractFloat64Exp(float64 a);
  73 inline int16 extractFloat32Exp(float32 a);
  74 inline flag extractFloat32Sign(float32 a);
  75 inline bits32 extractFloat32Frac(float32 a);
  76 inline float64 packFloat64(flag zSign, int16 zExp, bits64 zSig);
  77 inline void shift64RightJamming(bits64 a, int16 count, bits64 * zPtr);
  78 inline float32 packFloat32(flag zSign, int16 zExp, bits32 zSig);
  79 inline void shift32RightJamming(bits32 a, int16 count, bits32 * zPtr);
  80 float64 float64_sub(float64 a, float64 b);
  81 float32 float32_sub(float32 a, float32 b);
  82 float32 float32_add(float32 a, float32 b);
  83 float64 float64_add(float64 a, float64 b);
  84 float64 float64_div(float64 a, float64 b);
  85 float32 float32_div(float32 a, float32 b);
  86 float32 float32_mul(float32 a, float32 b);
  87 float64 float64_mul(float64 a, float64 b);
  88 inline void add128(bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 * z0Ptr,
  89                    bits64 * z1Ptr);
  90 inline void sub128(bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 * z0Ptr,
  91                    bits64 * z1Ptr);
  92 inline void mul64To128(bits64 a, bits64 b, bits64 * z0Ptr, bits64 * z1Ptr);
  93
  94 static int8 countLeadingZeros32(bits32 a);
  95 static int8 countLeadingZeros64(bits64 a);
  96 static float64 normalizeRoundAndPackFloat64(flag zSign, int16 zExp,
  97                                             bits64 zSig);
  98 static float64 subFloat64Sigs(float64 a, float64 b, flag zSign);
  99 static float64 addFloat64Sigs(float64 a, float64 b, flag zSign);
 100 static float32 roundAndPackFloat32(flag zSign, int16 zExp, bits32 zSig);
 101 static float32 normalizeRoundAndPackFloat32(flag zSign, int16 zExp,
 102                                             bits32 zSig);
 103 static float64 roundAndPackFloat64(flag zSign, int16 zExp, bits64 zSig);
 104 static float32 subFloat32Sigs(float32 a, float32 b, flag zSign);
 105 static float32 addFloat32Sigs(float32 a, float32 b, flag zSign);
 106 static void normalizeFloat64Subnormal(bits64 aSig, int16 * zExpPtr,
 107                                       bits64 * zSigPtr);
 108 static bits64 estimateDiv128To64(bits64 a0, bits64 a1, bits64 b);
 109 static void normalizeFloat32Subnormal(bits32 aSig, int16 * zExpPtr,
 110                                       bits32 * zSigPtr);
 111
 112 inline bits64 extractFloat64Frac(float64 a)
 113 {
 114         return a & LIT64(0x000FFFFFFFFFFFFF);
 115 }
 116
 117 inline flag extractFloat64Sign(float64 a)
 118 {
 119         return a >> 63;
 120 }
 121
 122 inline int16 extractFloat64Exp(float64 a)
 123 {
 124         return (a >> 52) & 0x7FF;
 125 }
 126
 127 inline int16 extractFloat32Exp(float32 a)
 128 {
 129         return (a >> 23) & 0xFF;
 130 }
 131
 132 inline flag extractFloat32Sign(float32 a)
 133 {
 134         return a >> 31;
 135 }
 136
 137 inline bits32 extractFloat32Frac(float32 a)
 138 {
 139         return a & 0x007FFFFF;
 140 }
 141
 142 inline float64 packFloat64(flag zSign, int16 zExp, bits64 zSig)
 143 {
 144         return (((bits64) zSign) << 63) + (((bits64) zExp) << 52) + zSig;
 145 }
 146
 147 inline void shift64RightJamming(bits64 a, int16 count, bits64 * zPtr)
 148 {
 149         bits64 z;
 150
 151         if (count == 0) {
 152                 z = a;
 153         } else if (count < 64) {
 154                 z = (a >> count) | ((a << ((-count) & 63)) != 0);
 155         } else {
 156                 z = (a != 0);
 157         }
 158         *zPtr = z;
 159 }
 160
 161 static int8 countLeadingZeros32(bits32 a)
 162 {
 163         static const int8 countLeadingZerosHigh[] = {
 164                 8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
 165                 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
 166                 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 167                 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 168                 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 169                 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 170                 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 171                 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 172                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 173                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 174                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 175                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 176                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 177                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 178                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 179                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
 180         };
 181         int8 shiftCount;
 182
 183         shiftCount = 0;
 184         if (a < 0x10000) {
 185                 shiftCount += 16;
 186                 a <<= 16;
 187         }
 188         if (a < 0x1000000) {
 189                 shiftCount += 8;
 190                 a <<= 8;
 191         }
 192         shiftCount += countLeadingZerosHigh[a >> 24];
 193         return shiftCount;
 194
 195 }
 196
 197 static int8 countLeadingZeros64(bits64 a)
 198 {
 199         int8 shiftCount;
 200
 201         shiftCount = 0;
 202         if (a < ((bits64) 1) << 32) {
 203                 shiftCount += 32;
 204         } else {
 205                 a >>= 32;
 206         }
 207         shiftCount += countLeadingZeros32(a);
 208         return shiftCount;
 209
 210 }
 211
 212 static float64 normalizeRoundAndPackFloat64(flag zSign, int16 zExp, bits64 zSig)
 213 {
 214         int8 shiftCount;
 215
 216         shiftCount = countLeadingZeros64(zSig) - 1;
 217         return roundAndPackFloat64(zSign, zExp - shiftCount,
 218                                    zSig << shiftCount);
 219
 220 }
 221
 222 static float64 subFloat64Sigs(float64 a, float64 b, flag zSign)
 223 {
 224         int16 aExp, bExp, zExp;
 225         bits64 aSig, bSig, zSig;
 226         int16 expDiff;
 227
 228         aSig = extractFloat64Frac(a);
 229         aExp = extractFloat64Exp(a);
 230         bSig = extractFloat64Frac(b);
 231         bExp = extractFloat64Exp(b);
 232         expDiff = aExp - bExp;
 233         aSig <<= 10;
 234         bSig <<= 10;
 235         if (0 < expDiff)
 236                 goto aExpBigger;
 237         if (expDiff < 0)
 238                 goto bExpBigger;
 239         if (aExp == 0) {
 240                 aExp = 1;
 241                 bExp = 1;
 242         }
 243         if (bSig < aSig)
 244                 goto aBigger;
 245         if (aSig < bSig)
 246                 goto bBigger;
 247         return packFloat64(float_rounding_mode() == FPSCR_RM_ZERO, 0, 0);
 248       bExpBigger:
 249         if (bExp == 0x7FF) {
 250                 return packFloat64(zSign ^ 1, 0x7FF, 0);
 251         }
 252         if (aExp == 0) {
 253                 ++expDiff;
 254         } else {
 255                 aSig |= LIT64(0x4000000000000000);
 256         }
 257         shift64RightJamming(aSig, -expDiff, &aSig);
 258         bSig |= LIT64(0x4000000000000000);
 259       bBigger:
 260         zSig = bSig - aSig;
 261         zExp = bExp;
 262         zSign ^= 1;
 263         goto normalizeRoundAndPack;
 264       aExpBigger:
 265         if (aExp == 0x7FF) {
 266                 return a;
 267         }
 268         if (bExp == 0) {
 269                 --expDiff;
 270         } else {
 271                 bSig |= LIT64(0x4000000000000000);
 272         }
 273         shift64RightJamming(bSig, expDiff, &bSig);
 274         aSig |= LIT64(0x4000000000000000);
 275       aBigger:
 276         zSig = aSig - bSig;
 277         zExp = aExp;
 278       normalizeRoundAndPack:
 279         --zExp;
 280         return normalizeRoundAndPackFloat64(zSign, zExp, zSig);
 281
 282 }
 283 static float64 addFloat64Sigs(float64 a, float64 b, flag zSign)
 284 {
 285         int16 aExp, bExp, zExp;
 286         bits64 aSig, bSig, zSig;
 287         int16 expDiff;
 288
 289         aSig = extractFloat64Frac(a);
 290         aExp = extractFloat64Exp(a);
 291         bSig = extractFloat64Frac(b);
 292         bExp = extractFloat64Exp(b);
 293         expDiff = aExp - bExp;
 294         aSig <<= 9;
 295         bSig <<= 9;
 296         if (0 < expDiff) {
 297                 if (aExp == 0x7FF) {
 298                         return a;
 299                 }
 300                 if (bExp == 0) {
 301                         --expDiff;
 302                 } else {
 303                         bSig |= LIT64(0x2000000000000000);
 304                 }
 305                 shift64RightJamming(bSig, expDiff, &bSig);
 306                 zExp = aExp;
 307         } else if (expDiff < 0) {
 308                 if (bExp == 0x7FF) {
 309                         return packFloat64(zSign, 0x7FF, 0);
 310                 }
 311                 if (aExp == 0) {
 312                         ++expDiff;
 313                 } else {
 314                         aSig |= LIT64(0x2000000000000000);
 315                 }
 316                 shift64RightJamming(aSig, -expDiff, &aSig);
 317                 zExp = bExp;
 318         } else {
 319                 if (aExp == 0x7FF) {
 320                         return a;
 321                 }
 322                 if (aExp == 0)
 323                         return packFloat64(zSign, 0, (aSig + bSig) >> 9);
 324                 zSig = LIT64(0x4000000000000000) + aSig + bSig;
 325                 zExp = aExp;
 326                 goto roundAndPack;
 327         }
 328         aSig |= LIT64(0x2000000000000000);
 329         zSig = (aSig + bSig) << 1;
 330         --zExp;
 331         if ((sbits64) zSig < 0) {
 332                 zSig = aSig + bSig;
 333                 ++zExp;
 334         }
 335       roundAndPack:
 336         return roundAndPackFloat64(zSign, zExp, zSig);
 337
 338 }
 339
 340 inline float32 packFloat32(flag zSign, int16 zExp, bits32 zSig)
 341 {
 342         return (((bits32) zSign) << 31) + (((bits32) zExp) << 23) + zSig;
 343 }
 344
 345 inline void shift32RightJamming(bits32 a, int16 count, bits32 * zPtr)
 346 {
 347         bits32 z;
 348         if (count == 0) {
 349                 z = a;
 350         } else if (count < 32) {
 351                 z = (a >> count) | ((a << ((-count) & 31)) != 0);
 352         } else {
 353                 z = (a != 0);
 354         }
 355         *zPtr = z;
 356 }
 357
 358 static float32 roundAndPackFloat32(flag zSign, int16 zExp, bits32 zSig)
 359 {
 360         flag roundNearestEven;
 361         int8 roundIncrement, roundBits;
 362         flag isTiny;
 363
 364         /* SH4 has only 2 rounding modes - round to nearest and round to zero */
 365         roundNearestEven = (float_rounding_mode() == FPSCR_RM_NEAREST);
 366         roundIncrement = 0x40;
 367         if (!roundNearestEven) {
 368                 roundIncrement = 0;
 369         }
 370         roundBits = zSig & 0x7F;
 371         if (0xFD <= (bits16) zExp) {
 372                 if ((0xFD < zExp)
 373                     || ((zExp == 0xFD)
 374                         && ((sbits32) (zSig + roundIncrement) < 0))
 375                     ) {
 376                         float_raise(FPSCR_CAUSE_OVERFLOW | FPSCR_CAUSE_INEXACT);
 377                         return packFloat32(zSign, 0xFF,
 378                                            0) - (roundIncrement == 0);
 379                 }
 380                 if (zExp < 0) {
 381                         isTiny = (zExp < -1)
 382                             || (zSig + roundIncrement < 0x80000000);
 383                         shift32RightJamming(zSig, -zExp, &zSig);
 384                         zExp = 0;
 385                         roundBits = zSig & 0x7F;
 386                         if (isTiny && roundBits)
 387                                 float_raise(FPSCR_CAUSE_UNDERFLOW);
 388                 }
 389         }
 390         if (roundBits)
 391                 float_raise(FPSCR_CAUSE_INEXACT);
 392         zSig = (zSig + roundIncrement) >> 7;
 393         zSig &= ~(((roundBits ^ 0x40) == 0) & roundNearestEven);
 394         if (zSig == 0)
 395                 zExp = 0;
 396         return packFloat32(zSign, zExp, zSig);
 397
 398 }
 399
 400 static float32 normalizeRoundAndPackFloat32(flag zSign, int16 zExp, bits32 zSig)
 401 {
 402         int8 shiftCount;
 403
 404         shiftCount = countLeadingZeros32(zSig) - 1;
 405         return roundAndPackFloat32(zSign, zExp - shiftCount,
 406                                    zSig << shiftCount);
 407 }
 408
 409 static float64 roundAndPackFloat64(flag zSign, int16 zExp, bits64 zSig)
 410 {
 411         flag roundNearestEven;
 412         int16 roundIncrement, roundBits;
 413         flag isTiny;
 414
 415         /* SH4 has only 2 rounding modes - round to nearest and round to zero */
 416         roundNearestEven = (float_rounding_mode() == FPSCR_RM_NEAREST);
 417         roundIncrement = 0x200;
 418         if (!roundNearestEven) {
 419                 roundIncrement = 0;
 420         }
 421         roundBits = zSig & 0x3FF;
 422         if (0x7FD <= (bits16) zExp) {
 423                 if ((0x7FD < zExp)
 424                     || ((zExp == 0x7FD)
 425                         && ((sbits64) (zSig + roundIncrement) < 0))
 426                     ) {
 427                         float_raise(FPSCR_CAUSE_OVERFLOW | FPSCR_CAUSE_INEXACT);
 428                         return packFloat64(zSign, 0x7FF,
 429                                            0) - (roundIncrement == 0);
 430                 }
 431                 if (zExp < 0) {
 432                         isTiny = (zExp < -1)
 433                             || (zSig + roundIncrement <
 434                                 LIT64(0x8000000000000000));
 435                         shift64RightJamming(zSig, -zExp, &zSig);
 436                         zExp = 0;
 437                         roundBits = zSig & 0x3FF;
 438                         if (isTiny && roundBits)
 439                                 float_raise(FPSCR_CAUSE_UNDERFLOW);
 440                 }
 441         }
 442         if (roundBits)
 443                 float_raise(FPSCR_CAUSE_INEXACT);
 444         zSig = (zSig + roundIncrement) >> 10;
 445         zSig &= ~(((roundBits ^ 0x200) == 0) & roundNearestEven);
 446         if (zSig == 0)
 447                 zExp = 0;
 448         return packFloat64(zSign, zExp, zSig);
 449
 450 }
 451
 452 static float32 subFloat32Sigs(float32 a, float32 b, flag zSign)
 453 {
 454         int16 aExp, bExp, zExp;
 455         bits32 aSig, bSig, zSig;
 456         int16 expDiff;
 457
 458         aSig = extractFloat32Frac(a);
 459         aExp = extractFloat32Exp(a);
 460         bSig = extractFloat32Frac(b);
 461         bExp = extractFloat32Exp(b);
 462         expDiff = aExp - bExp;
 463         aSig <<= 7;
 464         bSig <<= 7;
 465         if (0 < expDiff)
 466                 goto aExpBigger;
 467         if (expDiff < 0)
 468                 goto bExpBigger;
 469         if (aExp == 0) {
 470                 aExp = 1;
 471                 bExp = 1;
 472         }
 473         if (bSig < aSig)
 474                 goto aBigger;
 475         if (aSig < bSig)
 476                 goto bBigger;
 477         return packFloat32(float_rounding_mode() == FPSCR_RM_ZERO, 0, 0);
 478       bExpBigger:
 479         if (bExp == 0xFF) {
 480                 return packFloat32(zSign ^ 1, 0xFF, 0);
 481         }
 482         if (aExp == 0) {
 483                 ++expDiff;
 484         } else {
 485                 aSig |= 0x40000000;
 486         }
 487         shift32RightJamming(aSig, -expDiff, &aSig);
 488         bSig |= 0x40000000;
 489       bBigger:
 490         zSig = bSig - aSig;
 491         zExp = bExp;
 492         zSign ^= 1;
 493         goto normalizeRoundAndPack;
 494       aExpBigger:
 495         if (aExp == 0xFF) {
 496                 return a;
 497         }
 498         if (bExp == 0) {
 499                 --expDiff;
 500         } else {
 501                 bSig |= 0x40000000;
 502         }
 503         shift32RightJamming(bSig, expDiff, &bSig);
 504         aSig |= 0x40000000;
 505       aBigger:
 506         zSig = aSig - bSig;
 507         zExp = aExp;
 508       normalizeRoundAndPack:
 509         --zExp;
 510         return normalizeRoundAndPackFloat32(zSign, zExp, zSig);
 511
 512 }
 513
 514 static float32 addFloat32Sigs(float32 a, float32 b, flag zSign)
 515 {
 516         int16 aExp, bExp, zExp;
 517         bits32 aSig, bSig, zSig;
 518         int16 expDiff;
 519
 520         aSig = extractFloat32Frac(a);
 521         aExp = extractFloat32Exp(a);
 522         bSig = extractFloat32Frac(b);
 523         bExp = extractFloat32Exp(b);
 524         expDiff = aExp - bExp;
 525         aSig <<= 6;
 526         bSig <<= 6;
 527         if (0 < expDiff) {
 528                 if (aExp == 0xFF) {
 529                         return a;
 530                 }
 531                 if (bExp == 0) {
 532                         --expDiff;
 533                 } else {
 534                         bSig |= 0x20000000;
 535                 }
 536                 shift32RightJamming(bSig, expDiff, &bSig);
 537                 zExp = aExp;
 538         } else if (expDiff < 0) {
 539                 if (bExp == 0xFF) {
 540                         return packFloat32(zSign, 0xFF, 0);
 541                 }
 542                 if (aExp == 0) {
 543                         ++expDiff;
 544                 } else {
 545                         aSig |= 0x20000000;
 546                 }
 547                 shift32RightJamming(aSig, -expDiff, &aSig);
 548                 zExp = bExp;
 549         } else {
 550                 if (aExp == 0xFF) {
 551                         return a;
 552                 }
 553                 if (aExp == 0)
 554                         return packFloat32(zSign, 0, (aSig + bSig) >> 6);
 555                 zSig = 0x40000000 + aSig + bSig;
 556                 zExp = aExp;
 557                 goto roundAndPack;
 558         }
 559         aSig |= 0x20000000;
 560         zSig = (aSig + bSig) << 1;
 561         --zExp;
 562         if ((sbits32) zSig < 0) {
 563                 zSig = aSig + bSig;
 564                 ++zExp;
 565         }
 566       roundAndPack:
 567         return roundAndPackFloat32(zSign, zExp, zSig);
 568
 569 }
 570
 571 float64 float64_sub(float64 a, float64 b)
 572 {
 573         flag aSign, bSign;
 574
 575         aSign = extractFloat64Sign(a);
 576         bSign = extractFloat64Sign(b);
 577         if (aSign == bSign) {
 578                 return subFloat64Sigs(a, b, aSign);
 579         } else {
 580                 return addFloat64Sigs(a, b, aSign);
 581         }
 582
 583 }
 584
 585 float32 float32_sub(float32 a, float32 b)
 586 {
 587         flag aSign, bSign;
 588
 589         aSign = extractFloat32Sign(a);
 590         bSign = extractFloat32Sign(b);
 591         if (aSign == bSign) {
 592                 return subFloat32Sigs(a, b, aSign);
 593         } else {
 594                 return addFloat32Sigs(a, b, aSign);
 595         }
 596
 597 }
 598
 599 float32 float32_add(float32 a, float32 b)
 600 {
 601         flag aSign, bSign;
 602
 603         aSign = extractFloat32Sign(a);
 604         bSign = extractFloat32Sign(b);
 605         if (aSign == bSign) {
 606                 return addFloat32Sigs(a, b, aSign);
 607         } else {
 608                 return subFloat32Sigs(a, b, aSign);
 609         }
 610
 611 }
 612
 613 float64 float64_add(float64 a, float64 b)
 614 {
 615         flag aSign, bSign;
 616
 617         aSign = extractFloat64Sign(a);
 618         bSign = extractFloat64Sign(b);
 619         if (aSign == bSign) {
 620                 return addFloat64Sigs(a, b, aSign);
 621         } else {
 622                 return subFloat64Sigs(a, b, aSign);
 623         }
 624 }
 625
 626 static void
 627 normalizeFloat64Subnormal(bits64 aSig, int16 * zExpPtr, bits64 * zSigPtr)
 628 {
 629         int8 shiftCount;
 630
 631         shiftCount = countLeadingZeros64(aSig) - 11;
 632         *zSigPtr = aSig << shiftCount;
 633         *zExpPtr = 1 - shiftCount;
 634 }
 635
 636 inline void add128(bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 * z0Ptr,
 637                    bits64 * z1Ptr)
 638 {
 639         bits64 z1;
 640
 641         z1 = a1 + b1;
 642         *z1Ptr = z1;
 643         *z0Ptr = a0 + b0 + (z1 < a1);
 644 }
 645
 646 inline void
 647 sub128(bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 * z0Ptr,
 648        bits64 * z1Ptr)
 649 {
 650         *z1Ptr = a1 - b1;
 651         *z0Ptr = a0 - b0 - (a1 < b1);
 652 }
 653
 654 static bits64 estimateDiv128To64(bits64 a0, bits64 a1, bits64 b)
 655 {
 656         bits64 b0, b1;
 657         bits64 rem0, rem1, term0, term1;
 658         bits64 z;
 659         if (b <= a0)
 660                 return LIT64(0xFFFFFFFFFFFFFFFF);
 661         b0 = b >> 32;
 662         z = (b0 << 32 <= a0) ? LIT64(0xFFFFFFFF00000000) : (a0 / b0) << 32;
 663         mul64To128(b, z, &term0, &term1);
 664         sub128(a0, a1, term0, term1, &rem0, &rem1);
 665         while (((sbits64) rem0) < 0) {
 666                 z -= LIT64(0x100000000);
 667                 b1 = b << 32;
 668                 add128(rem0, rem1, b0, b1, &rem0, &rem1);
 669         }
 670         rem0 = (rem0 << 32) | (rem1 >> 32);
 671         z |= (b0 << 32 <= rem0) ? 0xFFFFFFFF : rem0 / b0;
 672         return z;
 673 }
 674
 675 inline void mul64To128(bits64 a, bits64 b, bits64 * z0Ptr, bits64 * z1Ptr)
 676 {
 677         bits32 aHigh, aLow, bHigh, bLow;
 678         bits64 z0, zMiddleA, zMiddleB, z1;
 679
 680         aLow = a;
 681         aHigh = a >> 32;
 682         bLow = b;
 683         bHigh = b >> 32;
 684         z1 = ((bits64) aLow) * bLow;
 685         zMiddleA = ((bits64) aLow) * bHigh;
 686         zMiddleB = ((bits64) aHigh) * bLow;
 687         z0 = ((bits64) aHigh) * bHigh;
 688         zMiddleA += zMiddleB;
 689         z0 += (((bits64) (zMiddleA < zMiddleB)) << 32) + (zMiddleA >> 32);
 690         zMiddleA <<= 32;
 691         z1 += zMiddleA;
 692         z0 += (z1 < zMiddleA);
 693         *z1Ptr = z1;
 694         *z0Ptr = z0;
 695
 696 }
 697
 698 static void normalizeFloat32Subnormal(bits32 aSig, int16 * zExpPtr,
 699                                       bits32 * zSigPtr)
 700 {
 701         int8 shiftCount;
 702
 703         shiftCount = countLeadingZeros32(aSig) - 8;
 704         *zSigPtr = aSig << shiftCount;
 705         *zExpPtr = 1 - shiftCount;
 706
 707 }
 708
 709 float64 float64_div(float64 a, float64 b)
 710 {
 711         flag aSign, bSign, zSign;
 712         int16 aExp, bExp, zExp;
 713         bits64 aSig, bSig, zSig;
 714         bits64 rem0, rem1;
 715         bits64 term0, term1;
 716
 717         aSig = extractFloat64Frac(a);
 718         aExp = extractFloat64Exp(a);
 719         aSign = extractFloat64Sign(a);
 720         bSig = extractFloat64Frac(b);
 721         bExp = extractFloat64Exp(b);
 722         bSign = extractFloat64Sign(b);
 723         zSign = aSign ^ bSign;
 724         if (aExp == 0x7FF) {
 725                 if (bExp == 0x7FF) {
 726                 }
 727                 return packFloat64(zSign, 0x7FF, 0);
 728         }
 729         if (bExp == 0x7FF) {
 730                 return packFloat64(zSign, 0, 0);
 731         }
 732         if (bExp == 0) {
 733                 if (bSig == 0) {
 734                         if ((aExp | aSig) == 0) {
 735                                 float_raise(FPSCR_CAUSE_INVALID);
 736                         }
 737                         return packFloat64(zSign, 0x7FF, 0);
 738                 }
 739                 normalizeFloat64Subnormal(bSig, &bExp, &bSig);
 740         }
 741         if (aExp == 0) {
 742                 if (aSig == 0)
 743                         return packFloat64(zSign, 0, 0);
 744                 normalizeFloat64Subnormal(aSig, &aExp, &aSig);
 745         }
 746         zExp = aExp - bExp + 0x3FD;
 747         aSig = (aSig | LIT64(0x0010000000000000)) << 10;
 748         bSig = (bSig | LIT64(0x0010000000000000)) << 11;
 749         if (bSig <= (aSig + aSig)) {
 750                 aSig >>= 1;
 751                 ++zExp;
 752         }
 753         zSig = estimateDiv128To64(aSig, 0, bSig);
 754         if ((zSig & 0x1FF) <= 2) {
 755                 mul64To128(bSig, zSig, &term0, &term1);
 756                 sub128(aSig, 0, term0, term1, &rem0, &rem1);
 757                 while ((sbits64) rem0 < 0) {
 758                         --zSig;
 759                         add128(rem0, rem1, 0, bSig, &rem0, &rem1);
 760                 }
 761                 zSig |= (rem1 != 0);
 762         }
 763         return roundAndPackFloat64(zSign, zExp, zSig);
 764
 765 }
 766
 767 float32 float32_div(float32 a, float32 b)
 768 {
 769         flag aSign, bSign, zSign;
 770         int16 aExp, bExp, zExp;
 771         bits32 aSig, bSig, zSig;
 772
 773         aSig = extractFloat32Frac(a);
 774         aExp = extractFloat32Exp(a);
 775         aSign = extractFloat32Sign(a);
 776         bSig = extractFloat32Frac(b);
 777         bExp = extractFloat32Exp(b);
 778         bSign = extractFloat32Sign(b);
 779         zSign = aSign ^ bSign;
 780         if (aExp == 0xFF) {
 781                 if (bExp == 0xFF) {
 782                 }
 783                 return packFloat32(zSign, 0xFF, 0);
 784         }
 785         if (bExp == 0xFF) {
 786                 return packFloat32(zSign, 0, 0);
 787         }
 788         if (bExp == 0) {
 789                 if (bSig == 0) {
 790                         return packFloat32(zSign, 0xFF, 0);
 791                 }
 792                 normalizeFloat32Subnormal(bSig, &bExp, &bSig);
 793         }
 794         if (aExp == 0) {
 795                 if (aSig == 0)
 796                         return packFloat32(zSign, 0, 0);
 797                 normalizeFloat32Subnormal(aSig, &aExp, &aSig);
 798         }
 799         zExp = aExp - bExp + 0x7D;
 800         aSig = (aSig | 0x00800000) << 7;
 801         bSig = (bSig | 0x00800000) << 8;
 802         if (bSig <= (aSig + aSig)) {
 803                 aSig >>= 1;
 804                 ++zExp;
 805         }
 806         zSig = (((bits64) aSig) << 32) / bSig;
 807         if ((zSig & 0x3F) == 0) {
 808                 zSig |= (((bits64) bSig) * zSig != ((bits64) aSig) << 32);
 809         }
 810         return roundAndPackFloat32(zSign, zExp, zSig);
 811
 812 }
 813
 814 float32 float32_mul(float32 a, float32 b)
 815 {
 816         char aSign, bSign, zSign;
 817         int aExp, bExp, zExp;
 818         unsigned int aSig, bSig;
 819         unsigned long long zSig64;
 820         unsigned int zSig;
 821
 822         aSig = extractFloat32Frac(a);
 823         aExp = extractFloat32Exp(a);
 824         aSign = extractFloat32Sign(a);
 825         bSig = extractFloat32Frac(b);
 826         bExp = extractFloat32Exp(b);
 827         bSign = extractFloat32Sign(b);
 828         zSign = aSign ^ bSign;
 829         if (aExp == 0) {
 830                 if (aSig == 0)
 831                         return packFloat32(zSign, 0, 0);
 832                 normalizeFloat32Subnormal(aSig, &aExp, &aSig);
 833         }
 834         if (bExp == 0) {
 835                 if (bSig == 0)
 836                         return packFloat32(zSign, 0, 0);
 837                 normalizeFloat32Subnormal(bSig, &bExp, &bSig);
 838         }
 839         if ((bExp == 0xff && bSig == 0) || (aExp == 0xff && aSig == 0))
 840                 return roundAndPackFloat32(zSign, 0xff, 0);
 841
 842         zExp = aExp + bExp - 0x7F;
 843         aSig = (aSig | 0x00800000) << 7;
 844         bSig = (bSig | 0x00800000) << 8;
 845         shift64RightJamming(((unsigned long long)aSig) * bSig, 32, &zSig64);
 846         zSig = zSig64;
 847         if (0 <= (signed int)(zSig << 1)) {
 848                 zSig <<= 1;
 849                 --zExp;
 850         }
 851         return roundAndPackFloat32(zSign, zExp, zSig);
 852
 853 }
 854
 855 float64 float64_mul(float64 a, float64 b)
 856 {
 857         char aSign, bSign, zSign;
 858         int aExp, bExp, zExp;
 859         unsigned long long int aSig, bSig, zSig0, zSig1;
 860
 861         aSig = extractFloat64Frac(a);
 862         aExp = extractFloat64Exp(a);
 863         aSign = extractFloat64Sign(a);
 864         bSig = extractFloat64Frac(b);
 865         bExp = extractFloat64Exp(b);
 866         bSign = extractFloat64Sign(b);
 867         zSign = aSign ^ bSign;
 868
 869         if (aExp == 0) {
 870                 if (aSig == 0)
 871                         return packFloat64(zSign, 0, 0);
 872                 normalizeFloat64Subnormal(aSig, &aExp, &aSig);
 873         }
 874         if (bExp == 0) {
 875                 if (bSig == 0)
 876                         return packFloat64(zSign, 0, 0);
 877                 normalizeFloat64Subnormal(bSig, &bExp, &bSig);
 878         }
 879         if ((aExp == 0x7ff && aSig == 0) || (bExp == 0x7ff && bSig == 0))
 880                 return roundAndPackFloat64(zSign, 0x7ff, 0);
 881
 882         zExp = aExp + bExp - 0x3FF;
 883         aSig = (aSig | 0x0010000000000000LL) << 10;
 884         bSig = (bSig | 0x0010000000000000LL) << 11;
 885         mul64To128(aSig, bSig, &zSig0, &zSig1);
 886         zSig0 |= (zSig1 != 0);
 887         if (0 <= (signed long long int)(zSig0 << 1)) {
 888                 zSig0 <<= 1;
 889                 --zExp;
 890         }
 891         return roundAndPackFloat64(zSign, zExp, zSig0);
 892 }