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