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