Return to real.c CVS log

Up to [Apple XNU] / GNUtools / cc

Annotation of GNUtools/cc/real.c, revision 1.1

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.