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