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