gcc/real.c

   1 /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
   2    and support for XFmode IEEE extended real floating point arithmetic.
   3    Copyright (C) 1993, 1994, 1995 Free Software Foundation, Inc.
   4    Contributed by Stephen L. Moshier (moshier@world.std.com).
   5
   6 This file is part of GNU CC.
   7
   8 GNU CC 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 2, or (at your option)
  11 any later version.
  12
  13 GNU CC 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 GNU CC; see the file COPYING.  If not, write to
  20 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  21
  22 #include <stdio.h>
  23 #include <errno.h>
  24 #include "config.h"
  25 #include "tree.h"
  26
  27 #ifndef errno
  28 extern int errno;
  29 #endif
  30
  31 /* To enable support of XFmode extended real floating point, define
  32 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
  33
  34 To support cross compilation between IEEE, VAX and IBM floating
  35 point formats, define REAL_ARITHMETIC in the tm.h file.
  36
  37 In either case the machine files (tm.h) must not contain any code
  38 that tries to use host floating point arithmetic to convert
  39 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
  40 etc.  In cross-compile situations a REAL_VALUE_TYPE may not
  41 be intelligible to the host computer's native arithmetic.
  42
  43 The emulator defaults to the host's floating point format so that
  44 its decimal conversion functions can be used if desired (see
  45 real.h).
  46
  47 The first part of this file interfaces gcc to a floating point
  48 arithmetic suite that was not written with gcc in mind.  Avoid
  49 changing the low-level arithmetic routines unless you have suitable
  50 test programs available.  A special version of the PARANOIA floating
  51 point arithmetic tester, modified for this purpose, can be found on
  52 usc.edu: /pub/C-numanal/ieeetest.zoo.  Other tests, and libraries of
  53 XFmode and TFmode transcendental functions, can be obtained by ftp from
  54 netlib.att.com: netlib/cephes.   */
  55 \f
  56 /* Type of computer arithmetic.
  57    Only one of DEC, IBM, IEEE, or UNK should get defined.
  58
  59    `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
  60    to big-endian IEEE floating-point data structure.  This definition
  61    should work in SFmode `float' type and DFmode `double' type on
  62    virtually all big-endian IEEE machines.  If LONG_DOUBLE_TYPE_SIZE
  63    has been defined to be 96, then IEEE also invokes the particular
  64    XFmode (`long double' type) data structure used by the Motorola
  65    680x0 series processors.
  66
  67    `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
  68    little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
  69    has been defined to be 96, then IEEE also invokes the particular
  70    XFmode `long double' data structure used by the Intel 80x86 series
  71    processors.
  72
  73    `DEC' refers specifically to the Digital Equipment Corp PDP-11
  74    and VAX floating point data structure.  This model currently
  75    supports no type wider than DFmode.
  76
  77    `IBM' refers specifically to the IBM System/370 and compatible
  78    floating point data structure.  This model currently supports
  79    no type wider than DFmode.  The IBM conversions were contributed by
  80    frank@atom.ansto.gov.au (Frank Crawford).
  81
  82    If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
  83    then `long double' and `double' are both implemented, but they
  84    both mean DFmode.  In this case, the software floating-point
  85    support available here is activated by writing
  86       #define REAL_ARITHMETIC
  87    in tm.h.
  88
  89    The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
  90    and may deactivate XFmode since `long double' is used to refer
  91    to both modes.
  92
  93    The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
  94    contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
  95    separate the floating point unit's endian-ness from that of
  96    the integer addressing.  This permits one to define a big-endian
  97    FPU on a little-endian machine (e.g., ARM).  An extension to
  98    BYTES_BIG_ENDIAN may be required for some machines in the future.
  99    These optional macros may be defined in tm.h.  In real.h, they
 100    default to WORDS_BIG_ENDIAN, etc., so there is no need to define
 101    them for any normal host or target machine on which the floats
 102    and the integers have the same endian-ness.   */
 103
 104
 105 /* The following converts gcc macros into the ones used by this file.  */
 106
 107 /* REAL_ARITHMETIC defined means that macros in real.h are
 108    defined to call emulator functions.  */
 109 #ifdef REAL_ARITHMETIC
 110
 111 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
 112 /* PDP-11, Pro350, VAX: */
 113 #define DEC 1
 114 #else /* it's not VAX */
 115 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
 116 /* IBM System/370 style */
 117 #define IBM 1
 118 #else /* it's also not an IBM */
 119 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
 120 #define IEEE
 121 #else /* it's not IEEE either */
 122 /* UNKnown arithmetic.  We don't support this and can't go on. */
 123 unknown arithmetic type
 124 #define UNK 1
 125 #endif /* not IEEE */
 126 #endif /* not IBM */
 127 #endif /* not VAX */
 128
 129 #define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
 130
 131 #else
 132 /* REAL_ARITHMETIC not defined means that the *host's* data
 133    structure will be used.  It may differ by endian-ness from the
 134    target machine's structure and will get its ends swapped
 135    accordingly (but not here).  Probably only the decimal <-> binary
 136    functions in this file will actually be used in this case.  */
 137
 138 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
 139 #define DEC 1
 140 #else /* it's not VAX */
 141 #if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
 142 /* IBM System/370 style */
 143 #define IBM 1
 144 #else /* it's also not an IBM */
 145 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
 146 #define IEEE
 147 #else /* it's not IEEE either */
 148 unknown arithmetic type
 149 #define UNK 1
 150 #endif /* not IEEE */
 151 #endif /* not IBM */
 152 #endif /* not VAX */
 153
 154 #define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
 155
 156 #endif /* REAL_ARITHMETIC not defined */
 157
 158 /* Define INFINITY for support of infinity.
 159    Define NANS for support of Not-a-Number's (NaN's).  */
 160 #if !defined(DEC) && !defined(IBM)
 161 #define INFINITY
 162 #define NANS
 163 #endif
 164
 165 /* Support of NaNs requires support of infinity. */
 166 #ifdef NANS
 167 #ifndef INFINITY
 168 #define INFINITY
 169 #endif
 170 #endif
 171 \f
 172 /* Find a host integer type that is at least 16 bits wide,
 173    and another type at least twice whatever that size is. */
 174
 175 #if HOST_BITS_PER_CHAR >= 16
 176 #define EMUSHORT char
 177 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
 178 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
 179 #else
 180 #if HOST_BITS_PER_SHORT >= 16
 181 #define EMUSHORT short
 182 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
 183 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
 184 #else
 185 #if HOST_BITS_PER_INT >= 16
 186 #define EMUSHORT int
 187 #define EMUSHORT_SIZE HOST_BITS_PER_INT
 188 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
 189 #else
 190 #if HOST_BITS_PER_LONG >= 16
 191 #define EMUSHORT long
 192 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
 193 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
 194 #else
 195 /*  You will have to modify this program to have a smaller unit size. */
 196 #define EMU_NON_COMPILE
 197 #endif
 198 #endif
 199 #endif
 200 #endif
 201
 202 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
 203 #define EMULONG short
 204 #else
 205 #if HOST_BITS_PER_INT >= EMULONG_SIZE
 206 #define EMULONG int
 207 #else
 208 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
 209 #define EMULONG long
 210 #else
 211 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
 212 #define EMULONG long long int
 213 #else
 214 /*  You will have to modify this program to have a smaller unit size. */
 215 #define EMU_NON_COMPILE
 216 #endif
 217 #endif
 218 #endif
 219 #endif
 220
 221
 222 /* The host interface doesn't work if no 16-bit size exists. */
 223 #if EMUSHORT_SIZE != 16
 224 #define EMU_NON_COMPILE
 225 #endif
 226
 227 /* OK to continue compilation. */
 228 #ifndef EMU_NON_COMPILE
 229
 230 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
 231    In GET_REAL and PUT_REAL, r and e are pointers.
 232    A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
 233    in memory, with no holes.  */
 234
 235 #if LONG_DOUBLE_TYPE_SIZE == 96
 236 /* Number of 16 bit words in external e type format */
 237 #define NE 6
 238 #define MAXDECEXP 4932
 239 #define MINDECEXP -4956
 240 #define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
 241 #define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
 242 #else /* no XFmode */
 243 #if LONG_DOUBLE_TYPE_SIZE == 128
 244 #define NE 10
 245 #define MAXDECEXP 4932
 246 #define MINDECEXP -4977
 247 #define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
 248 #define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
 249 #else
 250 #define NE 6
 251 #define MAXDECEXP 4932
 252 #define MINDECEXP -4956
 253 #ifdef REAL_ARITHMETIC
 254 /* Emulator uses target format internally
 255    but host stores it in host endian-ness. */
 256
 257 #define GET_REAL(r,e)                                           \
 258 do {                                                            \
 259      if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN)  \
 260        e53toe ((unsigned EMUSHORT*) (r), (e));                  \
 261      else                                                       \
 262        {                                                        \
 263          unsigned EMUSHORT w[4];                                \
 264          w[3] = ((EMUSHORT *) r)[0];                            \
 265          w[2] = ((EMUSHORT *) r)[1];                            \
 266          w[1] = ((EMUSHORT *) r)[2];                            \
 267          w[0] = ((EMUSHORT *) r)[3];                            \
 268          e53toe (w, (e));                                       \
 269        }                                                        \
 270    } while (0)
 271
 272 #define PUT_REAL(e,r)                                           \
 273 do {                                                            \
 274      if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN)  \
 275        etoe53 ((e), (unsigned EMUSHORT *) (r));                 \
 276      else                                                       \
 277        {                                                        \
 278          unsigned EMUSHORT w[4];                                \
 279          etoe53 ((e), w);                                       \
 280          *((EMUSHORT *) r) = w[3];                              \
 281          *((EMUSHORT *) r + 1) = w[2];                          \
 282          *((EMUSHORT *) r + 2) = w[1];                          \
 283          *((EMUSHORT *) r + 3) = w[0];                          \
 284        }                                                        \
 285    } while (0)
 286
 287 #else /* not REAL_ARITHMETIC */
 288
 289 /* emulator uses host format */
 290 #define GET_REAL(r,e) e53toe ((unsigned EMUSHORT *) (r), (e))
 291 #define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
 292
 293 #endif /* not REAL_ARITHMETIC */
 294 #endif /* not TFmode */
 295 #endif /* no XFmode */
 296
 297
 298 /* Number of 16 bit words in internal format */
 299 #define NI (NE+3)
 300
 301 /* Array offset to exponent */
 302 #define E 1
 303
 304 /* Array offset to high guard word */
 305 #define M 2
 306
 307 /* Number of bits of precision */
 308 #define NBITS ((NI-4)*16)
 309
 310 /* Maximum number of decimal digits in ASCII conversion
 311  * = NBITS*log10(2)
 312  */
 313 #define NDEC (NBITS*8/27)
 314
 315 /* The exponent of 1.0 */
 316 #define EXONE (0x3fff)
 317
 318 extern int extra_warnings;
 319 extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
 320 extern unsigned EMUSHORT elog2[], esqrt2[];
 321
 322 static void endian      PROTO((unsigned EMUSHORT *, long *,
 323                                enum machine_mode));
 324 static void eclear      PROTO((unsigned EMUSHORT *));
 325 static void emov        PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 326 static void eabs        PROTO((unsigned EMUSHORT *));
 327 static void eneg        PROTO((unsigned EMUSHORT *));
 328 static int eisneg       PROTO((unsigned EMUSHORT *));
 329 static int eisinf       PROTO((unsigned EMUSHORT *));
 330 static int eisnan       PROTO((unsigned EMUSHORT *));
 331 static void einfin      PROTO((unsigned EMUSHORT *));
 332 static void enan        PROTO((unsigned EMUSHORT *, int));
 333 static void emovi       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 334 static void emovo       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 335 static void ecleaz      PROTO((unsigned EMUSHORT *));
 336 static void ecleazs     PROTO((unsigned EMUSHORT *));
 337 static void emovz       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 338 static void einan       PROTO((unsigned EMUSHORT *));
 339 static int eiisnan      PROTO((unsigned EMUSHORT *));
 340 static int eiisneg      PROTO((unsigned EMUSHORT *));
 341 static void eiinfin     PROTO((unsigned EMUSHORT *));
 342 static int eiisinf      PROTO((unsigned EMUSHORT *));
 343 static int ecmpm        PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 344 static void eshdn1      PROTO((unsigned EMUSHORT *));
 345 static void eshup1      PROTO((unsigned EMUSHORT *));
 346 static void eshdn8      PROTO((unsigned EMUSHORT *));
 347 static void eshup8      PROTO((unsigned EMUSHORT *));
 348 static void eshup6      PROTO((unsigned EMUSHORT *));
 349 static void eshdn6      PROTO((unsigned EMUSHORT *));
 350 static void eaddm       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));\f
 351 static void esubm       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 352 static void m16m        PROTO((unsigned int, unsigned short *,
 353                                unsigned short *));
 354 static int edivm        PROTO((unsigned short *, unsigned short *));
 355 static int emulm        PROTO((unsigned short *, unsigned short *));
 356 static void emdnorm     PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
 357 static void esub        PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 358                                unsigned EMUSHORT *));
 359 static void eadd        PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 360                                unsigned EMUSHORT *));
 361 static void eadd1       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 362                                unsigned EMUSHORT *));
 363 static void ediv        PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 364                                unsigned EMUSHORT *));
 365 static void emul        PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 366                                unsigned EMUSHORT *));
 367 static void e53toe      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 368 static void e64toe      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 369 static void e113toe     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 370 static void e24toe      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 371 static void etoe113     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 372 static void toe113      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 373 static void etoe64      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 374 static void toe64       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 375 static void etoe53      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 376 static void toe53       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 377 static void etoe24      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 378 static void toe24       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 379 static int ecmp         PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 380 static void eround      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 381 static void ltoe        PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
 382 static void ultoe       PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
 383 static void eifrac      PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
 384                                unsigned EMUSHORT *));
 385 static void euifrac     PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
 386                                unsigned EMUSHORT *));
 387 static int eshift       PROTO((unsigned EMUSHORT *, int));
 388 static int enormlz      PROTO((unsigned EMUSHORT *));
 389 static void e24toasc    PROTO((unsigned EMUSHORT *, char *, int));
 390 static void e53toasc    PROTO((unsigned EMUSHORT *, char *, int));
 391 static void e64toasc    PROTO((unsigned EMUSHORT *, char *, int));
 392 static void e113toasc   PROTO((unsigned EMUSHORT *, char *, int));
 393 static void etoasc      PROTO((unsigned EMUSHORT *, char *, int));
 394 static void asctoe24    PROTO((char *, unsigned EMUSHORT *));
 395 static void asctoe53    PROTO((char *, unsigned EMUSHORT *));
 396 static void asctoe64    PROTO((char *, unsigned EMUSHORT *));
 397 static void asctoe113   PROTO((char *, unsigned EMUSHORT *));
 398 static void asctoe      PROTO((char *, unsigned EMUSHORT *));
 399 static void asctoeg     PROTO((char *, unsigned EMUSHORT *, int));
 400 static void efloor      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 401 static void efrexp      PROTO((unsigned EMUSHORT *, int *,
 402                                unsigned EMUSHORT *));
 403 static void eldexp      PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
 404 static void eremain     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 405                                unsigned EMUSHORT *));
 406 static void eiremain    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 407 static void mtherr      PROTO((char *, int));
 408 static void dectoe      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 409 static void etodec      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 410 static void todec       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 411 static void ibmtoe      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 412                                enum machine_mode));
 413 static void etoibm      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 414                                enum machine_mode));
 415 static void toibm       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
 416                                enum machine_mode));
 417 static void make_nan    PROTO((unsigned EMUSHORT *, int, enum machine_mode));
 418 static void uditoe      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 419 static void ditoe       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 420 static void etoudi      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 421 static void etodi       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 422 static void esqrt       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
 423 \f
 424 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
 425    swapping ends if required, into output array of longs.  The
 426    result is normally passed to fprintf by the ASM_OUTPUT_ macros.   */
 427
 428 static void
 429 endian (e, x, mode)
 430      unsigned EMUSHORT e[];
 431      long x[];
 432      enum machine_mode mode;
 433 {
 434   unsigned long th, t;
 435
 436   if (REAL_WORDS_BIG_ENDIAN)
 437     {
 438       switch (mode)
 439         {
 440
 441         case TFmode:
 442           /* Swap halfwords in the fourth long. */
 443           th = (unsigned long) e[6] & 0xffff;
 444           t = (unsigned long) e[7] & 0xffff;
 445           t |= th << 16;
 446           x[3] = (long) t;
 447
 448         case XFmode:
 449
 450           /* Swap halfwords in the third long. */
 451           th = (unsigned long) e[4] & 0xffff;
 452           t = (unsigned long) e[5] & 0xffff;
 453           t |= th << 16;
 454           x[2] = (long) t;
 455           /* fall into the double case */
 456
 457         case DFmode:
 458
 459           /* swap halfwords in the second word */
 460           th = (unsigned long) e[2] & 0xffff;
 461           t = (unsigned long) e[3] & 0xffff;
 462           t |= th << 16;
 463           x[1] = (long) t;
 464           /* fall into the float case */
 465
 466         case HFmode:
 467         case SFmode:
 468
 469           /* swap halfwords in the first word */
 470           th = (unsigned long) e[0] & 0xffff;
 471           t = (unsigned long) e[1] & 0xffff;
 472           t |= th << 16;
 473           x[0] = t;
 474           break;
 475
 476         default:
 477           abort ();
 478         }
 479     }
 480   else
 481     {
 482       /* Pack the output array without swapping. */
 483
 484       switch (mode)
 485         {
 486
 487         case TFmode:
 488
 489           /* Pack the fourth long. */
 490           th = (unsigned long) e[7] & 0xffff;
 491           t = (unsigned long) e[6] & 0xffff;
 492           t |= th << 16;
 493           x[3] = (long) t;
 494
 495         case XFmode:
 496
 497           /* Pack the third long.
 498              Each element of the input REAL_VALUE_TYPE array has 16 useful bits
 499              in it.  */
 500           th = (unsigned long) e[5] & 0xffff;
 501           t = (unsigned long) e[4] & 0xffff;
 502           t |= th << 16;
 503           x[2] = (long) t;
 504           /* fall into the double case */
 505
 506         case DFmode:
 507
 508           /* pack the second long */
 509           th = (unsigned long) e[3] & 0xffff;
 510           t = (unsigned long) e[2] & 0xffff;
 511           t |= th << 16;
 512           x[1] = (long) t;
 513           /* fall into the float case */
 514
 515         case HFmode:
 516         case SFmode:
 517
 518           /* pack the first long */
 519           th = (unsigned long) e[1] & 0xffff;
 520           t = (unsigned long) e[0] & 0xffff;
 521           t |= th << 16;
 522           x[0] = t;
 523           break;
 524
 525         default:
 526           abort ();
 527         }
 528     }
 529 }
 530
 531
 532 /* This is the implementation of the REAL_ARITHMETIC macro.  */
 533
 534 void
 535 earith (value, icode, r1, r2)
 536      REAL_VALUE_TYPE *value;
 537      int icode;
 538      REAL_VALUE_TYPE *r1;
 539      REAL_VALUE_TYPE *r2;
 540 {
 541   unsigned EMUSHORT d1[NE], d2[NE], v[NE];
 542   enum tree_code code;
 543
 544   GET_REAL (r1, d1);
 545   GET_REAL (r2, d2);
 546 #ifdef NANS
 547 /*  Return NaN input back to the caller. */
 548   if (eisnan (d1))
 549     {
 550       PUT_REAL (d1, value);
 551       return;
 552     }
 553   if (eisnan (d2))
 554     {
 555       PUT_REAL (d2, value);
 556       return;
 557     }
 558 #endif
 559   code = (enum tree_code) icode;
 560   switch (code)
 561     {
 562     case PLUS_EXPR:
 563       eadd (d2, d1, v);
 564       break;
 565
 566     case MINUS_EXPR:
 567       esub (d2, d1, v);         /* d1 - d2 */
 568       break;
 569
 570     case MULT_EXPR:
 571       emul (d2, d1, v);
 572       break;
 573
 574     case RDIV_EXPR:
 575 #ifndef REAL_INFINITY
 576       if (ecmp (d2, ezero) == 0)
 577         {
 578 #ifdef NANS
 579         enan (v, eisneg (d1) ^ eisneg (d2));
 580         break;
 581 #else
 582         abort ();
 583 #endif
 584         }
 585 #endif
 586       ediv (d2, d1, v); /* d1/d2 */
 587       break;
 588
 589     case MIN_EXPR:              /* min (d1,d2) */
 590       if (ecmp (d1, d2) < 0)
 591         emov (d1, v);
 592       else
 593         emov (d2, v);
 594       break;
 595
 596     case MAX_EXPR:              /* max (d1,d2) */
 597       if (ecmp (d1, d2) > 0)
 598         emov (d1, v);
 599       else
 600         emov (d2, v);
 601       break;
 602     default:
 603       emov (ezero, v);
 604       break;
 605     }
 606 PUT_REAL (v, value);
 607 }
 608
 609
 610 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
 611    implements REAL_VALUE_RNDZINT (x) (etrunci (x)).  */
 612
 613 REAL_VALUE_TYPE
 614 etrunci (x)
 615      REAL_VALUE_TYPE x;
 616 {
 617   unsigned EMUSHORT f[NE], g[NE];
 618   REAL_VALUE_TYPE r;
 619   HOST_WIDE_INT l;
 620
 621   GET_REAL (&x, g);
 622 #ifdef NANS
 623   if (eisnan (g))
 624     return (x);
 625 #endif
 626   eifrac (g, &l, f);
 627   ltoe (&l, g);
 628   PUT_REAL (g, &r);
 629   return (r);
 630 }
 631
 632
 633 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
 634    implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)).  */
 635
 636 REAL_VALUE_TYPE
 637 etruncui (x)
 638      REAL_VALUE_TYPE x;
 639 {
 640   unsigned EMUSHORT f[NE], g[NE];
 641   REAL_VALUE_TYPE r;
 642   unsigned HOST_WIDE_INT l;
 643
 644   GET_REAL (&x, g);
 645 #ifdef NANS
 646   if (eisnan (g))
 647     return (x);
 648 #endif
 649   euifrac (g, &l, f);
 650   ultoe (&l, g);
 651   PUT_REAL (g, &r);
 652   return (r);
 653 }
 654
 655
 656 /* This is the REAL_VALUE_ATOF function.  It converts a decimal string to
 657    binary, rounding off as indicated by the machine_mode argument.  Then it
 658    promotes the rounded value to REAL_VALUE_TYPE.  */
 659
 660 REAL_VALUE_TYPE
 661 ereal_atof (s, t)
 662      char *s;
 663      enum machine_mode t;
 664 {
 665   unsigned EMUSHORT tem[NE], e[NE];
 666   REAL_VALUE_TYPE r;
 667
 668   switch (t)
 669     {
 670     case HFmode:
 671     case SFmode:
 672       asctoe24 (s, tem);
 673       e24toe (tem, e);
 674       break;
 675     case DFmode:
 676       asctoe53 (s, tem);
 677       e53toe (tem, e);
 678       break;
 679     case XFmode:
 680       asctoe64 (s, tem);
 681       e64toe (tem, e);
 682       break;
 683     case TFmode:
 684       asctoe113 (s, tem);
 685       e113toe (tem, e);
 686       break;
 687     default:
 688       asctoe (s, e);
 689     }
 690   PUT_REAL (e, &r);
 691   return (r);
 692 }
 693
 694
 695 /* Expansion of REAL_NEGATE.  */
 696
 697 REAL_VALUE_TYPE
 698 ereal_negate (x)
 699      REAL_VALUE_TYPE x;
 700 {
 701   unsigned EMUSHORT e[NE];
 702   REAL_VALUE_TYPE r;
 703
 704   GET_REAL (&x, e);
 705   eneg (e);
 706   PUT_REAL (e, &r);
 707   return (r);
 708 }
 709
 710
 711 /* Round real toward zero to HOST_WIDE_INT;
 712    implements REAL_VALUE_FIX (x).  */
 713
 714 HOST_WIDE_INT
 715 efixi (x)
 716      REAL_VALUE_TYPE x;
 717 {
 718   unsigned EMUSHORT f[NE], g[NE];
 719   HOST_WIDE_INT l;
 720
 721   GET_REAL (&x, f);
 722 #ifdef NANS
 723   if (eisnan (f))
 724     {
 725       warning ("conversion from NaN to int");
 726       return (-1);
 727     }
 728 #endif
 729   eifrac (f, &l, g);
 730   return l;
 731 }
 732
 733 /* Round real toward zero to unsigned HOST_WIDE_INT
 734    implements  REAL_VALUE_UNSIGNED_FIX (x).
 735    Negative input returns zero.  */
 736
 737 unsigned HOST_WIDE_INT
 738 efixui (x)
 739      REAL_VALUE_TYPE x;
 740 {
 741   unsigned EMUSHORT f[NE], g[NE];
 742   unsigned HOST_WIDE_INT l;
 743
 744   GET_REAL (&x, f);
 745 #ifdef NANS
 746   if (eisnan (f))
 747     {
 748       warning ("conversion from NaN to unsigned int");
 749       return (-1);
 750     }
 751 #endif
 752   euifrac (f, &l, g);
 753   return l;
 754 }
 755
 756
 757 /* REAL_VALUE_FROM_INT macro.  */
 758
 759 void
 760 ereal_from_int (d, i, j)
 761      REAL_VALUE_TYPE *d;
 762      HOST_WIDE_INT i, j;
 763 {
 764   unsigned EMUSHORT df[NE], dg[NE];
 765   HOST_WIDE_INT low, high;
 766   int sign;
 767
 768   sign = 0;
 769   low = i;
 770   if ((high = j) < 0)
 771     {
 772       sign = 1;
 773       /* complement and add 1 */
 774       high = ~high;
 775       if (low)
 776         low = -low;
 777       else
 778         high += 1;
 779     }
 780   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
 781   ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
 782   emul (dg, df, dg);
 783   ultoe ((unsigned HOST_WIDE_INT *) &low, df);
 784   eadd (df, dg, dg);
 785   if (sign)
 786     eneg (dg);
 787   PUT_REAL (dg, d);
 788 }
 789
 790
 791 /* REAL_VALUE_FROM_UNSIGNED_INT macro.   */
 792
 793 void
 794 ereal_from_uint (d, i, j)
 795      REAL_VALUE_TYPE *d;
 796      unsigned HOST_WIDE_INT i, j;
 797 {
 798   unsigned EMUSHORT df[NE], dg[NE];
 799   unsigned HOST_WIDE_INT low, high;
 800
 801   low = i;
 802   high = j;
 803   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
 804   ultoe (&high, dg);
 805   emul (dg, df, dg);
 806   ultoe (&low, df);
 807   eadd (df, dg, dg);
 808   PUT_REAL (dg, d);
 809 }
 810
 811
 812 /* REAL_VALUE_TO_INT macro.  */
 813
 814 void
 815 ereal_to_int (low, high, rr)
 816      HOST_WIDE_INT *low, *high;
 817      REAL_VALUE_TYPE rr;
 818 {
 819   unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
 820   int s;
 821
 822   GET_REAL (&rr, d);
 823 #ifdef NANS
 824   if (eisnan (d))
 825     {
 826       warning ("conversion from NaN to int");
 827       *low = -1;
 828       *high = -1;
 829       return;
 830     }
 831 #endif
 832   /* convert positive value */
 833   s = 0;
 834   if (eisneg (d))
 835     {
 836       eneg (d);
 837       s = 1;
 838     }
 839   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
 840   ediv (df, d, dg);             /* dg = d / 2^32 is the high word */
 841   euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
 842   emul (df, dh, dg);            /* fractional part is the low word */
 843   euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);
 844   if (s)
 845     {
 846       /* complement and add 1 */
 847       *high = ~(*high);
 848       if (*low)
 849         *low = -(*low);
 850       else
 851         *high += 1;
 852     }
 853 }
 854
 855
 856 /* REAL_VALUE_LDEXP macro.  */
 857
 858 REAL_VALUE_TYPE
 859 ereal_ldexp (x, n)
 860      REAL_VALUE_TYPE x;
 861      int n;
 862 {
 863   unsigned EMUSHORT e[NE], y[NE];
 864   REAL_VALUE_TYPE r;
 865
 866   GET_REAL (&x, e);
 867 #ifdef NANS
 868   if (eisnan (e))
 869     return (x);
 870 #endif
 871   eldexp (e, n, y);
 872   PUT_REAL (y, &r);
 873   return (r);
 874 }
 875
 876 /* These routines are conditionally compiled because functions
 877    of the same names may be defined in fold-const.c.  */
 878
 879 #ifdef REAL_ARITHMETIC
 880
 881 /* Check for infinity in a REAL_VALUE_TYPE. */
 882
 883 int
 884 target_isinf (x)
 885      REAL_VALUE_TYPE x;
 886 {
 887   unsigned EMUSHORT e[NE];
 888
 889 #ifdef INFINITY
 890   GET_REAL (&x, e);
 891   return (eisinf (e));
 892 #else
 893   return 0;
 894 #endif
 895 }
 896
 897 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
 898
 899 int
 900 target_isnan (x)
 901      REAL_VALUE_TYPE x;
 902 {
 903   unsigned EMUSHORT e[NE];
 904
 905 #ifdef NANS
 906   GET_REAL (&x, e);
 907   return (eisnan (e));
 908 #else
 909   return (0);
 910 #endif
 911 }
 912
 913
 914 /* Check for a negative REAL_VALUE_TYPE number.
 915    This just checks the sign bit, so that -0 counts as negative. */
 916
 917 int
 918 target_negative (x)
 919      REAL_VALUE_TYPE x;
 920 {
 921   return ereal_isneg (x);
 922 }
 923
 924 /* Expansion of REAL_VALUE_TRUNCATE.
 925    The result is in floating point, rounded to nearest or even.  */
 926
 927 REAL_VALUE_TYPE
 928 real_value_truncate (mode, arg)
 929      enum machine_mode mode;
 930      REAL_VALUE_TYPE arg;
 931 {
 932   unsigned EMUSHORT e[NE], t[NE];
 933   REAL_VALUE_TYPE r;
 934
 935   GET_REAL (&arg, e);
 936 #ifdef NANS
 937   if (eisnan (e))
 938     return (arg);
 939 #endif
 940   eclear (t);
 941   switch (mode)
 942     {
 943     case TFmode:
 944       etoe113 (e, t);
 945       e113toe (t, t);
 946       break;
 947
 948     case XFmode:
 949       etoe64 (e, t);
 950       e64toe (t, t);
 951       break;
 952
 953     case DFmode:
 954       etoe53 (e, t);
 955       e53toe (t, t);
 956       break;
 957
 958     case HFmode:
 959     case SFmode:
 960       etoe24 (e, t);
 961       e24toe (t, t);
 962       break;
 963
 964     case SImode:
 965       r = etrunci (arg);
 966       return (r);
 967
 968     /* If an unsupported type was requested, presume that
 969        the machine files know something useful to do with
 970        the unmodified value.  */
 971
 972     default:
 973       return (arg);
 974     }
 975   PUT_REAL (t, &r);
 976   return (r);
 977 }
 978
 979 #endif /* REAL_ARITHMETIC defined */
 980
 981 /* Used for debugging--print the value of R in human-readable format
 982    on stderr.  */
 983
 984 void
 985 debug_real (r)
 986      REAL_VALUE_TYPE r;
 987 {
 988   char dstr[30];
 989
 990   REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
 991   fprintf (stderr, "%s", dstr);
 992 }
 993
 994 \f
 995 /* The following routines convert REAL_VALUE_TYPE to the various floating
 996    point formats that are meaningful to supported computers.
 997
 998    The results are returned in 32-bit pieces, each piece stored in a `long'.
 999    This is so they can be printed by statements like
1000
1001       fprintf (file, "%lx, %lx", L[0],  L[1]);
1002
1003    that will work on both narrow- and wide-word host computers.  */
1004
1005 /* Convert R to a 128-bit long double precision value.  The output array L
1006    contains four 32-bit pieces of the result, in the order they would appear
1007    in memory.  */
1008
1009 void
1010 etartdouble (r, l)
1011      REAL_VALUE_TYPE r;
1012      long l[];
1013 {
1014   unsigned EMUSHORT e[NE];
1015
1016   GET_REAL (&r, e);
1017   etoe113 (e, e);
1018   endian (e, l, TFmode);
1019 }
1020
1021 /* Convert R to a double extended precision value.  The output array L
1022    contains three 32-bit pieces of the result, in the order they would
1023    appear in memory.  */
1024
1025 void
1026 etarldouble (r, l)
1027      REAL_VALUE_TYPE r;
1028      long l[];
1029 {
1030   unsigned EMUSHORT e[NE];
1031
1032   GET_REAL (&r, e);
1033   etoe64 (e, e);
1034   endian (e, l, XFmode);
1035 }
1036
1037 /* Convert R to a double precision value.  The output array L contains two
1038    32-bit pieces of the result, in the order they would appear in memory.  */
1039
1040 void
1041 etardouble (r, l)
1042      REAL_VALUE_TYPE r;
1043      long l[];
1044 {
1045   unsigned EMUSHORT e[NE];
1046
1047   GET_REAL (&r, e);
1048   etoe53 (e, e);
1049   endian (e, l, DFmode);
1050 }
1051
1052 /* Convert R to a single precision float value stored in the least-significant
1053    bits of a `long'.  */
1054
1055 long
1056 etarsingle (r)
1057      REAL_VALUE_TYPE r;
1058 {
1059   unsigned EMUSHORT e[NE];
1060   long l;
1061
1062   GET_REAL (&r, e);
1063   etoe24 (e, e);
1064   endian (e, &l, SFmode);
1065   return ((long) l);
1066 }
1067
1068 /* Convert X to a decimal ASCII string S for output to an assembly
1069    language file.  Note, there is no standard way to spell infinity or
1070    a NaN, so these values may require special treatment in the tm.h
1071    macros.  */
1072
1073 void
1074 ereal_to_decimal (x, s)
1075      REAL_VALUE_TYPE x;
1076      char *s;
1077 {
1078   unsigned EMUSHORT e[NE];
1079
1080   GET_REAL (&x, e);
1081   etoasc (e, s, 20);
1082 }
1083
1084 /* Compare X and Y.  Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1085    or -2 if either is a NaN.   */
1086
1087 int
1088 ereal_cmp (x, y)
1089      REAL_VALUE_TYPE x, y;
1090 {
1091   unsigned EMUSHORT ex[NE], ey[NE];
1092
1093   GET_REAL (&x, ex);
1094   GET_REAL (&y, ey);
1095   return (ecmp (ex, ey));
1096 }
1097
1098 /*  Return 1 if the sign bit of X is set, else return 0.  */
1099
1100 int
1101 ereal_isneg (x)
1102      REAL_VALUE_TYPE x;
1103 {
1104   unsigned EMUSHORT ex[NE];
1105
1106   GET_REAL (&x, ex);
1107   return (eisneg (ex));
1108 }
1109
1110 /* End of REAL_ARITHMETIC interface */
1111 \f
1112 /*
1113   Extended precision IEEE binary floating point arithmetic routines
1114
1115   Numbers are stored in C language as arrays of 16-bit unsigned
1116   short integers.  The arguments of the routines are pointers to
1117   the arrays.
1118
1119   External e type data structure, similar to Intel 8087 chip
1120   temporary real format but possibly with a larger significand:
1121
1122         NE-1 significand words  (least significant word first,
1123                                  most significant bit is normally set)
1124         exponent                (value = EXONE for 1.0,
1125                                 top bit is the sign)
1126
1127
1128   Internal exploded e-type data structure of a number (a "word" is 16 bits):
1129
1130   ei[0] sign word       (0 for positive, 0xffff for negative)
1131   ei[1] biased exponent (value = EXONE for the number 1.0)
1132   ei[2] high guard word (always zero after normalization)
1133   ei[3]
1134   to ei[NI-2]   significand     (NI-4 significand words,
1135                                  most significant word first,
1136                                  most significant bit is set)
1137   ei[NI-1]      low guard word  (0x8000 bit is rounding place)
1138
1139
1140
1141                 Routines for external format e-type numbers
1142
1143         asctoe (string, e)      ASCII string to extended double e type
1144         asctoe64 (string, &d)   ASCII string to long double
1145         asctoe53 (string, &d)   ASCII string to double
1146         asctoe24 (string, &f)   ASCII string to single
1147         asctoeg (string, e, prec) ASCII string to specified precision
1148         e24toe (&f, e)          IEEE single precision to e type
1149         e53toe (&d, e)          IEEE double precision to e type
1150         e64toe (&d, e)          IEEE long double precision to e type
1151         e113toe (&d, e)         128-bit long double precision to e type
1152         eabs (e)                        absolute value
1153         eadd (a, b, c)          c = b + a
1154         eclear (e)              e = 0
1155         ecmp (a, b)             Returns 1 if a > b, 0 if a == b,
1156                                 -1 if a < b, -2 if either a or b is a NaN.
1157         ediv (a, b, c)          c = b / a
1158         efloor (a, b)           truncate to integer, toward -infinity
1159         efrexp (a, exp, s)      extract exponent and significand
1160         eifrac (e, &l, frac)    e to HOST_WIDE_INT and e type fraction
1161         euifrac (e, &l, frac)   e to unsigned HOST_WIDE_INT and e type fraction
1162         einfin (e)              set e to infinity, leaving its sign alone
1163         eldexp (a, n, b)        multiply by 2**n
1164         emov (a, b)             b = a
1165         emul (a, b, c)          c = b * a
1166         eneg (e)                        e = -e
1167         eround (a, b)           b = nearest integer value to a
1168         esub (a, b, c)          c = b - a
1169         e24toasc (&f, str, n)   single to ASCII string, n digits after decimal
1170         e53toasc (&d, str, n)   double to ASCII string, n digits after decimal
1171         e64toasc (&d, str, n)   80-bit long double to ASCII string
1172         e113toasc (&d, str, n)  128-bit long double to ASCII string
1173         etoasc (e, str, n)      e to ASCII string, n digits after decimal
1174         etoe24 (e, &f)          convert e type to IEEE single precision
1175         etoe53 (e, &d)          convert e type to IEEE double precision
1176         etoe64 (e, &d)          convert e type to IEEE long double precision
1177         ltoe (&l, e)            HOST_WIDE_INT to e type
1178         ultoe (&l, e)           unsigned HOST_WIDE_INT to e type
1179         eisneg (e)              1 if sign bit of e != 0, else 0
1180         eisinf (e)              1 if e has maximum exponent (non-IEEE)
1181                                 or is infinite (IEEE)
1182         eisnan (e)              1 if e is a NaN
1183
1184
1185                 Routines for internal format exploded e-type numbers
1186
1187         eaddm (ai, bi)          add significands, bi = bi + ai
1188         ecleaz (ei)             ei = 0
1189         ecleazs (ei)            set ei = 0 but leave its sign alone
1190         ecmpm (ai, bi)          compare significands, return 1, 0, or -1
1191         edivm (ai, bi)          divide  significands, bi = bi / ai
1192         emdnorm (ai,l,s,exp)    normalize and round off
1193         emovi (a, ai)           convert external a to internal ai
1194         emovo (ai, a)           convert internal ai to external a
1195         emovz (ai, bi)          bi = ai, low guard word of bi = 0
1196         emulm (ai, bi)          multiply significands, bi = bi * ai
1197         enormlz (ei)            left-justify the significand
1198         eshdn1 (ai)             shift significand and guards down 1 bit
1199         eshdn8 (ai)             shift down 8 bits
1200         eshdn6 (ai)             shift down 16 bits
1201         eshift (ai, n)          shift ai n bits up (or down if n < 0)
1202         eshup1 (ai)             shift significand and guards up 1 bit
1203         eshup8 (ai)             shift up 8 bits
1204         eshup6 (ai)             shift up 16 bits
1205         esubm (ai, bi)          subtract significands, bi = bi - ai
1206         eiisinf (ai)            1 if infinite
1207         eiisnan (ai)            1 if a NaN
1208         eiisneg (ai)            1 if sign bit of ai != 0, else 0
1209         einan (ai)              set ai = NaN
1210         eiinfin (ai)            set ai = infinity
1211
1212   The result is always normalized and rounded to NI-4 word precision
1213   after each arithmetic operation.
1214
1215   Exception flags are NOT fully supported.
1216
1217   Signaling NaN's are NOT supported; they are treated the same
1218   as quiet NaN's.
1219
1220   Define INFINITY for support of infinity; otherwise a
1221   saturation arithmetic is implemented.
1222
1223   Define NANS for support of Not-a-Number items; otherwise the
1224   arithmetic will never produce a NaN output, and might be confused
1225   by a NaN input.
1226   If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1227   either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1228   may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1229   if in doubt.
1230
1231   Denormals are always supported here where appropriate (e.g., not
1232   for conversion to DEC numbers).  */
1233
1234 /* Definitions for error codes that are passed to the common error handling
1235    routine mtherr.
1236
1237    For Digital Equipment PDP-11 and VAX computers, certain
1238   IBM systems, and others that use numbers with a 56-bit
1239   significand, the symbol DEC should be defined.  In this
1240   mode, most floating point constants are given as arrays
1241   of octal integers to eliminate decimal to binary conversion
1242   errors that might be introduced by the compiler.
1243
1244   For computers, such as IBM PC, that follow the IEEE
1245   Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1246   Std 754-1985), the symbol IEEE should be defined.
1247   These numbers have 53-bit significands.  In this mode, constants
1248   are provided as arrays of hexadecimal 16 bit integers.
1249   The endian-ness of generated values is controlled by
1250   REAL_WORDS_BIG_ENDIAN.
1251
1252   To accommodate other types of computer arithmetic, all
1253   constants are also provided in a normal decimal radix
1254   which one can hope are correctly converted to a suitable
1255   format by the available C language compiler.  To invoke
1256   this mode, the symbol UNK is defined.
1257
1258   An important difference among these modes is a predefined
1259   set of machine arithmetic constants for each.  The numbers
1260   MACHEP (the machine roundoff error), MAXNUM (largest number
1261   represented), and several other parameters are preset by
1262   the configuration symbol.  Check the file const.c to
1263   ensure that these values are correct for your computer.
1264
1265   For ANSI C compatibility, define ANSIC equal to 1.  Currently
1266   this affects only the atan2 function and others that use it. */
1267
1268 /* Constant definitions for math error conditions.  */
1269
1270 #define DOMAIN          1       /* argument domain error */
1271 #define SING            2       /* argument singularity */
1272 #define OVERFLOW        3       /* overflow range error */
1273 #define UNDERFLOW       4       /* underflow range error */
1274 #define TLOSS           5       /* total loss of precision */
1275 #define PLOSS           6       /* partial loss of precision */
1276 #define INVALID         7       /* NaN-producing operation */
1277
1278 /*  e type constants used by high precision check routines */
1279
1280 #if LONG_DOUBLE_TYPE_SIZE == 128
1281 /* 0.0 */
1282 unsigned EMUSHORT ezero[NE] =
1283  {0x0000, 0x0000, 0x0000, 0x0000,
1284   0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1285 extern unsigned EMUSHORT ezero[];
1286
1287 /* 5.0E-1 */
1288 unsigned EMUSHORT ehalf[NE] =
1289  {0x0000, 0x0000, 0x0000, 0x0000,
1290   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1291 extern unsigned EMUSHORT ehalf[];
1292
1293 /* 1.0E0 */
1294 unsigned EMUSHORT eone[NE] =
1295  {0x0000, 0x0000, 0x0000, 0x0000,
1296   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1297 extern unsigned EMUSHORT eone[];
1298
1299 /* 2.0E0 */
1300 unsigned EMUSHORT etwo[NE] =
1301  {0x0000, 0x0000, 0x0000, 0x0000,
1302   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1303 extern unsigned EMUSHORT etwo[];
1304
1305 /* 3.2E1 */
1306 unsigned EMUSHORT e32[NE] =
1307  {0x0000, 0x0000, 0x0000, 0x0000,
1308   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1309 extern unsigned EMUSHORT e32[];
1310
1311 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1312 unsigned EMUSHORT elog2[NE] =
1313  {0x40f3, 0xf6af, 0x03f2, 0xb398,
1314   0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1315 extern unsigned EMUSHORT elog2[];
1316
1317 /* 1.41421356237309504880168872420969807856967187537695E0 */
1318 unsigned EMUSHORT esqrt2[NE] =
1319  {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1320   0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1321 extern unsigned EMUSHORT esqrt2[];
1322
1323 /* 3.14159265358979323846264338327950288419716939937511E0 */
1324 unsigned EMUSHORT epi[NE] =
1325  {0x2902, 0x1cd1, 0x80dc, 0x628b,
1326   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1327 extern unsigned EMUSHORT epi[];
1328
1329 #else
1330 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1331 unsigned EMUSHORT ezero[NE] =
1332  {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1333 unsigned EMUSHORT ehalf[NE] =
1334  {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1335 unsigned EMUSHORT eone[NE] =
1336  {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1337 unsigned EMUSHORT etwo[NE] =
1338  {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1339 unsigned EMUSHORT e32[NE] =
1340  {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1341 unsigned EMUSHORT elog2[NE] =
1342  {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1343 unsigned EMUSHORT esqrt2[NE] =
1344  {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1345 unsigned EMUSHORT epi[NE] =
1346  {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1347 #endif
1348
1349 /* Control register for rounding precision.
1350    This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits.  */
1351
1352 int rndprc = NBITS;
1353 extern int rndprc;
1354
1355 /*  Clear out entire e-type number X.  */
1356
1357 static void
1358 eclear (x)
1359      register unsigned EMUSHORT *x;
1360 {
1361   register int i;
1362
1363   for (i = 0; i < NE; i++)
1364     *x++ = 0;
1365 }
1366
1367 /* Move e-type number from A to B.  */
1368
1369 static void
1370 emov (a, b)
1371      register unsigned EMUSHORT *a, *b;
1372 {
1373   register int i;
1374
1375   for (i = 0; i < NE; i++)
1376     *b++ = *a++;
1377 }
1378
1379
1380 /* Absolute value of e-type X.  */
1381
1382 static void
1383 eabs (x)
1384      unsigned EMUSHORT x[];
1385 {
1386   /* sign is top bit of last word of external format */
1387   x[NE - 1] &= 0x7fff;
1388 }
1389
1390 /* Negate the e-type number X.  */
1391
1392 static void
1393 eneg (x)
1394      unsigned EMUSHORT x[];
1395 {
1396
1397   x[NE - 1] ^= 0x8000;          /* Toggle the sign bit */
1398 }
1399
1400 /* Return 1 if sign bit of e-type number X is nonzero, else zero.  */
1401
1402 static int
1403 eisneg (x)
1404      unsigned EMUSHORT x[];
1405 {
1406
1407   if (x[NE - 1] & 0x8000)
1408     return (1);
1409   else
1410     return (0);
1411 }
1412
1413 /* Return 1 if e-type number X is infinity, else return zero.  */
1414
1415 static int
1416 eisinf (x)
1417      unsigned EMUSHORT x[];
1418 {
1419
1420 #ifdef NANS
1421   if (eisnan (x))
1422     return (0);
1423 #endif
1424   if ((x[NE - 1] & 0x7fff) == 0x7fff)
1425     return (1);
1426   else
1427     return (0);
1428 }
1429
1430 /* Check if e-type number is not a number.  The bit pattern is one that we
1431    defined, so we know for sure how to detect it.  */
1432
1433 static int
1434 eisnan (x)
1435      unsigned EMUSHORT x[];
1436 {
1437 #ifdef NANS
1438   int i;
1439
1440   /* NaN has maximum exponent */
1441   if ((x[NE - 1] & 0x7fff) != 0x7fff)
1442     return (0);
1443   /* ... and non-zero significand field. */
1444   for (i = 0; i < NE - 1; i++)
1445     {
1446       if (*x++ != 0)
1447         return (1);
1448     }
1449 #endif
1450
1451   return (0);
1452 }
1453
1454 /*  Fill e-type number X with infinity pattern (IEEE)
1455     or largest possible number (non-IEEE). */
1456
1457 static void
1458 einfin (x)
1459      register unsigned EMUSHORT *x;
1460 {
1461   register int i;
1462
1463 #ifdef INFINITY
1464   for (i = 0; i < NE - 1; i++)
1465     *x++ = 0;
1466   *x |= 32767;
1467 #else
1468   for (i = 0; i < NE - 1; i++)
1469     *x++ = 0xffff;
1470   *x |= 32766;
1471   if (rndprc < NBITS)
1472     {
1473       if (rndprc == 113)
1474         {
1475           *(x - 9) = 0;
1476           *(x - 8) = 0;
1477         }
1478       if (rndprc == 64)
1479         {
1480           *(x - 5) = 0;
1481         }
1482       if (rndprc == 53)
1483         {
1484           *(x - 4) = 0xf800;
1485         }
1486       else
1487         {
1488           *(x - 4) = 0;
1489           *(x - 3) = 0;
1490           *(x - 2) = 0xff00;
1491         }
1492     }
1493 #endif
1494 }
1495
1496 /* Output an e-type NaN.
1497    This generates Intel's quiet NaN pattern for extended real.
1498    The exponent is 7fff, the leading mantissa word is c000.  */
1499
1500 static void
1501 enan (x, sign)
1502      register unsigned EMUSHORT *x;
1503      int sign;
1504 {
1505   register int i;
1506
1507   for (i = 0; i < NE - 2; i++)
1508     *x++ = 0;
1509   *x++ = 0xc000;
1510   *x = (sign << 15) | 0x7fff;
1511 }
1512
1513 /* Move in an e-type number A, converting it to exploded e-type B.  */
1514
1515 static void
1516 emovi (a, b)
1517      unsigned EMUSHORT *a, *b;
1518 {
1519   register unsigned EMUSHORT *p, *q;
1520   int i;
1521
1522   q = b;
1523   p = a + (NE - 1);             /* point to last word of external number */
1524   /* get the sign bit */
1525   if (*p & 0x8000)
1526     *q++ = 0xffff;
1527   else
1528     *q++ = 0;
1529   /* get the exponent */
1530   *q = *p--;
1531   *q++ &= 0x7fff;               /* delete the sign bit */
1532 #ifdef INFINITY
1533   if ((*(q - 1) & 0x7fff) == 0x7fff)
1534     {
1535 #ifdef NANS
1536       if (eisnan (a))
1537         {
1538           *q++ = 0;
1539           for (i = 3; i < NI; i++)
1540             *q++ = *p--;
1541           return;
1542         }
1543 #endif
1544
1545       for (i = 2; i < NI; i++)
1546         *q++ = 0;
1547       return;
1548     }
1549 #endif
1550
1551   /* clear high guard word */
1552   *q++ = 0;
1553   /* move in the significand */
1554   for (i = 0; i < NE - 1; i++)
1555     *q++ = *p--;
1556   /* clear low guard word */
1557   *q = 0;
1558 }
1559
1560 /* Move out exploded e-type number A, converting it to e type B.  */
1561
1562 static void
1563 emovo (a, b)
1564      unsigned EMUSHORT *a, *b;
1565 {
1566   register unsigned EMUSHORT *p, *q;
1567   unsigned EMUSHORT i;
1568   int j;
1569
1570   p = a;
1571   q = b + (NE - 1);             /* point to output exponent */
1572   /* combine sign and exponent */
1573   i = *p++;
1574   if (i)
1575     *q-- = *p++ | 0x8000;
1576   else
1577     *q-- = *p++;
1578 #ifdef INFINITY
1579   if (*(p - 1) == 0x7fff)
1580     {
1581 #ifdef NANS
1582       if (eiisnan (a))
1583         {
1584           enan (b, eiisneg (a));
1585           return;
1586         }
1587 #endif
1588       einfin (b);
1589         return;
1590     }
1591 #endif
1592   /* skip over guard word */
1593   ++p;
1594   /* move the significand */
1595   for (j = 0; j < NE - 1; j++)
1596     *q-- = *p++;
1597 }
1598
1599 /* Clear out exploded e-type number XI.  */
1600
1601 static void
1602 ecleaz (xi)
1603      register unsigned EMUSHORT *xi;
1604 {
1605   register int i;
1606
1607   for (i = 0; i < NI; i++)
1608     *xi++ = 0;
1609 }
1610
1611 /* Clear out exploded e-type XI, but don't touch the sign. */
1612
1613 static void
1614 ecleazs (xi)
1615      register unsigned EMUSHORT *xi;
1616 {
1617   register int i;
1618
1619   ++xi;
1620   for (i = 0; i < NI - 1; i++)
1621     *xi++ = 0;
1622 }
1623
1624 /* Move exploded e-type number from A to B.  */
1625
1626 static void
1627 emovz (a, b)
1628      register unsigned EMUSHORT *a, *b;
1629 {
1630   register int i;
1631
1632   for (i = 0; i < NI - 1; i++)
1633     *b++ = *a++;
1634   /* clear low guard word */
1635   *b = 0;
1636 }
1637
1638 /* Generate exploded e-type NaN.
1639    The explicit pattern for this is maximum exponent and
1640    top two significant bits set.  */
1641
1642 static void
1643 einan (x)
1644      unsigned EMUSHORT x[];
1645 {
1646
1647   ecleaz (x);
1648   x[E] = 0x7fff;
1649   x[M + 1] = 0xc000;
1650 }
1651
1652 /* Return nonzero if exploded e-type X is a NaN. */
1653
1654 static int
1655 eiisnan (x)
1656      unsigned EMUSHORT x[];
1657 {
1658   int i;
1659
1660   if ((x[E] & 0x7fff) == 0x7fff)
1661     {
1662       for (i = M + 1; i < NI; i++)
1663         {
1664           if (x[i] != 0)
1665             return (1);
1666         }
1667     }
1668   return (0);
1669 }
1670
1671 /* Return nonzero if sign of exploded e-type X is nonzero.  */
1672
1673 static int
1674 eiisneg (x)
1675      unsigned EMUSHORT x[];
1676 {
1677
1678   return x[0] != 0;
1679 }
1680
1681 /* Fill exploded e-type X with infinity pattern.
1682    This has maximum exponent and significand all zeros.  */
1683
1684 static void
1685 eiinfin (x)
1686      unsigned EMUSHORT x[];
1687 {
1688
1689   ecleaz (x);
1690   x[E] = 0x7fff;
1691 }
1692
1693 /* Return nonzero if exploded e-type X is infinite. */
1694
1695 static int
1696 eiisinf (x)
1697      unsigned EMUSHORT x[];
1698 {
1699
1700 #ifdef NANS
1701   if (eiisnan (x))
1702     return (0);
1703 #endif
1704   if ((x[E] & 0x7fff) == 0x7fff)
1705     return (1);
1706   return (0);
1707 }
1708
1709
1710 /* Compare significands of numbers in internal exploded e-type format.
1711    Guard words are included in the comparison.
1712
1713    Returns      +1 if a > b
1714                  0 if a == b
1715                 -1 if a < b   */
1716
1717 static int
1718 ecmpm (a, b)
1719      register unsigned EMUSHORT *a, *b;
1720 {
1721   int i;
1722
1723   a += M;                       /* skip up to significand area */
1724   b += M;
1725   for (i = M; i < NI; i++)
1726     {
1727       if (*a++ != *b++)
1728         goto difrnt;
1729     }
1730   return (0);
1731
1732  difrnt:
1733   if (*(--a) > *(--b))
1734     return (1);
1735   else
1736     return (-1);
1737 }
1738
1739 /* Shift significand of exploded e-type X down by 1 bit.  */
1740
1741 static void
1742 eshdn1 (x)
1743      register unsigned EMUSHORT *x;
1744 {
1745   register unsigned EMUSHORT bits;
1746   int i;
1747
1748   x += M;                       /* point to significand area */
1749
1750   bits = 0;
1751   for (i = M; i < NI; i++)
1752     {
1753       if (*x & 1)
1754         bits |= 1;
1755       *x >>= 1;
1756       if (bits & 2)
1757         *x |= 0x8000;
1758       bits <<= 1;
1759       ++x;
1760     }
1761 }
1762
1763 /* Shift significand of exploded e-type X up by 1 bit.  */
1764
1765 static void
1766 eshup1 (x)
1767      register unsigned EMUSHORT *x;
1768 {
1769   register unsigned EMUSHORT bits;
1770   int i;
1771
1772   x += NI - 1;
1773   bits = 0;
1774
1775   for (i = M; i < NI; i++)
1776     {
1777       if (*x & 0x8000)
1778         bits |= 1;
1779       *x <<= 1;
1780       if (bits & 2)
1781         *x |= 1;
1782       bits <<= 1;
1783       --x;
1784     }
1785 }
1786
1787
1788 /* Shift significand of exploded e-type X down by 8 bits.  */
1789
1790 static void
1791 eshdn8 (x)
1792      register unsigned EMUSHORT *x;
1793 {
1794   register unsigned EMUSHORT newbyt, oldbyt;
1795   int i;
1796
1797   x += M;
1798   oldbyt = 0;
1799   for (i = M; i < NI; i++)
1800     {
1801       newbyt = *x << 8;
1802       *x >>= 8;
1803       *x |= oldbyt;
1804       oldbyt = newbyt;
1805       ++x;
1806     }
1807 }
1808
1809 /* Shift significand of exploded e-type X up by 8 bits.  */
1810
1811 static void
1812 eshup8 (x)
1813      register unsigned EMUSHORT *x;
1814 {
1815   int i;
1816   register unsigned EMUSHORT newbyt, oldbyt;
1817
1818   x += NI - 1;
1819   oldbyt = 0;
1820
1821   for (i = M; i < NI; i++)
1822     {
1823       newbyt = *x >> 8;
1824       *x <<= 8;
1825       *x |= oldbyt;
1826       oldbyt = newbyt;
1827       --x;
1828     }
1829 }
1830
1831 /* Shift significand of exploded e-type X up by 16 bits.  */
1832
1833 static void
1834 eshup6 (x)
1835      register unsigned EMUSHORT *x;
1836 {
1837   int i;
1838   register unsigned EMUSHORT *p;
1839
1840   p = x + M;
1841   x += M + 1;
1842
1843   for (i = M; i < NI - 1; i++)
1844     *p++ = *x++;
1845
1846   *p = 0;
1847 }
1848
1849 /* Shift significand of exploded e-type X down by 16 bits.  */
1850
1851 static void
1852 eshdn6 (x)
1853      register unsigned EMUSHORT *x;
1854 {
1855   int i;
1856   register unsigned EMUSHORT *p;
1857
1858   x += NI - 1;
1859   p = x + 1;
1860
1861   for (i = M; i < NI - 1; i++)
1862     *(--p) = *(--x);
1863
1864   *(--p) = 0;
1865 }
1866
1867 /* Add significands of exploded e-type X and Y.  X + Y replaces Y.  */
1868
1869 static void
1870 eaddm (x, y)
1871      unsigned EMUSHORT *x, *y;
1872 {
1873   register unsigned EMULONG a;
1874   int i;
1875   unsigned int carry;
1876
1877   x += NI - 1;
1878   y += NI - 1;
1879   carry = 0;
1880   for (i = M; i < NI; i++)
1881     {
1882       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1883       if (a & 0x10000)
1884         carry = 1;
1885       else
1886         carry = 0;
1887       *y = (unsigned EMUSHORT) a;
1888       --x;
1889       --y;
1890     }
1891 }
1892
1893 /* Subtract significands of exploded e-type X and Y.  Y - X replaces Y.  */
1894
1895 static void
1896 esubm (x, y)
1897      unsigned EMUSHORT *x, *y;
1898 {
1899   unsigned EMULONG a;
1900   int i;
1901   unsigned int carry;
1902
1903   x += NI - 1;
1904   y += NI - 1;
1905   carry = 0;
1906   for (i = M; i < NI; i++)
1907     {
1908       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1909       if (a & 0x10000)
1910         carry = 1;
1911       else
1912         carry = 0;
1913       *y = (unsigned EMUSHORT) a;
1914       --x;
1915       --y;
1916     }
1917 }
1918
1919
1920 static unsigned EMUSHORT equot[NI];
1921
1922
1923 #if 0
1924 /* Radix 2 shift-and-add versions of multiply and divide  */
1925
1926
1927 /* Divide significands */
1928
1929 int
1930 edivm (den, num)
1931      unsigned EMUSHORT den[], num[];
1932 {
1933   int i;
1934   register unsigned EMUSHORT *p, *q;
1935   unsigned EMUSHORT j;
1936
1937   p = &equot[0];
1938   *p++ = num[0];
1939   *p++ = num[1];
1940
1941   for (i = M; i < NI; i++)
1942     {
1943       *p++ = 0;
1944     }
1945
1946   /* Use faster compare and subtraction if denominator has only 15 bits of
1947      significance.  */
1948
1949   p = &den[M + 2];
1950   if (*p++ == 0)
1951     {
1952       for (i = M + 3; i < NI; i++)
1953         {
1954           if (*p++ != 0)
1955             goto fulldiv;
1956         }
1957       if ((den[M + 1] & 1) != 0)
1958         goto fulldiv;
1959       eshdn1 (num);
1960       eshdn1 (den);
1961
1962       p = &den[M + 1];
1963       q = &num[M + 1];
1964
1965       for (i = 0; i < NBITS + 2; i++)
1966         {
1967           if (*p <= *q)
1968             {
1969               *q -= *p;
1970               j = 1;
1971             }
1972           else
1973             {
1974               j = 0;
1975             }
1976           eshup1 (equot);
1977           equot[NI - 2] |= j;
1978           eshup1 (num);
1979         }
1980       goto divdon;
1981     }
1982
1983   /* The number of quotient bits to calculate is NBITS + 1 scaling guard
1984      bit + 1 roundoff bit.  */
1985
1986  fulldiv:
1987
1988   p = &equot[NI - 2];
1989   for (i = 0; i < NBITS + 2; i++)
1990     {
1991       if (ecmpm (den, num) <= 0)
1992         {
1993           esubm (den, num);
1994           j = 1;                /* quotient bit = 1 */
1995         }
1996       else
1997         j = 0;
1998       eshup1 (equot);
1999       *p |= j;
2000       eshup1 (num);
2001     }
2002
2003  divdon:
2004
2005   eshdn1 (equot);
2006   eshdn1 (equot);
2007
2008   /* test for nonzero remainder after roundoff bit */
2009   p = &num[M];
2010   j = 0;
2011   for (i = M; i < NI; i++)
2012     {
2013       j |= *p++;
2014     }
2015   if (j)
2016     j = 1;
2017
2018
2019   for (i = 0; i < NI; i++)
2020     num[i] = equot[i];
2021   return ((int) j);
2022 }
2023
2024
2025 /* Multiply significands */
2026 int
2027 emulm (a, b)
2028      unsigned EMUSHORT a[], b[];
2029 {
2030   unsigned EMUSHORT *p, *q;
2031   int i, j, k;
2032
2033   equot[0] = b[0];
2034   equot[1] = b[1];
2035   for (i = M; i < NI; i++)
2036     equot[i] = 0;
2037
2038   p = &a[NI - 2];
2039   k = NBITS;
2040   while (*p == 0)               /* significand is not supposed to be zero */
2041     {
2042       eshdn6 (a);
2043       k -= 16;
2044     }
2045   if ((*p & 0xff) == 0)
2046     {
2047       eshdn8 (a);
2048       k -= 8;
2049     }
2050
2051   q = &equot[NI - 1];
2052   j = 0;
2053   for (i = 0; i < k; i++)
2054     {
2055       if (*p & 1)
2056         eaddm (b, equot);
2057       /* remember if there were any nonzero bits shifted out */
2058       if (*q & 1)
2059         j |= 1;
2060       eshdn1 (a);
2061       eshdn1 (equot);
2062     }
2063
2064   for (i = 0; i < NI; i++)
2065     b[i] = equot[i];
2066
2067   /* return flag for lost nonzero bits */
2068   return (j);
2069 }
2070
2071 #else
2072
2073 /* Radix 65536 versions of multiply and divide.  */
2074
2075 /* Multiply significand of e-type number B
2076    by 16-bit quantity A, return e-type result to C. */
2077
2078 static void
2079 m16m (a, b, c)
2080      unsigned int a;
2081      unsigned EMUSHORT b[], c[];
2082 {
2083   register unsigned EMUSHORT *pp;
2084   register unsigned EMULONG carry;
2085   unsigned EMUSHORT *ps;
2086   unsigned EMUSHORT p[NI];
2087   unsigned EMULONG aa, m;
2088   int i;
2089
2090   aa = a;
2091   pp = &p[NI-2];
2092   *pp++ = 0;
2093   *pp = 0;
2094   ps = &b[NI-1];
2095
2096   for (i=M+1; i<NI; i++)
2097     {
2098       if (*ps == 0)
2099         {
2100           --ps;
2101           --pp;
2102           *(pp-1) = 0;
2103         }
2104       else
2105         {
2106           m = (unsigned EMULONG) aa * *ps--;
2107           carry = (m & 0xffff) + *pp;
2108           *pp-- = (unsigned EMUSHORT)carry;
2109           carry = (carry >> 16) + (m >> 16) + *pp;
2110           *pp = (unsigned EMUSHORT)carry;
2111           *(pp-1) = carry >> 16;
2112         }
2113     }
2114   for (i=M; i<NI; i++)
2115     c[i] = p[i];
2116 }
2117
2118 /* Divide significands of exploded e-types NUM / DEN.  Neither the
2119    numerator NUM nor the denominator DEN is permitted to have its high guard
2120    word nonzero.  */
2121
2122 static int
2123 edivm (den, num)
2124      unsigned EMUSHORT den[], num[];
2125 {
2126   int i;
2127   register unsigned EMUSHORT *p;
2128   unsigned EMULONG tnum;
2129   unsigned EMUSHORT j, tdenm, tquot;
2130   unsigned EMUSHORT tprod[NI+1];
2131
2132   p = &equot[0];
2133   *p++ = num[0];
2134   *p++ = num[1];
2135
2136   for (i=M; i<NI; i++)
2137     {
2138       *p++ = 0;
2139     }
2140   eshdn1 (num);
2141   tdenm = den[M+1];
2142   for (i=M; i<NI; i++)
2143     {
2144       /* Find trial quotient digit (the radix is 65536). */
2145       tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
2146
2147       /* Do not execute the divide instruction if it will overflow. */
2148       if ((tdenm * 0xffffL) < tnum)
2149         tquot = 0xffff;
2150       else
2151         tquot = tnum / tdenm;
2152       /* Multiply denominator by trial quotient digit. */
2153       m16m ((unsigned int)tquot, den, tprod);
2154       /* The quotient digit may have been overestimated. */
2155       if (ecmpm (tprod, num) > 0)
2156         {
2157           tquot -= 1;
2158           esubm (den, tprod);
2159           if (ecmpm (tprod, num) > 0)
2160             {
2161               tquot -= 1;
2162               esubm (den, tprod);
2163             }
2164         }
2165       esubm (tprod, num);
2166       equot[i] = tquot;
2167       eshup6(num);
2168     }
2169   /* test for nonzero remainder after roundoff bit */
2170   p = &num[M];
2171   j = 0;
2172   for (i=M; i<NI; i++)
2173     {
2174       j |= *p++;
2175     }
2176   if (j)
2177     j = 1;
2178
2179   for (i=0; i<NI; i++)
2180     num[i] = equot[i];
2181
2182   return ((int)j);
2183 }
2184
2185 /* Multiply significands of exploded e-type A and B, result in B.  */
2186
2187 static int
2188 emulm (a, b)
2189      unsigned EMUSHORT a[], b[];
2190 {
2191   unsigned EMUSHORT *p, *q;
2192   unsigned EMUSHORT pprod[NI];
2193   unsigned EMUSHORT j;
2194   int i;
2195
2196   equot[0] = b[0];
2197   equot[1] = b[1];
2198   for (i=M; i<NI; i++)
2199     equot[i] = 0;
2200
2201   j = 0;
2202   p = &a[NI-1];
2203   q = &equot[NI-1];
2204   for (i=M+1; i<NI; i++)
2205     {
2206       if (*p == 0)
2207         {
2208           --p;
2209         }
2210       else
2211         {
2212           m16m ((unsigned int) *p--, b, pprod);
2213           eaddm(pprod, equot);
2214         }
2215       j |= *q;
2216       eshdn6(equot);
2217     }
2218
2219   for (i=0; i<NI; i++)
2220     b[i] = equot[i];
2221
2222   /* return flag for lost nonzero bits */
2223   return ((int)j);
2224 }
2225 #endif
2226
2227
2228 /* Normalize and round off.
2229
2230   The internal format number to be rounded is S.
2231   Input LOST is 0 if the value is exact.  This is the so-called sticky bit.
2232
2233   Input SUBFLG indicates whether the number was obtained
2234   by a subtraction operation.  In that case if LOST is nonzero
2235   then the number is slightly smaller than indicated.
2236
2237   Input EXP is the biased exponent, which may be negative.
2238   the exponent field of S is ignored but is replaced by
2239   EXP as adjusted by normalization and rounding.
2240
2241   Input RCNTRL is the rounding control.  If it is nonzero, the
2242   returned value will be rounded to RNDPRC bits.
2243
2244   For future reference:  In order for emdnorm to round off denormal
2245    significands at the right point, the input exponent must be
2246    adjusted to be the actual value it would have after conversion to
2247    the final floating point type.  This adjustment has been
2248    implemented for all type conversions (etoe53, etc.) and decimal
2249    conversions, but not for the arithmetic functions (eadd, etc.).
2250    Data types having standard 15-bit exponents are not affected by
2251    this, but SFmode and DFmode are affected. For example, ediv with
2252    rndprc = 24 will not round correctly to 24-bit precision if the
2253    result is denormal.   */
2254
2255 static int rlast = -1;
2256 static int rw = 0;
2257 static unsigned EMUSHORT rmsk = 0;
2258 static unsigned EMUSHORT rmbit = 0;
2259 static unsigned EMUSHORT rebit = 0;
2260 static int re = 0;
2261 static unsigned EMUSHORT rbit[NI];
2262
2263 static void
2264 emdnorm (s, lost, subflg, exp, rcntrl)
2265      unsigned EMUSHORT s[];
2266      int lost;
2267      int subflg;
2268      EMULONG exp;
2269      int rcntrl;
2270 {
2271   int i, j;
2272   unsigned EMUSHORT r;
2273
2274   /* Normalize */
2275   j = enormlz (s);
2276
2277   /* a blank significand could mean either zero or infinity. */
2278 #ifndef INFINITY
2279   if (j > NBITS)
2280     {
2281       ecleazs (s);
2282       return;
2283     }
2284 #endif
2285   exp -= j;
2286 #ifndef INFINITY
2287   if (exp >= 32767L)
2288     goto overf;
2289 #else
2290   if ((j > NBITS) && (exp < 32767))
2291     {
2292       ecleazs (s);
2293       return;
2294     }
2295 #endif
2296   if (exp < 0L)
2297     {
2298       if (exp > (EMULONG) (-NBITS - 1))
2299         {
2300           j = (int) exp;
2301           i = eshift (s, j);
2302           if (i)
2303             lost = 1;
2304         }
2305       else
2306         {
2307           ecleazs (s);
2308           return;
2309         }
2310     }
2311   /* Round off, unless told not to by rcntrl. */
2312   if (rcntrl == 0)
2313     goto mdfin;
2314   /* Set up rounding parameters if the control register changed. */
2315   if (rndprc != rlast)
2316     {
2317       ecleaz (rbit);
2318       switch (rndprc)
2319         {
2320         default:
2321         case NBITS:
2322           rw = NI - 1;          /* low guard word */
2323           rmsk = 0xffff;
2324           rmbit = 0x8000;
2325           re = rw - 1;
2326           rebit = 1;
2327           break;
2328         case 113:
2329           rw = 10;
2330           rmsk = 0x7fff;
2331           rmbit = 0x4000;
2332           rebit = 0x8000;
2333           re = rw;
2334           break;
2335         case 64:
2336           rw = 7;
2337           rmsk = 0xffff;
2338           rmbit = 0x8000;
2339           re = rw - 1;
2340           rebit = 1;
2341           break;
2342           /* For DEC or IBM arithmetic */
2343         case 56:
2344           rw = 6;
2345           rmsk = 0xff;
2346           rmbit = 0x80;
2347           rebit = 0x100;
2348           re = rw;
2349           break;
2350         case 53:
2351           rw = 6;
2352           rmsk = 0x7ff;
2353           rmbit = 0x0400;
2354           rebit = 0x800;
2355           re = rw;
2356           break;
2357         case 24:
2358           rw = 4;
2359           rmsk = 0xff;
2360           rmbit = 0x80;
2361           rebit = 0x100;
2362           re = rw;
2363           break;
2364         }
2365       rbit[re] = rebit;
2366       rlast = rndprc;
2367     }
2368
2369   /* Shift down 1 temporarily if the data structure has an implied
2370      most significant bit and the number is denormal.  */
2371   if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
2372     {
2373       lost |= s[NI - 1] & 1;
2374       eshdn1 (s);
2375     }
2376   /* Clear out all bits below the rounding bit,
2377      remembering in r if any were nonzero.  */
2378   r = s[rw] & rmsk;
2379   if (rndprc < NBITS)
2380     {
2381       i = rw + 1;
2382       while (i < NI)
2383         {
2384           if (s[i])
2385             r |= 1;
2386           s[i] = 0;
2387           ++i;
2388         }
2389     }
2390   s[rw] &= ~rmsk;
2391   if ((r & rmbit) != 0)
2392     {
2393       if (r == rmbit)
2394         {
2395           if (lost == 0)
2396             {                   /* round to even */
2397               if ((s[re] & rebit) == 0)
2398                 goto mddone;
2399             }
2400           else
2401             {
2402               if (subflg != 0)
2403                 goto mddone;
2404             }
2405         }
2406       eaddm (rbit, s);
2407     }
2408  mddone:
2409   if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
2410     {
2411       eshup1 (s);
2412     }
2413   if (s[2] != 0)
2414     {                           /* overflow on roundoff */
2415       eshdn1 (s);
2416       exp += 1;
2417     }
2418  mdfin:
2419   s[NI - 1] = 0;
2420   if (exp >= 32767L)
2421     {
2422 #ifndef INFINITY
2423     overf:
2424 #endif
2425 #ifdef INFINITY
2426       s[1] = 32767;
2427       for (i = 2; i < NI - 1; i++)
2428         s[i] = 0;
2429       if (extra_warnings)
2430         warning ("floating point overflow");
2431 #else
2432       s[1] = 32766;
2433       s[2] = 0;
2434       for (i = M + 1; i < NI - 1; i++)
2435         s[i] = 0xffff;
2436       s[NI - 1] = 0;
2437       if ((rndprc < 64) || (rndprc == 113))
2438         {
2439           s[rw] &= ~rmsk;
2440           if (rndprc == 24)
2441             {
2442               s[5] = 0;
2443               s[6] = 0;
2444             }
2445         }
2446 #endif
2447       return;
2448     }
2449   if (exp < 0)
2450     s[1] = 0;
2451   else
2452     s[1] = (unsigned EMUSHORT) exp;
2453 }
2454
2455 /*  Subtract.  C = B - A, all e type numbers.  */
2456
2457 static int subflg = 0;
2458
2459 static void
2460 esub (a, b, c)
2461      unsigned EMUSHORT *a, *b, *c;
2462 {
2463
2464 #ifdef NANS
2465   if (eisnan (a))
2466     {
2467       emov (a, c);
2468       return;
2469     }
2470   if (eisnan (b))
2471     {
2472       emov (b, c);
2473       return;
2474     }
2475 /* Infinity minus infinity is a NaN.
2476    Test for subtracting infinities of the same sign. */
2477   if (eisinf (a) && eisinf (b)
2478       && ((eisneg (a) ^ eisneg (b)) == 0))
2479     {
2480       mtherr ("esub", INVALID);
2481       enan (c, 0);
2482       return;
2483     }
2484 #endif
2485   subflg = 1;
2486   eadd1 (a, b, c);
2487 }
2488
2489 /* Add.  C = A + B, all e type. */
2490
2491 static void
2492 eadd (a, b, c)
2493      unsigned EMUSHORT *a, *b, *c;
2494 {
2495
2496 #ifdef NANS
2497 /* NaN plus anything is a NaN. */
2498   if (eisnan (a))
2499     {
2500       emov (a, c);
2501       return;
2502     }
2503   if (eisnan (b))
2504     {
2505       emov (b, c);
2506       return;
2507     }
2508 /* Infinity minus infinity is a NaN.
2509    Test for adding infinities of opposite signs. */
2510   if (eisinf (a) && eisinf (b)
2511       && ((eisneg (a) ^ eisneg (b)) != 0))
2512     {
2513       mtherr ("esub", INVALID);
2514       enan (c, 0);
2515       return;
2516     }
2517 #endif
2518   subflg = 0;
2519   eadd1 (a, b, c);
2520 }
2521
2522 /* Arithmetic common to both addition and subtraction.  */
2523
2524 static void
2525 eadd1 (a, b, c)
2526      unsigned EMUSHORT *a, *b, *c;
2527 {
2528   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2529   int i, lost, j, k;
2530   EMULONG lt, lta, ltb;
2531
2532 #ifdef INFINITY
2533   if (eisinf (a))
2534     {
2535       emov (a, c);
2536       if (subflg)
2537         eneg (c);
2538       return;
2539     }
2540   if (eisinf (b))
2541     {
2542       emov (b, c);
2543       return;
2544     }
2545 #endif
2546   emovi (a, ai);
2547   emovi (b, bi);
2548   if (subflg)
2549     ai[0] = ~ai[0];
2550
2551   /* compare exponents */
2552   lta = ai[E];
2553   ltb = bi[E];
2554   lt = lta - ltb;
2555   if (lt > 0L)
2556     {                           /* put the larger number in bi */
2557       emovz (bi, ci);
2558       emovz (ai, bi);
2559       emovz (ci, ai);
2560       ltb = bi[E];
2561       lt = -lt;
2562     }
2563   lost = 0;
2564   if (lt != 0L)
2565     {
2566       if (lt < (EMULONG) (-NBITS - 1))
2567         goto done;              /* answer same as larger addend */
2568       k = (int) lt;
2569       lost = eshift (ai, k);    /* shift the smaller number down */
2570     }
2571   else
2572     {
2573       /* exponents were the same, so must compare significands */
2574       i = ecmpm (ai, bi);
2575       if (i == 0)
2576         {                       /* the numbers are identical in magnitude */
2577           /* if different signs, result is zero */
2578           if (ai[0] != bi[0])
2579             {
2580               eclear (c);
2581               return;
2582             }
2583           /* if same sign, result is double */
2584           /* double denomalized tiny number */
2585           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2586             {
2587               eshup1 (bi);
2588               goto done;
2589             }
2590           /* add 1 to exponent unless both are zero! */
2591           for (j = 1; j < NI - 1; j++)
2592             {
2593               if (bi[j] != 0)
2594                 {
2595                   /* This could overflow, but let emovo take care of that. */
2596                   ltb += 1;
2597                   break;
2598                 }
2599             }
2600           bi[E] = (unsigned EMUSHORT) ltb;
2601           goto done;
2602         }
2603       if (i > 0)
2604         {                       /* put the larger number in bi */
2605           emovz (bi, ci);
2606           emovz (ai, bi);
2607           emovz (ci, ai);
2608         }
2609     }
2610   if (ai[0] == bi[0])
2611     {
2612       eaddm (ai, bi);
2613       subflg = 0;
2614     }
2615   else
2616     {
2617       esubm (ai, bi);
2618       subflg = 1;
2619     }
2620   emdnorm (bi, lost, subflg, ltb, 64);
2621
2622  done:
2623   emovo (bi, c);
2624 }
2625
2626 /* Divide: C = B/A, all e type.  */
2627
2628 static void
2629 ediv (a, b, c)
2630      unsigned EMUSHORT *a, *b, *c;
2631 {
2632   unsigned EMUSHORT ai[NI], bi[NI];
2633   int i;
2634   EMULONG lt, lta, ltb;
2635
2636 #ifdef NANS
2637 /* Return any NaN input. */
2638   if (eisnan (a))
2639     {
2640     emov (a, c);
2641     return;
2642     }
2643   if (eisnan (b))
2644     {
2645     emov (b, c);
2646     return;
2647     }
2648 /* Zero over zero, or infinity over infinity, is a NaN. */
2649   if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
2650       || (eisinf (a) && eisinf (b)))
2651     {
2652     mtherr ("ediv", INVALID);
2653     enan (c, eisneg (a) ^ eisneg (b));
2654     return;
2655     }
2656 #endif
2657 /* Infinity over anything else is infinity. */
2658 #ifdef INFINITY
2659   if (eisinf (b))
2660     {
2661       if (eisneg (a) ^ eisneg (b))
2662         *(c + (NE - 1)) = 0x8000;
2663       else
2664         *(c + (NE - 1)) = 0;
2665       einfin (c);
2666       return;
2667     }
2668 /* Anything else over infinity is zero. */
2669   if (eisinf (a))
2670     {
2671       eclear (c);
2672       return;
2673     }
2674 #endif
2675   emovi (a, ai);
2676   emovi (b, bi);
2677   lta = ai[E];
2678   ltb = bi[E];
2679   if (bi[E] == 0)
2680     {                           /* See if numerator is zero. */
2681       for (i = 1; i < NI - 1; i++)
2682         {
2683           if (bi[i] != 0)
2684             {
2685               ltb -= enormlz (bi);
2686               goto dnzro1;
2687             }
2688         }
2689       eclear (c);
2690       return;
2691     }
2692  dnzro1:
2693
2694   if (ai[E] == 0)
2695     {                           /* possible divide by zero */
2696       for (i = 1; i < NI - 1; i++)
2697         {
2698           if (ai[i] != 0)
2699             {
2700               lta -= enormlz (ai);
2701               goto dnzro2;
2702             }
2703         }
2704       if (ai[0] == bi[0])
2705         *(c + (NE - 1)) = 0;
2706       else
2707         *(c + (NE - 1)) = 0x8000;
2708 /* Divide by zero is not an invalid operation.
2709    It is a divide-by-zero operation!   */
2710       einfin (c);
2711       mtherr ("ediv", SING);
2712       return;
2713     }
2714  dnzro2:
2715
2716   i = edivm (ai, bi);
2717   /* calculate exponent */
2718   lt = ltb - lta + EXONE;
2719   emdnorm (bi, i, 0, lt, 64);
2720   /* set the sign */
2721   if (ai[0] == bi[0])
2722     bi[0] = 0;
2723   else
2724     bi[0] = 0Xffff;
2725   emovo (bi, c);
2726 }
2727
2728 /* Multiply e-types A and B, return e-type product C.   */
2729
2730 static void
2731 emul (a, b, c)
2732      unsigned EMUSHORT *a, *b, *c;
2733 {
2734   unsigned EMUSHORT ai[NI], bi[NI];
2735   int i, j;
2736   EMULONG lt, lta, ltb;
2737
2738 #ifdef NANS
2739 /* NaN times anything is the same NaN. */
2740   if (eisnan (a))
2741     {
2742     emov (a, c);
2743     return;
2744     }
2745   if (eisnan (b))
2746     {
2747     emov (b, c);
2748     return;
2749     }
2750 /* Zero times infinity is a NaN. */
2751   if ((eisinf (a) && (ecmp (b, ezero) == 0))
2752       || (eisinf (b) && (ecmp (a, ezero) == 0)))
2753     {
2754     mtherr ("emul", INVALID);
2755     enan (c, eisneg (a) ^ eisneg (b));
2756     return;
2757     }
2758 #endif
2759 /* Infinity times anything else is infinity. */
2760 #ifdef INFINITY
2761   if (eisinf (a) || eisinf (b))
2762     {
2763       if (eisneg (a) ^ eisneg (b))
2764         *(c + (NE - 1)) = 0x8000;
2765       else
2766         *(c + (NE - 1)) = 0;
2767       einfin (c);
2768       return;
2769     }
2770 #endif
2771   emovi (a, ai);
2772   emovi (b, bi);
2773   lta = ai[E];
2774   ltb = bi[E];
2775   if (ai[E] == 0)
2776     {
2777       for (i = 1; i < NI - 1; i++)
2778         {
2779           if (ai[i] != 0)
2780             {
2781               lta -= enormlz (ai);
2782               goto mnzer1;
2783             }
2784         }
2785       eclear (c);
2786       return;
2787     }
2788  mnzer1:
2789
2790   if (bi[E] == 0)
2791     {
2792       for (i = 1; i < NI - 1; i++)
2793         {
2794           if (bi[i] != 0)
2795             {
2796               ltb -= enormlz (bi);
2797               goto mnzer2;
2798             }
2799         }
2800       eclear (c);
2801       return;
2802     }
2803  mnzer2:
2804
2805   /* Multiply significands */
2806   j = emulm (ai, bi);
2807   /* calculate exponent */
2808   lt = lta + ltb - (EXONE - 1);
2809   emdnorm (bi, j, 0, lt, 64);
2810   /* calculate sign of product */
2811   if (ai[0] == bi[0])
2812     bi[0] = 0;
2813   else
2814     bi[0] = 0xffff;
2815   emovo (bi, c);
2816 }
2817
2818 /* Convert double precision PE to e-type Y.  */
2819
2820 static void
2821 e53toe (pe, y)
2822      unsigned EMUSHORT *pe, *y;
2823 {
2824 #ifdef DEC
2825
2826   dectoe (pe, y);
2827
2828 #else
2829 #ifdef IBM
2830
2831   ibmtoe (pe, y, DFmode);
2832
2833 #else
2834   register unsigned EMUSHORT r;
2835   register unsigned EMUSHORT *e, *p;
2836   unsigned EMUSHORT yy[NI];
2837   int denorm, k;
2838
2839   e = pe;
2840   denorm = 0;                   /* flag if denormalized number */
2841   ecleaz (yy);
2842   if (! REAL_WORDS_BIG_ENDIAN)
2843     e += 3;
2844   r = *e;
2845   yy[0] = 0;
2846   if (r & 0x8000)
2847     yy[0] = 0xffff;
2848   yy[M] = (r & 0x0f) | 0x10;
2849   r &= ~0x800f;                 /* strip sign and 4 significand bits */
2850 #ifdef INFINITY
2851   if (r == 0x7ff0)
2852     {
2853 #ifdef NANS
2854       if (! REAL_WORDS_BIG_ENDIAN)
2855         {
2856           if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
2857               || (pe[1] != 0) || (pe[0] != 0))
2858             {
2859               enan (y, yy[0] != 0);
2860               return;
2861             }
2862         }
2863       else
2864         {
2865           if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
2866               || (pe[2] != 0) || (pe[3] != 0))
2867             {
2868               enan (y, yy[0] != 0);
2869               return;
2870             }
2871         }
2872 #endif  /* NANS */
2873       eclear (y);
2874       einfin (y);
2875       if (yy[0])
2876         eneg (y);
2877       return;
2878     }
2879 #endif  /* INFINITY */
2880   r >>= 4;
2881   /* If zero exponent, then the significand is denormalized.
2882      So take back the understood high significand bit. */
2883
2884   if (r == 0)
2885     {
2886       denorm = 1;
2887       yy[M] &= ~0x10;
2888     }
2889   r += EXONE - 01777;
2890   yy[E] = r;
2891   p = &yy[M + 1];
2892 #ifdef IEEE
2893   if (! REAL_WORDS_BIG_ENDIAN)
2894     {
2895       *p++ = *(--e);
2896       *p++ = *(--e);
2897       *p++ = *(--e);
2898     }
2899   else
2900     {
2901       ++e;
2902       *p++ = *e++;
2903       *p++ = *e++;
2904       *p++ = *e++;
2905     }
2906 #endif
2907   eshift (yy, -5);
2908   if (denorm)
2909     {                           /* if zero exponent, then normalize the significand */
2910       if ((k = enormlz (yy)) > NBITS)
2911         ecleazs (yy);
2912       else
2913         yy[E] -= (unsigned EMUSHORT) (k - 1);
2914     }
2915   emovo (yy, y);
2916 #endif /* not IBM */
2917 #endif /* not DEC */
2918 }
2919
2920 /* Convert double extended precision float PE to e type Y.  */
2921
2922 static void
2923 e64toe (pe, y)
2924      unsigned EMUSHORT *pe, *y;
2925 {
2926   unsigned EMUSHORT yy[NI];
2927   unsigned EMUSHORT *e, *p, *q;
2928   int i;
2929
2930   e = pe;
2931   p = yy;
2932   for (i = 0; i < NE - 5; i++)
2933     *p++ = 0;
2934 /* This precision is not ordinarily supported on DEC or IBM. */
2935 #ifdef DEC
2936   for (i = 0; i < 5; i++)
2937     *p++ = *e++;
2938 #endif
2939 #ifdef IBM
2940   p = &yy[0] + (NE - 1);
2941   *p-- = *e++;
2942   ++e;
2943   for (i = 0; i < 5; i++)
2944     *p-- = *e++;
2945 #endif
2946 #ifdef IEEE
2947   if (! REAL_WORDS_BIG_ENDIAN)
2948     {
2949       for (i = 0; i < 5; i++)
2950         *p++ = *e++;
2951     }
2952   else
2953     {
2954       p = &yy[0] + (NE - 1);
2955       *p-- = *e++;
2956       ++e;
2957       for (i = 0; i < 4; i++)
2958         *p-- = *e++;
2959     }
2960 #endif
2961 #ifdef INFINITY
2962   /* Point to the exponent field and check max exponent cases.  */
2963   p = &yy[NE - 1];
2964   if (*p == 0x7fff)
2965     {
2966 #ifdef NANS
2967       if (! REAL_WORDS_BIG_ENDIAN)
2968         {
2969           for (i = 0; i < 4; i++)
2970             {
2971               if ((i != 3 && pe[i] != 0)
2972                   /* Anything but 0x8000 here, including 0, is a NaN.  */
2973                   || (i == 3 && pe[i] != 0x8000))
2974                 {
2975                   enan (y, (*p & 0x8000) != 0);
2976                   return;
2977                 }
2978             }
2979         }
2980       else
2981         {
2982           for (i = 1; i <= 4; i++)
2983             {
2984               if (pe[i] != 0)
2985                 {
2986                   enan (y, (*p & 0x8000) != 0);
2987                   return;
2988                 }
2989             }
2990         }
2991 #endif /* NANS */
2992       eclear (y);
2993       einfin (y);
2994       if (*p & 0x8000)
2995         eneg (y);
2996       return;
2997     }
2998 #endif  /* INFINITY */
2999   p = yy;
3000   q = y;
3001   for (i = 0; i < NE; i++)
3002     *q++ = *p++;
3003 }
3004
3005 /* Convert 128-bit long double precision float PE to e type Y.  */
3006
3007 static void
3008 e113toe (pe, y)
3009      unsigned EMUSHORT *pe, *y;
3010 {
3011   register unsigned EMUSHORT r;
3012   unsigned EMUSHORT *e, *p;
3013   unsigned EMUSHORT yy[NI];
3014   int denorm, i;
3015
3016   e = pe;
3017   denorm = 0;
3018   ecleaz (yy);
3019 #ifdef IEEE
3020   if (! REAL_WORDS_BIG_ENDIAN)
3021     e += 7;
3022 #endif
3023   r = *e;
3024   yy[0] = 0;
3025   if (r & 0x8000)
3026     yy[0] = 0xffff;
3027   r &= 0x7fff;
3028 #ifdef INFINITY
3029   if (r == 0x7fff)
3030     {
3031 #ifdef NANS
3032       if (! REAL_WORDS_BIG_ENDIAN)
3033         {
3034           for (i = 0; i < 7; i++)
3035             {
3036               if (pe[i] != 0)
3037                 {
3038                   enan (y, yy[0] != 0);
3039                   return;
3040                 }
3041             }
3042         }
3043       else
3044         {
3045           for (i = 1; i < 8; i++)
3046             {
3047               if (pe[i] != 0)
3048                 {
3049                   enan (y, yy[0] != 0);
3050                   return;
3051                 }
3052             }
3053         }
3054 #endif /* NANS */
3055       eclear (y);
3056       einfin (y);
3057       if (yy[0])
3058         eneg (y);
3059       return;
3060     }
3061 #endif  /* INFINITY */
3062   yy[E] = r;
3063   p = &yy[M + 1];
3064 #ifdef IEEE
3065   if (! REAL_WORDS_BIG_ENDIAN)
3066     {
3067       for (i = 0; i < 7; i++)
3068         *p++ = *(--e);
3069     }
3070   else
3071     {
3072       ++e;
3073       for (i = 0; i < 7; i++)
3074         *p++ = *e++;
3075     }
3076 #endif
3077 /* If denormal, remove the implied bit; else shift down 1. */
3078   if (r == 0)
3079     {
3080       yy[M] = 0;
3081     }
3082   else
3083     {
3084       yy[M] = 1;
3085       eshift (yy, -1);
3086     }
3087   emovo (yy, y);
3088 }
3089
3090 /* Convert single precision float PE to e type Y.  */
3091
3092 static void
3093 e24toe (pe, y)
3094      unsigned EMUSHORT *pe, *y;
3095 {
3096 #ifdef IBM
3097
3098   ibmtoe (pe, y, SFmode);
3099
3100 #else
3101   register unsigned EMUSHORT r;
3102   register unsigned EMUSHORT *e, *p;
3103   unsigned EMUSHORT yy[NI];
3104   int denorm, k;
3105
3106   e = pe;
3107   denorm = 0;                   /* flag if denormalized number */
3108   ecleaz (yy);
3109 #ifdef IEEE
3110   if (! REAL_WORDS_BIG_ENDIAN)
3111     e += 1;
3112 #endif
3113 #ifdef DEC
3114   e += 1;
3115 #endif
3116   r = *e;
3117   yy[0] = 0;
3118   if (r & 0x8000)
3119     yy[0] = 0xffff;
3120   yy[M] = (r & 0x7f) | 0200;
3121   r &= ~0x807f;                 /* strip sign and 7 significand bits */
3122 #ifdef INFINITY
3123   if (r == 0x7f80)
3124     {
3125 #ifdef NANS
3126       if (REAL_WORDS_BIG_ENDIAN)
3127         {
3128           if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
3129             {
3130               enan (y, yy[0] != 0);
3131               return;
3132             }
3133         }
3134       else
3135         {
3136           if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
3137             {
3138               enan (y, yy[0] != 0);
3139               return;
3140             }
3141         }
3142 #endif  /* NANS */
3143       eclear (y);
3144       einfin (y);
3145       if (yy[0])
3146         eneg (y);
3147       return;
3148     }
3149 #endif  /* INFINITY */
3150   r >>= 7;
3151   /* If zero exponent, then the significand is denormalized.
3152      So take back the understood high significand bit. */
3153   if (r == 0)
3154     {
3155       denorm = 1;
3156       yy[M] &= ~0200;
3157     }
3158   r += EXONE - 0177;
3159   yy[E] = r;
3160   p = &yy[M + 1];
3161 #ifdef DEC
3162   *p++ = *(--e);
3163 #endif
3164 #ifdef IEEE
3165   if (! REAL_WORDS_BIG_ENDIAN)
3166     *p++ = *(--e);
3167   else
3168     {
3169       ++e;
3170       *p++ = *e++;
3171     }
3172 #endif
3173   eshift (yy, -8);
3174   if (denorm)
3175     {                           /* if zero exponent, then normalize the significand */
3176       if ((k = enormlz (yy)) > NBITS)
3177         ecleazs (yy);
3178       else
3179         yy[E] -= (unsigned EMUSHORT) (k - 1);
3180     }
3181   emovo (yy, y);
3182 #endif /* not IBM */
3183 }
3184
3185 /* Convert e-type X to IEEE 128-bit long double format E.  */
3186
3187 static void
3188 etoe113 (x, e)
3189      unsigned EMUSHORT *x, *e;
3190 {
3191   unsigned EMUSHORT xi[NI];
3192   EMULONG exp;
3193   int rndsav;
3194
3195 #ifdef NANS
3196   if (eisnan (x))
3197     {
3198       make_nan (e, eisneg (x), TFmode);
3199       return;
3200     }
3201 #endif
3202   emovi (x, xi);
3203   exp = (EMULONG) xi[E];
3204 #ifdef INFINITY
3205   if (eisinf (x))
3206     goto nonorm;
3207 #endif
3208   /* round off to nearest or even */
3209   rndsav = rndprc;
3210   rndprc = 113;
3211   emdnorm (xi, 0, 0, exp, 64);
3212   rndprc = rndsav;
3213  nonorm:
3214   toe113 (xi, e);
3215 }
3216
3217 /* Convert exploded e-type X, that has already been rounded to
3218    113-bit precision, to IEEE 128-bit long double format Y.  */
3219
3220 static void
3221 toe113 (a, b)
3222      unsigned EMUSHORT *a, *b;
3223 {
3224   register unsigned EMUSHORT *p, *q;
3225   unsigned EMUSHORT i;
3226
3227 #ifdef NANS
3228   if (eiisnan (a))
3229     {
3230       make_nan (b, eiisneg (a), TFmode);
3231       return;
3232     }
3233 #endif
3234   p = a;
3235   if (REAL_WORDS_BIG_ENDIAN)
3236     q = b;
3237   else
3238     q = b + 7;                  /* point to output exponent */
3239
3240   /* If not denormal, delete the implied bit. */
3241   if (a[E] != 0)
3242     {
3243       eshup1 (a);
3244     }
3245   /* combine sign and exponent */
3246   i = *p++;
3247   if (REAL_WORDS_BIG_ENDIAN)
3248     {
3249       if (i)
3250         *q++ = *p++ | 0x8000;
3251       else
3252         *q++ = *p++;
3253     }
3254   else
3255     {
3256       if (i)
3257         *q-- = *p++ | 0x8000;
3258       else
3259         *q-- = *p++;
3260     }
3261   /* skip over guard word */
3262   ++p;
3263   /* move the significand */
3264   if (REAL_WORDS_BIG_ENDIAN)
3265     {
3266       for (i = 0; i < 7; i++)
3267         *q++ = *p++;
3268     }
3269   else
3270     {
3271       for (i = 0; i < 7; i++)
3272         *q-- = *p++;
3273     }
3274 }
3275
3276 /* Convert e-type X to IEEE double extended format E.  */
3277
3278 static void
3279 etoe64 (x, e)
3280      unsigned EMUSHORT *x, *e;
3281 {
3282   unsigned EMUSHORT xi[NI];
3283   EMULONG exp;
3284   int rndsav;
3285
3286 #ifdef NANS
3287   if (eisnan (x))
3288     {
3289       make_nan (e, eisneg (x), XFmode);
3290       return;
3291     }
3292 #endif
3293   emovi (x, xi);
3294   /* adjust exponent for offset */
3295   exp = (EMULONG) xi[E];
3296 #ifdef INFINITY
3297   if (eisinf (x))
3298     goto nonorm;
3299 #endif
3300   /* round off to nearest or even */
3301   rndsav = rndprc;
3302   rndprc = 64;
3303   emdnorm (xi, 0, 0, exp, 64);
3304   rndprc = rndsav;
3305  nonorm:
3306   toe64 (xi, e);
3307 }
3308
3309 /* Convert exploded e-type X, that has already been rounded to
3310    64-bit precision, to IEEE double extended format Y.  */
3311
3312 static void
3313 toe64 (a, b)
3314      unsigned EMUSHORT *a, *b;
3315 {
3316   register unsigned EMUSHORT *p, *q;
3317   unsigned EMUSHORT i;
3318
3319 #ifdef NANS
3320   if (eiisnan (a))
3321     {
3322       make_nan (b, eiisneg (a), XFmode);
3323       return;
3324     }
3325 #endif
3326   p = a;
3327 #ifdef IBM
3328   q = b;
3329 #endif
3330 #ifdef DEC
3331   q = b + 4;
3332 #endif
3333 #ifdef IEEE
3334   if (REAL_WORDS_BIG_ENDIAN)
3335     q = b;
3336   else
3337     {
3338       q = b + 4;                        /* point to output exponent */
3339 #if LONG_DOUBLE_TYPE_SIZE == 96
3340       /* Clear the last two bytes of 12-byte Intel format */
3341       *(q+1) = 0;
3342 #endif
3343     }
3344 #endif
3345
3346   /* combine sign and exponent */
3347   i = *p++;
3348 #ifdef IBM
3349   if (i)
3350     *q++ = *p++ | 0x8000;
3351   else
3352     *q++ = *p++;
3353   *q++ = 0;
3354 #endif
3355 #ifdef DEC
3356   if (i)
3357     *q-- = *p++ | 0x8000;
3358   else
3359     *q-- = *p++;
3360 #endif
3361 #ifdef IEEE
3362   if (REAL_WORDS_BIG_ENDIAN)
3363     {
3364       if (i)
3365         *q++ = *p++ | 0x8000;
3366       else
3367         *q++ = *p++;
3368       *q++ = 0;
3369     }
3370   else
3371     {
3372       if (i)
3373         *q-- = *p++ | 0x8000;
3374       else
3375         *q-- = *p++;
3376     }
3377 #endif
3378   /* skip over guard word */
3379   ++p;
3380   /* move the significand */
3381 #ifdef IBM
3382   for (i = 0; i < 4; i++)
3383     *q++ = *p++;
3384 #endif
3385 #ifdef DEC
3386   for (i = 0; i < 4; i++)
3387     *q-- = *p++;
3388 #endif
3389 #ifdef IEEE
3390   if (REAL_WORDS_BIG_ENDIAN)
3391     {
3392       for (i = 0; i < 4; i++)
3393         *q++ = *p++;
3394     }
3395   else
3396     {
3397 #ifdef INFINITY
3398       if (eiisinf (a))
3399         {
3400           /* Intel long double infinity significand.  */
3401           *q-- = 0x8000;
3402           *q-- = 0;
3403           *q-- = 0;
3404           *q = 0;
3405           return;
3406         }
3407 #endif
3408       for (i = 0; i < 4; i++)
3409         *q-- = *p++;
3410     }
3411 #endif
3412 }
3413
3414 /* e type to double precision.  */
3415
3416 #ifdef DEC
3417 /* Convert e-type X to DEC-format double E.  */
3418
3419 static void
3420 etoe53 (x, e)
3421      unsigned EMUSHORT *x, *e;
3422 {
3423   etodec (x, e);                /* see etodec.c */
3424 }
3425
3426 /* Convert exploded e-type X, that has already been rounded to
3427    56-bit double precision, to DEC double Y.  */
3428
3429 static void
3430 toe53 (x, y)
3431      unsigned EMUSHORT *x, *y;
3432 {
3433   todec (x, y);
3434 }
3435
3436 #else
3437 #ifdef IBM
3438 /* Convert e-type X to IBM 370-format double E.  */
3439
3440 static void
3441 etoe53 (x, e)
3442      unsigned EMUSHORT *x, *e;
3443 {
3444   etoibm (x, e, DFmode);
3445 }
3446
3447 /* Convert exploded e-type X, that has already been rounded to
3448    56-bit precision, to IBM 370 double Y.  */
3449
3450 static void
3451 toe53 (x, y)
3452      unsigned EMUSHORT *x, *y;
3453 {
3454   toibm (x, y, DFmode);
3455 }
3456
3457 #else  /* it's neither DEC nor IBM */
3458
3459 /* Convert e-type X to IEEE double E.  */
3460
3461 static void
3462 etoe53 (x, e)
3463      unsigned EMUSHORT *x, *e;
3464 {
3465   unsigned EMUSHORT xi[NI];
3466   EMULONG exp;
3467   int rndsav;
3468
3469 #ifdef NANS
3470   if (eisnan (x))
3471     {
3472       make_nan (e, eisneg (x), DFmode);
3473       return;
3474     }
3475 #endif
3476   emovi (x, xi);
3477   /* adjust exponent for offsets */
3478   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
3479 #ifdef INFINITY
3480   if (eisinf (x))
3481     goto nonorm;
3482 #endif
3483   /* round off to nearest or even */
3484   rndsav = rndprc;
3485   rndprc = 53;
3486   emdnorm (xi, 0, 0, exp, 64);
3487   rndprc = rndsav;
3488  nonorm:
3489   toe53 (xi, e);
3490 }
3491
3492 /* Convert exploded e-type X, that has already been rounded to
3493    53-bit precision, to IEEE double Y.  */
3494
3495 static void
3496 toe53 (x, y)
3497      unsigned EMUSHORT *x, *y;
3498 {
3499   unsigned EMUSHORT i;
3500   unsigned EMUSHORT *p;
3501
3502 #ifdef NANS
3503   if (eiisnan (x))
3504     {
3505       make_nan (y, eiisneg (x), DFmode);
3506       return;
3507     }
3508 #endif
3509   p = &x[0];
3510 #ifdef IEEE
3511   if (! REAL_WORDS_BIG_ENDIAN)
3512     y += 3;
3513 #endif
3514   *y = 0;                       /* output high order */
3515   if (*p++)
3516     *y = 0x8000;                /* output sign bit */
3517
3518   i = *p++;
3519   if (i >= (unsigned int) 2047)
3520     {                           /* Saturate at largest number less than infinity. */
3521 #ifdef INFINITY
3522       *y |= 0x7ff0;
3523       if (! REAL_WORDS_BIG_ENDIAN)
3524         {
3525           *(--y) = 0;
3526           *(--y) = 0;
3527           *(--y) = 0;
3528         }
3529       else
3530         {
3531           ++y;
3532           *y++ = 0;
3533           *y++ = 0;
3534           *y++ = 0;
3535         }
3536 #else
3537       *y |= (unsigned EMUSHORT) 0x7fef;
3538       if (! REAL_WORDS_BIG_ENDIAN)
3539         {
3540           *(--y) = 0xffff;
3541           *(--y) = 0xffff;
3542           *(--y) = 0xffff;
3543         }
3544       else
3545         {
3546           ++y;
3547           *y++ = 0xffff;
3548           *y++ = 0xffff;
3549           *y++ = 0xffff;
3550         }
3551 #endif
3552       return;
3553     }
3554   if (i == 0)
3555     {
3556       eshift (x, 4);
3557     }
3558   else
3559     {
3560       i <<= 4;
3561       eshift (x, 5);
3562     }
3563   i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
3564   *y |= (unsigned EMUSHORT) i;  /* high order output already has sign bit set */
3565   if (! REAL_WORDS_BIG_ENDIAN)
3566     {
3567       *(--y) = *p++;
3568       *(--y) = *p++;
3569       *(--y) = *p;
3570     }
3571   else
3572     {
3573       ++y;
3574       *y++ = *p++;
3575       *y++ = *p++;
3576       *y++ = *p++;
3577     }
3578 }
3579
3580 #endif /* not IBM */
3581 #endif /* not DEC */
3582
3583
3584
3585 /* e type to single precision.  */
3586
3587 #ifdef IBM
3588 /* Convert e-type X to IBM 370 float E.  */
3589
3590 static void
3591 etoe24 (x, e)
3592      unsigned EMUSHORT *x, *e;
3593 {
3594   etoibm (x, e, SFmode);
3595 }
3596
3597 /* Convert exploded e-type X, that has already been rounded to
3598    float precision, to IBM 370 float Y.  */
3599
3600 static void
3601 toe24 (x, y)
3602      unsigned EMUSHORT *x, *y;
3603 {
3604   toibm (x, y, SFmode);
3605 }
3606
3607 #else
3608 /* Convert e-type X to IEEE float E.  DEC float is the same as IEEE float.  */
3609
3610 static void
3611 etoe24 (x, e)
3612      unsigned EMUSHORT *x, *e;
3613 {
3614   EMULONG exp;
3615   unsigned EMUSHORT xi[NI];
3616   int rndsav;
3617
3618 #ifdef NANS
3619   if (eisnan (x))
3620     {
3621       make_nan (e, eisneg (x), SFmode);
3622       return;
3623     }
3624 #endif
3625   emovi (x, xi);
3626   /* adjust exponent for offsets */
3627   exp = (EMULONG) xi[E] - (EXONE - 0177);
3628 #ifdef INFINITY
3629   if (eisinf (x))
3630     goto nonorm;
3631 #endif
3632   /* round off to nearest or even */
3633   rndsav = rndprc;
3634   rndprc = 24;
3635   emdnorm (xi, 0, 0, exp, 64);
3636   rndprc = rndsav;
3637  nonorm:
3638   toe24 (xi, e);
3639 }
3640
3641 /* Convert exploded e-type X, that has already been rounded to
3642    float precision, to IEEE float Y.  */
3643
3644 static void
3645 toe24 (x, y)
3646      unsigned EMUSHORT *x, *y;
3647 {
3648   unsigned EMUSHORT i;
3649   unsigned EMUSHORT *p;
3650
3651 #ifdef NANS
3652   if (eiisnan (x))
3653     {
3654       make_nan (y, eiisneg (x), SFmode);
3655       return;
3656     }
3657 #endif
3658   p = &x[0];
3659 #ifdef IEEE
3660   if (! REAL_WORDS_BIG_ENDIAN)
3661     y += 1;
3662 #endif
3663 #ifdef DEC
3664   y += 1;
3665 #endif
3666   *y = 0;                       /* output high order */
3667   if (*p++)
3668     *y = 0x8000;                /* output sign bit */
3669
3670   i = *p++;
3671 /* Handle overflow cases. */
3672   if (i >= 255)
3673     {
3674 #ifdef INFINITY
3675       *y |= (unsigned EMUSHORT) 0x7f80;
3676 #ifdef DEC
3677       *(--y) = 0;
3678 #endif
3679 #ifdef IEEE
3680       if (! REAL_WORDS_BIG_ENDIAN)
3681         *(--y) = 0;
3682       else
3683         {
3684           ++y;
3685           *y = 0;
3686         }
3687 #endif
3688 #else  /* no INFINITY */
3689       *y |= (unsigned EMUSHORT) 0x7f7f;
3690 #ifdef DEC
3691       *(--y) = 0xffff;
3692 #endif
3693 #ifdef IEEE
3694       if (! REAL_WORDS_BIG_ENDIAN)
3695         *(--y) = 0xffff;
3696       else
3697         {
3698           ++y;
3699           *y = 0xffff;
3700         }
3701 #endif
3702 #ifdef ERANGE
3703       errno = ERANGE;
3704 #endif
3705 #endif  /* no INFINITY */
3706       return;
3707     }
3708   if (i == 0)
3709     {
3710       eshift (x, 7);
3711     }
3712   else
3713     {
3714       i <<= 7;
3715       eshift (x, 8);
3716     }
3717   i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
3718   /* High order output already has sign bit set.  */
3719   *y |= i;
3720 #ifdef DEC
3721   *(--y) = *p;
3722 #endif
3723 #ifdef IEEE
3724   if (! REAL_WORDS_BIG_ENDIAN)
3725     *(--y) = *p;
3726   else
3727     {
3728       ++y;
3729       *y = *p;
3730     }
3731 #endif
3732 }
3733 #endif  /* not IBM */
3734
3735 /* Compare two e type numbers.
3736    Return +1 if a > b
3737            0 if a == b
3738           -1 if a < b
3739           -2 if either a or b is a NaN.  */
3740
3741 static int
3742 ecmp (a, b)
3743      unsigned EMUSHORT *a, *b;
3744 {
3745   unsigned EMUSHORT ai[NI], bi[NI];
3746   register unsigned EMUSHORT *p, *q;
3747   register int i;
3748   int msign;
3749
3750 #ifdef NANS
3751   if (eisnan (a)  || eisnan (b))
3752       return (-2);
3753 #endif
3754   emovi (a, ai);
3755   p = ai;
3756   emovi (b, bi);
3757   q = bi;
3758
3759   if (*p != *q)
3760     {                           /* the signs are different */
3761       /* -0 equals + 0 */
3762       for (i = 1; i < NI - 1; i++)
3763         {
3764           if (ai[i] != 0)
3765             goto nzro;
3766           if (bi[i] != 0)
3767             goto nzro;
3768         }
3769       return (0);
3770     nzro:
3771       if (*p == 0)
3772         return (1);
3773       else
3774         return (-1);
3775     }
3776   /* both are the same sign */
3777   if (*p == 0)
3778     msign = 1;
3779   else
3780     msign = -1;
3781   i = NI - 1;
3782   do
3783     {
3784       if (*p++ != *q++)
3785         {
3786           goto diff;
3787         }
3788     }
3789   while (--i > 0);
3790
3791   return (0);                   /* equality */
3792
3793  diff:
3794
3795   if (*(--p) > *(--q))
3796     return (msign);             /* p is bigger */
3797   else
3798     return (-msign);            /* p is littler */
3799 }
3800
3801 /* Find e-type nearest integer to X, as floor (X + 0.5).  */
3802
3803 static void
3804 eround (x, y)
3805      unsigned EMUSHORT *x, *y;
3806 {
3807   eadd (ehalf, x, y);
3808   efloor (y, y);
3809 }
3810
3811 /* Convert HOST_WIDE_INT LP to e type Y.  */
3812
3813 static void
3814 ltoe (lp, y)
3815      HOST_WIDE_INT *lp;
3816      unsigned EMUSHORT *y;
3817 {
3818   unsigned EMUSHORT yi[NI];
3819   unsigned HOST_WIDE_INT ll;
3820   int k;
3821
3822   ecleaz (yi);
3823   if (*lp < 0)
3824     {
3825       /* make it positive */
3826       ll = (unsigned HOST_WIDE_INT) (-(*lp));
3827       yi[0] = 0xffff;           /* put correct sign in the e type number */
3828     }
3829   else
3830     {
3831       ll = (unsigned HOST_WIDE_INT) (*lp);
3832     }
3833   /* move the long integer to yi significand area */
3834 #if HOST_BITS_PER_WIDE_INT == 64
3835   yi[M] = (unsigned EMUSHORT) (ll >> 48);
3836   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3837   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3838   yi[M + 3] = (unsigned EMUSHORT) ll;
3839   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
3840 #else
3841   yi[M] = (unsigned EMUSHORT) (ll >> 16);
3842   yi[M + 1] = (unsigned EMUSHORT) ll;
3843   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
3844 #endif
3845
3846   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3847     ecleaz (yi);                /* it was zero */
3848   else
3849     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
3850   emovo (yi, y);                /* output the answer */
3851 }
3852
3853 /* Convert unsigned HOST_WIDE_INT LP to e type Y.  */
3854
3855 static void
3856 ultoe (lp, y)
3857      unsigned HOST_WIDE_INT *lp;
3858      unsigned EMUSHORT *y;
3859 {
3860   unsigned EMUSHORT yi[NI];
3861   unsigned HOST_WIDE_INT ll;
3862   int k;
3863
3864   ecleaz (yi);
3865   ll = *lp;
3866
3867   /* move the long integer to ayi significand area */
3868 #if HOST_BITS_PER_WIDE_INT == 64
3869   yi[M] = (unsigned EMUSHORT) (ll >> 48);
3870   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3871   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3872   yi[M + 3] = (unsigned EMUSHORT) ll;
3873   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
3874 #else
3875   yi[M] = (unsigned EMUSHORT) (ll >> 16);
3876   yi[M + 1] = (unsigned EMUSHORT) ll;
3877   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
3878 #endif
3879
3880   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3881     ecleaz (yi);                /* it was zero */
3882   else
3883     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
3884   emovo (yi, y);                /* output the answer */
3885 }
3886
3887
3888 /* Find signed HOST_WIDE_INT integer I and floating point fractional
3889    part FRAC of e-type (packed internal format) floating point input X.
3890    The integer output I has the sign of the input, except that
3891    positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
3892    The output e-type fraction FRAC is the positive fractional
3893    part of abs (X).  */
3894
3895 static void
3896 eifrac (x, i, frac)
3897      unsigned EMUSHORT *x;
3898      HOST_WIDE_INT *i;
3899      unsigned EMUSHORT *frac;
3900 {
3901   unsigned EMUSHORT xi[NI];
3902   int j, k;
3903   unsigned HOST_WIDE_INT ll;
3904
3905   emovi (x, xi);
3906   k = (int) xi[E] - (EXONE - 1);
3907   if (k <= 0)
3908     {
3909       /* if exponent <= 0, integer = 0 and real output is fraction */
3910       *i = 0L;
3911       emovo (xi, frac);
3912       return;
3913     }
3914   if (k > (HOST_BITS_PER_WIDE_INT - 1))
3915     {
3916       /* long integer overflow: output large integer
3917          and correct fraction  */
3918       if (xi[0])
3919         *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
3920       else
3921         {
3922 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
3923           /* In this case, let it overflow and convert as if unsigned.  */
3924           euifrac (x, &ll, frac);
3925           *i = (HOST_WIDE_INT) ll;
3926           return;
3927 #else
3928           /* In other cases, return the largest positive integer.  */
3929           *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
3930 #endif
3931         }
3932       eshift (xi, k);
3933       if (extra_warnings)
3934         warning ("overflow on truncation to integer");
3935     }
3936   else if (k > 16)
3937     {
3938       /* Shift more than 16 bits: first shift up k-16 mod 16,
3939          then shift up by 16's.  */
3940       j = k - ((k >> 4) << 4);
3941       eshift (xi, j);
3942       ll = xi[M];
3943       k -= j;
3944       do
3945         {
3946           eshup6 (xi);
3947           ll = (ll << 16) | xi[M];
3948         }
3949       while ((k -= 16) > 0);
3950       *i = ll;
3951       if (xi[0])
3952         *i = -(*i);
3953     }
3954   else
3955       {
3956         /* shift not more than 16 bits */
3957           eshift (xi, k);
3958         *i = (HOST_WIDE_INT) xi[M] & 0xffff;
3959         if (xi[0])
3960           *i = -(*i);
3961       }
3962   xi[0] = 0;
3963   xi[E] = EXONE - 1;
3964   xi[M] = 0;
3965   if ((k = enormlz (xi)) > NBITS)
3966     ecleaz (xi);
3967   else
3968     xi[E] -= (unsigned EMUSHORT) k;
3969
3970   emovo (xi, frac);
3971 }
3972
3973
3974 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
3975    FRAC of e-type X.  A negative input yields integer output = 0 but
3976    correct fraction.  */
3977
3978 static void
3979 euifrac (x, i, frac)
3980      unsigned EMUSHORT *x;
3981      unsigned HOST_WIDE_INT *i;
3982      unsigned EMUSHORT *frac;
3983 {
3984   unsigned HOST_WIDE_INT ll;
3985   unsigned EMUSHORT xi[NI];
3986   int j, k;
3987
3988   emovi (x, xi);
3989   k = (int) xi[E] - (EXONE - 1);
3990   if (k <= 0)
3991     {
3992       /* if exponent <= 0, integer = 0 and argument is fraction */
3993       *i = 0L;
3994       emovo (xi, frac);
3995       return;
3996     }
3997   if (k > HOST_BITS_PER_WIDE_INT)
3998     {
3999       /* Long integer overflow: output large integer
4000          and correct fraction.
4001          Note, the BSD microvax compiler says that ~(0UL)
4002          is a syntax error.  */
4003       *i = ~(0L);
4004       eshift (xi, k);
4005       if (extra_warnings)
4006         warning ("overflow on truncation to unsigned integer");
4007     }
4008   else if (k > 16)
4009     {
4010       /* Shift more than 16 bits: first shift up k-16 mod 16,
4011          then shift up by 16's.  */
4012       j = k - ((k >> 4) << 4);
4013       eshift (xi, j);
4014       ll = xi[M];
4015       k -= j;
4016       do
4017         {
4018           eshup6 (xi);
4019           ll = (ll << 16) | xi[M];
4020         }
4021       while ((k -= 16) > 0);
4022       *i = ll;
4023     }
4024   else
4025     {
4026       /* shift not more than 16 bits */
4027       eshift (xi, k);
4028       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4029     }
4030
4031   if (xi[0])  /* A negative value yields unsigned integer 0. */
4032     *i = 0L;
4033
4034   xi[0] = 0;
4035   xi[E] = EXONE - 1;
4036   xi[M] = 0;
4037   if ((k = enormlz (xi)) > NBITS)
4038     ecleaz (xi);
4039   else
4040     xi[E] -= (unsigned EMUSHORT) k;
4041
4042   emovo (xi, frac);
4043 }
4044
4045 /* Shift the significand of exploded e-type X up or down by SC bits.  */
4046
4047 static int
4048 eshift (x, sc)
4049      unsigned EMUSHORT *x;
4050      int sc;
4051 {
4052   unsigned EMUSHORT lost;
4053   unsigned EMUSHORT *p;
4054
4055   if (sc == 0)
4056     return (0);
4057
4058   lost = 0;
4059   p = x + NI - 1;
4060
4061   if (sc < 0)
4062     {
4063       sc = -sc;
4064       while (sc >= 16)
4065         {
4066           lost |= *p;           /* remember lost bits */
4067           eshdn6 (x);
4068           sc -= 16;
4069         }
4070
4071       while (sc >= 8)
4072         {
4073           lost |= *p & 0xff;
4074           eshdn8 (x);
4075           sc -= 8;
4076         }
4077
4078       while (sc > 0)
4079         {
4080           lost |= *p & 1;
4081           eshdn1 (x);
4082           sc -= 1;
4083         }
4084     }
4085   else
4086     {
4087       while (sc >= 16)
4088         {
4089           eshup6 (x);
4090           sc -= 16;
4091         }
4092
4093       while (sc >= 8)
4094         {
4095           eshup8 (x);
4096           sc -= 8;
4097         }
4098
4099       while (sc > 0)
4100         {
4101           eshup1 (x);
4102           sc -= 1;
4103         }
4104     }
4105   if (lost)
4106     lost = 1;
4107   return ((int) lost);
4108 }
4109
4110 /* Shift normalize the significand area of exploded e-type X.
4111    Return the shift count (up = positive).  */
4112
4113 static int
4114 enormlz (x)
4115      unsigned EMUSHORT x[];
4116 {
4117   register unsigned EMUSHORT *p;
4118   int sc;
4119
4120   sc = 0;
4121   p = &x[M];
4122   if (*p != 0)
4123     goto normdn;
4124   ++p;
4125   if (*p & 0x8000)
4126     return (0);                 /* already normalized */
4127   while (*p == 0)
4128     {
4129       eshup6 (x);
4130       sc += 16;
4131
4132       /* With guard word, there are NBITS+16 bits available.
4133        Return true if all are zero.  */
4134       if (sc > NBITS)
4135         return (sc);
4136     }
4137   /* see if high byte is zero */
4138   while ((*p & 0xff00) == 0)
4139     {
4140       eshup8 (x);
4141       sc += 8;
4142     }
4143   /* now shift 1 bit at a time */
4144   while ((*p & 0x8000) == 0)
4145     {
4146       eshup1 (x);
4147       sc += 1;
4148       if (sc > NBITS)
4149         {
4150           mtherr ("enormlz", UNDERFLOW);
4151           return (sc);
4152         }
4153     }
4154   return (sc);
4155
4156   /* Normalize by shifting down out of the high guard word
4157      of the significand */
4158  normdn:
4159
4160   if (*p & 0xff00)
4161     {
4162       eshdn8 (x);
4163       sc -= 8;
4164     }
4165   while (*p != 0)
4166     {
4167       eshdn1 (x);
4168       sc -= 1;
4169
4170       if (sc < -NBITS)
4171         {
4172           mtherr ("enormlz", OVERFLOW);
4173           return (sc);
4174         }
4175     }
4176   return (sc);
4177 }
4178
4179 /* Powers of ten used in decimal <-> binary conversions.  */
4180
4181 #define NTEN 12
4182 #define MAXP 4096
4183
4184 #if LONG_DOUBLE_TYPE_SIZE == 128
4185 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4186 {
4187   {0x6576, 0x4a92, 0x804a, 0x153f,
4188    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4189   {0x6a32, 0xce52, 0x329a, 0x28ce,
4190    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4191   {0x526c, 0x50ce, 0xf18b, 0x3d28,
4192    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4193   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4194    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4195   {0x851e, 0xeab7, 0x98fe, 0x901b,
4196    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4197   {0x0235, 0x0137, 0x36b1, 0x336c,
4198    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4199   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4200    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4201   {0x0000, 0x0000, 0x0000, 0x0000,
4202    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4203   {0x0000, 0x0000, 0x0000, 0x0000,
4204    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4205   {0x0000, 0x0000, 0x0000, 0x0000,
4206    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4207   {0x0000, 0x0000, 0x0000, 0x0000,
4208    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4209   {0x0000, 0x0000, 0x0000, 0x0000,
4210    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4211   {0x0000, 0x0000, 0x0000, 0x0000,
4212    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4213 };
4214
4215 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4216 {
4217   {0x2030, 0xcffc, 0xa1c3, 0x8123,
4218    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4219   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4220    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4221   {0xf53f, 0xf698, 0x6bd3, 0x0158,
4222    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4223   {0xe731, 0x04d4, 0xe3f2, 0xd332,
4224    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4225   {0xa23e, 0x5308, 0xfefb, 0x1155,
4226    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4227   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4228    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4229   {0x2a20, 0x6224, 0x47b3, 0x98d7,
4230    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4231   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4232    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4233   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4234    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4235   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4236    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4237   {0xc155, 0xa4a8, 0x404e, 0x6113,
4238    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4239   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4240    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4241   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4242    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4243 };
4244 #else
4245 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4246 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4247 {
4248   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4249   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4250   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4251   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4252   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4253   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4254   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4255   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4256   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4257   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4258   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4259   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4260   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4261 };
4262
4263 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4264 {
4265   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4266   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4267   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4268   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4269   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4270   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4271   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4272   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4273   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4274   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4275   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4276   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4277   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4278 };
4279 #endif
4280
4281 /* Convert float value X to ASCII string STRING with NDIG digits after
4282    the decimal point.  */
4283
4284 static void
4285 e24toasc (x, string, ndigs)
4286      unsigned EMUSHORT x[];
4287      char *string;
4288      int ndigs;
4289 {
4290   unsigned EMUSHORT w[NI];
4291
4292   e24toe (x, w);
4293   etoasc (w, string, ndigs);
4294 }
4295
4296 /* Convert double value X to ASCII string STRING with NDIG digits after
4297    the decimal point.  */
4298
4299 static void
4300 e53toasc (x, string, ndigs)
4301      unsigned EMUSHORT x[];
4302      char *string;
4303      int ndigs;
4304 {
4305   unsigned EMUSHORT w[NI];
4306
4307   e53toe (x, w);
4308   etoasc (w, string, ndigs);
4309 }
4310
4311 /* Convert double extended value X to ASCII string STRING with NDIG digits
4312    after the decimal point.  */
4313
4314 static void
4315 e64toasc (x, string, ndigs)
4316      unsigned EMUSHORT x[];
4317      char *string;
4318      int ndigs;
4319 {
4320   unsigned EMUSHORT w[NI];
4321
4322   e64toe (x, w);
4323   etoasc (w, string, ndigs);
4324 }
4325
4326 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4327    after the decimal point.  */
4328
4329 static void
4330 e113toasc (x, string, ndigs)
4331      unsigned EMUSHORT x[];
4332      char *string;
4333      int ndigs;
4334 {
4335   unsigned EMUSHORT w[NI];
4336
4337   e113toe (x, w);
4338   etoasc (w, string, ndigs);
4339 }
4340
4341 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4342    the decimal point.  */
4343
4344 static char wstring[80];        /* working storage for ASCII output */
4345
4346 static void
4347 etoasc (x, string, ndigs)
4348      unsigned EMUSHORT x[];
4349      char *string;
4350      int ndigs;
4351 {
4352   EMUSHORT digit;
4353   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
4354   unsigned EMUSHORT *p, *r, *ten;
4355   unsigned EMUSHORT sign;
4356   int i, j, k, expon, rndsav;
4357   char *s, *ss;
4358   unsigned EMUSHORT m;
4359
4360
4361   rndsav = rndprc;
4362   ss = string;
4363   s = wstring;
4364   *ss = '\0';
4365   *s = '\0';
4366 #ifdef NANS
4367   if (eisnan (x))
4368     {
4369       sprintf (wstring, " NaN ");
4370       goto bxit;
4371     }
4372 #endif
4373   rndprc = NBITS;               /* set to full precision */
4374   emov (x, y);                  /* retain external format */
4375   if (y[NE - 1] & 0x8000)
4376     {
4377       sign = 0xffff;
4378       y[NE - 1] &= 0x7fff;
4379     }
4380   else
4381     {
4382       sign = 0;
4383     }
4384   expon = 0;
4385   ten = &etens[NTEN][0];
4386   emov (eone, t);
4387   /* Test for zero exponent */
4388   if (y[NE - 1] == 0)
4389     {
4390       for (k = 0; k < NE - 1; k++)
4391         {
4392           if (y[k] != 0)
4393             goto tnzro;         /* denormalized number */
4394         }
4395       goto isone;               /* valid all zeros */
4396     }
4397  tnzro:
4398
4399   /* Test for infinity. */
4400   if (y[NE - 1] == 0x7fff)
4401     {
4402       if (sign)
4403         sprintf (wstring, " -Infinity ");
4404       else
4405         sprintf (wstring, " Infinity ");
4406       goto bxit;
4407     }
4408
4409   /* Test for exponent nonzero but significand denormalized.
4410    * This is an error condition.
4411    */
4412   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4413     {
4414       mtherr ("etoasc", DOMAIN);
4415       sprintf (wstring, "NaN");
4416       goto bxit;
4417     }
4418
4419   /* Compare to 1.0 */
4420   i = ecmp (eone, y);
4421   if (i == 0)
4422     goto isone;
4423
4424   if (i == -2)
4425     abort ();
4426
4427   if (i < 0)
4428     {                           /* Number is greater than 1 */
4429       /* Convert significand to an integer and strip trailing decimal zeros. */
4430       emov (y, u);
4431       u[NE - 1] = EXONE + NBITS - 1;
4432
4433       p = &etens[NTEN - 4][0];
4434       m = 16;
4435       do
4436         {
4437           ediv (p, u, t);
4438           efloor (t, w);
4439           for (j = 0; j < NE - 1; j++)
4440             {
4441               if (t[j] != w[j])
4442                 goto noint;
4443             }
4444           emov (t, u);
4445           expon += (int) m;
4446         noint:
4447           p += NE;
4448           m >>= 1;
4449         }
4450       while (m != 0);
4451
4452       /* Rescale from integer significand */
4453       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4454       emov (u, y);
4455       /* Find power of 10 */
4456       emov (eone, t);
4457       m = MAXP;
4458       p = &etens[0][0];
4459       /* An unordered compare result shouldn't happen here. */
4460       while (ecmp (ten, u) <= 0)
4461         {
4462           if (ecmp (p, u) <= 0)
4463             {
4464               ediv (p, u, u);
4465               emul (p, t, t);
4466               expon += (int) m;
4467             }
4468           m >>= 1;
4469           if (m == 0)
4470             break;
4471           p += NE;
4472         }
4473     }
4474   else
4475     {                           /* Number is less than 1.0 */
4476       /* Pad significand with trailing decimal zeros. */
4477       if (y[NE - 1] == 0)
4478         {
4479           while ((y[NE - 2] & 0x8000) == 0)
4480             {
4481               emul (ten, y, y);
4482               expon -= 1;
4483             }
4484         }
4485       else
4486         {
4487           emovi (y, w);
4488           for (i = 0; i < NDEC + 1; i++)
4489             {
4490               if ((w[NI - 1] & 0x7) != 0)
4491                 break;
4492               /* multiply by 10 */
4493               emovz (w, u);
4494               eshdn1 (u);
4495               eshdn1 (u);
4496               eaddm (w, u);
4497               u[1] += 3;
4498               while (u[2] != 0)
4499                 {
4500                   eshdn1 (u);
4501                   u[1] += 1;
4502                 }
4503               if (u[NI - 1] != 0)
4504                 break;
4505               if (eone[NE - 1] <= u[1])
4506                 break;
4507               emovz (u, w);
4508               expon -= 1;
4509             }
4510           emovo (w, y);
4511         }
4512       k = -MAXP;
4513       p = &emtens[0][0];
4514       r = &etens[0][0];
4515       emov (y, w);
4516       emov (eone, t);
4517       while (ecmp (eone, w) > 0)
4518         {
4519           if (ecmp (p, w) >= 0)
4520             {
4521               emul (r, w, w);
4522               emul (r, t, t);
4523               expon += k;
4524             }
4525           k /= 2;
4526           if (k == 0)
4527             break;
4528           p += NE;
4529           r += NE;
4530         }
4531       ediv (t, eone, t);
4532     }
4533  isone:
4534   /* Find the first (leading) digit. */
4535   emovi (t, w);
4536   emovz (w, t);
4537   emovi (y, w);
4538   emovz (w, y);
4539   eiremain (t, y);
4540   digit = equot[NI - 1];
4541   while ((digit == 0) && (ecmp (y, ezero) != 0))
4542     {
4543       eshup1 (y);
4544       emovz (y, u);
4545       eshup1 (u);
4546       eshup1 (u);
4547       eaddm (u, y);
4548       eiremain (t, y);
4549       digit = equot[NI - 1];
4550       expon -= 1;
4551     }
4552   s = wstring;
4553   if (sign)
4554     *s++ = '-';
4555   else
4556     *s++ = ' ';
4557   /* Examine number of digits requested by caller. */
4558   if (ndigs < 0)
4559     ndigs = 0;
4560   if (ndigs > NDEC)
4561     ndigs = NDEC;
4562   if (digit == 10)
4563     {
4564       *s++ = '1';
4565       *s++ = '.';
4566       if (ndigs > 0)
4567         {
4568           *s++ = '0';
4569           ndigs -= 1;
4570         }
4571       expon += 1;
4572     }
4573   else
4574     {
4575       *s++ = (char)digit + '0';
4576       *s++ = '.';
4577     }
4578   /* Generate digits after the decimal point. */
4579   for (k = 0; k <= ndigs; k++)
4580     {
4581       /* multiply current number by 10, without normalizing */
4582       eshup1 (y);
4583       emovz (y, u);
4584       eshup1 (u);
4585       eshup1 (u);
4586       eaddm (u, y);
4587       eiremain (t, y);
4588       *s++ = (char) equot[NI - 1] + '0';
4589     }
4590   digit = equot[NI - 1];
4591   --s;
4592   ss = s;
4593   /* round off the ASCII string */
4594   if (digit > 4)
4595     {
4596       /* Test for critical rounding case in ASCII output. */
4597       if (digit == 5)
4598         {
4599           emovo (y, t);
4600           if (ecmp (t, ezero) != 0)
4601             goto roun;          /* round to nearest */
4602           if ((*(s - 1) & 1) == 0)
4603             goto doexp;         /* round to even */
4604         }
4605       /* Round up and propagate carry-outs */
4606     roun:
4607       --s;
4608       k = *s & 0x7f;
4609       /* Carry out to most significant digit? */
4610       if (k == '.')
4611         {
4612           --s;
4613           k = *s;
4614           k += 1;
4615           *s = (char) k;
4616           /* Most significant digit carries to 10? */
4617           if (k > '9')
4618             {
4619               expon += 1;
4620               *s = '1';
4621             }
4622           goto doexp;
4623         }
4624       /* Round up and carry out from less significant digits */
4625       k += 1;
4626       *s = (char) k;
4627       if (k > '9')
4628         {
4629           *s = '0';
4630           goto roun;
4631         }
4632     }
4633  doexp:
4634   /*
4635      if (expon >= 0)
4636      sprintf (ss, "e+%d", expon);
4637      else
4638      sprintf (ss, "e%d", expon);
4639      */
4640   sprintf (ss, "e%d", expon);
4641  bxit:
4642   rndprc = rndsav;
4643   /* copy out the working string */
4644   s = string;
4645   ss = wstring;
4646   while (*ss == ' ')            /* strip possible leading space */
4647     ++ss;
4648   while ((*s++ = *ss++) != '\0')
4649     ;
4650 }
4651
4652
4653 /* Convert ASCII string to floating point.
4654
4655    Numeric input is a free format decimal number of any length, with
4656    or without decimal point.  Entering E after the number followed by an
4657    integer number causes the second number to be interpreted as a power of
4658    10 to be multiplied by the first number (i.e., "scientific" notation).  */
4659
4660 /* Convert ASCII string S to single precision float value Y.  */
4661
4662 static void
4663 asctoe24 (s, y)
4664      char *s;
4665      unsigned EMUSHORT *y;
4666 {
4667   asctoeg (s, y, 24);
4668 }
4669
4670
4671 /* Convert ASCII string S to double precision value Y.  */
4672
4673 static void
4674 asctoe53 (s, y)
4675      char *s;
4676      unsigned EMUSHORT *y;
4677 {
4678 #if defined(DEC) || defined(IBM)
4679   asctoeg (s, y, 56);
4680 #else
4681   asctoeg (s, y, 53);
4682 #endif
4683 }
4684
4685
4686 /* Convert ASCII string S to double extended value Y.  */
4687
4688 static void
4689 asctoe64 (s, y)
4690      char *s;
4691      unsigned EMUSHORT *y;
4692 {
4693   asctoeg (s, y, 64);
4694 }
4695
4696 /* Convert ASCII string S to 128-bit long double Y.  */
4697
4698 static void
4699 asctoe113 (s, y)
4700      char *s;
4701      unsigned EMUSHORT *y;
4702 {
4703   asctoeg (s, y, 113);
4704 }
4705
4706 /* Convert ASCII string S to e type Y.  */
4707
4708 static void
4709 asctoe (s, y)
4710      char *s;
4711      unsigned EMUSHORT *y;
4712 {
4713   asctoeg (s, y, NBITS);
4714 }
4715
4716 /* Convert ASCII string SS to e type Y, with a specified rounding precision
4717    of OPREC bits. */
4718
4719 static void
4720 asctoeg (ss, y, oprec)
4721      char *ss;
4722      unsigned EMUSHORT *y;
4723      int oprec;
4724 {
4725   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
4726   int esign, decflg, sgnflg, nexp, exp, prec, lost;
4727   int k, trail, c, rndsav;
4728   EMULONG lexp;
4729   unsigned EMUSHORT nsign, *p;
4730   char *sp, *s, *lstr;
4731
4732   /* Copy the input string. */
4733   lstr = (char *) alloca (strlen (ss) + 1);
4734   s = ss;
4735   while (*s == ' ')             /* skip leading spaces */
4736     ++s;
4737   sp = lstr;
4738   while ((*sp++ = *s++) != '\0')
4739     ;
4740   s = lstr;
4741
4742   rndsav = rndprc;
4743   rndprc = NBITS;               /* Set to full precision */
4744   lost = 0;
4745   nsign = 0;
4746   decflg = 0;
4747   sgnflg = 0;
4748   nexp = 0;
4749   exp = 0;
4750   prec = 0;
4751   ecleaz (yy);
4752   trail = 0;
4753
4754  nxtcom:
4755   k = *s - '0';
4756   if ((k >= 0) && (k <= 9))
4757     {
4758       /* Ignore leading zeros */
4759       if ((prec == 0) && (decflg == 0) && (k == 0))
4760         goto donchr;
4761       /* Identify and strip trailing zeros after the decimal point. */
4762       if ((trail == 0) && (decflg != 0))
4763         {
4764           sp = s;
4765           while ((*sp >= '0') && (*sp <= '9'))
4766             ++sp;
4767           /* Check for syntax error */
4768           c = *sp & 0x7f;
4769           if ((c != 'e') && (c != 'E') && (c != '\0')
4770               && (c != '\n') && (c != '\r') && (c != ' ')
4771               && (c != ','))
4772             goto error;
4773           --sp;
4774           while (*sp == '0')
4775             *sp-- = 'z';
4776           trail = 1;
4777           if (*s == 'z')
4778             goto donchr;
4779         }
4780
4781       /* If enough digits were given to more than fill up the yy register,
4782          continuing until overflow into the high guard word yy[2]
4783          guarantees that there will be a roundoff bit at the top
4784          of the low guard word after normalization.  */
4785
4786       if (yy[2] == 0)
4787         {
4788           if (decflg)
4789             nexp += 1;          /* count digits after decimal point */
4790           eshup1 (yy);          /* multiply current number by 10 */
4791           emovz (yy, xt);
4792           eshup1 (xt);
4793           eshup1 (xt);
4794           eaddm (xt, yy);
4795           ecleaz (xt);
4796           xt[NI - 2] = (unsigned EMUSHORT) k;
4797           eaddm (xt, yy);
4798         }
4799       else
4800         {
4801           /* Mark any lost non-zero digit.  */
4802           lost |= k;
4803           /* Count lost digits before the decimal point.  */
4804           if (decflg == 0)
4805             nexp -= 1;
4806         }
4807       prec += 1;
4808       goto donchr;
4809     }
4810
4811   switch (*s)
4812     {
4813     case 'z':
4814       break;
4815     case 'E':
4816     case 'e':
4817       goto expnt;
4818     case '.':                   /* decimal point */
4819       if (decflg)
4820         goto error;
4821       ++decflg;
4822       break;
4823     case '-':
4824       nsign = 0xffff;
4825       if (sgnflg)
4826         goto error;
4827       ++sgnflg;
4828       break;
4829     case '+':
4830       if (sgnflg)
4831         goto error;
4832       ++sgnflg;
4833       break;
4834     case ',':
4835     case ' ':
4836     case '\0':
4837     case '\n':
4838     case '\r':
4839       goto daldone;
4840     case 'i':
4841     case 'I':
4842       goto infinite;
4843     default:
4844     error:
4845 #ifdef NANS
4846       einan (yy);
4847 #else
4848       mtherr ("asctoe", DOMAIN);
4849       eclear (yy);
4850 #endif
4851       goto aexit;
4852     }
4853  donchr:
4854   ++s;
4855   goto nxtcom;
4856
4857   /* Exponent interpretation */
4858  expnt:
4859
4860   esign = 1;
4861   exp = 0;
4862   ++s;
4863   /* check for + or - */
4864   if (*s == '-')
4865     {
4866       esign = -1;
4867       ++s;
4868     }
4869   if (*s == '+')
4870     ++s;
4871   while ((*s >= '0') && (*s <= '9'))
4872     {
4873       exp *= 10;
4874       exp += *s++ - '0';
4875       if (exp > -(MINDECEXP))
4876         {
4877           if (esign < 0)
4878             goto zero;
4879           else
4880             goto infinite;
4881         }
4882     }
4883   if (esign < 0)
4884     exp = -exp;
4885   if (exp > MAXDECEXP)
4886     {
4887  infinite:
4888       ecleaz (yy);
4889       yy[E] = 0x7fff;           /* infinity */
4890       goto aexit;
4891     }
4892   if (exp < MINDECEXP)
4893     {
4894  zero:
4895       ecleaz (yy);
4896       goto aexit;
4897     }
4898
4899  daldone:
4900   nexp = exp - nexp;
4901   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
4902   while ((nexp > 0) && (yy[2] == 0))
4903     {
4904       emovz (yy, xt);
4905       eshup1 (xt);
4906       eshup1 (xt);
4907       eaddm (yy, xt);
4908       eshup1 (xt);
4909       if (xt[2] != 0)
4910         break;
4911       nexp -= 1;
4912       emovz (xt, yy);
4913     }
4914   if ((k = enormlz (yy)) > NBITS)
4915     {
4916       ecleaz (yy);
4917       goto aexit;
4918     }
4919   lexp = (EXONE - 1 + NBITS) - k;
4920   emdnorm (yy, lost, 0, lexp, 64);
4921
4922   /* Convert to external format:
4923
4924      Multiply by 10**nexp.  If precision is 64 bits,
4925      the maximum relative error incurred in forming 10**n
4926      for 0 <= n <= 324 is 8.2e-20, at 10**180.
4927      For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
4928      For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */
4929
4930   lexp = yy[E];
4931   if (nexp == 0)
4932     {
4933       k = 0;
4934       goto expdon;
4935     }
4936   esign = 1;
4937   if (nexp < 0)
4938     {
4939       nexp = -nexp;
4940       esign = -1;
4941       if (nexp > 4096)
4942         {
4943           /* Punt.  Can't handle this without 2 divides. */
4944           emovi (etens[0], tt);
4945           lexp -= tt[E];
4946           k = edivm (tt, yy);
4947           lexp += EXONE;
4948           nexp -= 4096;
4949         }
4950     }
4951   p = &etens[NTEN][0];
4952   emov (eone, xt);
4953   exp = 1;
4954   do
4955     {
4956       if (exp & nexp)
4957         emul (p, xt, xt);
4958       p -= NE;
4959       exp = exp + exp;
4960     }
4961   while (exp <= MAXP);
4962
4963   emovi (xt, tt);
4964   if (esign < 0)
4965     {
4966       lexp -= tt[E];
4967       k = edivm (tt, yy);
4968       lexp += EXONE;
4969     }
4970   else
4971     {
4972       lexp += tt[E];
4973       k = emulm (tt, yy);
4974       lexp -= EXONE - 1;
4975     }
4976
4977  expdon:
4978
4979   /* Round and convert directly to the destination type */
4980   if (oprec == 53)
4981     lexp -= EXONE - 0x3ff;
4982 #ifdef IBM
4983   else if (oprec == 24 || oprec == 56)
4984     lexp -= EXONE - (0x41 << 2);
4985 #else
4986   else if (oprec == 24)
4987     lexp -= EXONE - 0177;
4988 #endif
4989 #ifdef DEC
4990   else if (oprec == 56)
4991     lexp -= EXONE - 0201;
4992 #endif
4993   rndprc = oprec;
4994   emdnorm (yy, k, 0, lexp, 64);
4995
4996  aexit:
4997
4998   rndprc = rndsav;
4999   yy[0] = nsign;
5000   switch (oprec)
5001     {
5002 #ifdef DEC
5003     case 56:
5004       todec (yy, y);            /* see etodec.c */
5005       break;
5006 #endif
5007 #ifdef IBM
5008     case 56:
5009       toibm (yy, y, DFmode);
5010       break;
5011 #endif
5012     case 53:
5013       toe53 (yy, y);
5014       break;
5015     case 24:
5016       toe24 (yy, y);
5017       break;
5018     case 64:
5019       toe64 (yy, y);
5020       break;
5021     case 113:
5022       toe113 (yy, y);
5023       break;
5024     case NBITS:
5025       emovo (yy, y);
5026       break;
5027     }
5028 }
5029
5030
5031
5032 /* Return Y = largest integer not greater than X (truncated toward minus
5033    infinity).  */
5034
5035 static unsigned EMUSHORT bmask[] =
5036 {
5037   0xffff,
5038   0xfffe,
5039   0xfffc,
5040   0xfff8,
5041   0xfff0,
5042   0xffe0,
5043   0xffc0,
5044   0xff80,
5045   0xff00,
5046   0xfe00,
5047   0xfc00,
5048   0xf800,
5049   0xf000,
5050   0xe000,
5051   0xc000,
5052   0x8000,
5053   0x0000,
5054 };
5055
5056 static void
5057 efloor (x, y)
5058      unsigned EMUSHORT x[], y[];
5059 {
5060   register unsigned EMUSHORT *p;
5061   int e, expon, i;
5062   unsigned EMUSHORT f[NE];
5063
5064   emov (x, f);                  /* leave in external format */
5065   expon = (int) f[NE - 1];
5066   e = (expon & 0x7fff) - (EXONE - 1);
5067   if (e <= 0)
5068     {
5069       eclear (y);
5070       goto isitneg;
5071     }
5072   /* number of bits to clear out */
5073   e = NBITS - e;
5074   emov (f, y);
5075   if (e <= 0)
5076     return;
5077
5078   p = &y[0];
5079   while (e >= 16)
5080     {
5081       *p++ = 0;
5082       e -= 16;
5083     }
5084   /* clear the remaining bits */
5085   *p &= bmask[e];
5086   /* truncate negatives toward minus infinity */
5087  isitneg:
5088
5089   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
5090     {
5091       for (i = 0; i < NE - 1; i++)
5092         {
5093           if (f[i] != y[i])
5094             {
5095               esub (eone, y, y);
5096               break;
5097             }
5098         }
5099     }
5100 }
5101
5102
5103 /* Return S and EXP such that  S * 2^EXP = X and .5 <= S < 1.
5104    For example, 1.1 = 0.55 * 2^1.  */
5105
5106 static void
5107 efrexp (x, exp, s)
5108      unsigned EMUSHORT x[];
5109      int *exp;
5110      unsigned EMUSHORT s[];
5111 {
5112   unsigned EMUSHORT xi[NI];
5113   EMULONG li;
5114
5115   emovi (x, xi);
5116   /*  Handle denormalized numbers properly using long integer exponent.  */
5117   li = (EMULONG) ((EMUSHORT) xi[1]);
5118
5119   if (li == 0)
5120     {
5121       li -= enormlz (xi);
5122     }
5123   xi[1] = 0x3ffe;
5124   emovo (xi, s);
5125   *exp = (int) (li - 0x3ffe);
5126 }
5127
5128 /* Return e type Y = X * 2^PWR2.  */
5129
5130 static void
5131 eldexp (x, pwr2, y)
5132      unsigned EMUSHORT x[];
5133      int pwr2;
5134      unsigned EMUSHORT y[];
5135 {
5136   unsigned EMUSHORT xi[NI];
5137   EMULONG li;
5138   int i;
5139
5140   emovi (x, xi);
5141   li = xi[1];
5142   li += pwr2;
5143   i = 0;
5144   emdnorm (xi, i, i, li, 64);
5145   emovo (xi, y);
5146 }
5147
5148
5149 /* C = remainder after dividing B by A, all e type values.
5150    Least significant integer quotient bits left in EQUOT.  */
5151
5152 static void
5153 eremain (a, b, c)
5154      unsigned EMUSHORT a[], b[], c[];
5155 {
5156   unsigned EMUSHORT den[NI], num[NI];
5157
5158 #ifdef NANS
5159   if (eisinf (b)
5160       || (ecmp (a, ezero) == 0)
5161       || eisnan (a)
5162       || eisnan (b))
5163     {
5164       enan (c, 0);
5165       return;
5166     }
5167 #endif
5168   if (ecmp (a, ezero) == 0)
5169     {
5170       mtherr ("eremain", SING);
5171       eclear (c);
5172       return;
5173     }
5174   emovi (a, den);
5175   emovi (b, num);
5176   eiremain (den, num);
5177   /* Sign of remainder = sign of quotient */
5178   if (a[0] == b[0])
5179     num[0] = 0;
5180   else
5181     num[0] = 0xffff;
5182   emovo (num, c);
5183 }
5184
5185 /*  Return quotient of exploded e-types NUM / DEN in EQUOT,
5186     remainder in NUM.  */
5187
5188 static void
5189 eiremain (den, num)
5190      unsigned EMUSHORT den[], num[];
5191 {
5192   EMULONG ld, ln;
5193   unsigned EMUSHORT j;
5194
5195   ld = den[E];
5196   ld -= enormlz (den);
5197   ln = num[E];
5198   ln -= enormlz (num);
5199   ecleaz (equot);
5200   while (ln >= ld)
5201     {
5202       if (ecmpm (den, num) <= 0)
5203         {
5204           esubm (den, num);
5205           j = 1;
5206         }
5207       else
5208           j = 0;
5209       eshup1 (equot);
5210       equot[NI - 1] |= j;
5211       eshup1 (num);
5212       ln -= 1;
5213     }
5214   emdnorm (num, 0, 0, ln, 0);
5215 }
5216
5217 /* Report an error condition CODE encountered in function NAME.
5218    CODE is one of the following:
5219
5220     Mnemonic        Value          Significance
5221
5222      DOMAIN            1       argument domain error
5223      SING              2       function singularity
5224      OVERFLOW          3       overflow range error
5225      UNDERFLOW         4       underflow range error
5226      TLOSS             5       total loss of precision
5227      PLOSS             6       partial loss of precision
5228      INVALID           7       NaN - producing operation
5229      EDOM             33       Unix domain error code
5230      ERANGE           34       Unix range error code
5231
5232    The order of appearance of the following messages is bound to the
5233    error codes defined above.  */
5234
5235 #define NMSGS 8
5236 static char *ermsg[NMSGS] =
5237 {
5238   "unknown",                    /* error code 0 */
5239   "domain",                     /* error code 1 */
5240   "singularity",                /* et seq.      */
5241   "overflow",
5242   "underflow",
5243   "total loss of precision",
5244   "partial loss of precision",
5245   "invalid operation"
5246 };
5247
5248 int merror = 0;
5249 extern int merror;
5250
5251 static void
5252 mtherr (name, code)
5253      char *name;
5254      int code;
5255 {
5256   char errstr[80];
5257
5258   /* The string passed by the calling program is supposed to be the
5259      name of the function in which the error occurred.
5260      The code argument selects which error message string will be printed.  */
5261
5262   if ((code <= 0) || (code >= NMSGS))
5263     code = 0;
5264   sprintf (errstr, " %s %s error", name, ermsg[code]);
5265   if (extra_warnings)
5266     warning (errstr);
5267   /* Set global error message word */
5268   merror = code + 1;
5269 }
5270
5271 #ifdef DEC
5272 /* Convert DEC double precision D to e type E.  */
5273
5274 static void
5275 dectoe (d, e)
5276      unsigned EMUSHORT *d;
5277      unsigned EMUSHORT *e;
5278 {
5279   unsigned EMUSHORT y[NI];
5280   register unsigned EMUSHORT r, *p;
5281
5282   ecleaz (y);                   /* start with a zero */
5283   p = y;                        /* point to our number */
5284   r = *d;                       /* get DEC exponent word */
5285   if (*d & (unsigned int) 0x8000)
5286     *p = 0xffff;                /* fill in our sign */
5287   ++p;                          /* bump pointer to our exponent word */
5288   r &= 0x7fff;                  /* strip the sign bit */
5289   if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
5290     goto done;
5291
5292
5293   r >>= 7;                      /* shift exponent word down 7 bits */
5294   r += EXONE - 0201;            /* subtract DEC exponent offset */
5295   /* add our e type exponent offset */
5296   *p++ = r;                     /* to form our exponent */
5297
5298   r = *d++;                     /* now do the high order mantissa */
5299   r &= 0177;                    /* strip off the DEC exponent and sign bits */
5300   r |= 0200;                    /* the DEC understood high order mantissa bit */
5301   *p++ = r;                     /* put result in our high guard word */
5302
5303   *p++ = *d++;                  /* fill in the rest of our mantissa */
5304   *p++ = *d++;
5305   *p = *d;
5306
5307   eshdn8 (y);                   /* shift our mantissa down 8 bits */
5308  done:
5309   emovo (y, e);
5310 }
5311
5312 /* Convert e type X to DEC double precision D.  */
5313
5314 static void
5315 etodec (x, d)
5316      unsigned EMUSHORT *x, *d;
5317 {
5318   unsigned EMUSHORT xi[NI];
5319   EMULONG exp;
5320   int rndsav;
5321
5322   emovi (x, xi);
5323   /* Adjust exponent for offsets.  */
5324   exp = (EMULONG) xi[E] - (EXONE - 0201);
5325   /* Round off to nearest or even.  */
5326   rndsav = rndprc;
5327   rndprc = 56;
5328   emdnorm (xi, 0, 0, exp, 64);
5329   rndprc = rndsav;
5330   todec (xi, d);
5331 }
5332
5333 /* Convert exploded e-type X, that has already been rounded to
5334    56-bit precision, to DEC format double Y.  */
5335
5336 static void
5337 todec (x, y)
5338      unsigned EMUSHORT *x, *y;
5339 {
5340   unsigned EMUSHORT i;
5341   unsigned EMUSHORT *p;
5342
5343   p = x;
5344   *y = 0;
5345   if (*p++)
5346     *y = 0100000;
5347   i = *p++;
5348   if (i == 0)
5349     {
5350       *y++ = 0;
5351       *y++ = 0;
5352       *y++ = 0;
5353       *y++ = 0;
5354       return;
5355     }
5356   if (i > 0377)
5357     {
5358       *y++ |= 077777;
5359       *y++ = 0xffff;
5360       *y++ = 0xffff;
5361       *y++ = 0xffff;
5362 #ifdef ERANGE
5363       errno = ERANGE;
5364 #endif
5365       return;
5366     }
5367   i &= 0377;
5368   i <<= 7;
5369   eshup8 (x);
5370   x[M] &= 0177;
5371   i |= x[M];
5372   *y++ |= i;
5373   *y++ = x[M + 1];
5374   *y++ = x[M + 2];
5375   *y++ = x[M + 3];
5376 }
5377 #endif /* DEC */
5378
5379 #ifdef IBM
5380 /* Convert IBM single/double precision to e type.  */
5381
5382 static void
5383 ibmtoe (d, e, mode)
5384      unsigned EMUSHORT *d;
5385      unsigned EMUSHORT *e;
5386      enum machine_mode mode;
5387 {
5388   unsigned EMUSHORT y[NI];
5389   register unsigned EMUSHORT r, *p;
5390   int rndsav;
5391
5392   ecleaz (y);                   /* start with a zero */
5393   p = y;                        /* point to our number */
5394   r = *d;                       /* get IBM exponent word */
5395   if (*d & (unsigned int) 0x8000)
5396     *p = 0xffff;                /* fill in our sign */
5397   ++p;                          /* bump pointer to our exponent word */
5398   r &= 0x7f00;                  /* strip the sign bit */
5399   r >>= 6;                      /* shift exponent word down 6 bits */
5400                                 /* in fact shift by 8 right and 2 left */
5401   r += EXONE - (0x41 << 2);     /* subtract IBM exponent offset */
5402                                 /* add our e type exponent offset */
5403   *p++ = r;                     /* to form our exponent */
5404
5405   *p++ = *d++ & 0xff;           /* now do the high order mantissa */
5406                                 /* strip off the IBM exponent and sign bits */
5407   if (mode != SFmode)           /* there are only 2 words in SFmode */
5408     {
5409       *p++ = *d++;              /* fill in the rest of our mantissa */
5410       *p++ = *d++;
5411     }
5412   *p = *d;
5413
5414   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5415     y[0] = y[E] = 0;
5416   else
5417     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
5418                               /* handle change in RADIX */
5419   emovo (y, e);
5420 }
5421
5422
5423
5424 /* Convert e type to IBM single/double precision.  */
5425
5426 static void
5427 etoibm (x, d, mode)
5428      unsigned EMUSHORT *x, *d;
5429      enum machine_mode mode;
5430 {
5431   unsigned EMUSHORT xi[NI];
5432   EMULONG exp;
5433   int rndsav;
5434
5435   emovi (x, xi);
5436   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));        /* adjust exponent for offsets */
5437                                                         /* round off to nearest or even */
5438   rndsav = rndprc;
5439   rndprc = 56;
5440   emdnorm (xi, 0, 0, exp, 64);
5441   rndprc = rndsav;
5442   toibm (xi, d, mode);
5443 }
5444
5445 static void
5446 toibm (x, y, mode)
5447      unsigned EMUSHORT *x, *y;
5448      enum machine_mode mode;
5449 {
5450   unsigned EMUSHORT i;
5451   unsigned EMUSHORT *p;
5452   int r;
5453
5454   p = x;
5455   *y = 0;
5456   if (*p++)
5457     *y = 0x8000;
5458   i = *p++;
5459   if (i == 0)
5460     {
5461       *y++ = 0;
5462       *y++ = 0;
5463       if (mode != SFmode)
5464         {
5465           *y++ = 0;
5466           *y++ = 0;
5467         }
5468       return;
5469     }
5470   r = i & 0x3;
5471   i >>= 2;
5472   if (i > 0x7f)
5473     {
5474       *y++ |= 0x7fff;
5475       *y++ = 0xffff;
5476       if (mode != SFmode)
5477         {
5478           *y++ = 0xffff;
5479           *y++ = 0xffff;
5480         }
5481 #ifdef ERANGE
5482       errno = ERANGE;
5483 #endif
5484       return;
5485     }
5486   i &= 0x7f;
5487   *y |= (i << 8);
5488   eshift (x, r + 5);
5489   *y++ |= x[M];
5490   *y++ = x[M + 1];
5491   if (mode != SFmode)
5492     {
5493       *y++ = x[M + 2];
5494       *y++ = x[M + 3];
5495     }
5496 }
5497 #endif /* IBM */
5498
5499 /* Output a binary NaN bit pattern in the target machine's format.  */
5500
5501 /* If special NaN bit patterns are required, define them in tm.h
5502    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
5503    patterns here. */
5504 #ifdef TFMODE_NAN
5505 TFMODE_NAN;
5506 #else
5507 #ifdef IEEE
5508 unsigned EMUSHORT TFbignan[8] =
5509  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5510 unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
5511 #endif
5512 #endif
5513
5514 #ifdef XFMODE_NAN
5515 XFMODE_NAN;
5516 #else
5517 #ifdef IEEE
5518 unsigned EMUSHORT XFbignan[6] =
5519  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5520 unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
5521 #endif
5522 #endif
5523
5524 #ifdef DFMODE_NAN
5525 DFMODE_NAN;
5526 #else
5527 #ifdef IEEE
5528 unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5529 unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
5530 #endif
5531 #endif
5532
5533 #ifdef SFMODE_NAN
5534 SFMODE_NAN;
5535 #else
5536 #ifdef IEEE
5537 unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
5538 unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
5539 #endif
5540 #endif
5541
5542
5543 static void
5544 make_nan (nan, sign, mode)
5545      unsigned EMUSHORT *nan;
5546      int sign;
5547      enum machine_mode mode;
5548 {
5549   int n;
5550   unsigned EMUSHORT *p;
5551
5552   switch (mode)
5553     {
5554 /* Possibly the `reserved operand' patterns on a VAX can be
5555    used like NaN's, but probably not in the same way as IEEE. */
5556 #if !defined(DEC) && !defined(IBM)
5557     case TFmode:
5558       n = 8;
5559       if (REAL_WORDS_BIG_ENDIAN)
5560         p = TFbignan;
5561       else
5562         p = TFlittlenan;
5563       break;
5564     case XFmode:
5565       n = 6;
5566       if (REAL_WORDS_BIG_ENDIAN)
5567         p = XFbignan;
5568       else
5569         p = XFlittlenan;
5570       break;
5571     case DFmode:
5572       n = 4;
5573       if (REAL_WORDS_BIG_ENDIAN)
5574         p = DFbignan;
5575       else
5576         p = DFlittlenan;
5577       break;
5578     case HFmode:
5579     case SFmode:
5580       n = 2;
5581       if (REAL_WORDS_BIG_ENDIAN)
5582         p = SFbignan;
5583       else
5584         p = SFlittlenan;
5585       break;
5586 #endif
5587     default:
5588       abort ();
5589     }
5590   if (REAL_WORDS_BIG_ENDIAN)
5591     *nan++ = (sign << 15) | *p++;
5592   while (--n != 0)
5593     *nan++ = *p++;
5594   if (! REAL_WORDS_BIG_ENDIAN)
5595     *nan = (sign << 15) | *p;
5596 }
5597
5598 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
5599    This is the inverse of the function `etarsingle' invoked by
5600    REAL_VALUE_TO_TARGET_SINGLE.  */
5601
5602 REAL_VALUE_TYPE
5603 ereal_from_float (f)
5604      HOST_WIDE_INT f;
5605 {
5606   REAL_VALUE_TYPE r;
5607   unsigned EMUSHORT s[2];
5608   unsigned EMUSHORT e[NE];
5609
5610   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
5611    This is the inverse operation to what the function `endian' does.  */
5612   if (REAL_WORDS_BIG_ENDIAN)
5613     {
5614       s[0] = (unsigned EMUSHORT) (f >> 16);
5615       s[1] = (unsigned EMUSHORT) f;
5616     }
5617   else
5618     {
5619       s[0] = (unsigned EMUSHORT) f;
5620       s[1] = (unsigned EMUSHORT) (f >> 16);
5621     }
5622   /* Convert and promote the target float to E-type. */
5623   e24toe (s, e);
5624   /* Output E-type to REAL_VALUE_TYPE. */
5625   PUT_REAL (e, &r);
5626   return r;
5627 }
5628
5629
5630 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
5631    This is the inverse of the function `etardouble' invoked by
5632    REAL_VALUE_TO_TARGET_DOUBLE.
5633
5634    The DFmode is stored as an array of HOST_WIDE_INT in the target's
5635    data format, with no holes in the bit packing.  The first element
5636    of the input array holds the bits that would come first in the
5637    target computer's memory.  */
5638
5639 REAL_VALUE_TYPE
5640 ereal_from_double (d)
5641      HOST_WIDE_INT d[];
5642 {
5643   REAL_VALUE_TYPE r;
5644   unsigned EMUSHORT s[4];
5645   unsigned EMUSHORT e[NE];
5646
5647   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
5648   if (REAL_WORDS_BIG_ENDIAN)
5649     {
5650       s[0] = (unsigned EMUSHORT) (d[0] >> 16);
5651       s[1] = (unsigned EMUSHORT) d[0];
5652 #if HOST_BITS_PER_WIDE_INT == 32
5653       s[2] = (unsigned EMUSHORT) (d[1] >> 16);
5654       s[3] = (unsigned EMUSHORT) d[1];
5655 #else
5656       /* In this case the entire target double is contained in the
5657          first array element.  The second element of the input is
5658          ignored.  */
5659       s[2] = (unsigned EMUSHORT) (d[0] >> 48);
5660       s[3] = (unsigned EMUSHORT) (d[0] >> 32);
5661 #endif
5662     }
5663   else
5664     {
5665       /* Target float words are little-endian.  */
5666       s[0] = (unsigned EMUSHORT) d[0];
5667       s[1] = (unsigned EMUSHORT) (d[0] >> 16);
5668 #if HOST_BITS_PER_WIDE_INT == 32
5669       s[2] = (unsigned EMUSHORT) d[1];
5670       s[3] = (unsigned EMUSHORT) (d[1] >> 16);
5671 #else
5672       s[2] = (unsigned EMUSHORT) (d[0] >> 32);
5673       s[3] = (unsigned EMUSHORT) (d[0] >> 48);
5674 #endif
5675     }
5676   /* Convert target double to E-type. */
5677   e53toe (s, e);
5678   /* Output E-type to REAL_VALUE_TYPE. */
5679   PUT_REAL (e, &r);
5680   return r;
5681 }
5682
5683
5684 /* Convert target computer unsigned 64-bit integer to e-type.
5685    The endian-ness of DImode follows the convention for integers,
5686    so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN.  */
5687
5688 static void
5689 uditoe (di, e)
5690      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
5691      unsigned EMUSHORT *e;
5692 {
5693   unsigned EMUSHORT yi[NI];
5694   int k;
5695
5696   ecleaz (yi);
5697   if (WORDS_BIG_ENDIAN)
5698     {
5699       for (k = M; k < M + 4; k++)
5700         yi[k] = *di++;
5701     }
5702   else
5703     {
5704       for (k = M + 3; k >= M; k--)
5705         yi[k] = *di++;
5706     }
5707   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
5708   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5709     ecleaz (yi);                /* it was zero */
5710   else
5711     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5712   emovo (yi, e);
5713 }
5714
5715 /* Convert target computer signed 64-bit integer to e-type. */
5716
5717 static void
5718 ditoe (di, e)
5719      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
5720      unsigned EMUSHORT *e;
5721 {
5722   unsigned EMULONG acc;
5723   unsigned EMUSHORT yi[NI];
5724   unsigned EMUSHORT carry;
5725   int k, sign;
5726
5727   ecleaz (yi);
5728   if (WORDS_BIG_ENDIAN)
5729     {
5730       for (k = M; k < M + 4; k++)
5731         yi[k] = *di++;
5732     }
5733   else
5734     {
5735       for (k = M + 3; k >= M; k--)
5736         yi[k] = *di++;
5737     }
5738   /* Take absolute value */
5739   sign = 0;
5740   if (yi[M] & 0x8000)
5741     {
5742       sign = 1;
5743       carry = 0;
5744       for (k = M + 3; k >= M; k--)
5745         {
5746           acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
5747           yi[k] = acc;
5748           carry = 0;
5749           if (acc & 0x10000)
5750             carry = 1;
5751         }
5752     }
5753   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
5754   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5755     ecleaz (yi);                /* it was zero */
5756   else
5757     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5758   emovo (yi, e);
5759   if (sign)
5760         eneg (e);
5761 }
5762
5763
5764 /* Convert e-type to unsigned 64-bit int. */
5765
5766 static void
5767 etoudi (x, i)
5768      unsigned EMUSHORT *x;
5769      unsigned EMUSHORT *i;
5770 {
5771   unsigned EMUSHORT xi[NI];
5772   int j, k;
5773
5774   emovi (x, xi);
5775   if (xi[0])
5776     {
5777       xi[M] = 0;
5778       goto noshift;
5779     }
5780   k = (int) xi[E] - (EXONE - 1);
5781   if (k <= 0)
5782     {
5783       for (j = 0; j < 4; j++)
5784         *i++ = 0;
5785       return;
5786     }
5787   if (k > 64)
5788     {
5789       for (j = 0; j < 4; j++)
5790         *i++ = 0xffff;
5791       if (extra_warnings)
5792         warning ("overflow on truncation to integer");
5793       return;
5794     }
5795   if (k > 16)
5796     {
5797       /* Shift more than 16 bits: first shift up k-16 mod 16,
5798          then shift up by 16's.  */
5799       j = k - ((k >> 4) << 4);
5800       if (j == 0)
5801         j = 16;
5802       eshift (xi, j);
5803       if (WORDS_BIG_ENDIAN)
5804         *i++ = xi[M];
5805       else
5806         {
5807           i += 3;
5808           *i-- = xi[M];
5809         }
5810       k -= j;
5811       do
5812         {
5813           eshup6 (xi);
5814           if (WORDS_BIG_ENDIAN)
5815             *i++ = xi[M];
5816           else
5817             *i-- = xi[M];
5818         }
5819       while ((k -= 16) > 0);
5820     }
5821   else
5822     {
5823         /* shift not more than 16 bits */
5824       eshift (xi, k);
5825
5826 noshift:
5827
5828       if (WORDS_BIG_ENDIAN)
5829         {
5830           i += 3;
5831           *i-- = xi[M];
5832           *i-- = 0;
5833           *i-- = 0;
5834           *i = 0;
5835         }
5836       else
5837         {
5838           *i++ = xi[M];
5839           *i++ = 0;
5840           *i++ = 0;
5841           *i = 0;
5842         }
5843     }
5844 }
5845
5846
5847 /* Convert e-type to signed 64-bit int. */
5848
5849 static void
5850 etodi (x, i)
5851      unsigned EMUSHORT *x;
5852      unsigned EMUSHORT *i;
5853 {
5854   unsigned EMULONG acc;
5855   unsigned EMUSHORT xi[NI];
5856   unsigned EMUSHORT carry;
5857   unsigned EMUSHORT *isave;
5858   int j, k;
5859
5860   emovi (x, xi);
5861   k = (int) xi[E] - (EXONE - 1);
5862   if (k <= 0)
5863     {
5864       for (j = 0; j < 4; j++)
5865         *i++ = 0;
5866       return;
5867     }
5868   if (k > 64)
5869     {
5870       for (j = 0; j < 4; j++)
5871         *i++ = 0xffff;
5872       if (extra_warnings)
5873         warning ("overflow on truncation to integer");
5874       return;
5875     }
5876   isave = i;
5877   if (k > 16)
5878     {
5879       /* Shift more than 16 bits: first shift up k-16 mod 16,
5880          then shift up by 16's.  */
5881       j = k - ((k >> 4) << 4);
5882       if (j == 0)
5883         j = 16;
5884       eshift (xi, j);
5885       if (WORDS_BIG_ENDIAN)
5886         *i++ = xi[M];
5887       else
5888         {
5889           i += 3;
5890           *i-- = xi[M];
5891         }
5892       k -= j;
5893       do
5894         {
5895           eshup6 (xi);
5896           if (WORDS_BIG_ENDIAN)
5897             *i++ = xi[M];
5898           else
5899             *i-- = xi[M];
5900         }
5901       while ((k -= 16) > 0);
5902     }
5903   else
5904     {
5905         /* shift not more than 16 bits */
5906       eshift (xi, k);
5907
5908       if (WORDS_BIG_ENDIAN)
5909         {
5910           i += 3;
5911           *i = xi[M];
5912           *i-- = 0;
5913           *i-- = 0;
5914           *i = 0;
5915         }
5916       else
5917         {
5918           *i++ = xi[M];
5919           *i++ = 0;
5920           *i++ = 0;
5921           *i = 0;
5922         }
5923     }
5924   /* Negate if negative */
5925   if (xi[0])
5926     {
5927       carry = 0;
5928       if (WORDS_BIG_ENDIAN)
5929         isave += 3;
5930       for (k = 0; k < 4; k++)
5931         {
5932           acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
5933           if (WORDS_BIG_ENDIAN)
5934             *isave-- = acc;
5935           else
5936             *isave++ = acc;
5937           carry = 0;
5938           if (acc & 0x10000)
5939             carry = 1;
5940         }
5941     }
5942 }
5943
5944
5945 /* Longhand square root routine. */
5946
5947
5948 static int esqinited = 0;
5949 static unsigned short sqrndbit[NI];
5950
5951 static void
5952 esqrt (x, y)
5953      unsigned EMUSHORT *x, *y;
5954 {
5955   unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
5956   EMULONG m, exp;
5957   int i, j, k, n, nlups;
5958
5959   if (esqinited == 0)
5960     {
5961       ecleaz (sqrndbit);
5962       sqrndbit[NI - 2] = 1;
5963       esqinited = 1;
5964     }
5965   /* Check for arg <= 0 */
5966   i = ecmp (x, ezero);
5967   if (i <= 0)
5968     {
5969       if (i == -1)
5970         {
5971           mtherr ("esqrt", DOMAIN);
5972           eclear (y);
5973         }
5974       else
5975         emov (x, y);
5976       return;
5977     }
5978
5979 #ifdef INFINITY
5980   if (eisinf (x))
5981     {
5982       eclear (y);
5983       einfin (y);
5984       return;
5985     }
5986 #endif
5987   /* Bring in the arg and renormalize if it is denormal. */
5988   emovi (x, xx);
5989   m = (EMULONG) xx[1];          /* local long word exponent */
5990   if (m == 0)
5991     m -= enormlz (xx);
5992
5993   /* Divide exponent by 2 */
5994   m -= 0x3ffe;
5995   exp = (unsigned short) ((m / 2) + 0x3ffe);
5996
5997   /* Adjust if exponent odd */
5998   if ((m & 1) != 0)
5999     {
6000       if (m > 0)
6001         exp += 1;
6002       eshdn1 (xx);
6003     }
6004
6005   ecleaz (sq);
6006   ecleaz (num);
6007   n = 8;                        /* get 8 bits of result per inner loop */
6008   nlups = rndprc;
6009   j = 0;
6010
6011   while (nlups > 0)
6012     {
6013       /* bring in next word of arg */
6014       if (j < NE)
6015         num[NI - 1] = xx[j + 3];
6016       /* Do additional bit on last outer loop, for roundoff. */
6017       if (nlups <= 8)
6018         n = nlups + 1;
6019       for (i = 0; i < n; i++)
6020         {
6021           /* Next 2 bits of arg */
6022           eshup1 (num);
6023           eshup1 (num);
6024           /* Shift up answer */
6025           eshup1 (sq);
6026           /* Make trial divisor */
6027           for (k = 0; k < NI; k++)
6028             temp[k] = sq[k];
6029           eshup1 (temp);
6030           eaddm (sqrndbit, temp);
6031           /* Subtract and insert answer bit if it goes in */
6032           if (ecmpm (temp, num) <= 0)
6033             {
6034               esubm (temp, num);
6035               sq[NI - 2] |= 1;
6036             }
6037         }
6038       nlups -= n;
6039       j += 1;
6040     }
6041
6042   /* Adjust for extra, roundoff loop done. */
6043   exp += (NBITS - 1) - rndprc;
6044
6045   /* Sticky bit = 1 if the remainder is nonzero. */
6046   k = 0;
6047   for (i = 3; i < NI; i++)
6048     k |= (int) num[i];
6049
6050   /* Renormalize and round off. */
6051   emdnorm (sq, k, 0, exp, 64);
6052   emovo (sq, y);
6053 }
6054 #endif /* EMU_NON_COMPILE not defined */
6055 \f
6056 /* Return the binary precision of the significand for a given
6057    floating point mode.  The mode can hold an integer value
6058    that many bits wide, without losing any bits.  */
6059
6060 int
6061 significand_size (mode)
6062      enum machine_mode mode;
6063 {
6064
6065 switch (mode)
6066   {
6067   case SFmode:
6068     return 24;
6069
6070   case DFmode:
6071 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6072     return 53;
6073 #else
6074 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6075     return 56;
6076 #else
6077 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6078     return 56;
6079 #else
6080     abort ();
6081 #endif
6082 #endif
6083 #endif
6084
6085   case XFmode:
6086     return 64;
6087   case TFmode:
6088     return 113;
6089
6090   default:
6091     abort ();
6092   }
6093 }