gcc/floatlib.c

   1 /*
   2 ** libgcc support for software floating point.
   3 ** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
   4 ** Permission is granted to do *anything* you want with this file,
   5 ** commercial or otherwise, provided this message remains intact.  So there!
   6 ** I would appreciate receiving any updates/patches/changes that anyone
   7 ** makes, and am willing to be the repository for said changes (am I
   8 ** making a big mistake?).
   9
  10 Warning! Only single-precision is actually implemented.  This file
  11 won't really be much use until double-precision is supported.
  12
  13 However, once that is done, this file might eventually become a
  14 replacement for libgcc1.c.  It might also make possible
  15 cross-compilation for an IEEE target machine from a non-IEEE
  16 host such as a VAX.
  17
  18 If you'd like to work on completing this, please talk to rms@gnu.ai.mit.edu.
  19
  20
  21 **
  22 ** Pat Wood
  23 ** Pipeline Associates, Inc.
  24 ** pipeline!phw@motown.com or
  25 ** sun!pipeline!phw or
  26 ** uunet!motown!pipeline!phw
  27 **
  28 ** 05/01/91 -- V1.0 -- first release to gcc mailing lists
  29 ** 05/04/91 -- V1.1 -- added float and double prototypes and return values
  30 **                  -- fixed problems with adding and subtracting zero
  31 **                  -- fixed rounding in truncdfsf2
  32 **                  -- fixed SWAP define and tested on 386
  33 */
  34
  35 /*
  36 ** The following are routines that replace the libgcc soft floating point
  37 ** routines that are called automatically when -msoft-float is selected.
  38 ** The support single and double precision IEEE format, with provisions
  39 ** for byte-swapped machines (tested on 386).  Some of the double-precision
  40 ** routines work at full precision, but most of the hard ones simply punt
  41 ** and call the single precision routines, producing a loss of accuracy.
  42 ** long long support is not assumed or included.
  43 ** Overall accuracy is close to IEEE (actually 68882) for single-precision
  44 ** arithmetic.  I think there may still be a 1 in 1000 chance of a bit
  45 ** being rounded the wrong way during a multiply.  I'm not fussy enough to
  46 ** bother with it, but if anyone is, knock yourself out.
  47 **
  48 ** Efficiency has only been addressed where it was obvious that something
  49 ** would make a big difference.  Anyone who wants to do this right for
  50 ** best speed should go in and rewrite in assembler.
  51 **
  52 ** I have tested this only on a 68030 workstation and 386/ix integrated
  53 ** in with -msoft-float.
  54 */
  55
  56 /* the following deal with IEEE single-precision numbers */
  57 #define D_PHANTOM_BIT   0x00100000
  58 #define EXCESS          126
  59 #define SIGNBIT         0x80000000
  60 #define HIDDEN          (1 << 23)
  61 #define SIGN(fp)        ((fp) & SIGNBIT)
  62 #define EXP(fp)         (((fp) >> 23) & 0xFF)
  63 #define MANT(fp)        (((fp) & 0x7FFFFF) | HIDDEN)
  64 #define PACK(s,e,m)     ((s) | ((e) << 23) | (m))
  65
  66 /* the following deal with IEEE double-precision numbers */
  67 #define EXCESSD         1022
  68 #define HIDDEND         (1 << 20)
  69 #define EXPD(fp)        (((fp.l.upper) >> 20) & 0x7FF)
  70 #define SIGND(fp)       ((fp.l.upper) & SIGNBIT)
  71 #define MANTD(fp)       (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
  72                                 (fp.l.lower >> 22))
  73
  74 /* define SWAP for 386/960 reverse-byte-order brain-damaged CPUs */
  75 union double_long
  76   {
  77     double d;
  78 #ifdef SWAP
  79     struct {
  80       unsigned long lower;
  81       long upper;
  82     } l;
  83 #else
  84     struct {
  85       long upper;
  86       unsigned long lower;
  87     } l;
  88 #endif
  89   };
  90
  91 union float_long
  92   {
  93     float f;
  94     long l;
  95   };
  96
  97    struct _ieee {
  98 #ifdef SWAP
  99       unsigned mantissa2 : 32;
 100       unsigned mantissa1 : 20;
 101       unsigned exponent  : 11;
 102       unsigned sign      : 1;
 103 #else
 104       unsigned exponent  : 11;
 105       unsigned sign      : 1;
 106       unsigned mantissa2 : 32;
 107       unsigned mantissa1 : 20;
 108 #endif
 109    };
 110
 111    union _doubleu {
 112       double d;
 113       struct _ieee ieee;
 114 #ifdef SWAP
 115       struct {
 116          unsigned long lower;
 117          long upper;
 118       } l;
 119 #else
 120       struct {
 121          long upper;
 122          unsigned long lower;
 123       } l;
 124 #endif
 125    };
 126
 127 /* add two floats */
 128
 129 float
 130 __addsf3 (float a1, float a2)
 131 {
 132   register long mant1, mant2;
 133   register union float_long fl1, fl2;
 134   register int exp1, exp2;
 135   int sign = 0;
 136
 137   fl1.f = a1;
 138   fl2.f = a2;
 139
 140   /* check for zero args */
 141   if (!fl1.l)
 142     return (fl2.f);
 143   if (!fl2.l)
 144     return (fl1.f);
 145
 146   exp1 = EXP (fl1.l);
 147   exp2 = EXP (fl2.l);
 148
 149   if (exp1 > exp2 + 25)
 150     return (fl1.l);
 151   if (exp2 > exp1 + 25)
 152     return (fl2.l);
 153
 154   /* do everything in excess precision so's we can round later */
 155   mant1 = MANT (fl1.l) << 6;
 156   mant2 = MANT (fl2.l) << 6;
 157
 158   if (SIGN (fl1.l))
 159     mant1 = -mant1;
 160   if (SIGN (fl2.l))
 161     mant2 = -mant2;
 162
 163   if (exp1 > exp2)
 164     {
 165       mant2 >>= exp1 - exp2;
 166     }
 167   else
 168     {
 169       mant1 >>= exp2 - exp1;
 170       exp1 = exp2;
 171     }
 172   mant1 += mant2;
 173
 174   if (mant1 < 0)
 175     {
 176       mant1 = -mant1;
 177       sign = SIGNBIT;
 178     }
 179   else if (!mant1)
 180     return (0);
 181
 182   /* normalize up */
 183   while (!(mant1 & 0xE0000000))
 184     {
 185       mant1 <<= 1;
 186       exp1--;
 187     }
 188
 189   /* normalize down? */
 190   if (mant1 & (1 << 30))
 191     {
 192       mant1 >>= 1;
 193       exp1++;
 194     }
 195
 196   /* round to even */
 197   mant1 += (mant1 & 0x40) ? 0x20 : 0x1F;
 198
 199   /* normalize down? */
 200   if (mant1 & (1 << 30))
 201     {
 202       mant1 >>= 1;
 203       exp1++;
 204     }
 205
 206   /* lose extra precision */
 207   mant1 >>= 6;
 208
 209   /* turn off hidden bit */
 210   mant1 &= ~HIDDEN;
 211
 212   /* pack up and go home */
 213   fl1.l = PACK (sign, exp1, mant1);
 214   return (fl1.f);
 215 }
 216
 217 /* subtract two floats */
 218
 219 float
 220 __subsf3 (float a1, float a2)
 221 {
 222   register union float_long fl1, fl2;
 223
 224   fl1.f = a1;
 225   fl2.f = a2;
 226
 227   /* check for zero args */
 228   if (!fl2.l)
 229     return (fl1.f);
 230   if (!fl1.l)
 231     return (-fl2.f);
 232
 233   /* twiddle sign bit and add */
 234   fl2.l ^= SIGNBIT;
 235   return __addsf3 (a1, fl2.f);
 236 }
 237
 238 /* compare two floats */
 239
 240 long
 241 __cmpsf2 (float a1, float a2)
 242 {
 243   register union float_long fl1, fl2;
 244
 245   fl1.f = a1;
 246   fl2.f = a2;
 247
 248   if (SIGN (fl1.l) && SIGN (fl2.l))
 249     {
 250       fl1.l ^= SIGNBIT;
 251       fl2.l ^= SIGNBIT;
 252     }
 253   if (fl1.l < fl2.l)
 254     return (-1);
 255   if (fl1.l > fl2.l)
 256     return (1);
 257   return (0);
 258 }
 259
 260 /* multiply two floats */
 261
 262 float
 263 __mulsf3 (float a1, float a2)
 264 {
 265   register union float_long fl1, fl2;
 266   register unsigned long result;
 267   register int exp;
 268   int sign;
 269
 270   fl1.f = a1;
 271   fl2.f = a2;
 272
 273   if (!fl1.l || !fl2.l)
 274     return (0);
 275
 276   /* compute sign and exponent */
 277   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
 278   exp = EXP (fl1.l) - EXCESS;
 279   exp += EXP (fl2.l);
 280
 281   fl1.l = MANT (fl1.l);
 282   fl2.l = MANT (fl2.l);
 283
 284   /* the multiply is done as one 16x16 multiply and two 16x8 multiples */
 285   result = (fl1.l >> 8) * (fl2.l >> 8);
 286   result += ((fl1.l & 0xFF) * (fl2.l >> 8)) >> 8;
 287   result += ((fl2.l & 0xFF) * (fl1.l >> 8)) >> 8;
 288
 289   if (result & 0x80000000)
 290     {
 291       /* round */
 292       result += 0x80;
 293       result >>= 8;
 294     }
 295   else
 296     {
 297       /* round */
 298       result += 0x40;
 299       result >>= 7;
 300       exp--;
 301     }
 302
 303   result &= ~HIDDEN;
 304
 305   /* pack up and go home */
 306   fl1.l = PACK (sign, exp, result);
 307   return (fl1.f);
 308 }
 309
 310 /* divide two floats */
 311
 312 float
 313 __divsf3 (float a1, float a2)
 314 {
 315   register union float_long fl1, fl2;
 316   register int result;
 317   register int mask;
 318   register int exp, sign;
 319
 320   fl1.f = a1;
 321   fl2.f = a2;
 322
 323   /* subtract exponents */
 324   exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS;
 325
 326   /* compute sign */
 327   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
 328
 329   /* divide by zero??? */
 330   if (!fl2.l)
 331     /* return NaN or -NaN */
 332     return (sign ? 0xFFFFFFFF : 0x7FFFFFFF);
 333
 334   /* numerator zero??? */
 335   if (!fl1.l)
 336     return (0);
 337
 338   /* now get mantissas */
 339   fl1.l = MANT (fl1.l);
 340   fl2.l = MANT (fl2.l);
 341
 342   /* this assures we have 25 bits of precision in the end */
 343   if (fl1.l < fl2.l)
 344     {
 345       fl1.l <<= 1;
 346       exp--;
 347     }
 348
 349   /* now we perform repeated subtraction of fl2.l from fl1.l */
 350   mask = 0x1000000;
 351   result = 0;
 352   while (mask)
 353     {
 354       if (fl1.l >= fl2.l)
 355         {
 356           result |= mask;
 357           fl1.l -= fl2.l;
 358         }
 359       fl1.l <<= 1;
 360       mask >>= 1;
 361     }
 362
 363   /* round */
 364   result += 1;
 365
 366   /* normalize down */
 367   exp++;
 368   result >>= 1;
 369
 370   result &= ~HIDDEN;
 371
 372   /* pack up and go home */
 373   fl1.l = PACK (sign, exp, result);
 374   return (fl1.f);
 375 }
 376
 377 /* convert int to double */
 378
 379 double
 380 __floatsidf (register long a1)
 381 {
 382   register int sign = 0, exp = 31 + EXCESSD;
 383   union double_long dl;
 384
 385   if (!a1)
 386     {
 387       dl.l.upper = dl.l.lower = 0;
 388       return (dl.d);
 389     }
 390
 391   if (a1 < 0)
 392     {
 393       sign = SIGNBIT;
 394       a1 = -a1;
 395     }
 396
 397   while (a1 < 0x1000000)
 398     {
 399       a1 <<= 4;
 400       exp -= 4;
 401     }
 402
 403   while (a1 < 0x40000000)
 404     {
 405       a1 <<= 1;
 406       exp--;
 407     }
 408
 409   /* pack up and go home */
 410   dl.l.upper = sign;
 411   dl.l.upper |= exp << 20;
 412   dl.l.upper |= (a1 >> 10) & ~HIDDEND;
 413   dl.l.lower = a1 << 22;
 414
 415   return (dl.d);
 416 }
 417
 418 /* negate a float */
 419
 420 float
 421 __negsf2 (float a1)
 422 {
 423   register union float_long fl1;
 424
 425   fl1.f = a1;
 426   if (!fl1.l)
 427     return (0);
 428
 429   fl1.l ^= SIGNBIT;
 430   return (fl1.f);
 431 }
 432
 433 /* negate a double */
 434
 435 double
 436 __negdf2 (double a1)
 437 {
 438   register union double_long dl1;
 439
 440   dl1.d = a1;
 441
 442   if (!dl1.l.upper && !dl1.l.lower)
 443       return (dl1.d);
 444
 445   dl1.l.upper ^= SIGNBIT;
 446   return (dl1.d);
 447 }
 448
 449 /* convert float to double */
 450
 451 double
 452 __extendsfdf2 (float a1)
 453 {
 454   register union float_long fl1;
 455   register union double_long dl;
 456   register int exp;
 457
 458   fl1.f = a1;
 459
 460   if (!fl1.l)
 461     {
 462       dl.l.upper = dl.l.lower = 0;
 463       return (dl.d);
 464     }
 465
 466   dl.l.upper = SIGN (fl1.l);
 467   exp = EXP (fl1.l) - EXCESS + EXCESSD;
 468   dl.l.upper |= exp << 20;
 469   dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
 470   dl.l.lower = MANT (fl1.l) << 29;
 471
 472   return (dl.d);
 473 }
 474
 475 /* convert double to float */
 476
 477 float
 478 __truncdfsf2 (double a1)
 479 {
 480   register int exp;
 481   register long mant;
 482   register union float_long fl;
 483   register union double_long dl1;
 484
 485   dl1.d = a1;
 486
 487   if (!dl1.l.upper && !dl1.l.lower)
 488     return (0);
 489
 490   exp = EXPD (dl1) - EXCESSD + EXCESS;
 491
 492   /* shift double mantissa 6 bits so we can round */
 493   mant = MANTD (dl1) >> 6;
 494
 495   /* now round and shift down */
 496   mant += 1;
 497   mant >>= 1;
 498
 499   /* did the round overflow? */
 500   if (mant & 0xFF000000)
 501     {
 502       mant >>= 1;
 503       exp++;
 504     }
 505
 506   mant &= ~HIDDEN;
 507
 508   /* pack up and go home */
 509   fl.l = PACK (SIGND (dl1), exp, mant);
 510   return (fl.f);
 511 }
 512
 513 /* compare two doubles */
 514
 515 long
 516 __cmpdf2 (double a1, double a2)
 517 {
 518   register union double_long dl1, dl2;
 519
 520   dl1.d = a1;
 521   dl2.d = a2;
 522
 523   if (SIGND (dl1) && SIGND (dl2))
 524     {
 525       dl1.l.upper ^= SIGNBIT;
 526       dl2.l.upper ^= SIGNBIT;
 527     }
 528   if (dl1.l.upper < dl2.l.upper)
 529     return (-1);
 530   if (dl1.l.upper > dl2.l.upper)
 531     return (1);
 532   if (dl1.l.lower < dl2.l.lower)
 533     return (-1);
 534   if (dl1.l.lower > dl2.l.lower)
 535     return (1);
 536   return (0);
 537 }
 538
 539 /* convert double to int */
 540
 541 long
 542 __fixdfsi (double a1)
 543 {
 544   register union double_long dl1;
 545   register int exp;
 546   register long l;
 547
 548   dl1.d = a1;
 549
 550   if (!dl1.l.upper && !dl1.l.lower)
 551     return (0);
 552
 553   exp = EXPD (dl1) - EXCESSD - 31;
 554   l = MANTD (dl1);
 555
 556   if (exp > 0)
 557     return (0x7FFFFFFF | SIGND (dl1)); /* largest integer */
 558
 559   /* shift down until exp = 0 or l = 0 */
 560   if (exp < 0 && exp > -32 && l)
 561     l >>= -exp;
 562   else
 563     return (0);
 564
 565   return (SIGND (dl1) ? -l : l);
 566 }
 567
 568 /* convert double to unsigned int */
 569
 570 unsigned
 571 long __fixunsdfsi (double a1)
 572 {
 573   register union double_long dl1;
 574   register int exp;
 575   register unsigned long l;
 576
 577   dl1.d = a1;
 578
 579   if (!dl1.l.upper && !dl1.l.lower)
 580     return (0);
 581
 582   exp = EXPD (dl1) - EXCESSD - 32;
 583   l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21));
 584
 585   if (exp > 0)
 586     return (0xFFFFFFFF);        /* largest integer */
 587
 588   /* shift down until exp = 0 or l = 0 */
 589   if (exp < 0 && exp > -32 && l)
 590     l >>= -exp;
 591   else
 592     return (0);
 593
 594   return (l);
 595 }
 596
 597 /* For now, the hard double-precision routines simply
 598    punt and do it in single */
 599 /* addtwo doubles */
 600
 601 double
 602 __adddf3 (double a1, double a2)
 603 {
 604   return ((float) a1 + (float) a2);
 605 }
 606
 607 /* subtract two doubles */
 608
 609 double
 610 __subdf3 (double a1, double a2)
 611 {
 612   return ((float) a1 - (float) a2);
 613 }
 614
 615 /* multiply two doubles */
 616
 617 double
 618 __muldf3 (double a1, double a2)
 619 {
 620   return ((float) a1 * (float) a2);
 621 }
 622
 623 /*
 624  *
 625  * Name:   Barrett Richardson
 626  * E-mail: barrett@iglou.com
 627  * When:   Thu Dec 15 10:31:11 EST 1994
 628  *
 629  *    callable function:
 630  *
 631  *       double __divdf3(double a1, double a2);
 632  *
 633  *       Does software divide of a1 / a2.
 634  *
 635  *       Based largely on __divsf3() in floatlib.c in the gcc
 636  *       distribution.
 637  *
 638  *    Purpose: To be used in conjunction with the -msoft-float
 639  *             option of gcc. You should be able to tack it to the
 640  *             end of floatlib.c included in the gcc distribution,
 641  *             and delete the __divdf3() already there which just
 642  *             calls the single precision function (or may just
 643  *             use the floating point processor with some configurations).
 644  *
 645  *   You may use this code for whatever your heart desires.
 646  */
 647
 648
 649
 650
 651 /*
 652  * Compare the the mantissas of two doubles.
 653  * Each mantissa is in two longs.
 654  *
 655  *   return      1   if x1's mantissa is greater than x2's
 656  *              -1   if x1's mantissa is less than x2's
 657  *               0   if the two mantissa's are equal.
 658  *
 659  *   The Mantissas won't fit into a 4 byte word, so they are
 660  *   broken up into two parts.
 661  *
 662  *   This function is used internally by __divdf3()
 663  */
 664
 665 int
 666 __dcmp (long x1m1, long x1m2, long x2m1, long x2m2)
 667 {
 668    if (x1m1 > x2m1)
 669       return 1;
 670
 671    if (x1m1 < x2m1)
 672       return -1;
 673
 674  /*  If the first word in the two mantissas were equal check the second word */
 675
 676    if (x1m2 > x2m2)
 677       return 1;
 678
 679    if (x1m2 < x2m2)
 680       return -1;
 681
 682    return 0;
 683 }
 684
 685
 686 /* divide two doubles */
 687
 688 double
 689 __divdf3 (double a1, double a2)
 690 {
 691
 692    int  sign,
 693         exponent,
 694         bit_bucket;
 695
 696    register unsigned long mantissa1,
 697                           mantissa2,
 698                           x1m1,
 699                           x1m2,
 700                           x2m1,
 701                           x2m2,
 702                           mask;
 703
 704    union _doubleu x1,
 705                   x2,
 706                   result;
 707
 708
 709    x1.d = a1;
 710    x2.d = a2;
 711
 712    exponent = x1.ieee.exponent - x2.ieee.exponent + EXCESSD;
 713
 714    sign = x1.ieee.sign ^ x2.ieee.sign;
 715
 716    x2.ieee.sign = 0;  /* don't want the sign bit to affect any zero */
 717                       /* comparisons when checking for zero divide  */
 718
 719    if (!x2.l.lower && !x2.l.upper) { /* check for zero divide */
 720       result.l.lower = 0x0;
 721       if (sign)
 722          result.l.upper = 0xFFF00000;   /* negative infinity */
 723       else
 724          result.l.upper = 0x7FF00000;   /* positive infinity */
 725       return result.d;
 726    }
 727
 728    if (!x1.l.upper && !x1.l.lower)  /* check for 0.0 numerator */
 729       return (0.0);
 730
 731    x1m1 = x1.ieee.mantissa1 | D_PHANTOM_BIT;  /* turn on phantom bit */
 732    x1m2 = x1.ieee.mantissa2;
 733
 734    x2m1 = x2.ieee.mantissa1 | D_PHANTOM_BIT;  /* turn on phantom bit */
 735    x2m2 = x2.ieee.mantissa2;
 736
 737    if (__dcmp(x1m1,x1m2,x2m1,x2m2) < 0) {
 738
 739    /* if x1's mantissa is less than x2's shift it left one and decrement */
 740    /* the exponent to accomodate the change in the mantissa              */
 741
 742       x1m1 <<= 1;               /*                          */
 743       bit_bucket = x1m2 >> 31;  /*  Shift mantissa left one */
 744       x1m1 |= bit_bucket;       /*                          */
 745       x1m2 <<= 1;               /*                          */
 746
 747       exponent--;
 748    }
 749
 750
 751    mantissa1 = 0;
 752    mantissa2 = 0;
 753
 754
 755   /* Get the first part of the results mantissa using successive */
 756   /* subtraction.                                                */
 757
 758    mask = 0x00200000;
 759    while (mask) {
 760
 761       if (__dcmp(x1m1,x1m2,x2m1,x2m2) >= 0) {
 762
 763      /* subtract x2's mantissa from x1's */
 764
 765          mantissa1 |= mask;   /* turn on a bit in the result */
 766
 767          if (x2m2 > x1m2)
 768             x1m1--;
 769          x1m2 -= x2m2;
 770          x1m1 -= x2m1;
 771       }
 772
 773       x1m1 <<= 1;               /*                          */
 774       bit_bucket = x1m2 >> 31;  /*  Shift mantissa left one */
 775       x1m1 |= bit_bucket;       /*                          */
 776       x1m2 <<= 1;               /*                          */
 777
 778       mask >>= 1;
 779    }
 780
 781   /* Get the second part of the results mantissa using successive */
 782   /* subtraction.                                                 */
 783
 784    mask = 0x80000000;
 785    while (mask) {
 786
 787       if (__dcmp(x1m1,x1m2,x2m1,x2m2) >= 0) {
 788
 789      /* subtract x2's mantissa from x1's */
 790
 791          mantissa2 |= mask;   /* turn on a bit in the result */
 792
 793          if (x2m2 > x1m2)
 794             x1m1--;
 795          x1m2 -= x2m2;
 796          x1m1 -= x2m1;
 797       }
 798       x1m1 <<= 1;               /*                          */
 799       bit_bucket = x1m2 >> 31;  /*  Shift mantissa left one */
 800       x1m1 |= bit_bucket;       /*                          */
 801       x1m2 <<= 1;               /*                          */
 802
 803       mask >>= 1;
 804    }
 805
 806   /* round up by adding 1 to mantissa */
 807
 808    if (mantissa2 == 0xFFFFFFFF) {  /* check for over flow */
 809
 810    /* spill if overflow */
 811
 812       mantissa2 = 0;
 813       mantissa1++;
 814    }
 815    else
 816       mantissa2++;
 817
 818    exponent++;   /* increment exponent (mantissa must be shifted right */
 819                  /* also)                                              */
 820
 821  /* shift mantissa right one and assume a phantom bit (which really gives */
 822  /* 53 bits of precision in the mantissa)                                 */
 823
 824    mantissa2 >>= 1;
 825    bit_bucket = mantissa1 & 1;
 826    mantissa2 |= (bit_bucket << 31);
 827    mantissa1 >>= 1;
 828
 829   /* put all the info into the result */
 830
 831    result.ieee.exponent  = exponent;
 832    result.ieee.sign      = sign;
 833    result.ieee.mantissa1 = mantissa1;
 834    result.ieee.mantissa2 = mantissa2;
 835
 836
 837    return result.d;
 838 }