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