sim/rx/fpu.c

   1 /* fpu.c --- FPU emulator for stand-alone RX simulator.
   2
   3 Copyright (C) 2008-2024 Free Software Foundation, Inc.
   4 Contributed by Red Hat, Inc.
   5
   6 This file is part of the GNU simulators.
   7
   8 This program is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or
  11 (at your option) any later version.
  12
  13 This program is distributed in the hope that it will be useful,
  14 but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 GNU General Public License for more details.
  17
  18 You should have received a copy of the GNU General Public License
  19 along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
  20
  21 /* This must come before any other includes.  */
  22 #include "defs.h"
  23
  24 #include <stdio.h>
  25 #include <stdlib.h>
  26
  27 #include "cpu.h"
  28 #include "fpu.h"
  29
  30 /* FPU encodings are as follows:
  31
  32    S EXPONENT MANTISSA
  33    1 12345678 12345678901234567890123
  34
  35    0 00000000 00000000000000000000000   +0
  36    1 00000000 00000000000000000000000   -0
  37
  38    X 00000000 00000000000000000000001   Denormals
  39    X 00000000 11111111111111111111111
  40
  41    X 00000001 XXXXXXXXXXXXXXXXXXXXXXX   Normals
  42    X 11111110 XXXXXXXXXXXXXXXXXXXXXXX
  43
  44    0 11111111 00000000000000000000000   +Inf
  45    1 11111111 00000000000000000000000   -Inf
  46
  47    X 11111111 0XXXXXXXXXXXXXXXXXXXXXX   SNaN (X != 0)
  48    X 11111111 1XXXXXXXXXXXXXXXXXXXXXX   QNaN (X != 0)
  49
  50 */
  51
  52 #define trace 0
  53 #define tprintf if (trace) printf
  54
  55 /* Some magic numbers.  */
  56 #define PLUS_MAX   0x7f7fffffUL
  57 #define MINUS_MAX  0xff7fffffUL
  58 #define PLUS_INF   0x7f800000UL
  59 #define MINUS_INF  0xff800000UL
  60 #define PLUS_ZERO  0x00000000UL
  61 #define MINUS_ZERO 0x80000000UL
  62
  63 #define FP_RAISE(e) fp_raise(FPSWBITS_C##e)
  64 static void
  65 fp_raise (int mask)
  66 {
  67   regs.r_fpsw |= mask;
  68   if (mask != FPSWBITS_CE)
  69     {
  70       if (regs.r_fpsw & (mask << FPSW_CESH))
  71         regs.r_fpsw |= (mask << FPSW_CFSH);
  72       if (regs.r_fpsw & FPSWBITS_FMASK)
  73         regs.r_fpsw |= FPSWBITS_FSUM;
  74       else
  75         regs.r_fpsw &= ~FPSWBITS_FSUM;
  76     }
  77 }
  78
  79 /* We classify all numbers as one of these.  They correspond to the
  80    rows/colums in the exception tables.  */
  81 typedef enum {
  82   FP_NORMAL,
  83   FP_PZERO,
  84   FP_NZERO,
  85   FP_PINFINITY,
  86   FP_NINFINITY,
  87   FP_DENORMAL,
  88   FP_QNAN,
  89   FP_SNAN
  90 } FP_Type;
  91
  92 #if defined DEBUG0
  93 static const char *fpt_names[] = {
  94   "Normal", "+0", "-0", "+Inf", "-Inf", "Denormal", "QNaN", "SNaN"
  95 };
  96 #endif
  97
  98 #define EXP_BIAS  127
  99 #define EXP_ZERO -127
 100 #define EXP_INF   128
 101
 102 #define MANT_BIAS 0x00080000UL
 103
 104 typedef struct {
 105   int exp;
 106   unsigned int mant; /* 24 bits */
 107   char type;
 108   char sign;
 109   fp_t orig_value;
 110 } FP_Parts;
 111
 112 static void
 113 fp_explode (fp_t f, FP_Parts *p)
 114 {
 115   int exp, mant, sign;
 116
 117   exp = ((f & 0x7f800000UL) >> 23);
 118   mant = f & 0x007fffffUL;
 119   sign = f & 0x80000000UL;
 120   /*printf("explode: %08x %x %2x %6x\n", f, sign, exp, mant);*/
 121
 122   p->sign = sign ? -1 : 1;
 123   p->exp = exp - EXP_BIAS;
 124   p->orig_value = f;
 125   p->mant = mant | 0x00800000UL;
 126
 127   if (p->exp == EXP_ZERO)
 128     {
 129       if (regs.r_fpsw & FPSWBITS_DN)
 130         mant = 0;
 131       if (mant)
 132         p->type = FP_DENORMAL;
 133       else
 134         {
 135           p->mant = 0;
 136           p->type = sign ? FP_NZERO : FP_PZERO;
 137         }
 138     }
 139   else if (p->exp == EXP_INF)
 140     {
 141       if (mant == 0)
 142         p->type = sign ? FP_NINFINITY : FP_PINFINITY;
 143       else if (mant & 0x00400000UL)
 144         p->type = FP_QNAN;
 145       else
 146         p->type = FP_SNAN;
 147     }
 148   else
 149     p->type = FP_NORMAL;
 150 }
 151
 152 static fp_t
 153 fp_implode (FP_Parts *p)
 154 {
 155   int exp, mant;
 156
 157   exp = p->exp + EXP_BIAS;
 158   mant = p->mant;
 159   /*printf("implode: exp %d mant 0x%x\n", exp, mant);*/
 160   if (p->type == FP_NORMAL)
 161     {
 162       while (mant
 163              && exp > 0
 164              && mant < 0x00800000UL)
 165         {
 166           mant <<= 1;
 167           exp --;
 168         }
 169       while (mant > 0x00ffffffUL)
 170         {
 171           mant >>= 1;
 172           exp ++;
 173         }
 174       if (exp < 0)
 175         {
 176           /* underflow */
 177           exp = 0;
 178           mant = 0;
 179           FP_RAISE (E);
 180         }
 181       if (exp >= 255)
 182         {
 183           /* overflow */
 184           exp = 255;
 185           mant = 0;
 186           FP_RAISE (O);
 187         }
 188     }
 189   mant &= 0x007fffffUL;
 190   exp &= 0xff;
 191   mant |= exp << 23;
 192   if (p->sign < 0)
 193     mant |= 0x80000000UL;
 194
 195   return mant;
 196 }
 197
 198 typedef union {
 199   unsigned long long ll;
 200   double d;
 201 } U_d_ll;
 202
 203 static int checked_format = 0;
 204
 205 /* We assume a double format like this:
 206    S[1] E[11] M[52]
 207 */
 208
 209 static double
 210 fp_to_double (FP_Parts *p)
 211 {
 212   U_d_ll u;
 213
 214   if (!checked_format)
 215     {
 216       u.d = 1.5;
 217       if (u.ll != 0x3ff8000000000000ULL)
 218         abort ();
 219       u.d = -225;
 220       if (u.ll != 0xc06c200000000000ULL)
 221         abort ();
 222       u.d = 10.1;
 223       if (u.ll != 0x4024333333333333ULL)
 224         abort ();
 225       checked_format = 1;
 226     }
 227
 228   u.ll = 0;
 229   if (p->sign < 0)
 230     u.ll |= (1ULL << 63);
 231   /* Make sure a zero encoding stays a zero.  */
 232   if (p->exp != -EXP_BIAS)
 233     u.ll |= ((unsigned long long)p->exp + 1023ULL) << 52;
 234   u.ll |= (unsigned long long) (p->mant & 0x007fffffUL) << (52 - 23);
 235   return u.d;
 236 }
 237
 238 static void
 239 double_to_fp (double d, FP_Parts *p)
 240 {
 241   int exp;
 242   U_d_ll u;
 243   int sign;
 244
 245   u.d = d;
 246
 247   sign = (u.ll & 0x8000000000000000ULL) ? 1 : 0;
 248   exp = u.ll >> 52;
 249   exp = (exp & 0x7ff);
 250
 251   if (exp == 0)
 252     {
 253       /* A generated denormal should show up as an underflow, not
 254          here.  */
 255       if (sign)
 256         fp_explode (MINUS_ZERO, p);
 257       else
 258         fp_explode (PLUS_ZERO, p);
 259       return;
 260     }
 261
 262   exp = exp - 1023;
 263   if ((exp + EXP_BIAS) > 254)
 264     {
 265       FP_RAISE (O);
 266       switch (regs.r_fpsw & FPSWBITS_RM)
 267         {
 268         case FPRM_NEAREST:
 269           if (sign)
 270             fp_explode (MINUS_INF, p);
 271           else
 272             fp_explode (PLUS_INF, p);
 273           break;
 274         case FPRM_ZERO:
 275           if (sign)
 276             fp_explode (MINUS_MAX, p);
 277           else
 278             fp_explode (PLUS_MAX, p);
 279           break;
 280         case FPRM_PINF:
 281           if (sign)
 282             fp_explode (MINUS_MAX, p);
 283           else
 284             fp_explode (PLUS_INF, p);
 285           break;
 286         case FPRM_NINF:
 287           if (sign)
 288             fp_explode (MINUS_INF, p);
 289           else
 290             fp_explode (PLUS_MAX, p);
 291           break;
 292         }
 293       return;
 294     }
 295   if ((exp + EXP_BIAS) < 1)
 296     {
 297       if (sign)
 298         fp_explode (MINUS_ZERO, p);
 299       else
 300         fp_explode (PLUS_ZERO, p);
 301       FP_RAISE (U);
 302     }
 303
 304   p->sign = sign ? -1 : 1;
 305   p->exp = exp;
 306   p->mant = u.ll >> (52-23) & 0x007fffffUL;
 307   p->mant |= 0x00800000UL;
 308   p->type = FP_NORMAL;
 309
 310   if (u.ll & 0x1fffffffULL)
 311     {
 312       switch (regs.r_fpsw & FPSWBITS_RM)
 313         {
 314         case FPRM_NEAREST:
 315           if (u.ll & 0x10000000ULL)
 316             p->mant ++;
 317           break;
 318         case FPRM_ZERO:
 319           break;
 320         case FPRM_PINF:
 321           if (sign == 1)
 322             p->mant ++;
 323           break;
 324         case FPRM_NINF:
 325           if (sign == -1)
 326             p->mant ++;
 327           break;
 328         }
 329       FP_RAISE (X);
 330     }
 331
 332 }
 333
 334 typedef enum {
 335   eNR,          /* Use the normal result.  */
 336   ePZ, eNZ,     /* +- zero */
 337   eSZ,          /* signed zero - XOR signs of ops together.  */
 338   eRZ,          /* +- zero depending on rounding mode.  */
 339   ePI, eNI,     /* +- Infinity */
 340   eSI,          /* signed infinity - XOR signs of ops together.  */
 341   eQN, eSN,     /* Quiet/Signalling NANs */
 342   eIn,          /* Invalid.  */
 343   eUn,          /* Unimplemented.  */
 344   eDZ,          /* Divide-by-zero.  */
 345   eLT,          /* less than */
 346   eGT,          /* greater than */
 347   eEQ,          /* equal to */
 348 } FP_ExceptionCases;
 349
 350 #if defined DEBUG0
 351 static const char *ex_names[] = {
 352   "NR", "PZ", "NZ", "SZ", "RZ", "PI", "NI", "SI", "QN", "SN", "IN", "Un", "DZ", "LT", "GT", "EQ"
 353 };
 354 #endif
 355
 356 /* This checks for all exceptional cases (not all FP exceptions) and
 357    returns TRUE if it is providing the result in *c.  If it returns
 358    FALSE, the caller should do the "normal" operation.  */
 359 static int
 360 check_exceptions (FP_Parts *a, FP_Parts *b, fp_t *c,
 361                   FP_ExceptionCases ex_tab[5][5],
 362                   FP_ExceptionCases *case_ret)
 363 {
 364   FP_ExceptionCases fpec;
 365
 366   if (a->type == FP_SNAN
 367       || b->type == FP_SNAN)
 368     fpec = eIn;
 369   else if (a->type == FP_QNAN
 370            || b->type == FP_QNAN)
 371     fpec = eQN;
 372   else if (a->type == FP_DENORMAL
 373            || b->type == FP_DENORMAL)
 374     fpec = eUn;
 375   else
 376     fpec = ex_tab[(int)(a->type)][(int)(b->type)];
 377
 378   /*printf("%s %s -> %s\n", fpt_names[(int)(a->type)], fpt_names[(int)(b->type)], ex_names[(int)(fpec)]);*/
 379
 380   if (case_ret)
 381     *case_ret = fpec;
 382
 383   switch (fpec)
 384     {
 385     case eNR:   /* Use the normal result.  */
 386       return 0;
 387
 388     case ePZ:   /* + zero */
 389       *c = 0x00000000;
 390       return 1;
 391
 392     case eNZ:   /* - zero */
 393       *c = 0x80000000;
 394       return 1;
 395
 396     case eSZ:   /* signed zero */
 397       *c = (a->sign == b->sign) ? PLUS_ZERO : MINUS_ZERO;
 398       return 1;
 399
 400     case eRZ:   /* +- zero depending on rounding mode.  */
 401       if ((regs.r_fpsw & FPSWBITS_RM) == FPRM_NINF)
 402         *c = 0x80000000;
 403       else
 404         *c = 0x00000000;
 405       return 1;
 406
 407     case ePI:   /* + Infinity */
 408       *c = 0x7F800000;
 409       return 1;
 410
 411     case eNI:   /* - Infinity */
 412       *c = 0xFF800000;
 413       return 1;
 414
 415     case eSI:   /* sign Infinity */
 416       *c = (a->sign == b->sign) ? PLUS_INF : MINUS_INF;
 417       return 1;
 418
 419     case eQN:   /* Quiet NANs */
 420       if (a->type == FP_QNAN)
 421         *c = a->orig_value;
 422       else
 423         *c = b->orig_value;
 424       return 1;
 425
 426     case eSN:   /* Signalling NANs */
 427       if (a->type == FP_SNAN)
 428         *c = a->orig_value;
 429       else
 430         *c = b->orig_value;
 431       FP_RAISE (V);
 432       return 1;
 433
 434     case eIn:   /* Invalid.  */
 435       FP_RAISE (V);
 436       if (a->type == FP_SNAN)
 437         *c = a->orig_value | 0x00400000;
 438       else if  (b->type == FP_SNAN)
 439         *c = b->orig_value | 0x00400000;
 440       else
 441         *c = 0x7fc00000;
 442       return 1;
 443
 444     case eUn:   /* Unimplemented.  */
 445       FP_RAISE (E);
 446       return 1;
 447
 448     case eDZ:   /* Division-by-zero.  */
 449       *c = (a->sign == b->sign) ? PLUS_INF : MINUS_INF;
 450       FP_RAISE (Z);
 451       return 1;
 452
 453     default:
 454       return 0;
 455     }
 456 }
 457
 458 #define CHECK_EXCEPTIONS(FPPa, FPPb, fpc, ex_tab) \
 459   if (check_exceptions (&FPPa, &FPPb, &fpc, ex_tab, 0)) \
 460     return fpc;
 461
 462 /* For each operation, we have two tables of how nonnormal cases are
 463    handled.  The DN=0 case is first, followed by the DN=1 case, with
 464    each table using the following layout: */
 465
 466 static FP_ExceptionCases ex_add_tab[5][5] = {
 467   /* N   +0   -0   +In  -In */
 468   { eNR, eNR, eNR, ePI, eNI }, /* Normal */
 469   { eNR, ePZ, eRZ, ePI, eNI }, /* +0   */
 470   { eNR, eRZ, eNZ, ePI, eNI }, /* -0   */
 471   { ePI, ePI, ePI, ePI, eIn }, /* +Inf */
 472   { eNI, eNI, eNI, eIn, eNI }, /* -Inf */
 473 };
 474
 475 fp_t
 476 rxfp_add (fp_t fa, fp_t fb)
 477 {
 478   FP_Parts a, b, c;
 479   fp_t rv;
 480   double da, db;
 481
 482   fp_explode (fa, &a);
 483   fp_explode (fb, &b);
 484   CHECK_EXCEPTIONS (a, b, rv, ex_add_tab);
 485
 486   da = fp_to_double (&a);
 487   db = fp_to_double (&b);
 488   tprintf("%g + %g = %g\n", da, db, da+db);
 489
 490   double_to_fp (da+db, &c);
 491   rv = fp_implode (&c);
 492   return rv;
 493 }
 494
 495 static FP_ExceptionCases ex_sub_tab[5][5] = {
 496   /* N   +0   -0   +In  -In */
 497   { eNR, eNR, eNR, eNI, ePI }, /* Normal */
 498   { eNR, eRZ, ePZ, eNI, ePI }, /* +0   */
 499   { eNR, eNZ, eRZ, eNI, ePI }, /* -0   */
 500   { ePI, ePI, ePI, eIn, ePI }, /* +Inf */
 501   { eNI, eNI, eNI, eNI, eIn }, /* -Inf */
 502 };
 503
 504 fp_t
 505 rxfp_sub (fp_t fa, fp_t fb)
 506 {
 507   FP_Parts a, b, c;
 508   fp_t rv;
 509   double da, db;
 510
 511   fp_explode (fa, &a);
 512   fp_explode (fb, &b);
 513   CHECK_EXCEPTIONS (a, b, rv, ex_sub_tab);
 514
 515   da = fp_to_double (&a);
 516   db = fp_to_double (&b);
 517   tprintf("%g - %g = %g\n", da, db, da-db);
 518
 519   double_to_fp (da-db, &c);
 520   rv = fp_implode (&c);
 521
 522   return rv;
 523 }
 524
 525 static FP_ExceptionCases ex_mul_tab[5][5] = {
 526   /* N   +0   -0   +In  -In */
 527   { eNR, eNR, eNR, eSI, eSI }, /* Normal */
 528   { eNR, ePZ, eNZ, eIn, eIn }, /* +0   */
 529   { eNR, eNZ, ePZ, eIn, eIn }, /* -0   */
 530   { eSI, eIn, eIn, ePI, eNI }, /* +Inf */
 531   { eSI, eIn, eIn, eNI, ePI }, /* -Inf */
 532 };
 533
 534 fp_t
 535 rxfp_mul (fp_t fa, fp_t fb)
 536 {
 537   FP_Parts a, b, c;
 538   fp_t rv;
 539   double da, db;
 540
 541   fp_explode (fa, &a);
 542   fp_explode (fb, &b);
 543   CHECK_EXCEPTIONS (a, b, rv, ex_mul_tab);
 544
 545   da = fp_to_double (&a);
 546   db = fp_to_double (&b);
 547   tprintf("%g x %g = %g\n", da, db, da*db);
 548
 549   double_to_fp (da*db, &c);
 550   rv = fp_implode (&c);
 551
 552   return rv;
 553 }
 554
 555 static FP_ExceptionCases ex_div_tab[5][5] = {
 556   /* N   +0   -0   +In  -In */
 557   { eNR, eDZ, eDZ, eSZ, eSZ }, /* Normal */
 558   { eSZ, eIn, eIn, ePZ, eNZ }, /* +0   */
 559   { eSZ, eIn, eIn, eNZ, ePZ }, /* -0   */
 560   { eSI, ePI, eNI, eIn, eIn }, /* +Inf */
 561   { eSI, eNI, ePI, eIn, eIn }, /* -Inf */
 562 };
 563
 564 fp_t
 565 rxfp_div (fp_t fa, fp_t fb)
 566 {
 567   FP_Parts a, b, c;
 568   fp_t rv;
 569   double da, db;
 570
 571   fp_explode (fa, &a);
 572   fp_explode (fb, &b);
 573   CHECK_EXCEPTIONS (a, b, rv, ex_div_tab);
 574
 575   da = fp_to_double (&a);
 576   db = fp_to_double (&b);
 577   tprintf("%g / %g = %g\n", da, db, da/db);
 578
 579   double_to_fp (da/db, &c);
 580   rv = fp_implode (&c);
 581
 582   return rv;
 583 }
 584
 585 static FP_ExceptionCases ex_cmp_tab[5][5] = {
 586   /* N   +0   -0   +In  -In */
 587   { eNR, eNR, eNR, eLT, eGT }, /* Normal */
 588   { eNR, eEQ, eEQ, eLT, eGT }, /* +0   */
 589   { eNR, eEQ, eEQ, eLT, eGT }, /* -0   */
 590   { eGT, eGT, eGT, eEQ, eGT }, /* +Inf */
 591   { eLT, eLT, eLT, eLT, eEQ }, /* -Inf */
 592 };
 593
 594 void
 595 rxfp_cmp (fp_t fa, fp_t fb)
 596 {
 597   FP_Parts a, b;
 598   fp_t c;
 599   FP_ExceptionCases reason;
 600   int flags = 0;
 601   double da, db;
 602
 603   fp_explode (fa, &a);
 604   fp_explode (fb, &b);
 605
 606   if (check_exceptions (&a, &b, &c, ex_cmp_tab, &reason))
 607     {
 608       if (reason == eQN)
 609         {
 610           /* Special case - incomparable.  */
 611           set_flags (FLAGBIT_Z | FLAGBIT_S | FLAGBIT_O, FLAGBIT_O);
 612           return;
 613         }
 614       return;
 615     }
 616
 617   switch (reason)
 618     {
 619     case eEQ:
 620       flags = FLAGBIT_Z;
 621       break;
 622     case eLT:
 623       flags = FLAGBIT_S;
 624       break;
 625     case eGT:
 626       flags = 0;
 627       break;
 628     case eNR:
 629       da = fp_to_double (&a);
 630       db = fp_to_double (&b);
 631       tprintf("fcmp: %g cmp %g\n", da, db);
 632       if (da < db)
 633         flags = FLAGBIT_S;
 634       else if (da == db)
 635         flags = FLAGBIT_Z;
 636       else
 637         flags = 0;
 638       break;
 639     default:
 640       abort();
 641     }
 642
 643   set_flags (FLAGBIT_Z | FLAGBIT_S | FLAGBIT_O, flags);
 644 }
 645
 646 long
 647 rxfp_ftoi (fp_t fa, int round_mode)
 648 {
 649   FP_Parts a;
 650   fp_t rv;
 651   int sign;
 652   int whole_bits, frac_bits;
 653
 654   fp_explode (fa, &a);
 655   sign = fa & 0x80000000UL;
 656
 657   switch (a.type)
 658     {
 659     case FP_NORMAL:
 660       break;
 661     case FP_PZERO:
 662     case FP_NZERO:
 663       return 0;
 664     case FP_PINFINITY:
 665       FP_RAISE (V);
 666       return 0x7fffffffL;
 667     case FP_NINFINITY:
 668       FP_RAISE (V);
 669       return 0x80000000L;
 670     case FP_DENORMAL:
 671       FP_RAISE (E);
 672       return 0;
 673     case FP_QNAN:
 674     case FP_SNAN:
 675       FP_RAISE (V);
 676       return sign ? 0x80000000U : 0x7fffffff;
 677     }
 678
 679   if (a.exp >= 31)
 680     {
 681       FP_RAISE (V);
 682       return sign ? 0x80000000U : 0x7fffffff;
 683     }
 684
 685   a.exp -= 23;
 686
 687   if (a.exp <= -25)
 688     {
 689       /* Less than 0.49999 */
 690       frac_bits = a.mant;
 691       whole_bits = 0;
 692     }
 693   else if (a.exp < 0)
 694     {
 695       frac_bits = a.mant << (32 + a.exp);
 696       whole_bits = a.mant >> (-a.exp);
 697     }
 698   else
 699     {
 700       frac_bits = 0;
 701       whole_bits = a.mant << a.exp;
 702     }
 703
 704   if (frac_bits)
 705     {
 706       switch (round_mode & 3)
 707         {
 708         case FPRM_NEAREST:
 709           if (frac_bits & 0x80000000UL)
 710             whole_bits ++;
 711           break;
 712         case FPRM_ZERO:
 713           break;
 714         case FPRM_PINF:
 715           if (!sign)
 716             whole_bits ++;
 717           break;
 718         case FPRM_NINF:
 719           if (sign)
 720             whole_bits ++;
 721           break;
 722         }
 723     }
 724
 725   rv = sign ? -whole_bits : whole_bits;
 726
 727   return rv;
 728 }
 729
 730 fp_t
 731 rxfp_itof (long fa, int round_mode)
 732 {
 733   fp_t rv;
 734   int sign = 0;
 735   unsigned int frac_bits;
 736   volatile unsigned int whole_bits;
 737   FP_Parts a = {0};
 738
 739   if (fa == 0)
 740     return PLUS_ZERO;
 741
 742   if (fa < 0)
 743     {
 744       fa = -fa;
 745       sign = 1;
 746       a.sign = -1;
 747     }
 748   else
 749     a.sign = 1;
 750
 751   whole_bits = fa;
 752   a.exp = 31;
 753
 754   while (! (whole_bits & 0x80000000UL))
 755     {
 756       a.exp --;
 757       whole_bits <<= 1;
 758     }
 759   frac_bits = whole_bits & 0xff;
 760   whole_bits = whole_bits >> 8;
 761
 762   if (frac_bits)
 763     {
 764       /* We must round */
 765       switch (round_mode & 3)
 766         {
 767         case FPRM_NEAREST:
 768           if (frac_bits & 0x80)
 769             whole_bits ++;
 770           break;
 771         case FPRM_ZERO:
 772           break;
 773         case FPRM_PINF:
 774           if (!sign)
 775             whole_bits ++;
 776           break;
 777         case FPRM_NINF:
 778           if (sign)
 779             whole_bits ++;
 780           break;
 781         }
 782     }
 783
 784   a.mant = whole_bits;
 785   if (whole_bits & 0xff000000UL)
 786     {
 787       a.mant >>= 1;
 788       a.exp ++;
 789     }
 790
 791   rv = fp_implode (&a);
 792   return rv;
 793 }
 794