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      Intel long double denormals also lose one bit of precision.  */
2372   if ((exp <= 0) && (rndprc != NBITS)
2373       && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2374     {
2375       lost |= s[NI - 1] & 1;
2376       eshdn1 (s);
2377     }
2378   /* Clear out all bits below the rounding bit,
2379      remembering in r if any were nonzero.  */
2380   r = s[rw] & rmsk;
2381   if (rndprc < NBITS)
2382     {
2383       i = rw + 1;
2384       while (i < NI)
2385         {
2386           if (s[i])
2387             r |= 1;
2388           s[i] = 0;
2389           ++i;
2390         }
2391     }
2392   s[rw] &= ~rmsk;
2393   if ((r & rmbit) != 0)
2394     {
2395       if (r == rmbit)
2396         {
2397           if (lost == 0)
2398             {                   /* round to even */
2399               if ((s[re] & rebit) == 0)
2400                 goto mddone;
2401             }
2402           else
2403             {
2404               if (subflg != 0)
2405                 goto mddone;
2406             }
2407         }
2408       eaddm (rbit, s);
2409     }
2410  mddone:
2411 /* Undo the temporary shift for denormal values. */
2412   if ((exp <= 0) && (rndprc != NBITS)
2413       && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2414     {
2415       eshup1 (s);
2416     }
2417   if (s[2] != 0)
2418     {                           /* overflow on roundoff */
2419       eshdn1 (s);
2420       exp += 1;
2421     }
2422  mdfin:
2423   s[NI - 1] = 0;
2424   if (exp >= 32767L)
2425     {
2426 #ifndef INFINITY
2427     overf:
2428 #endif
2429 #ifdef INFINITY
2430       s[1] = 32767;
2431       for (i = 2; i < NI - 1; i++)
2432         s[i] = 0;
2433       if (extra_warnings)
2434         warning ("floating point overflow");
2435 #else
2436       s[1] = 32766;
2437       s[2] = 0;
2438       for (i = M + 1; i < NI - 1; i++)
2439         s[i] = 0xffff;
2440       s[NI - 1] = 0;
2441       if ((rndprc < 64) || (rndprc == 113))
2442         {
2443           s[rw] &= ~rmsk;
2444           if (rndprc == 24)
2445             {
2446               s[5] = 0;
2447               s[6] = 0;
2448             }
2449         }
2450 #endif
2451       return;
2452     }
2453   if (exp < 0)
2454     s[1] = 0;
2455   else
2456     s[1] = (unsigned EMUSHORT) exp;
2457 }
2458
2459 /*  Subtract.  C = B - A, all e type numbers.  */
2460
2461 static int subflg = 0;
2462
2463 static void
2464 esub (a, b, c)
2465      unsigned EMUSHORT *a, *b, *c;
2466 {
2467
2468 #ifdef NANS
2469   if (eisnan (a))
2470     {
2471       emov (a, c);
2472       return;
2473     }
2474   if (eisnan (b))
2475     {
2476       emov (b, c);
2477       return;
2478     }
2479 /* Infinity minus infinity is a NaN.
2480    Test for subtracting infinities of the same sign. */
2481   if (eisinf (a) && eisinf (b)
2482       && ((eisneg (a) ^ eisneg (b)) == 0))
2483     {
2484       mtherr ("esub", INVALID);
2485       enan (c, 0);
2486       return;
2487     }
2488 #endif
2489   subflg = 1;
2490   eadd1 (a, b, c);
2491 }
2492
2493 /* Add.  C = A + B, all e type. */
2494
2495 static void
2496 eadd (a, b, c)
2497      unsigned EMUSHORT *a, *b, *c;
2498 {
2499
2500 #ifdef NANS
2501 /* NaN plus anything is a NaN. */
2502   if (eisnan (a))
2503     {
2504       emov (a, c);
2505       return;
2506     }
2507   if (eisnan (b))
2508     {
2509       emov (b, c);
2510       return;
2511     }
2512 /* Infinity minus infinity is a NaN.
2513    Test for adding infinities of opposite signs. */
2514   if (eisinf (a) && eisinf (b)
2515       && ((eisneg (a) ^ eisneg (b)) != 0))
2516     {
2517       mtherr ("esub", INVALID);
2518       enan (c, 0);
2519       return;
2520     }
2521 #endif
2522   subflg = 0;
2523   eadd1 (a, b, c);
2524 }
2525
2526 /* Arithmetic common to both addition and subtraction.  */
2527
2528 static void
2529 eadd1 (a, b, c)
2530      unsigned EMUSHORT *a, *b, *c;
2531 {
2532   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2533   int i, lost, j, k;
2534   EMULONG lt, lta, ltb;
2535
2536 #ifdef INFINITY
2537   if (eisinf (a))
2538     {
2539       emov (a, c);
2540       if (subflg)
2541         eneg (c);
2542       return;
2543     }
2544   if (eisinf (b))
2545     {
2546       emov (b, c);
2547       return;
2548     }
2549 #endif
2550   emovi (a, ai);
2551   emovi (b, bi);
2552   if (subflg)
2553     ai[0] = ~ai[0];
2554
2555   /* compare exponents */
2556   lta = ai[E];
2557   ltb = bi[E];
2558   lt = lta - ltb;
2559   if (lt > 0L)
2560     {                           /* put the larger number in bi */
2561       emovz (bi, ci);
2562       emovz (ai, bi);
2563       emovz (ci, ai);
2564       ltb = bi[E];
2565       lt = -lt;
2566     }
2567   lost = 0;
2568   if (lt != 0L)
2569     {
2570       if (lt < (EMULONG) (-NBITS - 1))
2571         goto done;              /* answer same as larger addend */
2572       k = (int) lt;
2573       lost = eshift (ai, k);    /* shift the smaller number down */
2574     }
2575   else
2576     {
2577       /* exponents were the same, so must compare significands */
2578       i = ecmpm (ai, bi);
2579       if (i == 0)
2580         {                       /* the numbers are identical in magnitude */
2581           /* if different signs, result is zero */
2582           if (ai[0] != bi[0])
2583             {
2584               eclear (c);
2585               return;
2586             }
2587           /* if same sign, result is double */
2588           /* double denormalized tiny number */
2589           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2590             {
2591               eshup1 (bi);
2592               goto done;
2593             }
2594           /* add 1 to exponent unless both are zero! */
2595           for (j = 1; j < NI - 1; j++)
2596             {
2597               if (bi[j] != 0)
2598                 {
2599                   /* This could overflow, but let emovo take care of that. */
2600                   ltb += 1;
2601                   break;
2602                 }
2603             }
2604           bi[E] = (unsigned EMUSHORT) ltb;
2605           goto done;
2606         }
2607       if (i > 0)
2608         {                       /* put the larger number in bi */
2609           emovz (bi, ci);
2610           emovz (ai, bi);
2611           emovz (ci, ai);
2612         }
2613     }
2614   if (ai[0] == bi[0])
2615     {
2616       eaddm (ai, bi);
2617       subflg = 0;
2618     }
2619   else
2620     {
2621       esubm (ai, bi);
2622       subflg = 1;
2623     }
2624   emdnorm (bi, lost, subflg, ltb, 64);
2625
2626  done:
2627   emovo (bi, c);
2628 }
2629
2630 /* Divide: C = B/A, all e type.  */
2631
2632 static void
2633 ediv (a, b, c)
2634      unsigned EMUSHORT *a, *b, *c;
2635 {
2636   unsigned EMUSHORT ai[NI], bi[NI];
2637   int i;
2638   EMULONG lt, lta, ltb;
2639
2640 #ifdef NANS
2641 /* Return any NaN input. */
2642   if (eisnan (a))
2643     {
2644     emov (a, c);
2645     return;
2646     }
2647   if (eisnan (b))
2648     {
2649     emov (b, c);
2650     return;
2651     }
2652 /* Zero over zero, or infinity over infinity, is a NaN. */
2653   if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
2654       || (eisinf (a) && eisinf (b)))
2655     {
2656     mtherr ("ediv", INVALID);
2657     enan (c, eisneg (a) ^ eisneg (b));
2658     return;
2659     }
2660 #endif
2661 /* Infinity over anything else is infinity. */
2662 #ifdef INFINITY
2663   if (eisinf (b))
2664     {
2665       if (eisneg (a) ^ eisneg (b))
2666         *(c + (NE - 1)) = 0x8000;
2667       else
2668         *(c + (NE - 1)) = 0;
2669       einfin (c);
2670       return;
2671     }
2672 /* Anything else over infinity is zero. */
2673   if (eisinf (a))
2674     {
2675       eclear (c);
2676       return;
2677     }
2678 #endif
2679   emovi (a, ai);
2680   emovi (b, bi);
2681   lta = ai[E];
2682   ltb = bi[E];
2683   if (bi[E] == 0)
2684     {                           /* See if numerator is zero. */
2685       for (i = 1; i < NI - 1; i++)
2686         {
2687           if (bi[i] != 0)
2688             {
2689               ltb -= enormlz (bi);
2690               goto dnzro1;
2691             }
2692         }
2693       eclear (c);
2694       return;
2695     }
2696  dnzro1:
2697
2698   if (ai[E] == 0)
2699     {                           /* possible divide by zero */
2700       for (i = 1; i < NI - 1; i++)
2701         {
2702           if (ai[i] != 0)
2703             {
2704               lta -= enormlz (ai);
2705               goto dnzro2;
2706             }
2707         }
2708       if (ai[0] == bi[0])
2709         *(c + (NE - 1)) = 0;
2710       else
2711         *(c + (NE - 1)) = 0x8000;
2712 /* Divide by zero is not an invalid operation.
2713    It is a divide-by-zero operation!   */
2714       einfin (c);
2715       mtherr ("ediv", SING);
2716       return;
2717     }
2718  dnzro2:
2719
2720   i = edivm (ai, bi);
2721   /* calculate exponent */
2722   lt = ltb - lta + EXONE;
2723   emdnorm (bi, i, 0, lt, 64);
2724   /* set the sign */
2725   if (ai[0] == bi[0])
2726     bi[0] = 0;
2727   else
2728     bi[0] = 0Xffff;
2729   emovo (bi, c);
2730 }
2731
2732 /* Multiply e-types A and B, return e-type product C.   */
2733
2734 static void
2735 emul (a, b, c)
2736      unsigned EMUSHORT *a, *b, *c;
2737 {
2738   unsigned EMUSHORT ai[NI], bi[NI];
2739   int i, j;
2740   EMULONG lt, lta, ltb;
2741
2742 #ifdef NANS
2743 /* NaN times anything is the same NaN. */
2744   if (eisnan (a))
2745     {
2746     emov (a, c);
2747     return;
2748     }
2749   if (eisnan (b))
2750     {
2751     emov (b, c);
2752     return;
2753     }
2754 /* Zero times infinity is a NaN. */
2755   if ((eisinf (a) && (ecmp (b, ezero) == 0))
2756       || (eisinf (b) && (ecmp (a, ezero) == 0)))
2757     {
2758     mtherr ("emul", INVALID);
2759     enan (c, eisneg (a) ^ eisneg (b));
2760     return;
2761     }
2762 #endif
2763 /* Infinity times anything else is infinity. */
2764 #ifdef INFINITY
2765   if (eisinf (a) || eisinf (b))
2766     {
2767       if (eisneg (a) ^ eisneg (b))
2768         *(c + (NE - 1)) = 0x8000;
2769       else
2770         *(c + (NE - 1)) = 0;
2771       einfin (c);
2772       return;
2773     }
2774 #endif
2775   emovi (a, ai);
2776   emovi (b, bi);
2777   lta = ai[E];
2778   ltb = bi[E];
2779   if (ai[E] == 0)
2780     {
2781       for (i = 1; i < NI - 1; i++)
2782         {
2783           if (ai[i] != 0)
2784             {
2785               lta -= enormlz (ai);
2786               goto mnzer1;
2787             }
2788         }
2789       eclear (c);
2790       return;
2791     }
2792  mnzer1:
2793
2794   if (bi[E] == 0)
2795     {
2796       for (i = 1; i < NI - 1; i++)
2797         {
2798           if (bi[i] != 0)
2799             {
2800               ltb -= enormlz (bi);
2801               goto mnzer2;
2802             }
2803         }
2804       eclear (c);
2805       return;
2806     }
2807  mnzer2:
2808
2809   /* Multiply significands */
2810   j = emulm (ai, bi);
2811   /* calculate exponent */
2812   lt = lta + ltb - (EXONE - 1);
2813   emdnorm (bi, j, 0, lt, 64);
2814   /* calculate sign of product */
2815   if (ai[0] == bi[0])
2816     bi[0] = 0;
2817   else
2818     bi[0] = 0xffff;
2819   emovo (bi, c);
2820 }
2821
2822 /* Convert double precision PE to e-type Y.  */
2823
2824 static void
2825 e53toe (pe, y)
2826      unsigned EMUSHORT *pe, *y;
2827 {
2828 #ifdef DEC
2829
2830   dectoe (pe, y);
2831
2832 #else
2833 #ifdef IBM
2834
2835   ibmtoe (pe, y, DFmode);
2836
2837 #else
2838   register unsigned EMUSHORT r;
2839   register unsigned EMUSHORT *e, *p;
2840   unsigned EMUSHORT yy[NI];
2841   int denorm, k;
2842
2843   e = pe;
2844   denorm = 0;                   /* flag if denormalized number */
2845   ecleaz (yy);
2846   if (! REAL_WORDS_BIG_ENDIAN)
2847     e += 3;
2848   r = *e;
2849   yy[0] = 0;
2850   if (r & 0x8000)
2851     yy[0] = 0xffff;
2852   yy[M] = (r & 0x0f) | 0x10;
2853   r &= ~0x800f;                 /* strip sign and 4 significand bits */
2854 #ifdef INFINITY
2855   if (r == 0x7ff0)
2856     {
2857 #ifdef NANS
2858       if (! REAL_WORDS_BIG_ENDIAN)
2859         {
2860           if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
2861               || (pe[1] != 0) || (pe[0] != 0))
2862             {
2863               enan (y, yy[0] != 0);
2864               return;
2865             }
2866         }
2867       else
2868         {
2869           if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
2870               || (pe[2] != 0) || (pe[3] != 0))
2871             {
2872               enan (y, yy[0] != 0);
2873               return;
2874             }
2875         }
2876 #endif  /* NANS */
2877       eclear (y);
2878       einfin (y);
2879       if (yy[0])
2880         eneg (y);
2881       return;
2882     }
2883 #endif  /* INFINITY */
2884   r >>= 4;
2885   /* If zero exponent, then the significand is denormalized.
2886      So take back the understood high significand bit. */
2887
2888   if (r == 0)
2889     {
2890       denorm = 1;
2891       yy[M] &= ~0x10;
2892     }
2893   r += EXONE - 01777;
2894   yy[E] = r;
2895   p = &yy[M + 1];
2896 #ifdef IEEE
2897   if (! REAL_WORDS_BIG_ENDIAN)
2898     {
2899       *p++ = *(--e);
2900       *p++ = *(--e);
2901       *p++ = *(--e);
2902     }
2903   else
2904     {
2905       ++e;
2906       *p++ = *e++;
2907       *p++ = *e++;
2908       *p++ = *e++;
2909     }
2910 #endif
2911   eshift (yy, -5);
2912   if (denorm)
2913     {                           /* if zero exponent, then normalize the significand */
2914       if ((k = enormlz (yy)) > NBITS)
2915         ecleazs (yy);
2916       else
2917         yy[E] -= (unsigned EMUSHORT) (k - 1);
2918     }
2919   emovo (yy, y);
2920 #endif /* not IBM */
2921 #endif /* not DEC */
2922 }
2923
2924 /* Convert double extended precision float PE to e type Y.  */
2925
2926 static void
2927 e64toe (pe, y)
2928      unsigned EMUSHORT *pe, *y;
2929 {
2930   unsigned EMUSHORT yy[NI];
2931   unsigned EMUSHORT *e, *p, *q;
2932   int i;
2933
2934   e = pe;
2935   p = yy;
2936   for (i = 0; i < NE - 5; i++)
2937     *p++ = 0;
2938 /* This precision is not ordinarily supported on DEC or IBM. */
2939 #ifdef DEC
2940   for (i = 0; i < 5; i++)
2941     *p++ = *e++;
2942 #endif
2943 #ifdef IBM
2944   p = &yy[0] + (NE - 1);
2945   *p-- = *e++;
2946   ++e;
2947   for (i = 0; i < 5; i++)
2948     *p-- = *e++;
2949 #endif
2950 #ifdef IEEE
2951   if (! REAL_WORDS_BIG_ENDIAN)
2952     {
2953       for (i = 0; i < 5; i++)
2954         *p++ = *e++;
2955
2956       /* For denormal long double Intel format, shift significand up one
2957          -- but only if the top significand bit is zero.  A top bit of 1
2958          is "pseudodenormal" when the exponent is zero.  */
2959       if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
2960         {
2961           unsigned EMUSHORT temp[NI];
2962
2963           emovi(yy, temp);
2964           eshup1(temp);
2965           emovo(temp,y);
2966           return;
2967         }
2968     }
2969   else
2970     {
2971       p = &yy[0] + (NE - 1);
2972       *p-- = *e++;
2973       ++e;
2974       for (i = 0; i < 4; i++)
2975         *p-- = *e++;
2976     }
2977 #endif
2978 #ifdef INFINITY
2979   /* Point to the exponent field and check max exponent cases.  */
2980   p = &yy[NE - 1];
2981   if (*p == 0x7fff)
2982     {
2983 #ifdef NANS
2984       if (! REAL_WORDS_BIG_ENDIAN)
2985         {
2986           for (i = 0; i < 4; i++)
2987             {
2988               if ((i != 3 && pe[i] != 0)
2989                   /* Anything but 0x8000 here, including 0, is a NaN.  */
2990                   || (i == 3 && pe[i] != 0x8000))
2991                 {
2992                   enan (y, (*p & 0x8000) != 0);
2993                   return;
2994                 }
2995             }
2996         }
2997       else
2998         {
2999           for (i = 1; i <= 4; i++)
3000             {
3001               if (pe[i] != 0)
3002                 {
3003                   enan (y, (*p & 0x8000) != 0);
3004                   return;
3005                 }
3006             }
3007         }
3008 #endif /* NANS */
3009       eclear (y);
3010       einfin (y);
3011       if (*p & 0x8000)
3012         eneg (y);
3013       return;
3014     }
3015 #endif  /* INFINITY */
3016   p = yy;
3017   q = y;
3018   for (i = 0; i < NE; i++)
3019     *q++ = *p++;
3020 }
3021
3022 /* Convert 128-bit long double precision float PE to e type Y.  */
3023
3024 static void
3025 e113toe (pe, y)
3026      unsigned EMUSHORT *pe, *y;
3027 {
3028   register unsigned EMUSHORT r;
3029   unsigned EMUSHORT *e, *p;
3030   unsigned EMUSHORT yy[NI];
3031   int denorm, i;
3032
3033   e = pe;
3034   denorm = 0;
3035   ecleaz (yy);
3036 #ifdef IEEE
3037   if (! REAL_WORDS_BIG_ENDIAN)
3038     e += 7;
3039 #endif
3040   r = *e;
3041   yy[0] = 0;
3042   if (r & 0x8000)
3043     yy[0] = 0xffff;
3044   r &= 0x7fff;
3045 #ifdef INFINITY
3046   if (r == 0x7fff)
3047     {
3048 #ifdef NANS
3049       if (! REAL_WORDS_BIG_ENDIAN)
3050         {
3051           for (i = 0; i < 7; i++)
3052             {
3053               if (pe[i] != 0)
3054                 {
3055                   enan (y, yy[0] != 0);
3056                   return;
3057                 }
3058             }
3059         }
3060       else
3061         {
3062           for (i = 1; i < 8; i++)
3063             {
3064               if (pe[i] != 0)
3065                 {
3066                   enan (y, yy[0] != 0);
3067                   return;
3068                 }
3069             }
3070         }
3071 #endif /* NANS */
3072       eclear (y);
3073       einfin (y);
3074       if (yy[0])
3075         eneg (y);
3076       return;
3077     }
3078 #endif  /* INFINITY */
3079   yy[E] = r;
3080   p = &yy[M + 1];
3081 #ifdef IEEE
3082   if (! REAL_WORDS_BIG_ENDIAN)
3083     {
3084       for (i = 0; i < 7; i++)
3085         *p++ = *(--e);
3086     }
3087   else
3088     {
3089       ++e;
3090       for (i = 0; i < 7; i++)
3091         *p++ = *e++;
3092     }
3093 #endif
3094 /* If denormal, remove the implied bit; else shift down 1. */
3095   if (r == 0)
3096     {
3097       yy[M] = 0;
3098     }
3099   else
3100     {
3101       yy[M] = 1;
3102       eshift (yy, -1);
3103     }
3104   emovo (yy, y);
3105 }
3106
3107 /* Convert single precision float PE to e type Y.  */
3108
3109 static void
3110 e24toe (pe, y)
3111      unsigned EMUSHORT *pe, *y;
3112 {
3113 #ifdef IBM
3114
3115   ibmtoe (pe, y, SFmode);
3116
3117 #else
3118   register unsigned EMUSHORT r;
3119   register unsigned EMUSHORT *e, *p;
3120   unsigned EMUSHORT yy[NI];
3121   int denorm, k;
3122
3123   e = pe;
3124   denorm = 0;                   /* flag if denormalized number */
3125   ecleaz (yy);
3126 #ifdef IEEE
3127   if (! REAL_WORDS_BIG_ENDIAN)
3128     e += 1;
3129 #endif
3130 #ifdef DEC
3131   e += 1;
3132 #endif
3133   r = *e;
3134   yy[0] = 0;
3135   if (r & 0x8000)
3136     yy[0] = 0xffff;
3137   yy[M] = (r & 0x7f) | 0200;
3138   r &= ~0x807f;                 /* strip sign and 7 significand bits */
3139 #ifdef INFINITY
3140   if (r == 0x7f80)
3141     {
3142 #ifdef NANS
3143       if (REAL_WORDS_BIG_ENDIAN)
3144         {
3145           if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
3146             {
3147               enan (y, yy[0] != 0);
3148               return;
3149             }
3150         }
3151       else
3152         {
3153           if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
3154             {
3155               enan (y, yy[0] != 0);
3156               return;
3157             }
3158         }
3159 #endif  /* NANS */
3160       eclear (y);
3161       einfin (y);
3162       if (yy[0])
3163         eneg (y);
3164       return;
3165     }
3166 #endif  /* INFINITY */
3167   r >>= 7;
3168   /* If zero exponent, then the significand is denormalized.
3169      So take back the understood high significand bit. */
3170   if (r == 0)
3171     {
3172       denorm = 1;
3173       yy[M] &= ~0200;
3174     }
3175   r += EXONE - 0177;
3176   yy[E] = r;
3177   p = &yy[M + 1];
3178 #ifdef DEC
3179   *p++ = *(--e);
3180 #endif
3181 #ifdef IEEE
3182   if (! REAL_WORDS_BIG_ENDIAN)
3183     *p++ = *(--e);
3184   else
3185     {
3186       ++e;
3187       *p++ = *e++;
3188     }
3189 #endif
3190   eshift (yy, -8);
3191   if (denorm)
3192     {                           /* if zero exponent, then normalize the significand */
3193       if ((k = enormlz (yy)) > NBITS)
3194         ecleazs (yy);
3195       else
3196         yy[E] -= (unsigned EMUSHORT) (k - 1);
3197     }
3198   emovo (yy, y);
3199 #endif /* not IBM */
3200 }
3201
3202 /* Convert e-type X to IEEE 128-bit long double format E.  */
3203
3204 static void
3205 etoe113 (x, e)
3206      unsigned EMUSHORT *x, *e;
3207 {
3208   unsigned EMUSHORT xi[NI];
3209   EMULONG exp;
3210   int rndsav;
3211
3212 #ifdef NANS
3213   if (eisnan (x))
3214     {
3215       make_nan (e, eisneg (x), TFmode);
3216       return;
3217     }
3218 #endif
3219   emovi (x, xi);
3220   exp = (EMULONG) xi[E];
3221 #ifdef INFINITY
3222   if (eisinf (x))
3223     goto nonorm;
3224 #endif
3225   /* round off to nearest or even */
3226   rndsav = rndprc;
3227   rndprc = 113;
3228   emdnorm (xi, 0, 0, exp, 64);
3229   rndprc = rndsav;
3230  nonorm:
3231   toe113 (xi, e);
3232 }
3233
3234 /* Convert exploded e-type X, that has already been rounded to
3235    113-bit precision, to IEEE 128-bit long double format Y.  */
3236
3237 static void
3238 toe113 (a, b)
3239      unsigned EMUSHORT *a, *b;
3240 {
3241   register unsigned EMUSHORT *p, *q;
3242   unsigned EMUSHORT i;
3243
3244 #ifdef NANS
3245   if (eiisnan (a))
3246     {
3247       make_nan (b, eiisneg (a), TFmode);
3248       return;
3249     }
3250 #endif
3251   p = a;
3252   if (REAL_WORDS_BIG_ENDIAN)
3253     q = b;
3254   else
3255     q = b + 7;                  /* point to output exponent */
3256
3257   /* If not denormal, delete the implied bit. */
3258   if (a[E] != 0)
3259     {
3260       eshup1 (a);
3261     }
3262   /* combine sign and exponent */
3263   i = *p++;
3264   if (REAL_WORDS_BIG_ENDIAN)
3265     {
3266       if (i)
3267         *q++ = *p++ | 0x8000;
3268       else
3269         *q++ = *p++;
3270     }
3271   else
3272     {
3273       if (i)
3274         *q-- = *p++ | 0x8000;
3275       else
3276         *q-- = *p++;
3277     }
3278   /* skip over guard word */
3279   ++p;
3280   /* move the significand */
3281   if (REAL_WORDS_BIG_ENDIAN)
3282     {
3283       for (i = 0; i < 7; i++)
3284         *q++ = *p++;
3285     }
3286   else
3287     {
3288       for (i = 0; i < 7; i++)
3289         *q-- = *p++;
3290     }
3291 }
3292
3293 /* Convert e-type X to IEEE double extended format E.  */
3294
3295 static void
3296 etoe64 (x, e)
3297      unsigned EMUSHORT *x, *e;
3298 {
3299   unsigned EMUSHORT xi[NI];
3300   EMULONG exp;
3301   int rndsav;
3302
3303 #ifdef NANS
3304   if (eisnan (x))
3305     {
3306       make_nan (e, eisneg (x), XFmode);
3307       return;
3308     }
3309 #endif
3310   emovi (x, xi);
3311   /* adjust exponent for offset */
3312   exp = (EMULONG) xi[E];
3313 #ifdef INFINITY
3314   if (eisinf (x))
3315     goto nonorm;
3316 #endif
3317   /* round off to nearest or even */
3318   rndsav = rndprc;
3319   rndprc = 64;
3320   emdnorm (xi, 0, 0, exp, 64);
3321   rndprc = rndsav;
3322  nonorm:
3323   toe64 (xi, e);
3324 }
3325
3326 /* Convert exploded e-type X, that has already been rounded to
3327    64-bit precision, to IEEE double extended format Y.  */
3328
3329 static void
3330 toe64 (a, b)
3331      unsigned EMUSHORT *a, *b;
3332 {
3333   register unsigned EMUSHORT *p, *q;
3334   unsigned EMUSHORT i;
3335
3336 #ifdef NANS
3337   if (eiisnan (a))
3338     {
3339       make_nan (b, eiisneg (a), XFmode);
3340       return;
3341     }
3342 #endif
3343   /* Shift denormal long double Intel format significand down one bit.  */
3344   if ((a[E] == 0) && ! REAL_WORDS_BIG_ENDIAN)
3345     eshdn1 (a);
3346   p = a;
3347 #ifdef IBM
3348   q = b;
3349 #endif
3350 #ifdef DEC
3351   q = b + 4;
3352 #endif
3353 #ifdef IEEE
3354   if (REAL_WORDS_BIG_ENDIAN)
3355     q = b;
3356   else
3357     {
3358       q = b + 4;                        /* point to output exponent */
3359 #if LONG_DOUBLE_TYPE_SIZE == 96
3360       /* Clear the last two bytes of 12-byte Intel format */
3361       *(q+1) = 0;
3362 #endif
3363     }
3364 #endif
3365
3366   /* combine sign and exponent */
3367   i = *p++;
3368 #ifdef IBM
3369   if (i)
3370     *q++ = *p++ | 0x8000;
3371   else
3372     *q++ = *p++;
3373   *q++ = 0;
3374 #endif
3375 #ifdef DEC
3376   if (i)
3377     *q-- = *p++ | 0x8000;
3378   else
3379     *q-- = *p++;
3380 #endif
3381 #ifdef IEEE
3382   if (REAL_WORDS_BIG_ENDIAN)
3383     {
3384       if (i)
3385         *q++ = *p++ | 0x8000;
3386       else
3387         *q++ = *p++;
3388       *q++ = 0;
3389     }
3390   else
3391     {
3392       if (i)
3393         *q-- = *p++ | 0x8000;
3394       else
3395         *q-- = *p++;
3396     }
3397 #endif
3398   /* skip over guard word */
3399   ++p;
3400   /* move the significand */
3401 #ifdef IBM
3402   for (i = 0; i < 4; i++)
3403     *q++ = *p++;
3404 #endif
3405 #ifdef DEC
3406   for (i = 0; i < 4; i++)
3407     *q-- = *p++;
3408 #endif
3409 #ifdef IEEE
3410   if (REAL_WORDS_BIG_ENDIAN)
3411     {
3412       for (i = 0; i < 4; i++)
3413         *q++ = *p++;
3414     }
3415   else
3416     {
3417 #ifdef INFINITY
3418       if (eiisinf (a))
3419         {
3420           /* Intel long double infinity significand.  */
3421           *q-- = 0x8000;
3422           *q-- = 0;
3423           *q-- = 0;
3424           *q = 0;
3425           return;
3426         }
3427 #endif
3428       for (i = 0; i < 4; i++)
3429         *q-- = *p++;
3430     }
3431 #endif
3432 }
3433
3434 /* e type to double precision.  */
3435
3436 #ifdef DEC
3437 /* Convert e-type X to DEC-format double E.  */
3438
3439 static void
3440 etoe53 (x, e)
3441      unsigned EMUSHORT *x, *e;
3442 {
3443   etodec (x, e);                /* see etodec.c */
3444 }
3445
3446 /* Convert exploded e-type X, that has already been rounded to
3447    56-bit double precision, to DEC double Y.  */
3448
3449 static void
3450 toe53 (x, y)
3451      unsigned EMUSHORT *x, *y;
3452 {
3453   todec (x, y);
3454 }
3455
3456 #else
3457 #ifdef IBM
3458 /* Convert e-type X to IBM 370-format double E.  */
3459
3460 static void
3461 etoe53 (x, e)
3462      unsigned EMUSHORT *x, *e;
3463 {
3464   etoibm (x, e, DFmode);
3465 }
3466
3467 /* Convert exploded e-type X, that has already been rounded to
3468    56-bit precision, to IBM 370 double Y.  */
3469
3470 static void
3471 toe53 (x, y)
3472      unsigned EMUSHORT *x, *y;
3473 {
3474   toibm (x, y, DFmode);
3475 }
3476
3477 #else  /* it's neither DEC nor IBM */
3478
3479 /* Convert e-type X to IEEE double E.  */
3480
3481 static void
3482 etoe53 (x, e)
3483      unsigned EMUSHORT *x, *e;
3484 {
3485   unsigned EMUSHORT xi[NI];
3486   EMULONG exp;
3487   int rndsav;
3488
3489 #ifdef NANS
3490   if (eisnan (x))
3491     {
3492       make_nan (e, eisneg (x), DFmode);
3493       return;
3494     }
3495 #endif
3496   emovi (x, xi);
3497   /* adjust exponent for offsets */
3498   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
3499 #ifdef INFINITY
3500   if (eisinf (x))
3501     goto nonorm;
3502 #endif
3503   /* round off to nearest or even */
3504   rndsav = rndprc;
3505   rndprc = 53;
3506   emdnorm (xi, 0, 0, exp, 64);
3507   rndprc = rndsav;
3508  nonorm:
3509   toe53 (xi, e);
3510 }
3511
3512 /* Convert exploded e-type X, that has already been rounded to
3513    53-bit precision, to IEEE double Y.  */
3514
3515 static void
3516 toe53 (x, y)
3517      unsigned EMUSHORT *x, *y;
3518 {
3519   unsigned EMUSHORT i;
3520   unsigned EMUSHORT *p;
3521
3522 #ifdef NANS
3523   if (eiisnan (x))
3524     {
3525       make_nan (y, eiisneg (x), DFmode);
3526       return;
3527     }
3528 #endif
3529   p = &x[0];
3530 #ifdef IEEE
3531   if (! REAL_WORDS_BIG_ENDIAN)
3532     y += 3;
3533 #endif
3534   *y = 0;                       /* output high order */
3535   if (*p++)
3536     *y = 0x8000;                /* output sign bit */
3537
3538   i = *p++;
3539   if (i >= (unsigned int) 2047)
3540     {                           /* Saturate at largest number less than infinity. */
3541 #ifdef INFINITY
3542       *y |= 0x7ff0;
3543       if (! REAL_WORDS_BIG_ENDIAN)
3544         {
3545           *(--y) = 0;
3546           *(--y) = 0;
3547           *(--y) = 0;
3548         }
3549       else
3550         {
3551           ++y;
3552           *y++ = 0;
3553           *y++ = 0;
3554           *y++ = 0;
3555         }
3556 #else
3557       *y |= (unsigned EMUSHORT) 0x7fef;
3558       if (! REAL_WORDS_BIG_ENDIAN)
3559         {
3560           *(--y) = 0xffff;
3561           *(--y) = 0xffff;
3562           *(--y) = 0xffff;
3563         }
3564       else
3565         {
3566           ++y;
3567           *y++ = 0xffff;
3568           *y++ = 0xffff;
3569           *y++ = 0xffff;
3570         }
3571 #endif
3572       return;
3573     }
3574   if (i == 0)
3575     {
3576       eshift (x, 4);
3577     }
3578   else
3579     {
3580       i <<= 4;
3581       eshift (x, 5);
3582     }
3583   i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
3584   *y |= (unsigned EMUSHORT) i;  /* high order output already has sign bit set */
3585   if (! REAL_WORDS_BIG_ENDIAN)
3586     {
3587       *(--y) = *p++;
3588       *(--y) = *p++;
3589       *(--y) = *p;
3590     }
3591   else
3592     {
3593       ++y;
3594       *y++ = *p++;
3595       *y++ = *p++;
3596       *y++ = *p++;
3597     }
3598 }
3599
3600 #endif /* not IBM */
3601 #endif /* not DEC */
3602
3603
3604
3605 /* e type to single precision.  */
3606
3607 #ifdef IBM
3608 /* Convert e-type X to IBM 370 float E.  */
3609
3610 static void
3611 etoe24 (x, e)
3612      unsigned EMUSHORT *x, *e;
3613 {
3614   etoibm (x, e, SFmode);
3615 }
3616
3617 /* Convert exploded e-type X, that has already been rounded to
3618    float precision, to IBM 370 float Y.  */
3619
3620 static void
3621 toe24 (x, y)
3622      unsigned EMUSHORT *x, *y;
3623 {
3624   toibm (x, y, SFmode);
3625 }
3626
3627 #else
3628 /* Convert e-type X to IEEE float E.  DEC float is the same as IEEE float.  */
3629
3630 static void
3631 etoe24 (x, e)
3632      unsigned EMUSHORT *x, *e;
3633 {
3634   EMULONG exp;
3635   unsigned EMUSHORT xi[NI];
3636   int rndsav;
3637
3638 #ifdef NANS
3639   if (eisnan (x))
3640     {
3641       make_nan (e, eisneg (x), SFmode);
3642       return;
3643     }
3644 #endif
3645   emovi (x, xi);
3646   /* adjust exponent for offsets */
3647   exp = (EMULONG) xi[E] - (EXONE - 0177);
3648 #ifdef INFINITY
3649   if (eisinf (x))
3650     goto nonorm;
3651 #endif
3652   /* round off to nearest or even */
3653   rndsav = rndprc;
3654   rndprc = 24;
3655   emdnorm (xi, 0, 0, exp, 64);
3656   rndprc = rndsav;
3657  nonorm:
3658   toe24 (xi, e);
3659 }
3660
3661 /* Convert exploded e-type X, that has already been rounded to
3662    float precision, to IEEE float Y.  */
3663
3664 static void
3665 toe24 (x, y)
3666      unsigned EMUSHORT *x, *y;
3667 {
3668   unsigned EMUSHORT i;
3669   unsigned EMUSHORT *p;
3670
3671 #ifdef NANS
3672   if (eiisnan (x))
3673     {
3674       make_nan (y, eiisneg (x), SFmode);
3675       return;
3676     }
3677 #endif
3678   p = &x[0];
3679 #ifdef IEEE
3680   if (! REAL_WORDS_BIG_ENDIAN)
3681     y += 1;
3682 #endif
3683 #ifdef DEC
3684   y += 1;
3685 #endif
3686   *y = 0;                       /* output high order */
3687   if (*p++)
3688     *y = 0x8000;                /* output sign bit */
3689
3690   i = *p++;
3691 /* Handle overflow cases. */
3692   if (i >= 255)
3693     {
3694 #ifdef INFINITY
3695       *y |= (unsigned EMUSHORT) 0x7f80;
3696 #ifdef DEC
3697       *(--y) = 0;
3698 #endif
3699 #ifdef IEEE
3700       if (! REAL_WORDS_BIG_ENDIAN)
3701         *(--y) = 0;
3702       else
3703         {
3704           ++y;
3705           *y = 0;
3706         }
3707 #endif
3708 #else  /* no INFINITY */
3709       *y |= (unsigned EMUSHORT) 0x7f7f;
3710 #ifdef DEC
3711       *(--y) = 0xffff;
3712 #endif
3713 #ifdef IEEE
3714       if (! REAL_WORDS_BIG_ENDIAN)
3715         *(--y) = 0xffff;
3716       else
3717         {
3718           ++y;
3719           *y = 0xffff;
3720         }
3721 #endif
3722 #ifdef ERANGE
3723       errno = ERANGE;
3724 #endif
3725 #endif  /* no INFINITY */
3726       return;
3727     }
3728   if (i == 0)
3729     {
3730       eshift (x, 7);
3731     }
3732   else
3733     {
3734       i <<= 7;
3735       eshift (x, 8);
3736     }
3737   i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
3738   /* High order output already has sign bit set.  */
3739   *y |= i;
3740 #ifdef DEC
3741   *(--y) = *p;
3742 #endif
3743 #ifdef IEEE
3744   if (! REAL_WORDS_BIG_ENDIAN)
3745     *(--y) = *p;
3746   else
3747     {
3748       ++y;
3749       *y = *p;
3750     }
3751 #endif
3752 }
3753 #endif  /* not IBM */
3754
3755 /* Compare two e type numbers.
3756    Return +1 if a > b
3757            0 if a == b
3758           -1 if a < b
3759           -2 if either a or b is a NaN.  */
3760
3761 static int
3762 ecmp (a, b)
3763      unsigned EMUSHORT *a, *b;
3764 {
3765   unsigned EMUSHORT ai[NI], bi[NI];
3766   register unsigned EMUSHORT *p, *q;
3767   register int i;
3768   int msign;
3769
3770 #ifdef NANS
3771   if (eisnan (a)  || eisnan (b))
3772       return (-2);
3773 #endif
3774   emovi (a, ai);
3775   p = ai;
3776   emovi (b, bi);
3777   q = bi;
3778
3779   if (*p != *q)
3780     {                           /* the signs are different */
3781       /* -0 equals + 0 */
3782       for (i = 1; i < NI - 1; i++)
3783         {
3784           if (ai[i] != 0)
3785             goto nzro;
3786           if (bi[i] != 0)
3787             goto nzro;
3788         }
3789       return (0);
3790     nzro:
3791       if (*p == 0)
3792         return (1);
3793       else
3794         return (-1);
3795     }
3796   /* both are the same sign */
3797   if (*p == 0)
3798     msign = 1;
3799   else
3800     msign = -1;
3801   i = NI - 1;
3802   do
3803     {
3804       if (*p++ != *q++)
3805         {
3806           goto diff;
3807         }
3808     }
3809   while (--i > 0);
3810
3811   return (0);                   /* equality */
3812
3813  diff:
3814
3815   if (*(--p) > *(--q))
3816     return (msign);             /* p is bigger */
3817   else
3818     return (-msign);            /* p is littler */
3819 }
3820
3821 /* Find e-type nearest integer to X, as floor (X + 0.5).  */
3822
3823 static void
3824 eround (x, y)
3825      unsigned EMUSHORT *x, *y;
3826 {
3827   eadd (ehalf, x, y);
3828   efloor (y, y);
3829 }
3830
3831 /* Convert HOST_WIDE_INT LP to e type Y.  */
3832
3833 static void
3834 ltoe (lp, y)
3835      HOST_WIDE_INT *lp;
3836      unsigned EMUSHORT *y;
3837 {
3838   unsigned EMUSHORT yi[NI];
3839   unsigned HOST_WIDE_INT ll;
3840   int k;
3841
3842   ecleaz (yi);
3843   if (*lp < 0)
3844     {
3845       /* make it positive */
3846       ll = (unsigned HOST_WIDE_INT) (-(*lp));
3847       yi[0] = 0xffff;           /* put correct sign in the e type number */
3848     }
3849   else
3850     {
3851       ll = (unsigned HOST_WIDE_INT) (*lp);
3852     }
3853   /* move the long integer to yi significand area */
3854 #if HOST_BITS_PER_WIDE_INT == 64
3855   yi[M] = (unsigned EMUSHORT) (ll >> 48);
3856   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3857   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3858   yi[M + 3] = (unsigned EMUSHORT) ll;
3859   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
3860 #else
3861   yi[M] = (unsigned EMUSHORT) (ll >> 16);
3862   yi[M + 1] = (unsigned EMUSHORT) ll;
3863   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
3864 #endif
3865
3866   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3867     ecleaz (yi);                /* it was zero */
3868   else
3869     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
3870   emovo (yi, y);                /* output the answer */
3871 }
3872
3873 /* Convert unsigned HOST_WIDE_INT LP to e type Y.  */
3874
3875 static void
3876 ultoe (lp, y)
3877      unsigned HOST_WIDE_INT *lp;
3878      unsigned EMUSHORT *y;
3879 {
3880   unsigned EMUSHORT yi[NI];
3881   unsigned HOST_WIDE_INT ll;
3882   int k;
3883
3884   ecleaz (yi);
3885   ll = *lp;
3886
3887   /* move the long integer to ayi significand area */
3888 #if HOST_BITS_PER_WIDE_INT == 64
3889   yi[M] = (unsigned EMUSHORT) (ll >> 48);
3890   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3891   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3892   yi[M + 3] = (unsigned EMUSHORT) ll;
3893   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
3894 #else
3895   yi[M] = (unsigned EMUSHORT) (ll >> 16);
3896   yi[M + 1] = (unsigned EMUSHORT) ll;
3897   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
3898 #endif
3899
3900   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3901     ecleaz (yi);                /* it was zero */
3902   else
3903     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
3904   emovo (yi, y);                /* output the answer */
3905 }
3906
3907
3908 /* Find signed HOST_WIDE_INT integer I and floating point fractional
3909    part FRAC of e-type (packed internal format) floating point input X.
3910    The integer output I has the sign of the input, except that
3911    positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
3912    The output e-type fraction FRAC is the positive fractional
3913    part of abs (X).  */
3914
3915 static void
3916 eifrac (x, i, frac)
3917      unsigned EMUSHORT *x;
3918      HOST_WIDE_INT *i;
3919      unsigned EMUSHORT *frac;
3920 {
3921   unsigned EMUSHORT xi[NI];
3922   int j, k;
3923   unsigned HOST_WIDE_INT ll;
3924
3925   emovi (x, xi);
3926   k = (int) xi[E] - (EXONE - 1);
3927   if (k <= 0)
3928     {
3929       /* if exponent <= 0, integer = 0 and real output is fraction */
3930       *i = 0L;
3931       emovo (xi, frac);
3932       return;
3933     }
3934   if (k > (HOST_BITS_PER_WIDE_INT - 1))
3935     {
3936       /* long integer overflow: output large integer
3937          and correct fraction  */
3938       if (xi[0])
3939         *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
3940       else
3941         {
3942 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
3943           /* In this case, let it overflow and convert as if unsigned.  */
3944           euifrac (x, &ll, frac);
3945           *i = (HOST_WIDE_INT) ll;
3946           return;
3947 #else
3948           /* In other cases, return the largest positive integer.  */
3949           *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
3950 #endif
3951         }
3952       eshift (xi, k);
3953       if (extra_warnings)
3954         warning ("overflow on truncation to integer");
3955     }
3956   else if (k > 16)
3957     {
3958       /* Shift more than 16 bits: first shift up k-16 mod 16,
3959          then shift up by 16's.  */
3960       j = k - ((k >> 4) << 4);
3961       eshift (xi, j);
3962       ll = xi[M];
3963       k -= j;
3964       do
3965         {
3966           eshup6 (xi);
3967           ll = (ll << 16) | xi[M];
3968         }
3969       while ((k -= 16) > 0);
3970       *i = ll;
3971       if (xi[0])
3972         *i = -(*i);
3973     }
3974   else
3975       {
3976         /* shift not more than 16 bits */
3977           eshift (xi, k);
3978         *i = (HOST_WIDE_INT) xi[M] & 0xffff;
3979         if (xi[0])
3980           *i = -(*i);
3981       }
3982   xi[0] = 0;
3983   xi[E] = EXONE - 1;
3984   xi[M] = 0;
3985   if ((k = enormlz (xi)) > NBITS)
3986     ecleaz (xi);
3987   else
3988     xi[E] -= (unsigned EMUSHORT) k;
3989
3990   emovo (xi, frac);
3991 }
3992
3993
3994 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
3995    FRAC of e-type X.  A negative input yields integer output = 0 but
3996    correct fraction.  */
3997
3998 static void
3999 euifrac (x, i, frac)
4000      unsigned EMUSHORT *x;
4001      unsigned HOST_WIDE_INT *i;
4002      unsigned EMUSHORT *frac;
4003 {
4004   unsigned HOST_WIDE_INT ll;
4005   unsigned EMUSHORT xi[NI];
4006   int j, k;
4007
4008   emovi (x, xi);
4009   k = (int) xi[E] - (EXONE - 1);
4010   if (k <= 0)
4011     {
4012       /* if exponent <= 0, integer = 0 and argument is fraction */
4013       *i = 0L;
4014       emovo (xi, frac);
4015       return;
4016     }
4017   if (k > HOST_BITS_PER_WIDE_INT)
4018     {
4019       /* Long integer overflow: output large integer
4020          and correct fraction.
4021          Note, the BSD microvax compiler says that ~(0UL)
4022          is a syntax error.  */
4023       *i = ~(0L);
4024       eshift (xi, k);
4025       if (extra_warnings)
4026         warning ("overflow on truncation to unsigned integer");
4027     }
4028   else if (k > 16)
4029     {
4030       /* Shift more than 16 bits: first shift up k-16 mod 16,
4031          then shift up by 16's.  */
4032       j = k - ((k >> 4) << 4);
4033       eshift (xi, j);
4034       ll = xi[M];
4035       k -= j;
4036       do
4037         {
4038           eshup6 (xi);
4039           ll = (ll << 16) | xi[M];
4040         }
4041       while ((k -= 16) > 0);
4042       *i = ll;
4043     }
4044   else
4045     {
4046       /* shift not more than 16 bits */
4047       eshift (xi, k);
4048       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4049     }
4050
4051   if (xi[0])  /* A negative value yields unsigned integer 0. */
4052     *i = 0L;
4053
4054   xi[0] = 0;
4055   xi[E] = EXONE - 1;
4056   xi[M] = 0;
4057   if ((k = enormlz (xi)) > NBITS)
4058     ecleaz (xi);
4059   else
4060     xi[E] -= (unsigned EMUSHORT) k;
4061
4062   emovo (xi, frac);
4063 }
4064
4065 /* Shift the significand of exploded e-type X up or down by SC bits.  */
4066
4067 static int
4068 eshift (x, sc)
4069      unsigned EMUSHORT *x;
4070      int sc;
4071 {
4072   unsigned EMUSHORT lost;
4073   unsigned EMUSHORT *p;
4074
4075   if (sc == 0)
4076     return (0);
4077
4078   lost = 0;
4079   p = x + NI - 1;
4080
4081   if (sc < 0)
4082     {
4083       sc = -sc;
4084       while (sc >= 16)
4085         {
4086           lost |= *p;           /* remember lost bits */
4087           eshdn6 (x);
4088           sc -= 16;
4089         }
4090
4091       while (sc >= 8)
4092         {
4093           lost |= *p & 0xff;
4094           eshdn8 (x);
4095           sc -= 8;
4096         }
4097
4098       while (sc > 0)
4099         {
4100           lost |= *p & 1;
4101           eshdn1 (x);
4102           sc -= 1;
4103         }
4104     }
4105   else
4106     {
4107       while (sc >= 16)
4108         {
4109           eshup6 (x);
4110           sc -= 16;
4111         }
4112
4113       while (sc >= 8)
4114         {
4115           eshup8 (x);
4116           sc -= 8;
4117         }
4118
4119       while (sc > 0)
4120         {
4121           eshup1 (x);
4122           sc -= 1;
4123         }
4124     }
4125   if (lost)
4126     lost = 1;
4127   return ((int) lost);
4128 }
4129
4130 /* Shift normalize the significand area of exploded e-type X.
4131    Return the shift count (up = positive).  */
4132
4133 static int
4134 enormlz (x)
4135      unsigned EMUSHORT x[];
4136 {
4137   register unsigned EMUSHORT *p;
4138   int sc;
4139
4140   sc = 0;
4141   p = &x[M];
4142   if (*p != 0)
4143     goto normdn;
4144   ++p;
4145   if (*p & 0x8000)
4146     return (0);                 /* already normalized */
4147   while (*p == 0)
4148     {
4149       eshup6 (x);
4150       sc += 16;
4151
4152       /* With guard word, there are NBITS+16 bits available.
4153        Return true if all are zero.  */
4154       if (sc > NBITS)
4155         return (sc);
4156     }
4157   /* see if high byte is zero */
4158   while ((*p & 0xff00) == 0)
4159     {
4160       eshup8 (x);
4161       sc += 8;
4162     }
4163   /* now shift 1 bit at a time */
4164   while ((*p & 0x8000) == 0)
4165     {
4166       eshup1 (x);
4167       sc += 1;
4168       if (sc > NBITS)
4169         {
4170           mtherr ("enormlz", UNDERFLOW);
4171           return (sc);
4172         }
4173     }
4174   return (sc);
4175
4176   /* Normalize by shifting down out of the high guard word
4177      of the significand */
4178  normdn:
4179
4180   if (*p & 0xff00)
4181     {
4182       eshdn8 (x);
4183       sc -= 8;
4184     }
4185   while (*p != 0)
4186     {
4187       eshdn1 (x);
4188       sc -= 1;
4189
4190       if (sc < -NBITS)
4191         {
4192           mtherr ("enormlz", OVERFLOW);
4193           return (sc);
4194         }
4195     }
4196   return (sc);
4197 }
4198
4199 /* Powers of ten used in decimal <-> binary conversions.  */
4200
4201 #define NTEN 12
4202 #define MAXP 4096
4203
4204 #if LONG_DOUBLE_TYPE_SIZE == 128
4205 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4206 {
4207   {0x6576, 0x4a92, 0x804a, 0x153f,
4208    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4209   {0x6a32, 0xce52, 0x329a, 0x28ce,
4210    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4211   {0x526c, 0x50ce, 0xf18b, 0x3d28,
4212    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4213   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4214    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4215   {0x851e, 0xeab7, 0x98fe, 0x901b,
4216    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4217   {0x0235, 0x0137, 0x36b1, 0x336c,
4218    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4219   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4220    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4221   {0x0000, 0x0000, 0x0000, 0x0000,
4222    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4223   {0x0000, 0x0000, 0x0000, 0x0000,
4224    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4225   {0x0000, 0x0000, 0x0000, 0x0000,
4226    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4227   {0x0000, 0x0000, 0x0000, 0x0000,
4228    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4229   {0x0000, 0x0000, 0x0000, 0x0000,
4230    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4231   {0x0000, 0x0000, 0x0000, 0x0000,
4232    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4233 };
4234
4235 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4236 {
4237   {0x2030, 0xcffc, 0xa1c3, 0x8123,
4238    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4239   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4240    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4241   {0xf53f, 0xf698, 0x6bd3, 0x0158,
4242    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4243   {0xe731, 0x04d4, 0xe3f2, 0xd332,
4244    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4245   {0xa23e, 0x5308, 0xfefb, 0x1155,
4246    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4247   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4248    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4249   {0x2a20, 0x6224, 0x47b3, 0x98d7,
4250    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4251   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4252    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4253   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4254    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4255   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4256    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4257   {0xc155, 0xa4a8, 0x404e, 0x6113,
4258    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4259   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4260    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4261   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4262    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4263 };
4264 #else
4265 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4266 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4267 {
4268   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4269   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4270   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4271   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4272   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4273   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4274   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4275   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4276   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4277   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4278   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4279   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4280   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4281 };
4282
4283 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4284 {
4285   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4286   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4287   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4288   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4289   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4290   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4291   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4292   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4293   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4294   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4295   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4296   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4297   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4298 };
4299 #endif
4300
4301 /* Convert float value X to ASCII string STRING with NDIG digits after
4302    the decimal point.  */
4303
4304 static void
4305 e24toasc (x, string, ndigs)
4306      unsigned EMUSHORT x[];
4307      char *string;
4308      int ndigs;
4309 {
4310   unsigned EMUSHORT w[NI];
4311
4312   e24toe (x, w);
4313   etoasc (w, string, ndigs);
4314 }
4315
4316 /* Convert double value X to ASCII string STRING with NDIG digits after
4317    the decimal point.  */
4318
4319 static void
4320 e53toasc (x, string, ndigs)
4321      unsigned EMUSHORT x[];
4322      char *string;
4323      int ndigs;
4324 {
4325   unsigned EMUSHORT w[NI];
4326
4327   e53toe (x, w);
4328   etoasc (w, string, ndigs);
4329 }
4330
4331 /* Convert double extended value X to ASCII string STRING with NDIG digits
4332    after the decimal point.  */
4333
4334 static void
4335 e64toasc (x, string, ndigs)
4336      unsigned EMUSHORT x[];
4337      char *string;
4338      int ndigs;
4339 {
4340   unsigned EMUSHORT w[NI];
4341
4342   e64toe (x, w);
4343   etoasc (w, string, ndigs);
4344 }
4345
4346 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4347    after the decimal point.  */
4348
4349 static void
4350 e113toasc (x, string, ndigs)
4351      unsigned EMUSHORT x[];
4352      char *string;
4353      int ndigs;
4354 {
4355   unsigned EMUSHORT w[NI];
4356
4357   e113toe (x, w);
4358   etoasc (w, string, ndigs);
4359 }
4360
4361 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4362    the decimal point.  */
4363
4364 static char wstring[80];        /* working storage for ASCII output */
4365
4366 static void
4367 etoasc (x, string, ndigs)
4368      unsigned EMUSHORT x[];
4369      char *string;
4370      int ndigs;
4371 {
4372   EMUSHORT digit;
4373   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
4374   unsigned EMUSHORT *p, *r, *ten;
4375   unsigned EMUSHORT sign;
4376   int i, j, k, expon, rndsav;
4377   char *s, *ss;
4378   unsigned EMUSHORT m;
4379
4380
4381   rndsav = rndprc;
4382   ss = string;
4383   s = wstring;
4384   *ss = '\0';
4385   *s = '\0';
4386 #ifdef NANS
4387   if (eisnan (x))
4388     {
4389       sprintf (wstring, " NaN ");
4390       goto bxit;
4391     }
4392 #endif
4393   rndprc = NBITS;               /* set to full precision */
4394   emov (x, y);                  /* retain external format */
4395   if (y[NE - 1] & 0x8000)
4396     {
4397       sign = 0xffff;
4398       y[NE - 1] &= 0x7fff;
4399     }
4400   else
4401     {
4402       sign = 0;
4403     }
4404   expon = 0;
4405   ten = &etens[NTEN][0];
4406   emov (eone, t);
4407   /* Test for zero exponent */
4408   if (y[NE - 1] == 0)
4409     {
4410       for (k = 0; k < NE - 1; k++)
4411         {
4412           if (y[k] != 0)
4413             goto tnzro;         /* denormalized number */
4414         }
4415       goto isone;               /* valid all zeros */
4416     }
4417  tnzro:
4418
4419   /* Test for infinity. */
4420   if (y[NE - 1] == 0x7fff)
4421     {
4422       if (sign)
4423         sprintf (wstring, " -Infinity ");
4424       else
4425         sprintf (wstring, " Infinity ");
4426       goto bxit;
4427     }
4428
4429   /* Test for exponent nonzero but significand denormalized.
4430    * This is an error condition.
4431    */
4432   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4433     {
4434       mtherr ("etoasc", DOMAIN);
4435       sprintf (wstring, "NaN");
4436       goto bxit;
4437     }
4438
4439   /* Compare to 1.0 */
4440   i = ecmp (eone, y);
4441   if (i == 0)
4442     goto isone;
4443
4444   if (i == -2)
4445     abort ();
4446
4447   if (i < 0)
4448     {                           /* Number is greater than 1 */
4449       /* Convert significand to an integer and strip trailing decimal zeros. */
4450       emov (y, u);
4451       u[NE - 1] = EXONE + NBITS - 1;
4452
4453       p = &etens[NTEN - 4][0];
4454       m = 16;
4455       do
4456         {
4457           ediv (p, u, t);
4458           efloor (t, w);
4459           for (j = 0; j < NE - 1; j++)
4460             {
4461               if (t[j] != w[j])
4462                 goto noint;
4463             }
4464           emov (t, u);
4465           expon += (int) m;
4466         noint:
4467           p += NE;
4468           m >>= 1;
4469         }
4470       while (m != 0);
4471
4472       /* Rescale from integer significand */
4473       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4474       emov (u, y);
4475       /* Find power of 10 */
4476       emov (eone, t);
4477       m = MAXP;
4478       p = &etens[0][0];
4479       /* An unordered compare result shouldn't happen here. */
4480       while (ecmp (ten, u) <= 0)
4481         {
4482           if (ecmp (p, u) <= 0)
4483             {
4484               ediv (p, u, u);
4485               emul (p, t, t);
4486               expon += (int) m;
4487             }
4488           m >>= 1;
4489           if (m == 0)
4490             break;
4491           p += NE;
4492         }
4493     }
4494   else
4495     {                           /* Number is less than 1.0 */
4496       /* Pad significand with trailing decimal zeros. */
4497       if (y[NE - 1] == 0)
4498         {
4499           while ((y[NE - 2] & 0x8000) == 0)
4500             {
4501               emul (ten, y, y);
4502               expon -= 1;
4503             }
4504         }
4505       else
4506         {
4507           emovi (y, w);
4508           for (i = 0; i < NDEC + 1; i++)
4509             {
4510               if ((w[NI - 1] & 0x7) != 0)
4511                 break;
4512               /* multiply by 10 */
4513               emovz (w, u);
4514               eshdn1 (u);
4515               eshdn1 (u);
4516               eaddm (w, u);
4517               u[1] += 3;
4518               while (u[2] != 0)
4519                 {
4520                   eshdn1 (u);
4521                   u[1] += 1;
4522                 }
4523               if (u[NI - 1] != 0)
4524                 break;
4525               if (eone[NE - 1] <= u[1])
4526                 break;
4527               emovz (u, w);
4528               expon -= 1;
4529             }
4530           emovo (w, y);
4531         }
4532       k = -MAXP;
4533       p = &emtens[0][0];
4534       r = &etens[0][0];
4535       emov (y, w);
4536       emov (eone, t);
4537       while (ecmp (eone, w) > 0)
4538         {
4539           if (ecmp (p, w) >= 0)
4540             {
4541               emul (r, w, w);
4542               emul (r, t, t);
4543               expon += k;
4544             }
4545           k /= 2;
4546           if (k == 0)
4547             break;
4548           p += NE;
4549           r += NE;
4550         }
4551       ediv (t, eone, t);
4552     }
4553  isone:
4554   /* Find the first (leading) digit. */
4555   emovi (t, w);
4556   emovz (w, t);
4557   emovi (y, w);
4558   emovz (w, y);
4559   eiremain (t, y);
4560   digit = equot[NI - 1];
4561   while ((digit == 0) && (ecmp (y, ezero) != 0))
4562     {
4563       eshup1 (y);
4564       emovz (y, u);
4565       eshup1 (u);
4566       eshup1 (u);
4567       eaddm (u, y);
4568       eiremain (t, y);
4569       digit = equot[NI - 1];
4570       expon -= 1;
4571     }
4572   s = wstring;
4573   if (sign)
4574     *s++ = '-';
4575   else
4576     *s++ = ' ';
4577   /* Examine number of digits requested by caller. */
4578   if (ndigs < 0)
4579     ndigs = 0;
4580   if (ndigs > NDEC)
4581     ndigs = NDEC;
4582   if (digit == 10)
4583     {
4584       *s++ = '1';
4585       *s++ = '.';
4586       if (ndigs > 0)
4587         {
4588           *s++ = '0';
4589           ndigs -= 1;
4590         }
4591       expon += 1;
4592     }
4593   else
4594     {
4595       *s++ = (char)digit + '0';
4596       *s++ = '.';
4597     }
4598   /* Generate digits after the decimal point. */
4599   for (k = 0; k <= ndigs; k++)
4600     {
4601       /* multiply current number by 10, without normalizing */
4602       eshup1 (y);
4603       emovz (y, u);
4604       eshup1 (u);
4605       eshup1 (u);
4606       eaddm (u, y);
4607       eiremain (t, y);
4608       *s++ = (char) equot[NI - 1] + '0';
4609     }
4610   digit = equot[NI - 1];
4611   --s;
4612   ss = s;
4613   /* round off the ASCII string */
4614   if (digit > 4)
4615     {
4616       /* Test for critical rounding case in ASCII output. */
4617       if (digit == 5)
4618         {
4619           emovo (y, t);
4620           if (ecmp (t, ezero) != 0)
4621             goto roun;          /* round to nearest */
4622           if ((*(s - 1) & 1) == 0)
4623             goto doexp;         /* round to even */
4624         }
4625       /* Round up and propagate carry-outs */
4626     roun:
4627       --s;
4628       k = *s & 0x7f;
4629       /* Carry out to most significant digit? */
4630       if (k == '.')
4631         {
4632           --s;
4633           k = *s;
4634           k += 1;
4635           *s = (char) k;
4636           /* Most significant digit carries to 10? */
4637           if (k > '9')
4638             {
4639               expon += 1;
4640               *s = '1';
4641             }
4642           goto doexp;
4643         }
4644       /* Round up and carry out from less significant digits */
4645       k += 1;
4646       *s = (char) k;
4647       if (k > '9')
4648         {
4649           *s = '0';
4650           goto roun;
4651         }
4652     }
4653  doexp:
4654   /*
4655      if (expon >= 0)
4656      sprintf (ss, "e+%d", expon);
4657      else
4658      sprintf (ss, "e%d", expon);
4659      */
4660   sprintf (ss, "e%d", expon);
4661  bxit:
4662   rndprc = rndsav;
4663   /* copy out the working string */
4664   s = string;
4665   ss = wstring;
4666   while (*ss == ' ')            /* strip possible leading space */
4667     ++ss;
4668   while ((*s++ = *ss++) != '\0')
4669     ;
4670 }
4671
4672
4673 /* Convert ASCII string to floating point.
4674
4675    Numeric input is a free format decimal number of any length, with
4676    or without decimal point.  Entering E after the number followed by an
4677    integer number causes the second number to be interpreted as a power of
4678    10 to be multiplied by the first number (i.e., "scientific" notation).  */
4679
4680 /* Convert ASCII string S to single precision float value Y.  */
4681
4682 static void
4683 asctoe24 (s, y)
4684      char *s;
4685      unsigned EMUSHORT *y;
4686 {
4687   asctoeg (s, y, 24);
4688 }
4689
4690
4691 /* Convert ASCII string S to double precision value Y.  */
4692
4693 static void
4694 asctoe53 (s, y)
4695      char *s;
4696      unsigned EMUSHORT *y;
4697 {
4698 #if defined(DEC) || defined(IBM)
4699   asctoeg (s, y, 56);
4700 #else
4701   asctoeg (s, y, 53);
4702 #endif
4703 }
4704
4705
4706 /* Convert ASCII string S to double extended value Y.  */
4707
4708 static void
4709 asctoe64 (s, y)
4710      char *s;
4711      unsigned EMUSHORT *y;
4712 {
4713   asctoeg (s, y, 64);
4714 }
4715
4716 /* Convert ASCII string S to 128-bit long double Y.  */
4717
4718 static void
4719 asctoe113 (s, y)
4720      char *s;
4721      unsigned EMUSHORT *y;
4722 {
4723   asctoeg (s, y, 113);
4724 }
4725
4726 /* Convert ASCII string S to e type Y.  */
4727
4728 static void
4729 asctoe (s, y)
4730      char *s;
4731      unsigned EMUSHORT *y;
4732 {
4733   asctoeg (s, y, NBITS);
4734 }
4735
4736 /* Convert ASCII string SS to e type Y, with a specified rounding precision
4737    of OPREC bits. */
4738
4739 static void
4740 asctoeg (ss, y, oprec)
4741      char *ss;
4742      unsigned EMUSHORT *y;
4743      int oprec;
4744 {
4745   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
4746   int esign, decflg, sgnflg, nexp, exp, prec, lost;
4747   int k, trail, c, rndsav;
4748   EMULONG lexp;
4749   unsigned EMUSHORT nsign, *p;
4750   char *sp, *s, *lstr;
4751
4752   /* Copy the input string. */
4753   lstr = (char *) alloca (strlen (ss) + 1);
4754   s = ss;
4755   while (*s == ' ')             /* skip leading spaces */
4756     ++s;
4757   sp = lstr;
4758   while ((*sp++ = *s++) != '\0')
4759     ;
4760   s = lstr;
4761
4762   rndsav = rndprc;
4763   rndprc = NBITS;               /* Set to full precision */
4764   lost = 0;
4765   nsign = 0;
4766   decflg = 0;
4767   sgnflg = 0;
4768   nexp = 0;
4769   exp = 0;
4770   prec = 0;
4771   ecleaz (yy);
4772   trail = 0;
4773
4774  nxtcom:
4775   k = *s - '0';
4776   if ((k >= 0) && (k <= 9))
4777     {
4778       /* Ignore leading zeros */
4779       if ((prec == 0) && (decflg == 0) && (k == 0))
4780         goto donchr;
4781       /* Identify and strip trailing zeros after the decimal point. */
4782       if ((trail == 0) && (decflg != 0))
4783         {
4784           sp = s;
4785           while ((*sp >= '0') && (*sp <= '9'))
4786             ++sp;
4787           /* Check for syntax error */
4788           c = *sp & 0x7f;
4789           if ((c != 'e') && (c != 'E') && (c != '\0')
4790               && (c != '\n') && (c != '\r') && (c != ' ')
4791               && (c != ','))
4792             goto error;
4793           --sp;
4794           while (*sp == '0')
4795             *sp-- = 'z';
4796           trail = 1;
4797           if (*s == 'z')
4798             goto donchr;
4799         }
4800
4801       /* If enough digits were given to more than fill up the yy register,
4802          continuing until overflow into the high guard word yy[2]
4803          guarantees that there will be a roundoff bit at the top
4804          of the low guard word after normalization.  */
4805
4806       if (yy[2] == 0)
4807         {
4808           if (decflg)
4809             nexp += 1;          /* count digits after decimal point */
4810           eshup1 (yy);          /* multiply current number by 10 */
4811           emovz (yy, xt);
4812           eshup1 (xt);
4813           eshup1 (xt);
4814           eaddm (xt, yy);
4815           ecleaz (xt);
4816           xt[NI - 2] = (unsigned EMUSHORT) k;
4817           eaddm (xt, yy);
4818         }
4819       else
4820         {
4821           /* Mark any lost non-zero digit.  */
4822           lost |= k;
4823           /* Count lost digits before the decimal point.  */
4824           if (decflg == 0)
4825             nexp -= 1;
4826         }
4827       prec += 1;
4828       goto donchr;
4829     }
4830
4831   switch (*s)
4832     {
4833     case 'z':
4834       break;
4835     case 'E':
4836     case 'e':
4837       goto expnt;
4838     case '.':                   /* decimal point */
4839       if (decflg)
4840         goto error;
4841       ++decflg;
4842       break;
4843     case '-':
4844       nsign = 0xffff;
4845       if (sgnflg)
4846         goto error;
4847       ++sgnflg;
4848       break;
4849     case '+':
4850       if (sgnflg)
4851         goto error;
4852       ++sgnflg;
4853       break;
4854     case ',':
4855     case ' ':
4856     case '\0':
4857     case '\n':
4858     case '\r':
4859       goto daldone;
4860     case 'i':
4861     case 'I':
4862       goto infinite;
4863     default:
4864     error:
4865 #ifdef NANS
4866       einan (yy);
4867 #else
4868       mtherr ("asctoe", DOMAIN);
4869       eclear (yy);
4870 #endif
4871       goto aexit;
4872     }
4873  donchr:
4874   ++s;
4875   goto nxtcom;
4876
4877   /* Exponent interpretation */
4878  expnt:
4879
4880   esign = 1;
4881   exp = 0;
4882   ++s;
4883   /* check for + or - */
4884   if (*s == '-')
4885     {
4886       esign = -1;
4887       ++s;
4888     }
4889   if (*s == '+')
4890     ++s;
4891   while ((*s >= '0') && (*s <= '9'))
4892     {
4893       exp *= 10;
4894       exp += *s++ - '0';
4895       if (exp > -(MINDECEXP))
4896         {
4897           if (esign < 0)
4898             goto zero;
4899           else
4900             goto infinite;
4901         }
4902     }
4903   if (esign < 0)
4904     exp = -exp;
4905   if (exp > MAXDECEXP)
4906     {
4907  infinite:
4908       ecleaz (yy);
4909       yy[E] = 0x7fff;           /* infinity */
4910       goto aexit;
4911     }
4912   if (exp < MINDECEXP)
4913     {
4914  zero:
4915       ecleaz (yy);
4916       goto aexit;
4917     }
4918
4919  daldone:
4920   nexp = exp - nexp;
4921   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
4922   while ((nexp > 0) && (yy[2] == 0))
4923     {
4924       emovz (yy, xt);
4925       eshup1 (xt);
4926       eshup1 (xt);
4927       eaddm (yy, xt);
4928       eshup1 (xt);
4929       if (xt[2] != 0)
4930         break;
4931       nexp -= 1;
4932       emovz (xt, yy);
4933     }
4934   if ((k = enormlz (yy)) > NBITS)
4935     {
4936       ecleaz (yy);
4937       goto aexit;
4938     }
4939   lexp = (EXONE - 1 + NBITS) - k;
4940   emdnorm (yy, lost, 0, lexp, 64);
4941
4942   /* Convert to external format:
4943
4944      Multiply by 10**nexp.  If precision is 64 bits,
4945      the maximum relative error incurred in forming 10**n
4946      for 0 <= n <= 324 is 8.2e-20, at 10**180.
4947      For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
4948      For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */
4949
4950   lexp = yy[E];
4951   if (nexp == 0)
4952     {
4953       k = 0;
4954       goto expdon;
4955     }
4956   esign = 1;
4957   if (nexp < 0)
4958     {
4959       nexp = -nexp;
4960       esign = -1;
4961       if (nexp > 4096)
4962         {
4963           /* Punt.  Can't handle this without 2 divides. */
4964           emovi (etens[0], tt);
4965           lexp -= tt[E];
4966           k = edivm (tt, yy);
4967           lexp += EXONE;
4968           nexp -= 4096;
4969         }
4970     }
4971   p = &etens[NTEN][0];
4972   emov (eone, xt);
4973   exp = 1;
4974   do
4975     {
4976       if (exp & nexp)
4977         emul (p, xt, xt);
4978       p -= NE;
4979       exp = exp + exp;
4980     }
4981   while (exp <= MAXP);
4982
4983   emovi (xt, tt);
4984   if (esign < 0)
4985     {
4986       lexp -= tt[E];
4987       k = edivm (tt, yy);
4988       lexp += EXONE;
4989     }
4990   else
4991     {
4992       lexp += tt[E];
4993       k = emulm (tt, yy);
4994       lexp -= EXONE - 1;
4995     }
4996
4997  expdon:
4998
4999   /* Round and convert directly to the destination type */
5000   if (oprec == 53)
5001     lexp -= EXONE - 0x3ff;
5002 #ifdef IBM
5003   else if (oprec == 24 || oprec == 56)
5004     lexp -= EXONE - (0x41 << 2);
5005 #else
5006   else if (oprec == 24)
5007     lexp -= EXONE - 0177;
5008 #endif
5009 #ifdef DEC
5010   else if (oprec == 56)
5011     lexp -= EXONE - 0201;
5012 #endif
5013   rndprc = oprec;
5014   emdnorm (yy, k, 0, lexp, 64);
5015
5016  aexit:
5017
5018   rndprc = rndsav;
5019   yy[0] = nsign;
5020   switch (oprec)
5021     {
5022 #ifdef DEC
5023     case 56:
5024       todec (yy, y);            /* see etodec.c */
5025       break;
5026 #endif
5027 #ifdef IBM
5028     case 56:
5029       toibm (yy, y, DFmode);
5030       break;
5031 #endif
5032     case 53:
5033       toe53 (yy, y);
5034       break;
5035     case 24:
5036       toe24 (yy, y);
5037       break;
5038     case 64:
5039       toe64 (yy, y);
5040       break;
5041     case 113:
5042       toe113 (yy, y);
5043       break;
5044     case NBITS:
5045       emovo (yy, y);
5046       break;
5047     }
5048 }
5049
5050
5051
5052 /* Return Y = largest integer not greater than X (truncated toward minus
5053    infinity).  */
5054
5055 static unsigned EMUSHORT bmask[] =
5056 {
5057   0xffff,
5058   0xfffe,
5059   0xfffc,
5060   0xfff8,
5061   0xfff0,
5062   0xffe0,
5063   0xffc0,
5064   0xff80,
5065   0xff00,
5066   0xfe00,
5067   0xfc00,
5068   0xf800,
5069   0xf000,
5070   0xe000,
5071   0xc000,
5072   0x8000,
5073   0x0000,
5074 };
5075
5076 static void
5077 efloor (x, y)
5078      unsigned EMUSHORT x[], y[];
5079 {
5080   register unsigned EMUSHORT *p;
5081   int e, expon, i;
5082   unsigned EMUSHORT f[NE];
5083
5084   emov (x, f);                  /* leave in external format */
5085   expon = (int) f[NE - 1];
5086   e = (expon & 0x7fff) - (EXONE - 1);
5087   if (e <= 0)
5088     {
5089       eclear (y);
5090       goto isitneg;
5091     }
5092   /* number of bits to clear out */
5093   e = NBITS - e;
5094   emov (f, y);
5095   if (e <= 0)
5096     return;
5097
5098   p = &y[0];
5099   while (e >= 16)
5100     {
5101       *p++ = 0;
5102       e -= 16;
5103     }
5104   /* clear the remaining bits */
5105   *p &= bmask[e];
5106   /* truncate negatives toward minus infinity */
5107  isitneg:
5108
5109   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
5110     {
5111       for (i = 0; i < NE - 1; i++)
5112         {
5113           if (f[i] != y[i])
5114             {
5115               esub (eone, y, y);
5116               break;
5117             }
5118         }
5119     }
5120 }
5121
5122
5123 /* Return S and EXP such that  S * 2^EXP = X and .5 <= S < 1.
5124    For example, 1.1 = 0.55 * 2^1.  */
5125
5126 static void
5127 efrexp (x, exp, s)
5128      unsigned EMUSHORT x[];
5129      int *exp;
5130      unsigned EMUSHORT s[];
5131 {
5132   unsigned EMUSHORT xi[NI];
5133   EMULONG li;
5134
5135   emovi (x, xi);
5136   /*  Handle denormalized numbers properly using long integer exponent.  */
5137   li = (EMULONG) ((EMUSHORT) xi[1]);
5138
5139   if (li == 0)
5140     {
5141       li -= enormlz (xi);
5142     }
5143   xi[1] = 0x3ffe;
5144   emovo (xi, s);
5145   *exp = (int) (li - 0x3ffe);
5146 }
5147
5148 /* Return e type Y = X * 2^PWR2.  */
5149
5150 static void
5151 eldexp (x, pwr2, y)
5152      unsigned EMUSHORT x[];
5153      int pwr2;
5154      unsigned EMUSHORT y[];
5155 {
5156   unsigned EMUSHORT xi[NI];
5157   EMULONG li;
5158   int i;
5159
5160   emovi (x, xi);
5161   li = xi[1];
5162   li += pwr2;
5163   i = 0;
5164   emdnorm (xi, i, i, li, 64);
5165   emovo (xi, y);
5166 }
5167
5168
5169 /* C = remainder after dividing B by A, all e type values.
5170    Least significant integer quotient bits left in EQUOT.  */
5171
5172 static void
5173 eremain (a, b, c)
5174      unsigned EMUSHORT a[], b[], c[];
5175 {
5176   unsigned EMUSHORT den[NI], num[NI];
5177
5178 #ifdef NANS
5179   if (eisinf (b)
5180       || (ecmp (a, ezero) == 0)
5181       || eisnan (a)
5182       || eisnan (b))
5183     {
5184       enan (c, 0);
5185       return;
5186     }
5187 #endif
5188   if (ecmp (a, ezero) == 0)
5189     {
5190       mtherr ("eremain", SING);
5191       eclear (c);
5192       return;
5193     }
5194   emovi (a, den);
5195   emovi (b, num);
5196   eiremain (den, num);
5197   /* Sign of remainder = sign of quotient */
5198   if (a[0] == b[0])
5199     num[0] = 0;
5200   else
5201     num[0] = 0xffff;
5202   emovo (num, c);
5203 }
5204
5205 /*  Return quotient of exploded e-types NUM / DEN in EQUOT,
5206     remainder in NUM.  */
5207
5208 static void
5209 eiremain (den, num)
5210      unsigned EMUSHORT den[], num[];
5211 {
5212   EMULONG ld, ln;
5213   unsigned EMUSHORT j;
5214
5215   ld = den[E];
5216   ld -= enormlz (den);
5217   ln = num[E];
5218   ln -= enormlz (num);
5219   ecleaz (equot);
5220   while (ln >= ld)
5221     {
5222       if (ecmpm (den, num) <= 0)
5223         {
5224           esubm (den, num);
5225           j = 1;
5226         }
5227       else
5228           j = 0;
5229       eshup1 (equot);
5230       equot[NI - 1] |= j;
5231       eshup1 (num);
5232       ln -= 1;
5233     }
5234   emdnorm (num, 0, 0, ln, 0);
5235 }
5236
5237 /* Report an error condition CODE encountered in function NAME.
5238    CODE is one of the following:
5239
5240     Mnemonic        Value          Significance
5241
5242      DOMAIN            1       argument domain error
5243      SING              2       function singularity
5244      OVERFLOW          3       overflow range error
5245      UNDERFLOW         4       underflow range error
5246      TLOSS             5       total loss of precision
5247      PLOSS             6       partial loss of precision
5248      INVALID           7       NaN - producing operation
5249      EDOM             33       Unix domain error code
5250      ERANGE           34       Unix range error code
5251
5252    The order of appearance of the following messages is bound to the
5253    error codes defined above.  */
5254
5255 #define NMSGS 8
5256 static char *ermsg[NMSGS] =
5257 {
5258   "unknown",                    /* error code 0 */
5259   "domain",                     /* error code 1 */
5260   "singularity",                /* et seq.      */
5261   "overflow",
5262   "underflow",
5263   "total loss of precision",
5264   "partial loss of precision",
5265   "invalid operation"
5266 };
5267
5268 int merror = 0;
5269 extern int merror;
5270
5271 static void
5272 mtherr (name, code)
5273      char *name;
5274      int code;
5275 {
5276   char errstr[80];
5277
5278   /* The string passed by the calling program is supposed to be the
5279      name of the function in which the error occurred.
5280      The code argument selects which error message string will be printed.  */
5281
5282   if ((code <= 0) || (code >= NMSGS))
5283     code = 0;
5284   sprintf (errstr, " %s %s error", name, ermsg[code]);
5285   if (extra_warnings)
5286     warning (errstr);
5287   /* Set global error message word */
5288   merror = code + 1;
5289 }
5290
5291 #ifdef DEC
5292 /* Convert DEC double precision D to e type E.  */
5293
5294 static void
5295 dectoe (d, e)
5296      unsigned EMUSHORT *d;
5297      unsigned EMUSHORT *e;
5298 {
5299   unsigned EMUSHORT y[NI];
5300   register unsigned EMUSHORT r, *p;
5301
5302   ecleaz (y);                   /* start with a zero */
5303   p = y;                        /* point to our number */
5304   r = *d;                       /* get DEC exponent word */
5305   if (*d & (unsigned int) 0x8000)
5306     *p = 0xffff;                /* fill in our sign */
5307   ++p;                          /* bump pointer to our exponent word */
5308   r &= 0x7fff;                  /* strip the sign bit */
5309   if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
5310     goto done;
5311
5312
5313   r >>= 7;                      /* shift exponent word down 7 bits */
5314   r += EXONE - 0201;            /* subtract DEC exponent offset */
5315   /* add our e type exponent offset */
5316   *p++ = r;                     /* to form our exponent */
5317
5318   r = *d++;                     /* now do the high order mantissa */
5319   r &= 0177;                    /* strip off the DEC exponent and sign bits */
5320   r |= 0200;                    /* the DEC understood high order mantissa bit */
5321   *p++ = r;                     /* put result in our high guard word */
5322
5323   *p++ = *d++;                  /* fill in the rest of our mantissa */
5324   *p++ = *d++;
5325   *p = *d;
5326
5327   eshdn8 (y);                   /* shift our mantissa down 8 bits */
5328  done:
5329   emovo (y, e);
5330 }
5331
5332 /* Convert e type X to DEC double precision D.  */
5333
5334 static void
5335 etodec (x, d)
5336      unsigned EMUSHORT *x, *d;
5337 {
5338   unsigned EMUSHORT xi[NI];
5339   EMULONG exp;
5340   int rndsav;
5341
5342   emovi (x, xi);
5343   /* Adjust exponent for offsets.  */
5344   exp = (EMULONG) xi[E] - (EXONE - 0201);
5345   /* Round off to nearest or even.  */
5346   rndsav = rndprc;
5347   rndprc = 56;
5348   emdnorm (xi, 0, 0, exp, 64);
5349   rndprc = rndsav;
5350   todec (xi, d);
5351 }
5352
5353 /* Convert exploded e-type X, that has already been rounded to
5354    56-bit precision, to DEC format double Y.  */
5355
5356 static void
5357 todec (x, y)
5358      unsigned EMUSHORT *x, *y;
5359 {
5360   unsigned EMUSHORT i;
5361   unsigned EMUSHORT *p;
5362
5363   p = x;
5364   *y = 0;
5365   if (*p++)
5366     *y = 0100000;
5367   i = *p++;
5368   if (i == 0)
5369     {
5370       *y++ = 0;
5371       *y++ = 0;
5372       *y++ = 0;
5373       *y++ = 0;
5374       return;
5375     }
5376   if (i > 0377)
5377     {
5378       *y++ |= 077777;
5379       *y++ = 0xffff;
5380       *y++ = 0xffff;
5381       *y++ = 0xffff;
5382 #ifdef ERANGE
5383       errno = ERANGE;
5384 #endif
5385       return;
5386     }
5387   i &= 0377;
5388   i <<= 7;
5389   eshup8 (x);
5390   x[M] &= 0177;
5391   i |= x[M];
5392   *y++ |= i;
5393   *y++ = x[M + 1];
5394   *y++ = x[M + 2];
5395   *y++ = x[M + 3];
5396 }
5397 #endif /* DEC */
5398
5399 #ifdef IBM
5400 /* Convert IBM single/double precision to e type.  */
5401
5402 static void
5403 ibmtoe (d, e, mode)
5404      unsigned EMUSHORT *d;
5405      unsigned EMUSHORT *e;
5406      enum machine_mode mode;
5407 {
5408   unsigned EMUSHORT y[NI];
5409   register unsigned EMUSHORT r, *p;
5410   int rndsav;
5411
5412   ecleaz (y);                   /* start with a zero */
5413   p = y;                        /* point to our number */
5414   r = *d;                       /* get IBM exponent word */
5415   if (*d & (unsigned int) 0x8000)
5416     *p = 0xffff;                /* fill in our sign */
5417   ++p;                          /* bump pointer to our exponent word */
5418   r &= 0x7f00;                  /* strip the sign bit */
5419   r >>= 6;                      /* shift exponent word down 6 bits */
5420                                 /* in fact shift by 8 right and 2 left */
5421   r += EXONE - (0x41 << 2);     /* subtract IBM exponent offset */
5422                                 /* add our e type exponent offset */
5423   *p++ = r;                     /* to form our exponent */
5424
5425   *p++ = *d++ & 0xff;           /* now do the high order mantissa */
5426                                 /* strip off the IBM exponent and sign bits */
5427   if (mode != SFmode)           /* there are only 2 words in SFmode */
5428     {
5429       *p++ = *d++;              /* fill in the rest of our mantissa */
5430       *p++ = *d++;
5431     }
5432   *p = *d;
5433
5434   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5435     y[0] = y[E] = 0;
5436   else
5437     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
5438                               /* handle change in RADIX */
5439   emovo (y, e);
5440 }
5441
5442
5443
5444 /* Convert e type to IBM single/double precision.  */
5445
5446 static void
5447 etoibm (x, d, mode)
5448      unsigned EMUSHORT *x, *d;
5449      enum machine_mode mode;
5450 {
5451   unsigned EMUSHORT xi[NI];
5452   EMULONG exp;
5453   int rndsav;
5454
5455   emovi (x, xi);
5456   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));        /* adjust exponent for offsets */
5457                                                         /* round off to nearest or even */
5458   rndsav = rndprc;
5459   rndprc = 56;
5460   emdnorm (xi, 0, 0, exp, 64);
5461   rndprc = rndsav;
5462   toibm (xi, d, mode);
5463 }
5464
5465 static void
5466 toibm (x, y, mode)
5467      unsigned EMUSHORT *x, *y;
5468      enum machine_mode mode;
5469 {
5470   unsigned EMUSHORT i;
5471   unsigned EMUSHORT *p;
5472   int r;
5473
5474   p = x;
5475   *y = 0;
5476   if (*p++)
5477     *y = 0x8000;
5478   i = *p++;
5479   if (i == 0)
5480     {
5481       *y++ = 0;
5482       *y++ = 0;
5483       if (mode != SFmode)
5484         {
5485           *y++ = 0;
5486           *y++ = 0;
5487         }
5488       return;
5489     }
5490   r = i & 0x3;
5491   i >>= 2;
5492   if (i > 0x7f)
5493     {
5494       *y++ |= 0x7fff;
5495       *y++ = 0xffff;
5496       if (mode != SFmode)
5497         {
5498           *y++ = 0xffff;
5499           *y++ = 0xffff;
5500         }
5501 #ifdef ERANGE
5502       errno = ERANGE;
5503 #endif
5504       return;
5505     }
5506   i &= 0x7f;
5507   *y |= (i << 8);
5508   eshift (x, r + 5);
5509   *y++ |= x[M];
5510   *y++ = x[M + 1];
5511   if (mode != SFmode)
5512     {
5513       *y++ = x[M + 2];
5514       *y++ = x[M + 3];
5515     }
5516 }
5517 #endif /* IBM */
5518
5519 /* Output a binary NaN bit pattern in the target machine's format.  */
5520
5521 /* If special NaN bit patterns are required, define them in tm.h
5522    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
5523    patterns here. */
5524 #ifdef TFMODE_NAN
5525 TFMODE_NAN;
5526 #else
5527 #ifdef IEEE
5528 unsigned EMUSHORT TFbignan[8] =
5529  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5530 unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
5531 #endif
5532 #endif
5533
5534 #ifdef XFMODE_NAN
5535 XFMODE_NAN;
5536 #else
5537 #ifdef IEEE
5538 unsigned EMUSHORT XFbignan[6] =
5539  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5540 unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
5541 #endif
5542 #endif
5543
5544 #ifdef DFMODE_NAN
5545 DFMODE_NAN;
5546 #else
5547 #ifdef IEEE
5548 unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5549 unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
5550 #endif
5551 #endif
5552
5553 #ifdef SFMODE_NAN
5554 SFMODE_NAN;
5555 #else
5556 #ifdef IEEE
5557 unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
5558 unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
5559 #endif
5560 #endif
5561
5562
5563 static void
5564 make_nan (nan, sign, mode)
5565      unsigned EMUSHORT *nan;
5566      int sign;
5567      enum machine_mode mode;
5568 {
5569   int n;
5570   unsigned EMUSHORT *p;
5571
5572   switch (mode)
5573     {
5574 /* Possibly the `reserved operand' patterns on a VAX can be
5575    used like NaN's, but probably not in the same way as IEEE. */
5576 #if !defined(DEC) && !defined(IBM)
5577     case TFmode:
5578       n = 8;
5579       if (REAL_WORDS_BIG_ENDIAN)
5580         p = TFbignan;
5581       else
5582         p = TFlittlenan;
5583       break;
5584     case XFmode:
5585       n = 6;
5586       if (REAL_WORDS_BIG_ENDIAN)
5587         p = XFbignan;
5588       else
5589         p = XFlittlenan;
5590       break;
5591     case DFmode:
5592       n = 4;
5593       if (REAL_WORDS_BIG_ENDIAN)
5594         p = DFbignan;
5595       else
5596         p = DFlittlenan;
5597       break;
5598     case HFmode:
5599     case SFmode:
5600       n = 2;
5601       if (REAL_WORDS_BIG_ENDIAN)
5602         p = SFbignan;
5603       else
5604         p = SFlittlenan;
5605       break;
5606 #endif
5607     default:
5608       abort ();
5609     }
5610   if (REAL_WORDS_BIG_ENDIAN)
5611     *nan++ = (sign << 15) | *p++;
5612   while (--n != 0)
5613     *nan++ = *p++;
5614   if (! REAL_WORDS_BIG_ENDIAN)
5615     *nan = (sign << 15) | *p;
5616 }
5617
5618 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
5619    This is the inverse of the function `etarsingle' invoked by
5620    REAL_VALUE_TO_TARGET_SINGLE.  */
5621
5622 REAL_VALUE_TYPE
5623 ereal_from_float (f)
5624      HOST_WIDE_INT f;
5625 {
5626   REAL_VALUE_TYPE r;
5627   unsigned EMUSHORT s[2];
5628   unsigned EMUSHORT e[NE];
5629
5630   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
5631    This is the inverse operation to what the function `endian' does.  */
5632   if (REAL_WORDS_BIG_ENDIAN)
5633     {
5634       s[0] = (unsigned EMUSHORT) (f >> 16);
5635       s[1] = (unsigned EMUSHORT) f;
5636     }
5637   else
5638     {
5639       s[0] = (unsigned EMUSHORT) f;
5640       s[1] = (unsigned EMUSHORT) (f >> 16);
5641     }
5642   /* Convert and promote the target float to E-type. */
5643   e24toe (s, e);
5644   /* Output E-type to REAL_VALUE_TYPE. */
5645   PUT_REAL (e, &r);
5646   return r;
5647 }
5648
5649
5650 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
5651    This is the inverse of the function `etardouble' invoked by
5652    REAL_VALUE_TO_TARGET_DOUBLE.
5653
5654    The DFmode is stored as an array of HOST_WIDE_INT in the target's
5655    data format, with no holes in the bit packing.  The first element
5656    of the input array holds the bits that would come first in the
5657    target computer's memory.  */
5658
5659 REAL_VALUE_TYPE
5660 ereal_from_double (d)
5661      HOST_WIDE_INT d[];
5662 {
5663   REAL_VALUE_TYPE r;
5664   unsigned EMUSHORT s[4];
5665   unsigned EMUSHORT e[NE];
5666
5667   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
5668   if (REAL_WORDS_BIG_ENDIAN)
5669     {
5670       s[0] = (unsigned EMUSHORT) (d[0] >> 16);
5671       s[1] = (unsigned EMUSHORT) d[0];
5672 #if HOST_BITS_PER_WIDE_INT == 32
5673       s[2] = (unsigned EMUSHORT) (d[1] >> 16);
5674       s[3] = (unsigned EMUSHORT) d[1];
5675 #else
5676       /* In this case the entire target double is contained in the
5677          first array element.  The second element of the input is
5678          ignored.  */
5679       s[2] = (unsigned EMUSHORT) (d[0] >> 48);
5680       s[3] = (unsigned EMUSHORT) (d[0] >> 32);
5681 #endif
5682     }
5683   else
5684     {
5685       /* Target float words are little-endian.  */
5686       s[0] = (unsigned EMUSHORT) d[0];
5687       s[1] = (unsigned EMUSHORT) (d[0] >> 16);
5688 #if HOST_BITS_PER_WIDE_INT == 32
5689       s[2] = (unsigned EMUSHORT) d[1];
5690       s[3] = (unsigned EMUSHORT) (d[1] >> 16);
5691 #else
5692       s[2] = (unsigned EMUSHORT) (d[0] >> 32);
5693       s[3] = (unsigned EMUSHORT) (d[0] >> 48);
5694 #endif
5695     }
5696   /* Convert target double to E-type. */
5697   e53toe (s, e);
5698   /* Output E-type to REAL_VALUE_TYPE. */
5699   PUT_REAL (e, &r);
5700   return r;
5701 }
5702
5703
5704 /* Convert target computer unsigned 64-bit integer to e-type.
5705    The endian-ness of DImode follows the convention for integers,
5706    so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN.  */
5707
5708 static void
5709 uditoe (di, e)
5710      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
5711      unsigned EMUSHORT *e;
5712 {
5713   unsigned EMUSHORT yi[NI];
5714   int k;
5715
5716   ecleaz (yi);
5717   if (WORDS_BIG_ENDIAN)
5718     {
5719       for (k = M; k < M + 4; k++)
5720         yi[k] = *di++;
5721     }
5722   else
5723     {
5724       for (k = M + 3; k >= M; k--)
5725         yi[k] = *di++;
5726     }
5727   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
5728   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5729     ecleaz (yi);                /* it was zero */
5730   else
5731     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5732   emovo (yi, e);
5733 }
5734
5735 /* Convert target computer signed 64-bit integer to e-type. */
5736
5737 static void
5738 ditoe (di, e)
5739      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
5740      unsigned EMUSHORT *e;
5741 {
5742   unsigned EMULONG acc;
5743   unsigned EMUSHORT yi[NI];
5744   unsigned EMUSHORT carry;
5745   int k, sign;
5746
5747   ecleaz (yi);
5748   if (WORDS_BIG_ENDIAN)
5749     {
5750       for (k = M; k < M + 4; k++)
5751         yi[k] = *di++;
5752     }
5753   else
5754     {
5755       for (k = M + 3; k >= M; k--)
5756         yi[k] = *di++;
5757     }
5758   /* Take absolute value */
5759   sign = 0;
5760   if (yi[M] & 0x8000)
5761     {
5762       sign = 1;
5763       carry = 0;
5764       for (k = M + 3; k >= M; k--)
5765         {
5766           acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
5767           yi[k] = acc;
5768           carry = 0;
5769           if (acc & 0x10000)
5770             carry = 1;
5771         }
5772     }
5773   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
5774   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5775     ecleaz (yi);                /* it was zero */
5776   else
5777     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5778   emovo (yi, e);
5779   if (sign)
5780         eneg (e);
5781 }
5782
5783
5784 /* Convert e-type to unsigned 64-bit int. */
5785
5786 static void
5787 etoudi (x, i)
5788      unsigned EMUSHORT *x;
5789      unsigned EMUSHORT *i;
5790 {
5791   unsigned EMUSHORT xi[NI];
5792   int j, k;
5793
5794   emovi (x, xi);
5795   if (xi[0])
5796     {
5797       xi[M] = 0;
5798       goto noshift;
5799     }
5800   k = (int) xi[E] - (EXONE - 1);
5801   if (k <= 0)
5802     {
5803       for (j = 0; j < 4; j++)
5804         *i++ = 0;
5805       return;
5806     }
5807   if (k > 64)
5808     {
5809       for (j = 0; j < 4; j++)
5810         *i++ = 0xffff;
5811       if (extra_warnings)
5812         warning ("overflow on truncation to integer");
5813       return;
5814     }
5815   if (k > 16)
5816     {
5817       /* Shift more than 16 bits: first shift up k-16 mod 16,
5818          then shift up by 16's.  */
5819       j = k - ((k >> 4) << 4);
5820       if (j == 0)
5821         j = 16;
5822       eshift (xi, j);
5823       if (WORDS_BIG_ENDIAN)
5824         *i++ = xi[M];
5825       else
5826         {
5827           i += 3;
5828           *i-- = xi[M];
5829         }
5830       k -= j;
5831       do
5832         {
5833           eshup6 (xi);
5834           if (WORDS_BIG_ENDIAN)
5835             *i++ = xi[M];
5836           else
5837             *i-- = xi[M];
5838         }
5839       while ((k -= 16) > 0);
5840     }
5841   else
5842     {
5843         /* shift not more than 16 bits */
5844       eshift (xi, k);
5845
5846 noshift:
5847
5848       if (WORDS_BIG_ENDIAN)
5849         {
5850           i += 3;
5851           *i-- = xi[M];
5852           *i-- = 0;
5853           *i-- = 0;
5854           *i = 0;
5855         }
5856       else
5857         {
5858           *i++ = xi[M];
5859           *i++ = 0;
5860           *i++ = 0;
5861           *i = 0;
5862         }
5863     }
5864 }
5865
5866
5867 /* Convert e-type to signed 64-bit int. */
5868
5869 static void
5870 etodi (x, i)
5871      unsigned EMUSHORT *x;
5872      unsigned EMUSHORT *i;
5873 {
5874   unsigned EMULONG acc;
5875   unsigned EMUSHORT xi[NI];
5876   unsigned EMUSHORT carry;
5877   unsigned EMUSHORT *isave;
5878   int j, k;
5879
5880   emovi (x, xi);
5881   k = (int) xi[E] - (EXONE - 1);
5882   if (k <= 0)
5883     {
5884       for (j = 0; j < 4; j++)
5885         *i++ = 0;
5886       return;
5887     }
5888   if (k > 64)
5889     {
5890       for (j = 0; j < 4; j++)
5891         *i++ = 0xffff;
5892       if (extra_warnings)
5893         warning ("overflow on truncation to integer");
5894       return;
5895     }
5896   isave = i;
5897   if (k > 16)
5898     {
5899       /* Shift more than 16 bits: first shift up k-16 mod 16,
5900          then shift up by 16's.  */
5901       j = k - ((k >> 4) << 4);
5902       if (j == 0)
5903         j = 16;
5904       eshift (xi, j);
5905       if (WORDS_BIG_ENDIAN)
5906         *i++ = xi[M];
5907       else
5908         {
5909           i += 3;
5910           *i-- = xi[M];
5911         }
5912       k -= j;
5913       do
5914         {
5915           eshup6 (xi);
5916           if (WORDS_BIG_ENDIAN)
5917             *i++ = xi[M];
5918           else
5919             *i-- = xi[M];
5920         }
5921       while ((k -= 16) > 0);
5922     }
5923   else
5924     {
5925         /* shift not more than 16 bits */
5926       eshift (xi, k);
5927
5928       if (WORDS_BIG_ENDIAN)
5929         {
5930           i += 3;
5931           *i = xi[M];
5932           *i-- = 0;
5933           *i-- = 0;
5934           *i = 0;
5935         }
5936       else
5937         {
5938           *i++ = xi[M];
5939           *i++ = 0;
5940           *i++ = 0;
5941           *i = 0;
5942         }
5943     }
5944   /* Negate if negative */
5945   if (xi[0])
5946     {
5947       carry = 0;
5948       if (WORDS_BIG_ENDIAN)
5949         isave += 3;
5950       for (k = 0; k < 4; k++)
5951         {
5952           acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
5953           if (WORDS_BIG_ENDIAN)
5954             *isave-- = acc;
5955           else
5956             *isave++ = acc;
5957           carry = 0;
5958           if (acc & 0x10000)
5959             carry = 1;
5960         }
5961     }
5962 }
5963
5964
5965 /* Longhand square root routine. */
5966
5967
5968 static int esqinited = 0;
5969 static unsigned short sqrndbit[NI];
5970
5971 static void
5972 esqrt (x, y)
5973      unsigned EMUSHORT *x, *y;
5974 {
5975   unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
5976   EMULONG m, exp;
5977   int i, j, k, n, nlups;
5978
5979   if (esqinited == 0)
5980     {
5981       ecleaz (sqrndbit);
5982       sqrndbit[NI - 2] = 1;
5983       esqinited = 1;
5984     }
5985   /* Check for arg <= 0 */
5986   i = ecmp (x, ezero);
5987   if (i <= 0)
5988     {
5989       if (i == -1)
5990         {
5991           mtherr ("esqrt", DOMAIN);
5992           eclear (y);
5993         }
5994       else
5995         emov (x, y);
5996       return;
5997     }
5998
5999 #ifdef INFINITY
6000   if (eisinf (x))
6001     {
6002       eclear (y);
6003       einfin (y);
6004       return;
6005     }
6006 #endif
6007   /* Bring in the arg and renormalize if it is denormal. */
6008   emovi (x, xx);
6009   m = (EMULONG) xx[1];          /* local long word exponent */
6010   if (m == 0)
6011     m -= enormlz (xx);
6012
6013   /* Divide exponent by 2 */
6014   m -= 0x3ffe;
6015   exp = (unsigned short) ((m / 2) + 0x3ffe);
6016
6017   /* Adjust if exponent odd */
6018   if ((m & 1) != 0)
6019     {
6020       if (m > 0)
6021         exp += 1;
6022       eshdn1 (xx);
6023     }
6024
6025   ecleaz (sq);
6026   ecleaz (num);
6027   n = 8;                        /* get 8 bits of result per inner loop */
6028   nlups = rndprc;
6029   j = 0;
6030
6031   while (nlups > 0)
6032     {
6033       /* bring in next word of arg */
6034       if (j < NE)
6035         num[NI - 1] = xx[j + 3];
6036       /* Do additional bit on last outer loop, for roundoff. */
6037       if (nlups <= 8)
6038         n = nlups + 1;
6039       for (i = 0; i < n; i++)
6040         {
6041           /* Next 2 bits of arg */
6042           eshup1 (num);
6043           eshup1 (num);
6044           /* Shift up answer */
6045           eshup1 (sq);
6046           /* Make trial divisor */
6047           for (k = 0; k < NI; k++)
6048             temp[k] = sq[k];
6049           eshup1 (temp);
6050           eaddm (sqrndbit, temp);
6051           /* Subtract and insert answer bit if it goes in */
6052           if (ecmpm (temp, num) <= 0)
6053             {
6054               esubm (temp, num);
6055               sq[NI - 2] |= 1;
6056             }
6057         }
6058       nlups -= n;
6059       j += 1;
6060     }
6061
6062   /* Adjust for extra, roundoff loop done. */
6063   exp += (NBITS - 1) - rndprc;
6064
6065   /* Sticky bit = 1 if the remainder is nonzero. */
6066   k = 0;
6067   for (i = 3; i < NI; i++)
6068     k |= (int) num[i];
6069
6070   /* Renormalize and round off. */
6071   emdnorm (sq, k, 0, exp, 64);
6072   emovo (sq, y);
6073 }
6074 #endif /* EMU_NON_COMPILE not defined */
6075 \f
6076 /* Return the binary precision of the significand for a given
6077    floating point mode.  The mode can hold an integer value
6078    that many bits wide, without losing any bits.  */
6079
6080 int
6081 significand_size (mode)
6082      enum machine_mode mode;
6083 {
6084
6085 switch (mode)
6086   {
6087   case SFmode:
6088     return 24;
6089
6090   case DFmode:
6091 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6092     return 53;
6093 #else
6094 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6095     return 56;
6096 #else
6097 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6098     return 56;
6099 #else
6100     abort ();
6101 #endif
6102 #endif
6103 #endif
6104
6105   case XFmode:
6106     return 64;
6107   case TFmode:
6108     return 113;
6109
6110   default:
6111     abort ();
6112   }
6113 }