OpenCores
URL https://opencores.org/ocsvn/or1k/or1k/trunk

Subversion Repositories or1k

[/] [or1k/] [trunk/] [newlib-1.10.0/] [newlib/] [libc/] [stdlib/] [ldtoa.c] - Blame information for rev 1773

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 1010 ivang
 
2
 /* Extended precision arithmetic functions for long double I/O.
3
  * This program has been placed in the public domain.
4
  */
5
 
6
#include <_ansi.h>
7
#include <reent.h>
8
#include <string.h>
9
#include <stdlib.h>
10
#include "mprec.h"
11
 
12
/* These are the externally visible entries. */
13
/* linux name:  long double _IO_strtold (char *, char **); */
14
long double _strtold (char *, char **);
15
char * _ldtoa_r (struct _reent *, long double, int, int, int *, int *, char **);
16
int    _ldcheck (long double *);
17
#if 0
18
void _IO_ldtostr(long double *, char *, int, int, char);
19
#endif
20
 
21
 /* Number of 16 bit words in external x type format */
22
 #define NE 10
23
 
24
 /* Number of 16 bit words in internal format */
25
 #define NI (NE+3)
26
 
27
 /* Array offset to exponent */
28
 #define E 1
29
 
30
 /* Array offset to high guard word */
31
 #define M 2
32
 
33
 /* Number of bits of precision */
34
 #define NBITS ((NI-4)*16)
35
 
36
 /* Maximum number of decimal digits in ASCII conversion
37
  * = NBITS*log10(2)
38
  */
39
 #define NDEC (NBITS*8/27)
40
 
41
 /* The exponent of 1.0 */
42
 #define EXONE (0x3fff)
43
 
44
/* Control structure for long doublue conversion including rounding precision values.
45
 * rndprc can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
46
 */
47
typedef struct
48
{
49
  int rlast;
50
  int rndprc;
51
  int rw;
52
  int re;
53
  int outexpon;
54
  unsigned short rmsk;
55
  unsigned short rmbit;
56
  unsigned short rebit;
57
  unsigned short rbit[NI];
58
  unsigned short equot[NI];
59
} LDPARMS;
60
 
61
static void esub(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp);
62
static void emul(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp);
63
static void ediv(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp);
64
static int ecmp(short unsigned int *a, short unsigned int *b);
65
static int enormlz(short unsigned int *x);
66
static int eshift(short unsigned int *x, int sc);
67
static void eshup1(register short unsigned int *x);
68
static void eshup8(register short unsigned int *x);
69
static void eshup6(register short unsigned int *x);
70
static void eshdn1(register short unsigned int *x);
71
static void eshdn8(register short unsigned int *x);
72
static void eshdn6(register short unsigned int *x);
73
static void eneg(short unsigned int *x);
74
static void emov(register short unsigned int *a, register short unsigned int *b);
75
static void eclear(register short unsigned int *x);
76
static void einfin(register short unsigned int *x, register LDPARMS *ldp);
77
static void efloor(short unsigned int *x, short unsigned int *y, LDPARMS *ldp);
78
static void etoasc(short unsigned int *x, char *string, int ndigs, int outformat, LDPARMS *ldp);
79
 
80
#if LDBL_MANT_DIG == 24
81
static void e24toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
82
#elif LDBL_MANT_DIG == 53
83
static void e53toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
84
#elif LDBL_MANT_DIG == 64
85
static void e64toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
86
#else
87
static void e113toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
88
#endif
89
 
90
/*                                                      econst.c        */
91
/*  e type constants used by high precision check routines */
92
 
93
#if NE == 10
94
/* 0.0 */
95
static unsigned short ezero[NE] =
96
 {0x0000, 0x0000, 0x0000, 0x0000,
97
  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
98
 
99
/* 1.0E0 */
100
static unsigned short eone[NE] =
101
 {0x0000, 0x0000, 0x0000, 0x0000,
102
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
103
 
104
#else
105
 
106
/* 0.0 */
107
static unsigned short ezero[NE] = {
108
0, 0000000,0000000,0000000,0000000,0000000,};
109
/* 1.0E0 */
110
static unsigned short eone[NE] = {
111
0, 0000000,0000000,0000000,0100000,0x3fff,};
112
 
113
#endif
114
 
115
/* Debugging routine for displaying errors */
116
#ifdef DEBUG
117
/* Notice: the order of appearance of the following
118
 * messages is bound to the error codes defined
119
 * in mconf.h.
120
 */
121
static char *ermsg[7] = {
122
"unknown",      /* error code 0 */
123
"domain",       /* error code 1 */
124
"singularity",  /* et seq.      */
125
"overflow",
126
"underflow",
127
"total loss of precision",
128
"partial loss of precision"
129
};
130
#define mtherr(name, code) printf( "\n%s %s error\n", name, ermsg[code] );
131
#else
132
#define mtherr(name, code)
133
#endif
134
 
135
/*                                                      ieee.c
136
 *
137
 *    Extended precision IEEE binary floating point arithmetic routines
138
 *
139
 * Numbers are stored in C language as arrays of 16-bit unsigned
140
 * short integers.  The arguments of the routines are pointers to
141
 * the arrays.
142
 *
143
 *
144
 * External e type data structure, simulates Intel 8087 chip
145
 * temporary real format but possibly with a larger significand:
146
 *
147
 *      NE-1 significand words  (least significant word first,
148
 *                               most significant bit is normally set)
149
 *      exponent                (value = EXONE for 1.0,
150
 *                              top bit is the sign)
151
 *
152
 *
153
 * Internal data structure of a number (a "word" is 16 bits):
154
 *
155
 * ei[0]        sign word       (0 for positive, 0xffff for negative)
156
 * ei[1]        biased exponent (value = EXONE for the number 1.0)
157
 * ei[2]        high guard word (always zero after normalization)
158
 * ei[3]
159
 * to ei[NI-2]  significand     (NI-4 significand words,
160
 *                               most significant word first,
161
 *                               most significant bit is set)
162
 * ei[NI-1]     low guard word  (0x8000 bit is rounding place)
163
 *
164
 *
165
 *
166
 *              Routines for external format numbers
167
 *
168
 *      asctoe( string, e )     ASCII string to extended double e type
169
 *      asctoe64( string, &d )  ASCII string to long double
170
 *      asctoe53( string, &d )  ASCII string to double
171
 *      asctoe24( string, &f )  ASCII string to single
172
 *      asctoeg( string, e, prec, ldp ) ASCII string to specified precision
173
 *      e24toe( &f, e, ldp )    IEEE single precision to e type
174
 *      e53toe( &d, e, ldp )    IEEE double precision to e type
175
 *      e64toe( &d, e, ldp )    IEEE long double precision to e type
176
 *      e113toe( &d, e, ldp )   IEEE long double precision to e type
177
 *      eabs(e)                 absolute value
178
 *      eadd( a, b, c )         c = b + a
179
 *      eclear(e)               e = 0
180
 *      ecmp (a, b)             Returns 1 if a > b, 0 if a == b,
181
 *                              -1 if a < b, -2 if either a or b is a NaN.
182
 *      ediv( a, b, c, ldp )    c = b / a
183
 *      efloor( a, b, ldp )     truncate to integer, toward -infinity
184
 *      efrexp( a, exp, s )     extract exponent and significand
185
 *      eifrac( e, &l, frac )   e to long integer and e type fraction
186
 *      euifrac( e, &l, frac )  e to unsigned long integer and e type fraction
187
 *      einfin( e, ldp )        set e to infinity, leaving its sign alone
188
 *      eldexp( a, n, b )       multiply by 2**n
189
 *      emov( a, b )            b = a
190
 *      emul( a, b, c, ldp )    c = b * a
191
 *      eneg(e)                 e = -e
192
 *      eround( a, b )          b = nearest integer value to a
193
 *      esub( a, b, c, ldp )    c = b - a
194
 *      e24toasc( &f, str, n )  single to ASCII string, n digits after decimal
195
 *      e53toasc( &d, str, n )  double to ASCII string, n digits after decimal
196
 *      e64toasc( &d, str, n )  long double to ASCII string
197
 *      etoasc(e,str,n,fmt,ldp)e to ASCII string, n digits after decimal
198
 *      etoe24( e, &f )         convert e type to IEEE single precision
199
 *      etoe53( e, &d )         convert e type to IEEE double precision
200
 *      etoe64( e, &d )         convert e type to IEEE long double precision
201
 *      ltoe( &l, e )           long (32 bit) integer to e type
202
 *      ultoe( &l, e )          unsigned long (32 bit) integer to e type
203
 *      eisneg( e )             1 if sign bit of e != 0, else 0
204
 *      eisinf( e )             1 if e has maximum exponent (non-IEEE)
205
 *                              or is infinite (IEEE)
206
 *      eisnan( e )             1 if e is a NaN
207
 *      esqrt( a, b )           b = square root of a
208
 *
209
 *
210
 *              Routines for internal format numbers
211
 *
212
 *      eaddm( ai, bi )         add significands, bi = bi + ai
213
 *      ecleaz(ei)              ei = 0
214
 *      ecleazs(ei)             set ei = 0 but leave its sign alone
215
 *      ecmpm( ai, bi )         compare significands, return 1, 0, or -1
216
 *      edivm( ai, bi, ldp )    divide  significands, bi = bi / ai
217
 *      emdnorm(ai,l,s,exp,ldp) normalize and round off
218
 *      emovi( a, ai )          convert external a to internal ai
219
 *      emovo( ai, a, ldp )     convert internal ai to external a
220
 *      emovz( ai, bi )         bi = ai, low guard word of bi = 0
221
 *      emulm( ai, bi, ldp )    multiply significands, bi = bi * ai
222
 *      enormlz(ei)             left-justify the significand
223
 *      eshdn1( ai )            shift significand and guards down 1 bit
224
 *      eshdn8( ai )            shift down 8 bits
225
 *      eshdn6( ai )            shift down 16 bits
226
 *      eshift( ai, n )         shift ai n bits up (or down if n < 0)
227
 *      eshup1( ai )            shift significand and guards up 1 bit
228
 *      eshup8( ai )            shift up 8 bits
229
 *      eshup6( ai )            shift up 16 bits
230
 *      esubm( ai, bi )         subtract significands, bi = bi - ai
231
 *
232
 *
233
 * The result is always normalized and rounded to NI-4 word precision
234
 * after each arithmetic operation.
235
 *
236
 * Exception flags are NOT fully supported.
237
 *
238
 * Define INFINITY in mconf.h for support of infinity; otherwise a
239
 * saturation arithmetic is implemented.
240
 *
241
 * Define NANS for support of Not-a-Number items; otherwise the
242
 * arithmetic will never produce a NaN output, and might be confused
243
 * by a NaN input.
244
 * If NaN's are supported, the output of ecmp(a,b) is -2 if
245
 * either a or b is a NaN. This means asking if(ecmp(a,b) < 0)
246
 * may not be legitimate. Use if(ecmp(a,b) == -1) for less-than
247
 * if in doubt.
248
 * Signaling NaN's are NOT supported; they are treated the same
249
 * as quiet NaN's.
250
 *
251
 * Denormals are always supported here where appropriate (e.g., not
252
 * for conversion to DEC numbers).
253
 */
254
 
255
/*
256
 * Revision history:
257
 *
258
 *  5 Jan 84    PDP-11 assembly language version
259
 *  6 Dec 86    C language version
260
 * 30 Aug 88    100 digit version, improved rounding
261
 * 15 May 92    80-bit long double support
262
 * 22 Nov 00    Revised to fit into newlib by Jeff Johnston <jjohnstn@redhat.com>
263
 *
264
 * Author:  S. L. Moshier.
265
 *
266
 * Copyright (c) 1984,2000 S.L. Moshier
267
 *
268
 * Permission to use, copy, modify, and distribute this software for any
269
 * purpose without fee is hereby granted, provided that this entire notice
270
 * is included in all copies of any software which is or includes a copy
271
 * or modification of this software and in all copies of the supporting
272
 * documentation for such software.
273
 *
274
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
275
 * WARRANTY.  IN PARTICULAR,  THE AUTHOR MAKES NO REPRESENTATION
276
 * OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY OF THIS
277
 * SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
278
 *
279
 */
280
 
281
#include <stdio.h>
282
/* #include "\usr\include\stdio.h" */
283
/*#include "ehead.h"*/
284
/*#include "mconf.h"*/
285
/*                                                      mconf.h
286
 *
287
 *      Common include file for math routines
288
 *
289
 *
290
 *
291
 * SYNOPSIS:
292
 *
293
 * #include "mconf.h"
294
 *
295
 *
296
 *
297
 * DESCRIPTION:
298
 *
299
 * This file contains definitions for error codes that are
300
 * passed to the common error handling routine mtherr()
301
 * (which see).
302
 *
303
 * The file also includes a conditional assembly definition
304
 * for the type of computer arithmetic (IEEE, DEC, Motorola
305
 * IEEE, or UNKnown).
306
 *
307
 * For Digital Equipment PDP-11 and VAX computers, certain
308
 * IBM systems, and others that use numbers with a 56-bit
309
 * significand, the symbol DEC should be defined.  In this
310
 * mode, most floating point constants are given as arrays
311
 * of octal integers to eliminate decimal to binary conversion
312
 * errors that might be introduced by the compiler.
313
 *
314
 * For computers, such as IBM PC, that follow the IEEE
315
 * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
316
 * Std 754-1985), the symbol IBMPC should be defined.  These
317
 * numbers have 53-bit significands.  In this mode, constants
318
 * are provided as arrays of hexadecimal 16 bit integers.
319
 *
320
 * To accommodate other types of computer arithmetic, all
321
 * constants are also provided in a normal decimal radix
322
 * which one can hope are correctly converted to a suitable
323
 * format by the available C language compiler.  To invoke
324
 * this mode, the symbol UNK is defined.
325
 *
326
 * An important difference among these modes is a predefined
327
 * set of machine arithmetic constants for each.  The numbers
328
 * MACHEP (the machine roundoff error), MAXNUM (largest number
329
 * represented), and several other parameters are preset by
330
 * the configuration symbol.  Check the file const.c to
331
 * ensure that these values are correct for your computer.
332
 *
333
 * For ANSI C compatibility, define ANSIC equal to 1.  Currently
334
 * this affects only the atan2() function and others that use it.
335
 */
336
 
337
/* Constant definitions for math error conditions
338
 */
339
 
340
#define DOMAIN          1       /* argument domain error */
341
#define SING            2       /* argument singularity */
342
#define OVERFLOW        3       /* overflow range error */
343
#define UNDERFLOW       4       /* underflow range error */
344
#define TLOSS           5       /* total loss of precision */
345
#define PLOSS           6       /* partial loss of precision */
346
 
347
#define EDOM            33
348
#define ERANGE          34
349
 
350
typedef struct
351
        {
352
        double r;
353
        double i;
354
        }cmplx;
355
 
356
/* Type of computer arithmetic */
357
 
358
#ifndef DEC
359
#ifdef __IEEE_LITTLE_ENDIAN
360
#define IBMPC 1
361
#else  /* !__IEEE_LITTLE_ENDIAN */
362
#define MIEEE 1
363
#endif /* !__IEEE_LITTLE_ENDIAN */
364
#endif /* !DEC */
365
 
366
/* Define 1 for ANSI C atan2() function
367
 * See atan.c and clog.c.
368
 */
369
#define ANSIC 1
370
 
371
/*define VOLATILE volatile*/
372
#define VOLATILE 
373
 
374
#define NANS
375
#define INFINITY
376
 
377
/* NaN's require infinity support. */
378
#ifdef NANS
379
#ifndef INFINITY
380
#define INFINITY
381
#endif
382
#endif
383
 
384
/* This handles 64-bit long ints. */
385
#define LONGBITS (8 * sizeof(long))
386
 
387
 
388
static void eaddm(short unsigned int *x, short unsigned int *y);
389
static void esubm(short unsigned int *x, short unsigned int *y);
390
static void emdnorm(short unsigned int *s, int lost, int subflg, long int exp, int rcntrl, LDPARMS *ldp);
391
static int  asctoeg(char *ss, short unsigned int *y, int oprec, LDPARMS *ldp);
392
static void enan(short unsigned int *nan, int size);
393
#if LDBL_MANT_DIG == 24
394
static void toe24(short unsigned int *x, short unsigned int *y);
395
#elif LDBL_MANT_DIG == 53
396
static void toe53(short unsigned int *x, short unsigned int *y);
397
#elif LDBL_MANT_DIG == 64
398
static void toe64(short unsigned int *a, short unsigned int *b);
399
#else
400
static void toe113(short unsigned int *a, short unsigned int *b);
401
#endif
402
static void eiremain(short unsigned int *den, short unsigned int *num, LDPARMS *ldp);
403
static int ecmpm(register short unsigned int *a, register short unsigned int *b);
404
static int edivm(short unsigned int *den, short unsigned int *num, LDPARMS *ldp);
405
static int emulm(short unsigned int *a, short unsigned int *b, LDPARMS *ldp);
406
static int eisneg(short unsigned int *x);
407
static int eisinf(short unsigned int *x);
408
static void emovi(short unsigned int *a, short unsigned int *b);
409
static void emovo(short unsigned int *a, short unsigned int *b, LDPARMS *ldp);
410
static void emovz(register short unsigned int *a, register short unsigned int *b);
411
static void ecleaz(register short unsigned int *xi);
412
static void eadd1(short unsigned int *a, short unsigned int *b, short unsigned int *c, int subflg, LDPARMS *ldp);
413
static int eisnan(short unsigned int *x);
414
static int eiisnan(short unsigned int *x);
415
 
416
#ifdef DEC
417
static void etodec(), todec(), dectoe();
418
#endif
419
 
420
/*
421
; Clear out entire external format number.
422
;
423
; unsigned short x[];
424
; eclear( x );
425
*/
426
 
427
static void eclear(register short unsigned int *x)
428
{
429
register int i;
430
 
431
for( i=0; i<NE; i++ )
432
        *x++ = 0;
433
}
434
 
435
 
436
 
437
/* Move external format number from a to b.
438
 *
439
 * emov( a, b );
440
 */
441
 
442
static void emov(register short unsigned int *a, register short unsigned int *b)
443
{
444
register int i;
445
 
446
for( i=0; i<NE; i++ )
447
        *b++ = *a++;
448
}
449
 
450
 
451
/*
452
;       Negate external format number
453
;
454
;       unsigned short x[NE];
455
;       eneg( x );
456
*/
457
 
458
static void eneg(short unsigned int *x)
459
{
460
 
461
#ifdef NANS
462
if( eisnan(x) )
463
        return;
464
#endif
465
x[NE-1] ^= 0x8000; /* Toggle the sign bit */
466
}
467
 
468
 
469
 
470
/* Return 1 if external format number is negative,
471
 * else return zero.
472
 */
473
static int eisneg(short unsigned int *x)
474
{
475
 
476
#ifdef NANS
477
if( eisnan(x) )
478
        return( 0 );
479
#endif
480
if( x[NE-1] & 0x8000 )
481
        return( 1 );
482
else
483
        return( 0 );
484
}
485
 
486
 
487
/* Return 1 if external format number has maximum possible exponent,
488
 * else return zero.
489
 */
490
static int eisinf(short unsigned int *x)
491
{
492
 
493
if( (x[NE-1] & 0x7fff) == 0x7fff )
494
        {
495
#ifdef NANS
496
        if( eisnan(x) )
497
                return( 0 );
498
#endif
499
        return( 1 );
500
        }
501
else
502
        return( 0 );
503
}
504
 
505
/* Check if e-type number is not a number.
506
 */
507
static int eisnan(short unsigned int *x)
508
{
509
 
510
#ifdef NANS
511
int i;
512
/* NaN has maximum exponent */
513
if( (x[NE-1] & 0x7fff) != 0x7fff )
514
        return (0);
515
/* ... and non-zero significand field. */
516
for( i=0; i<NE-1; i++ )
517
        {
518
        if( *x++ != 0 )
519
                return (1);
520
        }
521
#endif
522
return (0);
523
}
524
 
525
/*
526
; Fill entire number, including exponent and significand, with
527
; largest possible number.  These programs implement a saturation
528
; value that is an ordinary, legal number.  A special value
529
; "infinity" may also be implemented; this would require tests
530
; for that value and implementation of special rules for arithmetic
531
; operations involving inifinity.
532
*/
533
 
534
static void einfin(register short unsigned int *x, register LDPARMS *ldp)
535
{
536
register int i;
537
 
538
#ifdef INFINITY
539
for( i=0; i<NE-1; i++ )
540
        *x++ = 0;
541
*x |= 32767;
542
ldp = ldp;
543
#else
544
for( i=0; i<NE-1; i++ )
545
        *x++ = 0xffff;
546
*x |= 32766;
547
if( ldp->rndprc < NBITS )
548
        {
549
        if (ldp->rndprc == 113)
550
                {
551
                *(x - 9) = 0;
552
                *(x - 8) = 0;
553
                }
554
        if( ldp->rndprc == 64 )
555
                {
556
                *(x-5) = 0;
557
                }
558
        if( ldp->rndprc == 53 )
559
                {
560
                *(x-4) = 0xf800;
561
                }
562
        else
563
                {
564
                *(x-4) = 0;
565
                *(x-3) = 0;
566
                *(x-2) = 0xff00;
567
                }
568
        }
569
#endif
570
}
571
 
572
/* Move in external format number,
573
 * converting it to internal format.
574
 */
575
static void emovi(short unsigned int *a, short unsigned int *b)
576
{
577
register unsigned short *p, *q;
578
int i;
579
 
580
q = b;
581
p = a + (NE-1); /* point to last word of external number */
582
/* get the sign bit */
583
if( *p & 0x8000 )
584
        *q++ = 0xffff;
585
else
586
        *q++ = 0;
587
/* get the exponent */
588
*q = *p--;
589
*q++ &= 0x7fff; /* delete the sign bit */
590
#ifdef INFINITY
591
if( (*(q-1) & 0x7fff) == 0x7fff )
592
        {
593
#ifdef NANS
594
        if( eisnan(a) )
595
                {
596
                *q++ = 0;
597
                for( i=3; i<NI; i++ )
598
                        *q++ = *p--;
599
                return;
600
                }
601
#endif
602
        for( i=2; i<NI; i++ )
603
                *q++ = 0;
604
        return;
605
        }
606
#endif
607
/* clear high guard word */
608
*q++ = 0;
609
/* move in the significand */
610
for( i=0; i<NE-1; i++ )
611
        *q++ = *p--;
612
/* clear low guard word */
613
*q = 0;
614
}
615
 
616
 
617
/* Move internal format number out,
618
 * converting it to external format.
619
 */
620
static void emovo(short unsigned int *a, short unsigned int *b, LDPARMS *ldp)
621
{
622
register unsigned short *p, *q;
623
unsigned short i;
624
 
625
p = a;
626
q = b + (NE-1); /* point to output exponent */
627
/* combine sign and exponent */
628
i = *p++;
629
if( i )
630
        *q-- = *p++ | 0x8000;
631
else
632
        *q-- = *p++;
633
#ifdef INFINITY
634
if( *(p-1) == 0x7fff )
635
        {
636
#ifdef NANS
637
        if( eiisnan(a) )
638
                {
639
                enan( b, NBITS );
640
                return;
641
                }
642
#endif
643
        einfin(b, ldp);
644
        return;
645
        }
646
#endif
647
/* skip over guard word */
648
++p;
649
/* move the significand */
650
for( i=0; i<NE-1; i++ )
651
        *q-- = *p++;
652
}
653
 
654
 
655
/* Clear out internal format number.
656
 */
657
 
658
static void ecleaz(register short unsigned int *xi)
659
{
660
register int i;
661
 
662
for( i=0; i<NI; i++ )
663
        *xi++ = 0;
664
}
665
 
666
/* same, but don't touch the sign. */
667
 
668
static void ecleazs(register short unsigned int *xi)
669
{
670
register int i;
671
 
672
++xi;
673
for(i=0; i<NI-1; i++)
674
        *xi++ = 0;
675
}
676
 
677
 
678
 
679
 
680
/* Move internal format number from a to b.
681
 */
682
static void emovz(register short unsigned int *a, register short unsigned int *b)
683
{
684
register int i;
685
 
686
for( i=0; i<NI-1; i++ )
687
        *b++ = *a++;
688
/* clear low guard word */
689
*b = 0;
690
}
691
 
692
/* Return nonzero if internal format number is a NaN.
693
 */
694
 
695
static int eiisnan (short unsigned int *x)
696
{
697
int i;
698
 
699
if( (x[E] & 0x7fff) == 0x7fff )
700
        {
701
        for( i=M+1; i<NI; i++ )
702
                {
703
                if( x[i] != 0 )
704
                        return(1);
705
                }
706
        }
707
return(0);
708
}
709
 
710
#if LDBL_MANT_DIG == 64
711
 
712
/* Return nonzero if internal format number is infinite. */
713
static int
714
eiisinf (x)
715
     unsigned short x[];
716
{
717
 
718
#ifdef NANS
719
  if (eiisnan (x))
720
    return (0);
721
#endif
722
  if ((x[E] & 0x7fff) == 0x7fff)
723
    return (1);
724
  return (0);
725
}
726
#endif /* LDBL_MANT_DIG == 64 */
727
 
728
/*
729
;       Compare significands of numbers in internal format.
730
;       Guard words are included in the comparison.
731
;
732
;       unsigned short a[NI], b[NI];
733
;       cmpm( a, b );
734
;
735
;       for the significands:
736
;       returns +1 if a > b
737
;                0 if a == b
738
;               -1 if a < b
739
*/
740
static int ecmpm(register short unsigned int *a, register short unsigned int *b)
741
{
742
int i;
743
 
744
a += M; /* skip up to significand area */
745
b += M;
746
for( i=M; i<NI; i++ )
747
        {
748
        if( *a++ != *b++ )
749
                goto difrnt;
750
        }
751
return(0);
752
 
753
difrnt:
754
if( *(--a) > *(--b) )
755
        return(1);
756
else
757
        return(-1);
758
}
759
 
760
 
761
/*
762
;       Shift significand down by 1 bit
763
*/
764
 
765
static void eshdn1(register short unsigned int *x)
766
{
767
register unsigned short bits;
768
int i;
769
 
770
x += M; /* point to significand area */
771
 
772
bits = 0;
773
for( i=M; i<NI; i++ )
774
        {
775
        if( *x & 1 )
776
                bits |= 1;
777
        *x >>= 1;
778
        if( bits & 2 )
779
                *x |= 0x8000;
780
        bits <<= 1;
781
        ++x;
782
        }
783
}
784
 
785
 
786
 
787
/*
788
;       Shift significand up by 1 bit
789
*/
790
 
791
static void eshup1(register short unsigned int *x)
792
{
793
register unsigned short bits;
794
int i;
795
 
796
x += NI-1;
797
bits = 0;
798
 
799
for( i=M; i<NI; i++ )
800
        {
801
        if( *x & 0x8000 )
802
                bits |= 1;
803
        *x <<= 1;
804
        if( bits & 2 )
805
                *x |= 1;
806
        bits <<= 1;
807
        --x;
808
        }
809
}
810
 
811
 
812
 
813
/*
814
;       Shift significand down by 8 bits
815
*/
816
 
817
static void eshdn8(register short unsigned int *x)
818
{
819
register unsigned short newbyt, oldbyt;
820
int i;
821
 
822
x += M;
823
oldbyt = 0;
824
for( i=M; i<NI; i++ )
825
        {
826
        newbyt = *x << 8;
827
        *x >>= 8;
828
        *x |= oldbyt;
829
        oldbyt = newbyt;
830
        ++x;
831
        }
832
}
833
 
834
/*
835
;       Shift significand up by 8 bits
836
*/
837
 
838
static void eshup8(register short unsigned int *x)
839
{
840
int i;
841
register unsigned short newbyt, oldbyt;
842
 
843
x += NI-1;
844
oldbyt = 0;
845
 
846
for( i=M; i<NI; i++ )
847
        {
848
        newbyt = *x >> 8;
849
        *x <<= 8;
850
        *x |= oldbyt;
851
        oldbyt = newbyt;
852
        --x;
853
        }
854
}
855
 
856
/*
857
;       Shift significand up by 16 bits
858
*/
859
 
860
static void eshup6(register short unsigned int *x)
861
{
862
int i;
863
register unsigned short *p;
864
 
865
p = x + M;
866
x += M + 1;
867
 
868
for( i=M; i<NI-1; i++ )
869
        *p++ = *x++;
870
 
871
*p = 0;
872
}
873
 
874
/*
875
;       Shift significand down by 16 bits
876
*/
877
 
878
static void eshdn6(register short unsigned int *x)
879
{
880
int i;
881
register unsigned short *p;
882
 
883
x += NI-1;
884
p = x + 1;
885
 
886
for( i=M; i<NI-1; i++ )
887
        *(--p) = *(--x);
888
 
889
*(--p) = 0;
890
}
891
 
892
/*
893
;       Add significands
894
;       x + y replaces y
895
*/
896
 
897
static void eaddm(short unsigned int *x, short unsigned int *y)
898
{
899
register unsigned long a;
900
int i;
901
unsigned int carry;
902
 
903
x += NI-1;
904
y += NI-1;
905
carry = 0;
906
for( i=M; i<NI; i++ )
907
        {
908
        a = (unsigned long )(*x) + (unsigned long )(*y) + carry;
909
        if( a & 0x10000 )
910
                carry = 1;
911
        else
912
                carry = 0;
913
        *y = (unsigned short )a;
914
        --x;
915
        --y;
916
        }
917
}
918
 
919
/*
920
;       Subtract significands
921
;       y - x replaces y
922
*/
923
 
924
static void esubm(short unsigned int *x, short unsigned int *y)
925
{
926
unsigned long a;
927
int i;
928
unsigned int carry;
929
 
930
x += NI-1;
931
y += NI-1;
932
carry = 0;
933
for( i=M; i<NI; i++ )
934
        {
935
        a = (unsigned long )(*y) - (unsigned long )(*x) - carry;
936
        if( a & 0x10000 )
937
                carry = 1;
938
        else
939
                carry = 0;
940
        *y = (unsigned short )a;
941
        --x;
942
        --y;
943
        }
944
}
945
 
946
 
947
/* Divide significands */
948
 
949
 
950
/* Multiply significand of e-type number b
951
by 16-bit quantity a, e-type result to c. */
952
 
953
static void m16m(short unsigned int a, short unsigned int *b, short unsigned int *c)
954
{
955
register unsigned short *pp;
956
register unsigned long carry;
957
unsigned short *ps;
958
unsigned short p[NI];
959
unsigned long aa, m;
960
int i;
961
 
962
aa = a;
963
pp = &p[NI-2];
964
*pp++ = 0;
965
*pp = 0;
966
ps = &b[NI-1];
967
 
968
for( i=M+1; i<NI; i++ )
969
        {
970
        if( *ps == 0 )
971
                {
972
                --ps;
973
                --pp;
974
                *(pp-1) = 0;
975
                }
976
        else
977
                {
978
                m = (unsigned long) aa * *ps--;
979
                carry = (m & 0xffff) + *pp;
980
                *pp-- = (unsigned short )carry;
981
                carry = (carry >> 16) + (m >> 16) + *pp;
982
                *pp = (unsigned short )carry;
983
                *(pp-1) = carry >> 16;
984
                }
985
        }
986
for( i=M; i<NI; i++ )
987
        c[i] = p[i];
988
}
989
 
990
 
991
/* Divide significands. Neither the numerator nor the denominator
992
is permitted to have its high guard word nonzero.  */
993
 
994
 
995
static int edivm(short unsigned int *den, short unsigned int *num, LDPARMS *ldp)
996
{
997
int i;
998
register unsigned short *p;
999
unsigned long tnum;
1000
unsigned short j, tdenm, tquot;
1001
unsigned short tprod[NI+1];
1002
unsigned short *equot = ldp->equot;
1003
 
1004
p = &equot[0];
1005
*p++ = num[0];
1006
*p++ = num[1];
1007
 
1008
for( i=M; i<NI; i++ )
1009
        {
1010
        *p++ = 0;
1011
        }
1012
eshdn1( num );
1013
tdenm = den[M+1];
1014
for( i=M; i<NI; i++ )
1015
        {
1016
        /* Find trial quotient digit (the radix is 65536). */
1017
        tnum = (((unsigned long) num[M]) << 16) + num[M+1];
1018
 
1019
        /* Do not execute the divide instruction if it will overflow. */
1020
        if( (tdenm * 0xffffUL) < tnum )
1021
                tquot = 0xffff;
1022
        else
1023
                tquot = tnum / tdenm;
1024
 
1025
                /* Prove that the divide worked. */
1026
/*
1027
        tcheck = (unsigned long )tquot * tdenm;
1028
        if( tnum - tcheck > tdenm )
1029
                tquot = 0xffff;
1030
*/
1031
        /* Multiply denominator by trial quotient digit. */
1032
        m16m( tquot, den, tprod );
1033
        /* The quotient digit may have been overestimated. */
1034
        if( ecmpm( tprod, num ) > 0 )
1035
                {
1036
                tquot -= 1;
1037
                esubm( den, tprod );
1038
                if( ecmpm( tprod, num ) > 0 )
1039
                        {
1040
                        tquot -= 1;
1041
                        esubm( den, tprod );
1042
                        }
1043
                }
1044
/*
1045
        if( ecmpm( tprod, num ) > 0 )
1046
                {
1047
                eshow( "tprod", tprod );
1048
                eshow( "num  ", num );
1049
                printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
1050
                         tnum, den[M+1], tquot );
1051
                }
1052
*/
1053
        esubm( tprod, num );
1054
/*
1055
        if( ecmpm( num, den ) >= 0 )
1056
                {
1057
                eshow( "num  ", num );
1058
                eshow( "den  ", den );
1059
                printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
1060
                         tnum, den[M+1], tquot );
1061
                }
1062
*/
1063
        equot[i] = tquot;
1064
        eshup6(num);
1065
        }
1066
/* test for nonzero remainder after roundoff bit */
1067
p = &num[M];
1068
j = 0;
1069
for( i=M; i<NI; i++ )
1070
        {
1071
        j |= *p++;
1072
        }
1073
if( j )
1074
        j = 1;
1075
 
1076
for( i=0; i<NI; i++ )
1077
        num[i] = equot[i];
1078
 
1079
return( (int )j );
1080
}
1081
 
1082
 
1083
 
1084
/* Multiply significands */
1085
static int emulm(short unsigned int *a, short unsigned int *b, LDPARMS *ldp)
1086
{
1087
unsigned short *p, *q;
1088
unsigned short pprod[NI];
1089
unsigned short j;
1090
int i;
1091
unsigned short *equot = ldp->equot;
1092
 
1093
equot[0] = b[0];
1094
equot[1] = b[1];
1095
for( i=M; i<NI; i++ )
1096
        equot[i] = 0;
1097
 
1098
j = 0;
1099
p = &a[NI-1];
1100
q = &equot[NI-1];
1101
for( i=M+1; i<NI; i++ )
1102
        {
1103
        if( *p == 0 )
1104
                {
1105
                --p;
1106
                }
1107
        else
1108
                {
1109
                m16m( *p--, b, pprod );
1110
                eaddm(pprod, equot);
1111
                }
1112
        j |= *q;
1113
        eshdn6(equot);
1114
        }
1115
 
1116
for( i=0; i<NI; i++ )
1117
        b[i] = equot[i];
1118
 
1119
/* return flag for lost nonzero bits */
1120
return( (int)j );
1121
}
1122
 
1123
 
1124
/*
1125
static void eshow(str, x)
1126
char *str;
1127
unsigned short *x;
1128
{
1129
int i;
1130
 
1131
printf( "%s ", str );
1132
for( i=0; i<NI; i++ )
1133
        printf( "%04x ", *x++ );
1134
printf( "\n" );
1135
}
1136
*/
1137
 
1138
 
1139
/*
1140
 * Normalize and round off.
1141
 *
1142
 * The internal format number to be rounded is "s".
1143
 * Input "lost" indicates whether the number is exact.
1144
 * This is the so-called sticky bit.
1145
 *
1146
 * Input "subflg" indicates whether the number was obtained
1147
 * by a subtraction operation.  In that case if lost is nonzero
1148
 * then the number is slightly smaller than indicated.
1149
 *
1150
 * Input "exp" is the biased exponent, which may be negative.
1151
 * the exponent field of "s" is ignored but is replaced by
1152
 * "exp" as adjusted by normalization and rounding.
1153
 *
1154
 * Input "rcntrl" is the rounding control.
1155
 */
1156
 
1157
 
1158
static void emdnorm(short unsigned int *s, int lost, int subflg, long int exp, int rcntrl, LDPARMS *ldp)
1159
{
1160
int i, j;
1161
unsigned short r;
1162
 
1163
/* Normalize */
1164
j = enormlz( s );
1165
 
1166
/* a blank significand could mean either zero or infinity. */
1167
#ifndef INFINITY
1168
if( j > NBITS )
1169
        {
1170
        ecleazs( s );
1171
        return;
1172
        }
1173
#endif
1174
exp -= j;
1175
#ifndef INFINITY
1176
if( exp >= 32767L )
1177
        goto overf;
1178
#else
1179
if( (j > NBITS) && (exp < 32767L) )
1180
        {
1181
        ecleazs( s );
1182
        return;
1183
        }
1184
#endif
1185
if( exp < 0L )
1186
        {
1187
        if( exp > (long )(-NBITS-1) )
1188
                {
1189
                j = (int )exp;
1190
                i = eshift( s, j );
1191
                if( i )
1192
                        lost = 1;
1193
                }
1194
        else
1195
                {
1196
                ecleazs( s );
1197
                return;
1198
                }
1199
        }
1200
/* Round off, unless told not to by rcntrl. */
1201
if( rcntrl == 0 )
1202
        goto mdfin;
1203
/* Set up rounding parameters if the control register changed. */
1204
if( ldp->rndprc != ldp->rlast )
1205
        {
1206
        ecleaz( ldp->rbit );
1207
        switch( ldp->rndprc )
1208
                {
1209
                default:
1210
                case NBITS:
1211
                        ldp->rw = NI-1; /* low guard word */
1212
                        ldp->rmsk = 0xffff;
1213
                        ldp->rmbit = 0x8000;
1214
                        ldp->rebit = 1;
1215
                        ldp->re = ldp->rw - 1;
1216
                        break;
1217
                case 113:
1218
                        ldp->rw = 10;
1219
                        ldp->rmsk = 0x7fff;
1220
                        ldp->rmbit = 0x4000;
1221
                        ldp->rebit = 0x8000;
1222
                        ldp->re = ldp->rw;
1223
                        break;
1224
                case 64:
1225
                        ldp->rw = 7;
1226
                        ldp->rmsk = 0xffff;
1227
                        ldp->rmbit = 0x8000;
1228
                        ldp->rebit = 1;
1229
                        ldp->re = ldp->rw-1;
1230
                        break;
1231
/* For DEC arithmetic */
1232
                case 56:
1233
                        ldp->rw = 6;
1234
                        ldp->rmsk = 0xff;
1235
                        ldp->rmbit = 0x80;
1236
                        ldp->rebit = 0x100;
1237
                        ldp->re = ldp->rw;
1238
                        break;
1239
                case 53:
1240
                        ldp->rw = 6;
1241
                        ldp->rmsk = 0x7ff;
1242
                        ldp->rmbit = 0x0400;
1243
                        ldp->rebit = 0x800;
1244
                        ldp->re = ldp->rw;
1245
                        break;
1246
                case 24:
1247
                        ldp->rw = 4;
1248
                        ldp->rmsk = 0xff;
1249
                        ldp->rmbit = 0x80;
1250
                        ldp->rebit = 0x100;
1251
                        ldp->re = ldp->rw;
1252
                        break;
1253
                }
1254
        ldp->rbit[ldp->re] = ldp->rebit;
1255
        ldp->rlast = ldp->rndprc;
1256
        }
1257
 
1258
/* Shift down 1 temporarily if the data structure has an implied
1259
 * most significant bit and the number is denormal.
1260
 * For rndprc = 64 or NBITS, there is no implied bit.
1261
 * But Intel long double denormals lose one bit of significance even so.
1262
 */
1263
#if IBMPC
1264
if( (exp <= 0) && (ldp->rndprc != NBITS) )
1265
#else
1266
if( (exp <= 0) && (ldp->rndprc != 64) && (ldp->rndprc != NBITS) )
1267
#endif
1268
        {
1269
        lost |= s[NI-1] & 1;
1270
        eshdn1(s);
1271
        }
1272
/* Clear out all bits below the rounding bit,
1273
 * remembering in r if any were nonzero.
1274
 */
1275
r = s[ldp->rw] & ldp->rmsk;
1276
if( ldp->rndprc < NBITS )
1277
        {
1278
        i = ldp->rw + 1;
1279
        while( i < NI )
1280
                {
1281
                if( s[i] )
1282
                        r |= 1;
1283
                s[i] = 0;
1284
                ++i;
1285
                }
1286
        }
1287
s[ldp->rw] &= ~ldp->rmsk;
1288
if( (r & ldp->rmbit) != 0 )
1289
        {
1290
        if( r == ldp->rmbit )
1291
                {
1292
                if( lost == 0 )
1293
                        { /* round to even */
1294
                        if( (s[ldp->re] & ldp->rebit) == 0 )
1295
                                goto mddone;
1296
                        }
1297
                else
1298
                        {
1299
                        if( subflg != 0 )
1300
                                goto mddone;
1301
                        }
1302
                }
1303
        eaddm( ldp->rbit, s );
1304
        }
1305
mddone:
1306
#if IBMPC
1307
if( (exp <= 0) && (ldp->rndprc != NBITS) )
1308
#else
1309
if( (exp <= 0) && (ldp->rndprc != 64) && (ldp->rndprc != NBITS) )
1310
#endif
1311
        {
1312
        eshup1(s);
1313
        }
1314
if( s[2] != 0 )
1315
        { /* overflow on roundoff */
1316
        eshdn1(s);
1317
        exp += 1;
1318
        }
1319
mdfin:
1320
s[NI-1] = 0;
1321
if( exp >= 32767L )
1322
        {
1323
#ifndef INFINITY
1324
overf:
1325
#endif
1326
#ifdef INFINITY
1327
        s[1] = 32767;
1328
        for( i=2; i<NI-1; i++ )
1329
                s[i] = 0;
1330
#else
1331
        s[1] = 32766;
1332
        s[2] = 0;
1333
        for( i=M+1; i<NI-1; i++ )
1334
                s[i] = 0xffff;
1335
        s[NI-1] = 0;
1336
        if( (ldp->rndprc < 64) || (ldp->rndprc == 113) )
1337
                {
1338
                s[ldp->rw] &= ~ldp->rmsk;
1339
                if( ldp->rndprc == 24 )
1340
                        {
1341
                        s[5] = 0;
1342
                        s[6] = 0;
1343
                        }
1344
                }
1345
#endif
1346
        return;
1347
        }
1348
if( exp < 0 )
1349
        s[1] = 0;
1350
else
1351
        s[1] = (unsigned short )exp;
1352
}
1353
 
1354
 
1355
 
1356
/*
1357
;       Subtract external format numbers.
1358
;
1359
;       unsigned short a[NE], b[NE], c[NE];
1360
;       LDPARMS *ldp;
1361
;       esub( a, b, c, ldp );    c = b - a
1362
*/
1363
 
1364
static void esub(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp)
1365
{
1366
 
1367
#ifdef NANS
1368
if( eisnan(a) )
1369
        {
1370
        emov (a, c);
1371
        return;
1372
        }
1373
if( eisnan(b) )
1374
        {
1375
        emov(b,c);
1376
        return;
1377
        }
1378
/* Infinity minus infinity is a NaN.
1379
 * Test for subtracting infinities of the same sign.
1380
 */
1381
if( eisinf(a) && eisinf(b) && ((eisneg (a) ^ eisneg (b)) == 0))
1382
        {
1383
        mtherr( "esub", DOMAIN );
1384
        enan( c, NBITS );
1385
        return;
1386
        }
1387
#endif
1388
eadd1( a, b, c, 1, ldp );
1389
}
1390
 
1391
 
1392
 
1393
static void eadd1(short unsigned int *a, short unsigned int *b, short unsigned int *c, int subflg, LDPARMS *ldp)
1394
{
1395
unsigned short ai[NI], bi[NI], ci[NI];
1396
int i, lost, j, k;
1397
long lt, lta, ltb;
1398
 
1399
#ifdef INFINITY
1400
if( eisinf(a) )
1401
        {
1402
        emov(a,c);
1403
        if( subflg )
1404
                eneg(c);
1405
        return;
1406
        }
1407
if( eisinf(b) )
1408
        {
1409
        emov(b,c);
1410
        return;
1411
        }
1412
#endif
1413
emovi( a, ai );
1414
emovi( b, bi );
1415
if( subflg )
1416
        ai[0] = ~ai[0];
1417
 
1418
/* compare exponents */
1419
lta = ai[E];
1420
ltb = bi[E];
1421
lt = lta - ltb;
1422
if( lt > 0L )
1423
        {       /* put the larger number in bi */
1424
        emovz( bi, ci );
1425
        emovz( ai, bi );
1426
        emovz( ci, ai );
1427
        ltb = bi[E];
1428
        lt = -lt;
1429
        }
1430
lost = 0;
1431
if( lt != 0L )
1432
        {
1433
        if( lt < (long )(-NBITS-1) )
1434
                goto done;      /* answer same as larger addend */
1435
        k = (int )lt;
1436
        lost = eshift( ai, k ); /* shift the smaller number down */
1437
        }
1438
else
1439
        {
1440
/* exponents were the same, so must compare significands */
1441
        i = ecmpm( ai, bi );
1442
        if( i == 0 )
1443
                { /* the numbers are identical in magnitude */
1444
                /* if different signs, result is zero */
1445
                if( ai[0] != bi[0] )
1446
                        {
1447
                        eclear(c);
1448
                        return;
1449
                        }
1450
                /* if same sign, result is double */
1451
                /* double denomalized tiny number */
1452
                if( (bi[E] == 0) && ((bi[3] & 0x8000) == 0) )
1453
                        {
1454
                        eshup1( bi );
1455
                        goto done;
1456
                        }
1457
                /* add 1 to exponent unless both are zero! */
1458
                for( j=1; j<NI-1; j++ )
1459
                        {
1460
                        if( bi[j] != 0 )
1461
                                {
1462
/* This could overflow, but let emovo take care of that. */
1463
                                ltb += 1;
1464
                                break;
1465
                                }
1466
                        }
1467
                bi[E] = (unsigned short )ltb;
1468
                goto done;
1469
                }
1470
        if( i > 0 )
1471
                {       /* put the larger number in bi */
1472
                emovz( bi, ci );
1473
                emovz( ai, bi );
1474
                emovz( ci, ai );
1475
                }
1476
        }
1477
if( ai[0] == bi[0] )
1478
        {
1479
        eaddm( ai, bi );
1480
        subflg = 0;
1481
        }
1482
else
1483
        {
1484
        esubm( ai, bi );
1485
        subflg = 1;
1486
        }
1487
emdnorm( bi, lost, subflg, ltb, 64, ldp );
1488
 
1489
done:
1490
emovo( bi, c, ldp );
1491
}
1492
 
1493
 
1494
 
1495
/*
1496
;       Divide.
1497
;
1498
;       unsigned short a[NE], b[NE], c[NE];
1499
;       LDPARMS *ldp;
1500
;       ediv( a, b, c, ldp );   c = b / a
1501
*/
1502
static void ediv(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp)
1503
{
1504
unsigned short ai[NI], bi[NI];
1505
int i;
1506
long lt, lta, ltb;
1507
 
1508
#ifdef NANS
1509
/* Return any NaN input. */
1510
if( eisnan(a) )
1511
        {
1512
        emov(a,c);
1513
        return;
1514
        }
1515
if( eisnan(b) )
1516
        {
1517
        emov(b,c);
1518
        return;
1519
        }
1520
/* Zero over zero, or infinity over infinity, is a NaN. */
1521
if( ((ecmp(a,ezero) == 0) && (ecmp(b,ezero) == 0))
1522
        || (eisinf (a) && eisinf (b)) )
1523
        {
1524
        mtherr( "ediv", DOMAIN );
1525
        enan( c, NBITS );
1526
        return;
1527
        }
1528
#endif
1529
/* Infinity over anything else is infinity. */
1530
#ifdef INFINITY
1531
if( eisinf(b) )
1532
        {
1533
        if( eisneg(a) ^ eisneg(b) )
1534
                *(c+(NE-1)) = 0x8000;
1535
        else
1536
                *(c+(NE-1)) = 0;
1537
        einfin(c, ldp);
1538
        return;
1539
        }
1540
if( eisinf(a) )
1541
        {
1542
        eclear(c);
1543
        return;
1544
        }
1545
#endif
1546
emovi( a, ai );
1547
emovi( b, bi );
1548
lta = ai[E];
1549
ltb = bi[E];
1550
if( bi[E] == 0 )
1551
        { /* See if numerator is zero. */
1552
        for( i=1; i<NI-1; i++ )
1553
                {
1554
                if( bi[i] != 0 )
1555
                        {
1556
                        ltb -= enormlz( bi );
1557
                        goto dnzro1;
1558
                        }
1559
                }
1560
        eclear(c);
1561
        return;
1562
        }
1563
dnzro1:
1564
 
1565
if( ai[E] == 0 )
1566
        {       /* possible divide by zero */
1567
        for( i=1; i<NI-1; i++ )
1568
                {
1569
                if( ai[i] != 0 )
1570
                        {
1571
                        lta -= enormlz( ai );
1572
                        goto dnzro2;
1573
                        }
1574
                }
1575
        if( ai[0] == bi[0] )
1576
                *(c+(NE-1)) = 0;
1577
        else
1578
                *(c+(NE-1)) = 0x8000;
1579
        einfin(c, ldp);
1580
        mtherr( "ediv", SING );
1581
        return;
1582
        }
1583
dnzro2:
1584
 
1585
i = edivm( ai, bi, ldp );
1586
/* calculate exponent */
1587
lt = ltb - lta + EXONE;
1588
emdnorm( bi, i, 0, lt, 64, ldp );
1589
/* set the sign */
1590
if( ai[0] == bi[0] )
1591
        bi[0] = 0;
1592
else
1593
        bi[0] = 0Xffff;
1594
emovo( bi, c, ldp );
1595
}
1596
 
1597
 
1598
 
1599
/*
1600
;       Multiply.
1601
;
1602
;       unsigned short a[NE], b[NE], c[NE];
1603
;       LDPARMS *ldp
1604
;       emul( a, b, c, ldp );   c = b * a
1605
*/
1606
static void emul(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp)
1607
{
1608
unsigned short ai[NI], bi[NI];
1609
int i, j;
1610
long lt, lta, ltb;
1611
 
1612
#ifdef NANS
1613
/* NaN times anything is the same NaN. */
1614
if( eisnan(a) )
1615
        {
1616
        emov(a,c);
1617
        return;
1618
        }
1619
if( eisnan(b) )
1620
        {
1621
        emov(b,c);
1622
        return;
1623
        }
1624
/* Zero times infinity is a NaN. */
1625
if( (eisinf(a) && (ecmp(b,ezero) == 0))
1626
        || (eisinf(b) && (ecmp(a,ezero) == 0)) )
1627
        {
1628
        mtherr( "emul", DOMAIN );
1629
        enan( c, NBITS );
1630
        return;
1631
        }
1632
#endif
1633
/* Infinity times anything else is infinity. */
1634
#ifdef INFINITY
1635
if( eisinf(a) || eisinf(b) )
1636
        {
1637
        if( eisneg(a) ^ eisneg(b) )
1638
                *(c+(NE-1)) = 0x8000;
1639
        else
1640
                *(c+(NE-1)) = 0;
1641
        einfin(c, ldp);
1642
        return;
1643
        }
1644
#endif
1645
emovi( a, ai );
1646
emovi( b, bi );
1647
lta = ai[E];
1648
ltb = bi[E];
1649
if( ai[E] == 0 )
1650
        {
1651
        for( i=1; i<NI-1; i++ )
1652
                {
1653
                if( ai[i] != 0 )
1654
                        {
1655
                        lta -= enormlz( ai );
1656
                        goto mnzer1;
1657
                        }
1658
                }
1659
        eclear(c);
1660
        return;
1661
        }
1662
mnzer1:
1663
 
1664
if( bi[E] == 0 )
1665
        {
1666
        for( i=1; i<NI-1; i++ )
1667
                {
1668
                if( bi[i] != 0 )
1669
                        {
1670
                        ltb -= enormlz( bi );
1671
                        goto mnzer2;
1672
                        }
1673
                }
1674
        eclear(c);
1675
        return;
1676
        }
1677
mnzer2:
1678
 
1679
/* Multiply significands */
1680
j = emulm( ai, bi, ldp );
1681
/* calculate exponent */
1682
lt = lta + ltb - (EXONE - 1);
1683
emdnorm( bi, j, 0, lt, 64, ldp );
1684
/* calculate sign of product */
1685
if( ai[0] == bi[0] )
1686
        bi[0] = 0;
1687
else
1688
        bi[0] = 0xffff;
1689
emovo( bi, c, ldp );
1690
}
1691
 
1692
 
1693
 
1694
#if LDBL_MANT_DIG > 64
1695
static void e113toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp)
1696
{
1697
register unsigned short r;
1698
unsigned short *e, *p;
1699
unsigned short yy[NI];
1700
int denorm, i;
1701
 
1702
e = pe;
1703
denorm = 0;
1704
ecleaz(yy);
1705
#ifdef IBMPC
1706
e += 7;
1707
#endif
1708
r = *e;
1709
yy[0] = 0;
1710
if( r & 0x8000 )
1711
        yy[0] = 0xffff;
1712
r &= 0x7fff;
1713
#ifdef INFINITY
1714
if( r == 0x7fff )
1715
        {
1716
#ifdef NANS
1717
#ifdef IBMPC
1718
        for( i=0; i<7; i++ )
1719
                {
1720
                if( pe[i] != 0 )
1721
                        {
1722
                        enan( y, NBITS );
1723
                        return;
1724
                        }
1725
                }
1726
#else  /* !IBMPC */
1727
        for( i=1; i<8; i++ )
1728
                {
1729
                if( pe[i] != 0 )
1730
                        {
1731
                        enan( y, NBITS );
1732
                        return;
1733
                        }
1734
                }
1735
#endif /* !IBMPC */
1736
#endif /* NANS */
1737
        eclear( y );
1738
        einfin( y, ldp );
1739
        if( *e & 0x8000 )
1740
                eneg(y);
1741
        return;
1742
        }
1743
#endif  /* INFINITY */
1744
yy[E] = r;
1745
p = &yy[M + 1];
1746
#ifdef IBMPC
1747
for( i=0; i<7; i++ )
1748
        *p++ = *(--e);
1749
#else  /* IBMPC */
1750
++e;
1751
for( i=0; i<7; i++ )
1752
        *p++ = *e++;
1753
#endif /* IBMPC */ 
1754
/* If denormal, remove the implied bit; else shift down 1. */
1755
if( r == 0 )
1756
        {
1757
        yy[M] = 0;
1758
        }
1759
else
1760
        {
1761
        yy[M] = 1;
1762
        eshift( yy, -1 );
1763
        }
1764
emovo(yy,y,ldp);
1765
}
1766
 
1767
/* move out internal format to ieee long double */
1768
static void toe113(short unsigned int *a, short unsigned int *b)
1769
{
1770
register unsigned short *p, *q;
1771
unsigned short i;
1772
 
1773
#ifdef NANS
1774
if( eiisnan(a) )
1775
        {
1776
        enan( b, 113 );
1777
        return;
1778
        }
1779
#endif
1780
p = a;
1781
#ifdef MIEEE
1782
q = b;
1783
#else
1784
q = b + 7;                      /* point to output exponent */
1785
#endif
1786
 
1787
/* If not denormal, delete the implied bit. */
1788
if( a[E] != 0 )
1789
        {
1790
        eshup1 (a);
1791
        }
1792
/* combine sign and exponent */
1793
i = *p++;
1794
#ifdef MIEEE
1795
if( i )
1796
        *q++ = *p++ | 0x8000;
1797
else
1798
        *q++ = *p++;
1799
#else
1800
if( i )
1801
        *q-- = *p++ | 0x8000;
1802
else
1803
        *q-- = *p++;
1804
#endif
1805
/* skip over guard word */
1806
++p;
1807
/* move the significand */
1808
#ifdef MIEEE
1809
for (i = 0; i < 7; i++)
1810
        *q++ = *p++;
1811
#else
1812
for (i = 0; i < 7; i++)
1813
        *q-- = *p++;
1814
#endif
1815
}
1816
#endif /* LDBL_MANT_DIG > 64 */
1817
 
1818
 
1819
#if LDBL_MANT_DIG == 64
1820
static void e64toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp)
1821
{
1822
unsigned short yy[NI];
1823
unsigned short *p, *q, *e;
1824
int i;
1825
 
1826
e = pe;
1827
p = yy;
1828
 
1829
for( i=0; i<NE-5; i++ )
1830
        *p++ = 0;
1831
#ifdef IBMPC
1832
for( i=0; i<5; i++ )
1833
        *p++ = *e++;
1834
#endif
1835
#ifdef DEC
1836
for( i=0; i<5; i++ )
1837
        *p++ = *e++;
1838
#endif
1839
#ifdef MIEEE
1840
p = &yy[0] + (NE-1);
1841
*p-- = *e++;
1842
++e;  /* MIEEE skips over 2nd short */
1843
for( i=0; i<4; i++ )
1844
        *p-- = *e++;
1845
#endif
1846
 
1847
#ifdef IBMPC
1848
/* For Intel long double, shift denormal significand up 1
1849
   -- but only if the top significand bit is zero.  */
1850
if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
1851
  {
1852
    unsigned short temp[NI+1];
1853
    emovi(yy, temp);
1854
    eshup1(temp);
1855
    emovo(temp,y,ldp);
1856
    return;
1857
  }
1858
#endif
1859
#ifdef INFINITY
1860
/* Point to the exponent field.  */
1861
p = &yy[NE-1];
1862
if( *p == 0x7fff )
1863
        {
1864
#ifdef NANS
1865
#ifdef IBMPC
1866
        for( i=0; i<4; i++ )
1867
                {
1868
                if((i != 3 && pe[i] != 0)
1869
                   /* Check for Intel long double infinity pattern.  */
1870
                   || (i == 3 && pe[i] != 0x8000))
1871
                        {
1872
                        enan( y, NBITS );
1873
                        return;
1874
                        }
1875
                }
1876
#endif
1877
#ifdef MIEEE
1878
        for( i=2; i<=5; i++ )
1879
                {
1880
                if( pe[i] != 0 )
1881
                        {
1882
                        enan( y, NBITS );
1883
                        return;
1884
                        }
1885
                }
1886
#endif
1887
#endif /* NANS */
1888
        eclear( y );
1889
        einfin( y, ldp );
1890
        if( *p & 0x8000 )
1891
                eneg(y);
1892
        return;
1893
        }
1894
#endif /* INFINITY */
1895
p = yy;
1896
q = y;
1897
for( i=0; i<NE; i++ )
1898
        *q++ = *p++;
1899
}
1900
 
1901
/* move out internal format to ieee long double */
1902
static void toe64(short unsigned int *a, short unsigned int *b)
1903
{
1904
register unsigned short *p, *q;
1905
unsigned short i;
1906
 
1907
#ifdef NANS
1908
if( eiisnan(a) )
1909
        {
1910
        enan( b, 64 );
1911
        return;
1912
        }
1913
#endif
1914
#ifdef IBMPC
1915
/* Shift Intel denormal significand down 1.  */
1916
if( a[E] == 0 )
1917
  eshdn1(a);
1918
#endif
1919
p = a;
1920
#ifdef MIEEE
1921
q = b;
1922
#else
1923
q = b + 4; /* point to output exponent */
1924
/* NOTE: Intel data type is 96 bits wide, clear the last word here. */
1925
*(q+1)= 0;
1926
#endif
1927
 
1928
/* combine sign and exponent */
1929
i = *p++;
1930
#ifdef MIEEE
1931
if( i )
1932
        *q++ = *p++ | 0x8000;
1933
else
1934
        *q++ = *p++;
1935
*q++ = 0; /* leave 2nd short blank */
1936
#else
1937
if( i )
1938
        *q-- = *p++ | 0x8000;
1939
else
1940
        *q-- = *p++;
1941
#endif
1942
/* skip over guard word */
1943
++p;
1944
/* move the significand */
1945
#ifdef MIEEE
1946
for( i=0; i<4; i++ )
1947
        *q++ = *p++;
1948
#else
1949
#ifdef INFINITY
1950
#ifdef IBMPC
1951
if (eiisinf (a))
1952
        {
1953
        /* Intel long double infinity.  */
1954
        *q-- = 0x8000;
1955
        *q-- = 0;
1956
        *q-- = 0;
1957
        *q = 0;
1958
        return;
1959
        }
1960
#endif /* IBMPC */
1961
#endif /* INFINITY */
1962
for( i=0; i<4; i++ )
1963
        *q-- = *p++;
1964
#endif
1965
}
1966
 
1967
#endif /* LDBL_MANT_DIG == 64 */
1968
 
1969
#if LDBL_MANT_DIG == 53
1970
/*
1971
; Convert IEEE double precision to e type
1972
;       double d;
1973
;       unsigned short x[N+2];
1974
;       e53toe( &d, x );
1975
*/
1976
static void e53toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp)
1977
{
1978
#ifdef DEC
1979
 
1980
dectoe( pe, y ); /* see etodec.c */
1981
 
1982
#else
1983
 
1984
register unsigned short r;
1985
register unsigned short *p, *e;
1986
unsigned short yy[NI];
1987
int denorm, k;
1988
 
1989
e = pe;
1990
denorm = 0;      /* flag if denormalized number */
1991
ecleaz(yy);
1992
#ifdef IBMPC
1993
e += 3;
1994
#endif
1995
#ifdef DEC
1996
e += 3;
1997
#endif 
1998
r = *e;
1999
yy[0] = 0;
2000
if( r & 0x8000 )
2001
        yy[0] = 0xffff;
2002
yy[M] = (r & 0x0f) | 0x10;
2003
r &= ~0x800f;   /* strip sign and 4 significand bits */
2004
#ifdef INFINITY
2005
if( r == 0x7ff0 )
2006
        {
2007
#ifdef NANS
2008
#ifdef IBMPC
2009
        if( ((pe[3] & 0xf) != 0) || (pe[2] != 0)
2010
                || (pe[1] != 0) || (pe[0] != 0) )
2011
                {
2012
                enan( y, NBITS );
2013
                return;
2014
                }
2015
#else  /* !IBMPC */
2016
        if( ((pe[0] & 0xf) != 0) || (pe[1] != 0)
2017
                 || (pe[2] != 0) || (pe[3] != 0) )
2018
                {
2019
                enan( y, NBITS );
2020
                return;
2021
                }
2022
#endif /* !IBMPC */
2023
#endif  /* NANS */
2024
        eclear( y );
2025
        einfin( y, ldp );
2026
        if( yy[0] )
2027
                eneg(y);
2028
        return;
2029
        }
2030
#endif
2031
r >>= 4;
2032
/* If zero exponent, then the significand is denormalized.
2033
 * So, take back the understood high significand bit. */
2034
if( r == 0 )
2035
        {
2036
        denorm = 1;
2037
        yy[M] &= ~0x10;
2038
        }
2039
r += EXONE - 01777;
2040
yy[E] = r;
2041
p = &yy[M+1];
2042
#ifdef IBMPC
2043
*p++ = *(--e);
2044
*p++ = *(--e);
2045
*p++ = *(--e);
2046
#else  /* !IBMPC */
2047
++e;
2048
*p++ = *e++;
2049
*p++ = *e++;
2050
*p++ = *e++;
2051
#endif /* !IBMPC */
2052
(void )eshift( yy, -5 );
2053
if( denorm )
2054
        { /* if zero exponent, then normalize the significand */
2055
        if( (k = enormlz(yy)) > NBITS )
2056
                ecleazs(yy);
2057
        else
2058
                yy[E] -= (unsigned short )(k-1);
2059
        }
2060
emovo( yy, y, ldp );
2061
#endif /* !DEC */
2062
}
2063
 
2064
/*
2065
; e type to IEEE double precision
2066
;       double d;
2067
;       unsigned short x[NE];
2068
;       etoe53( x, &d );
2069
*/
2070
 
2071
#ifdef DEC
2072
 
2073
static void etoe53( x, e )
2074
unsigned short *x, *e;
2075
{
2076
etodec( x, e ); /* see etodec.c */
2077
}
2078
 
2079
static void toe53( x, y )
2080
unsigned short *x, *y;
2081
{
2082
todec( x, y );
2083
}
2084
 
2085
#else
2086
 
2087
static void toe53(short unsigned int *x, short unsigned int *y)
2088
{
2089
unsigned short i;
2090
unsigned short *p;
2091
 
2092
 
2093
#ifdef NANS
2094
if( eiisnan(x) )
2095
        {
2096
        enan( y, 53 );
2097
        return;
2098
        }
2099
#endif
2100
p = &x[0];
2101
#ifdef IBMPC
2102
y += 3;
2103
#endif
2104
#ifdef DEC
2105
y += 3;
2106
#endif
2107
*y = 0;  /* output high order */
2108
if( *p++ )
2109
        *y = 0x8000;    /* output sign bit */
2110
 
2111
i = *p++;
2112
if( i >= (unsigned int )2047 )
2113
        {       /* Saturate at largest number less than infinity. */
2114
#ifdef INFINITY
2115
        *y |= 0x7ff0;
2116
#ifdef IBMPC
2117
        *(--y) = 0;
2118
        *(--y) = 0;
2119
        *(--y) = 0;
2120
#else /* !IBMPC */
2121
        ++y;
2122
        *y++ = 0;
2123
        *y++ = 0;
2124
        *y++ = 0;
2125
#endif /* IBMPC */
2126
#else /* !INFINITY */
2127
        *y |= (unsigned short )0x7fef;
2128
#ifdef IBMPC
2129
        *(--y) = 0xffff;
2130
        *(--y) = 0xffff;
2131
        *(--y) = 0xffff;
2132
#else /* !IBMPC */
2133
        ++y;
2134
        *y++ = 0xffff;
2135
        *y++ = 0xffff;
2136
        *y++ = 0xffff;
2137
#endif
2138
#endif /* !INFINITY */
2139
        return;
2140
        }
2141
if( i == 0 )
2142
        {
2143
        (void )eshift( x, 4 );
2144
        }
2145
else
2146
        {
2147
        i <<= 4;
2148
        (void )eshift( x, 5 );
2149
        }
2150
i |= *p++ & (unsigned short )0x0f;      /* *p = xi[M] */
2151
*y |= (unsigned short )i; /* high order output already has sign bit set */
2152
#ifdef IBMPC
2153
*(--y) = *p++;
2154
*(--y) = *p++;
2155
*(--y) = *p;
2156
#else /* !IBMPC */
2157
++y;
2158
*y++ = *p++;
2159
*y++ = *p++;
2160
*y++ = *p++;
2161
#endif /* !IBMPC */
2162
}
2163
 
2164
#endif /* not DEC */
2165
#endif /* LDBL_MANT_DIG == 53 */
2166
 
2167
#if LDBL_MANT_DIG == 24
2168
/*
2169
; Convert IEEE single precision to e type
2170
;       float d;
2171
;       unsigned short x[N+2];
2172
;       dtox( &d, x );
2173
*/
2174
void e24toe( short unsigned int *pe, short unsigned int *y, LDPARMS *ldp )
2175
{
2176
register unsigned short r;
2177
register unsigned short *p, *e;
2178
unsigned short yy[NI];
2179
int denorm, k;
2180
 
2181
e = pe;
2182
denorm = 0;      /* flag if denormalized number */
2183
ecleaz(yy);
2184
#ifdef IBMPC
2185
e += 1;
2186
#endif
2187
#ifdef DEC
2188
e += 1;
2189
#endif
2190
r = *e;
2191
yy[0] = 0;
2192
if( r & 0x8000 )
2193
        yy[0] = 0xffff;
2194
yy[M] = (r & 0x7f) | 0200;
2195
r &= ~0x807f;   /* strip sign and 7 significand bits */
2196
#ifdef INFINITY
2197
if( r == 0x7f80 )
2198
        {
2199
#ifdef NANS
2200
#ifdef MIEEE
2201
        if( ((pe[0] & 0x7f) != 0) || (pe[1] != 0) )
2202
                {
2203
                enan( y, NBITS );
2204
                return;
2205
                }
2206
#else  /* !MIEEE */
2207
        if( ((pe[1] & 0x7f) != 0) || (pe[0] != 0) )
2208
                {
2209
                enan( y, NBITS );
2210
                return;
2211
                }
2212
#endif /* !MIEEE */
2213
#endif  /* NANS */
2214
        eclear( y );
2215
        einfin( y, ldp );
2216
        if( yy[0] )
2217
                eneg(y);
2218
        return;
2219
        }
2220
#endif
2221
r >>= 7;
2222
/* If zero exponent, then the significand is denormalized.
2223
 * So, take back the understood high significand bit. */
2224
if( r == 0 )
2225
        {
2226
        denorm = 1;
2227
        yy[M] &= ~0200;
2228
        }
2229
r += EXONE - 0177;
2230
yy[E] = r;
2231
p = &yy[M+1];
2232
#ifdef IBMPC
2233
*p++ = *(--e);
2234
#endif
2235
#ifdef DEC
2236
*p++ = *(--e);
2237
#endif
2238
#ifdef MIEEE
2239
++e;
2240
*p++ = *e++;
2241
#endif
2242
(void )eshift( yy, -8 );
2243
if( denorm )
2244
        { /* if zero exponent, then normalize the significand */
2245
        if( (k = enormlz(yy)) > NBITS )
2246
                ecleazs(yy);
2247
        else
2248
                yy[E] -= (unsigned short )(k-1);
2249
        }
2250
emovo( yy, y, ldp );
2251
}
2252
 
2253
static void toe24(short unsigned int *x, short unsigned int *y)
2254
{
2255
unsigned short i;
2256
unsigned short *p;
2257
 
2258
#ifdef NANS
2259
if( eiisnan(x) )
2260
        {
2261
        enan( y, 24 );
2262
        return;
2263
        }
2264
#endif
2265
p = &x[0];
2266
#ifdef IBMPC
2267
y += 1;
2268
#endif
2269
#ifdef DEC
2270
y += 1;
2271
#endif
2272
*y = 0;  /* output high order */
2273
if( *p++ )
2274
        *y = 0x8000;    /* output sign bit */
2275
 
2276
i = *p++;
2277
if( i >= 255 )
2278
        {       /* Saturate at largest number less than infinity. */
2279
#ifdef INFINITY
2280
        *y |= (unsigned short )0x7f80;
2281
#ifdef IBMPC
2282
        *(--y) = 0;
2283
#endif
2284
#ifdef DEC
2285
        *(--y) = 0;
2286
#endif
2287
#ifdef MIEEE
2288
        ++y;
2289
        *y = 0;
2290
#endif
2291
#else /* !INFINITY */
2292
        *y |= (unsigned short )0x7f7f;
2293
#ifdef IBMPC
2294
        *(--y) = 0xffff;
2295
#endif
2296
#ifdef DEC
2297
        *(--y) = 0xffff;
2298
#endif
2299
#ifdef MIEEE
2300
        ++y;
2301
        *y = 0xffff;
2302
#endif
2303
#endif /* !INFINITY */
2304
        return;
2305
        }
2306
if( i == 0 )
2307
        {
2308
        (void )eshift( x, 7 );
2309
        }
2310
else
2311
        {
2312
        i <<= 7;
2313
        (void )eshift( x, 8 );
2314
        }
2315
i |= *p++ & (unsigned short )0x7f;      /* *p = xi[M] */
2316
*y |= i;        /* high order output already has sign bit set */
2317
#ifdef IBMPC
2318
*(--y) = *p;
2319
#endif
2320
#ifdef DEC
2321
*(--y) = *p;
2322
#endif
2323
#ifdef MIEEE
2324
++y;
2325
*y = *p;
2326
#endif
2327
}
2328
#endif /* LDBL_MANT_DIG == 24 */
2329
 
2330
/* Compare two e type numbers.
2331
 *
2332
 * unsigned short a[NE], b[NE];
2333
 * ecmp( a, b );
2334
 *
2335
 *  returns +1 if a > b
2336
 *           0 if a == b
2337
 *          -1 if a < b
2338
 *          -2 if either a or b is a NaN.
2339
 */
2340
static int ecmp(short unsigned int *a, short unsigned int *b)
2341
{
2342
unsigned short ai[NI], bi[NI];
2343
register unsigned short *p, *q;
2344
register int i;
2345
int msign;
2346
 
2347
#ifdef NANS
2348
if (eisnan (a)  || eisnan (b))
2349
        return( -2 );
2350
#endif
2351
emovi( a, ai );
2352
p = ai;
2353
emovi( b, bi );
2354
q = bi;
2355
 
2356
if( *p != *q )
2357
        { /* the signs are different */
2358
/* -0 equals + 0 */
2359
        for( i=1; i<NI-1; i++ )
2360
                {
2361
                if( ai[i] != 0 )
2362
                        goto nzro;
2363
                if( bi[i] != 0 )
2364
                        goto nzro;
2365
                }
2366
        return(0);
2367
nzro:
2368
        if( *p == 0 )
2369
                return( 1 );
2370
        else
2371
                return( -1 );
2372
        }
2373
/* both are the same sign */
2374
if( *p == 0 )
2375
        msign = 1;
2376
else
2377
        msign = -1;
2378
i = NI-1;
2379
do
2380
        {
2381
        if( *p++ != *q++ )
2382
                {
2383
                goto diff;
2384
                }
2385
        }
2386
while( --i > 0 );
2387
 
2388
return(0);       /* equality */
2389
 
2390
 
2391
 
2392
diff:
2393
 
2394
if( *(--p) > *(--q) )
2395
        return( msign );                /* p is bigger */
2396
else
2397
        return( -msign );       /* p is littler */
2398
}
2399
 
2400
 
2401
/*
2402
;       Shift significand
2403
;
2404
;       Shifts significand area up or down by the number of bits
2405
;       given by the variable sc.
2406
*/
2407
static int eshift(short unsigned int *x, int sc)
2408
{
2409
unsigned short lost;
2410
unsigned short *p;
2411
 
2412
if( sc == 0 )
2413
        return( 0 );
2414
 
2415
lost = 0;
2416
p = x + NI-1;
2417
 
2418
if( sc < 0 )
2419
        {
2420
        sc = -sc;
2421
        while( sc >= 16 )
2422
                {
2423
                lost |= *p;     /* remember lost bits */
2424
                eshdn6(x);
2425
                sc -= 16;
2426
                }
2427
 
2428
        while( sc >= 8 )
2429
                {
2430
                lost |= *p & 0xff;
2431
                eshdn8(x);
2432
                sc -= 8;
2433
                }
2434
 
2435
        while( sc > 0 )
2436
                {
2437
                lost |= *p & 1;
2438
                eshdn1(x);
2439
                sc -= 1;
2440
                }
2441
        }
2442
else
2443
        {
2444
        while( sc >= 16 )
2445
                {
2446
                eshup6(x);
2447
                sc -= 16;
2448
                }
2449
 
2450
        while( sc >= 8 )
2451
                {
2452
                eshup8(x);
2453
                sc -= 8;
2454
                }
2455
 
2456
        while( sc > 0 )
2457
                {
2458
                eshup1(x);
2459
                sc -= 1;
2460
                }
2461
        }
2462
if( lost )
2463
        lost = 1;
2464
return( (int )lost );
2465
}
2466
 
2467
 
2468
 
2469
/*
2470
;       normalize
2471
;
2472
; Shift normalizes the significand area pointed to by argument
2473
; shift count (up = positive) is returned.
2474
*/
2475
static int enormlz(short unsigned int *x)
2476
{
2477
register unsigned short *p;
2478
int sc;
2479
 
2480
sc = 0;
2481
p = &x[M];
2482
if( *p != 0 )
2483
        goto normdn;
2484
++p;
2485
if( *p & 0x8000 )
2486
        return( 0 );     /* already normalized */
2487
while( *p == 0 )
2488
        {
2489
        eshup6(x);
2490
        sc += 16;
2491
/* With guard word, there are NBITS+16 bits available.
2492
 * return true if all are zero.
2493
 */
2494
        if( sc > NBITS )
2495
                return( sc );
2496
        }
2497
/* see if high byte is zero */
2498
while( (*p & 0xff00) == 0 )
2499
        {
2500
        eshup8(x);
2501
        sc += 8;
2502
        }
2503
/* now shift 1 bit at a time */
2504
while( (*p  & 0x8000) == 0)
2505
        {
2506
        eshup1(x);
2507
        sc += 1;
2508
        if( sc > (NBITS+16) )
2509
                {
2510
                mtherr( "enormlz", UNDERFLOW );
2511
                return( sc );
2512
                }
2513
        }
2514
return( sc );
2515
 
2516
/* Normalize by shifting down out of the high guard word
2517
   of the significand */
2518
normdn:
2519
 
2520
if( *p & 0xff00 )
2521
        {
2522
        eshdn8(x);
2523
        sc -= 8;
2524
        }
2525
while( *p != 0 )
2526
        {
2527
        eshdn1(x);
2528
        sc -= 1;
2529
 
2530
        if( sc < -NBITS )
2531
                {
2532
                mtherr( "enormlz", OVERFLOW );
2533
                return( sc );
2534
                }
2535
        }
2536
return( sc );
2537
}
2538
 
2539
 
2540
 
2541
 
2542
/* Convert e type number to decimal format ASCII string.
2543
 * The constants are for 64 bit precision.
2544
 */
2545
 
2546
#define NTEN 12
2547
#define MAXP 4096
2548
 
2549
#if NE == 10
2550
static unsigned short etens[NTEN + 1][NE] =
2551
{
2552
  {0x6576, 0x4a92, 0x804a, 0x153f,
2553
   0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
2554
  {0x6a32, 0xce52, 0x329a, 0x28ce,
2555
   0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
2556
  {0x526c, 0x50ce, 0xf18b, 0x3d28,
2557
   0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
2558
  {0x9c66, 0x58f8, 0xbc50, 0x5c54,
2559
   0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
2560
  {0x851e, 0xeab7, 0x98fe, 0x901b,
2561
   0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
2562
  {0x0235, 0x0137, 0x36b1, 0x336c,
2563
   0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
2564
  {0x50f8, 0x25fb, 0xc76b, 0x6b71,
2565
   0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
2566
  {0x0000, 0x0000, 0x0000, 0x0000,
2567
   0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
2568
  {0x0000, 0x0000, 0x0000, 0x0000,
2569
   0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
2570
  {0x0000, 0x0000, 0x0000, 0x0000,
2571
   0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
2572
  {0x0000, 0x0000, 0x0000, 0x0000,
2573
   0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
2574
  {0x0000, 0x0000, 0x0000, 0x0000,
2575
   0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
2576
  {0x0000, 0x0000, 0x0000, 0x0000,
2577
   0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
2578
};
2579
 
2580
static unsigned short emtens[NTEN + 1][NE] =
2581
{
2582
  {0x2030, 0xcffc, 0xa1c3, 0x8123,
2583
   0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
2584
  {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
2585
   0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
2586
  {0xf53f, 0xf698, 0x6bd3, 0x0158,
2587
   0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
2588
  {0xe731, 0x04d4, 0xe3f2, 0xd332,
2589
   0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
2590
  {0xa23e, 0x5308, 0xfefb, 0x1155,
2591
   0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
2592
  {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
2593
   0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
2594
  {0x2a20, 0x6224, 0x47b3, 0x98d7,
2595
   0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
2596
  {0x0b5b, 0x4af2, 0xa581, 0x18ed,
2597
   0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
2598
  {0xbf71, 0xa9b3, 0x7989, 0xbe68,
2599
   0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
2600
  {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
2601
   0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
2602
  {0xc155, 0xa4a8, 0x404e, 0x6113,
2603
   0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
2604
  {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
2605
   0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
2606
  {0xcccd, 0xcccc, 0xcccc, 0xcccc,
2607
   0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
2608
};
2609
#else
2610
static unsigned short etens[NTEN+1][NE] = {
2611
{0xc94c,0x979a,0x8a20,0x5202,0xc460,0x7525,},/* 10**4096 */
2612
{0xa74d,0x5de4,0xc53d,0x3b5d,0x9e8b,0x5a92,},/* 10**2048 */
2613
{0x650d,0x0c17,0x8175,0x7586,0xc976,0x4d48,},
2614
{0xcc65,0x91c6,0xa60e,0xa0ae,0xe319,0x46a3,},
2615
{0xddbc,0xde8d,0x9df9,0xebfb,0xaa7e,0x4351,},
2616
{0xc66f,0x8cdf,0x80e9,0x47c9,0x93ba,0x41a8,},
2617
{0x3cbf,0xa6d5,0xffcf,0x1f49,0xc278,0x40d3,},
2618
{0xf020,0xb59d,0x2b70,0xada8,0x9dc5,0x4069,},
2619
{0x0000,0x0000,0x0400,0xc9bf,0x8e1b,0x4034,},
2620
{0x0000,0x0000,0x0000,0x2000,0xbebc,0x4019,},
2621
{0x0000,0x0000,0x0000,0x0000,0x9c40,0x400c,},
2622
{0x0000,0x0000,0x0000,0x0000,0xc800,0x4005,},
2623
{0x0000,0x0000,0x0000,0x0000,0xa000,0x4002,}, /* 10**1 */
2624
};
2625
 
2626
static unsigned short emtens[NTEN+1][NE] = {
2627
{0x2de4,0x9fde,0xd2ce,0x04c8,0xa6dd,0x0ad8,}, /* 10**-4096 */
2628
{0x4925,0x2de4,0x3436,0x534f,0xceae,0x256b,}, /* 10**-2048 */
2629
{0x87a6,0xc0bd,0xda57,0x82a5,0xa2a6,0x32b5,},
2630
{0x7133,0xd21c,0xdb23,0xee32,0x9049,0x395a,},
2631
{0xfa91,0x1939,0x637a,0x4325,0xc031,0x3cac,},
2632
{0xac7d,0xe4a0,0x64bc,0x467c,0xddd0,0x3e55,},
2633
{0x3f24,0xe9a5,0xa539,0xea27,0xa87f,0x3f2a,},
2634
{0x67de,0x94ba,0x4539,0x1ead,0xcfb1,0x3f94,},
2635
{0x4c2f,0xe15b,0xc44d,0x94be,0xe695,0x3fc9,},
2636
{0xfdc2,0xcefc,0x8461,0x7711,0xabcc,0x3fe4,},
2637
{0xd3c3,0x652b,0xe219,0x1758,0xd1b7,0x3ff1,},
2638
{0x3d71,0xd70a,0x70a3,0x0a3d,0xa3d7,0x3ff8,},
2639
{0xcccd,0xcccc,0xcccc,0xcccc,0xcccc,0x3ffb,}, /* 10**-1 */
2640
};
2641
#endif
2642
 
2643
 
2644
 
2645
/* ASCII string outputs for unix */
2646
 
2647
 
2648
#if 0
2649
void _IO_ldtostr(x, string, ndigs, flags, fmt)
2650
long double *x;
2651
char *string;
2652
int ndigs;
2653
int flags;
2654
char fmt;
2655
{
2656
unsigned short w[NI];
2657
char *t, *u;
2658
LDPARMS rnd;
2659
LDPARMS *ldp = &rnd;
2660
 
2661
rnd.rlast = -1;
2662
rnd.rndprc = NBITS;
2663
 
2664
if (sizeof(long double) == 16)
2665
  e113toe( (unsigned short *)x, w, ldp );
2666
else
2667
  e64toe( (unsigned short *)x, w, ldp );
2668
 
2669
etoasc( w, string, ndigs, -1, ldp );
2670
if( ndigs == 0 && flags == 0 )
2671
        {
2672
        /* Delete the decimal point unless alternate format.  */
2673
        t = string;
2674
        while( *t != '.' )
2675
                ++t;
2676
        u = t +  1;
2677
        while( *t != '\0' )
2678
                *t++ = *u++;
2679
        }
2680
if (*string == ' ')
2681
        {
2682
        t = string;
2683
        u = t + 1;
2684
        while( *t != '\0' )
2685
                *t++ = *u++;
2686
        }
2687
if (fmt == 'E')
2688
        {
2689
        t = string;
2690
        while( *t != 'e' )
2691
                ++t;
2692
        *t = 'E';
2693
        }
2694
}
2695
 
2696
#endif
2697
 
2698
/* This routine will not return more than NDEC+1 digits. */
2699
 
2700
char *
2701
_ldtoa_r (struct _reent *ptr, long double d, int mode, int ndigits, int *decpt,
2702
          int *sign, char **rve)
2703
{
2704
unsigned short e[NI];
2705
char *s, *p;
2706
int k;
2707
LDPARMS rnd;
2708
LDPARMS *ldp = &rnd;
2709
char *outstr;
2710
 
2711
rnd.rlast = -1;
2712
rnd.rndprc = NBITS;
2713
 
2714
/* reentrancy addition to use mprec storage pool */
2715
if (ptr->_result)
2716
  {
2717
    ptr->_result->_k = ptr->_result_k;
2718
    ptr->_result->_maxwds = 1 << ptr->_result_k;
2719
    Bfree (ptr, ptr->_result);
2720
    ptr->_result = 0;
2721
  }
2722
 
2723
#if LDBL_MANT_DIG == 24
2724
e24toe( (unsigned short *)&d, e, ldp );
2725
#elif LDBL_MANT_DIG == 53
2726
e53toe( (unsigned short *)&d, e, ldp );
2727
#elif LDBL_MANT_DIG == 64
2728
e64toe( (unsigned short *)&d, e, ldp );
2729
#else
2730
e113toe( (unsigned short *)&d, e, ldp );
2731
#endif
2732
 
2733
if( eisneg(e) )
2734
        *sign = 1;
2735
else
2736
        *sign = 0;
2737
/* Mode 3 is "f" format.  */
2738
if( mode != 3 )
2739
        ndigits -= 1;
2740
/* Mode 0 is for %.999 format, which is supposed to give a
2741
   minimum length string that will convert back to the same binary value.
2742
   For now, just ask for 20 digits which is enough but sometimes too many.  */
2743
if( mode == 0 )
2744
        ndigits = 20;
2745
/* This sanity limit must agree with the corresponding one in etoasc, to
2746
   keep straight the returned value of outexpon.  */
2747
if( ndigits > NDEC )
2748
        ndigits = NDEC;
2749
 
2750
/* reentrancy addition to use mprec storage pool */
2751
ptr->_result = Balloc (ptr, 3);
2752
ptr->_result_k = 3;
2753
outstr = (char *)ptr->_result;
2754
 
2755
etoasc( e, outstr, ndigits, mode, ldp );
2756
s =  outstr;
2757
if( eisinf(e) || eisnan(e) )
2758
        {
2759
        *decpt = 9999;
2760
        goto stripspaces;
2761
        }
2762
*decpt = ldp->outexpon + 1;
2763
 
2764
/* Transform the string returned by etoasc into what the caller wants.  */
2765
 
2766
/* Look for decimal point and delete it from the string. */
2767
s = outstr;
2768
while( *s != '\0' )
2769
        {
2770
        if( *s == '.' )
2771
               goto yesdecpt;
2772
        ++s;
2773
        }
2774
goto nodecpt;
2775
 
2776
yesdecpt:
2777
 
2778
/* Delete the decimal point.  */
2779
while( *s != '\0' )
2780
        {
2781
        *s = *(s+1);
2782
        ++s;
2783
        }
2784
 
2785
nodecpt:
2786
 
2787
/* Back up over the exponent field. */
2788
while( *s != 'E' && s > outstr)
2789
        --s;
2790
*s = '\0';
2791
 
2792
stripspaces:
2793
 
2794
/* Strip leading spaces and sign. */
2795
p = outstr;
2796
while( *p == ' ' || *p == '-')
2797
        ++p;
2798
 
2799
/* Find new end of string.  */
2800
s = outstr;
2801
while( (*s++ = *p++) != '\0' )
2802
        ;
2803
--s;
2804
 
2805
/* Strip trailing zeros.  */
2806
if( mode == 2 )
2807
        k = 1;
2808
else if( ndigits > ldp->outexpon )
2809
        k = ndigits;
2810
else
2811
        k = ldp->outexpon;
2812
 
2813
while( *(s-1) == '0' && ((s - outstr) > k))
2814
        *(--s) = '\0';
2815
 
2816
/* In f format, flush small off-scale values to zero.
2817
   Rounding has been taken care of by etoasc. */
2818
if( mode == 3 && ((ndigits + ldp->outexpon) < 0))
2819
        {
2820
        s = outstr;
2821
        *s = '\0';
2822
        *decpt = 0;
2823
        }
2824
 
2825
if( rve )
2826
        *rve = s;
2827
return outstr;
2828
}
2829
 
2830
/* Routine used to tell if long double is NaN or Infinity or regular number.
2831
   Returns:  0 = regular number
2832
             1 = Nan
2833
             2 = Infinity
2834
*/
2835
int
2836
_ldcheck (long double *d)
2837
{
2838
unsigned short e[NI];
2839
char *s, *p;
2840
int k;
2841
LDPARMS rnd;
2842
LDPARMS *ldp = &rnd;
2843
char *outstr;
2844
 
2845
rnd.rlast = -1;
2846
rnd.rndprc = NBITS;
2847
 
2848
#if LDBL_MANT_DIG == 24
2849
e24toe( (unsigned short *)d, e, ldp );
2850
#elif LDBL_MANT_DIG == 53
2851
e53toe( (unsigned short *)d, e, ldp );
2852
#elif LDBL_MANT_DIG == 64
2853
e64toe( (unsigned short *)d, e, ldp );
2854
#else
2855
e113toe( (unsigned short *)d, e, ldp );
2856
#endif
2857
 
2858
if( (e[NE-1] & 0x7fff) == 0x7fff )
2859
        {
2860
#ifdef NANS
2861
        if( eisnan(e) )
2862
                return( 1 );
2863
#endif
2864
        return( 2 );
2865
        }
2866
else
2867
        return( 0 );
2868
} /* _ldcheck */
2869
 
2870
static void etoasc(short unsigned int *x, char *string, int ndigits, int outformat, LDPARMS *ldp)
2871
{
2872
long digit;
2873
unsigned short y[NI], t[NI], u[NI], w[NI];
2874
unsigned short *p, *r, *ten;
2875
unsigned short sign;
2876
int i, j, k, expon, rndsav, ndigs;
2877
char *s, *ss;
2878
unsigned short m;
2879
unsigned short *equot = ldp->equot;
2880
 
2881
ndigs = ndigits;
2882
rndsav = ldp->rndprc;
2883
#ifdef NANS
2884
if( eisnan(x) )
2885
        {
2886
        sprintf( string, " NaN " );
2887
        expon = 9999;
2888
        goto bxit;
2889
        }
2890
#endif
2891
ldp->rndprc = NBITS;            /* set to full precision */
2892
emov( x, y ); /* retain external format */
2893
if( y[NE-1] & 0x8000 )
2894
        {
2895
        sign = 0xffff;
2896
        y[NE-1] &= 0x7fff;
2897
        }
2898
else
2899
        {
2900
        sign = 0;
2901
        }
2902
expon = 0;
2903
ten = &etens[NTEN][0];
2904
emov( eone, t );
2905
/* Test for zero exponent */
2906
if( y[NE-1] == 0 )
2907
        {
2908
        for( k=0; k<NE-1; k++ )
2909
                {
2910
                if( y[k] != 0 )
2911
                        goto tnzro; /* denormalized number */
2912
                }
2913
        goto isone; /* legal all zeros */
2914
        }
2915
tnzro:
2916
 
2917
/* Test for infinity.
2918
 */
2919
if( y[NE-1] == 0x7fff )
2920
        {
2921
        if( sign )
2922
                sprintf( string, " -Infinity " );
2923
        else
2924
                sprintf( string, " Infinity " );
2925
        expon = 9999;
2926
        goto bxit;
2927
        }
2928
 
2929
/* Test for exponent nonzero but significand denormalized.
2930
 * This is an error condition.
2931
 */
2932
if( (y[NE-1] != 0) && ((y[NE-2] & 0x8000) == 0) )
2933
        {
2934
        mtherr( "etoasc", DOMAIN );
2935
        sprintf( string, "NaN" );
2936
        expon = 9999;
2937
        goto bxit;
2938
        }
2939
 
2940
/* Compare to 1.0 */
2941
i = ecmp( eone, y );
2942
if( i == 0 )
2943
        goto isone;
2944
 
2945
if( i < 0 )
2946
        { /* Number is greater than 1 */
2947
/* Convert significand to an integer and strip trailing decimal zeros. */
2948
        emov( y, u );
2949
        u[NE-1] = EXONE + NBITS - 1;
2950
 
2951
        p = &etens[NTEN-4][0];
2952
        m = 16;
2953
do
2954
        {
2955
        ediv( p, u, t, ldp );
2956
        efloor( t, w, ldp );
2957
        for( j=0; j<NE-1; j++ )
2958
                {
2959
                if( t[j] != w[j] )
2960
                        goto noint;
2961
                }
2962
        emov( t, u );
2963
        expon += (int )m;
2964
noint:
2965
        p += NE;
2966
        m >>= 1;
2967
        }
2968
while( m != 0 );
2969
 
2970
/* Rescale from integer significand */
2971
        u[NE-1] += y[NE-1] - (unsigned int )(EXONE + NBITS - 1);
2972
        emov( u, y );
2973
/* Find power of 10 */
2974
        emov( eone, t );
2975
        m = MAXP;
2976
        p = &etens[0][0];
2977
        while( ecmp( ten, u ) <= 0 )
2978
                {
2979
                if( ecmp( p, u ) <= 0 )
2980
                        {
2981
                        ediv( p, u, u, ldp );
2982
                        emul( p, t, t, ldp );
2983
                        expon += (int )m;
2984
                        }
2985
                m >>= 1;
2986
                if( m == 0 )
2987
                        break;
2988
                p += NE;
2989
                }
2990
        }
2991
else
2992
        { /* Number is less than 1.0 */
2993
/* Pad significand with trailing decimal zeros. */
2994
        if( y[NE-1] == 0 )
2995
                {
2996
                while( (y[NE-2] & 0x8000) == 0 )
2997
                        {
2998
                        emul( ten, y, y, ldp );
2999
                        expon -= 1;
3000
                        }
3001
                }
3002
        else
3003
                {
3004
                emovi( y, w );
3005
                for( i=0; i<NDEC+1; i++ )
3006
                        {
3007
                        if( (w[NI-1] & 0x7) != 0 )
3008
                                break;
3009
/* multiply by 10 */
3010
                        emovz( w, u );
3011
                        eshdn1( u );
3012
                        eshdn1( u );
3013
                        eaddm( w, u );
3014
                        u[1] += 3;
3015
                        while( u[2] != 0 )
3016
                                {
3017
                                eshdn1(u);
3018
                                u[1] += 1;
3019
                                }
3020
                        if( u[NI-1] != 0 )
3021
                                break;
3022
                        if( eone[NE-1] <= u[1] )
3023
                                break;
3024
                        emovz( u, w );
3025
                        expon -= 1;
3026
                        }
3027
                emovo( w, y, ldp );
3028
                }
3029
        k = -MAXP;
3030
        p = &emtens[0][0];
3031
        r = &etens[0][0];
3032
        emov( y, w );
3033
        emov( eone, t );
3034
        while( ecmp( eone, w ) > 0 )
3035
                {
3036
                if( ecmp( p, w ) >= 0 )
3037
                        {
3038
                        emul( r, w, w, ldp );
3039
                        emul( r, t, t, ldp );
3040
                        expon += k;
3041
                        }
3042
                k /= 2;
3043
                if( k == 0 )
3044
                        break;
3045
                p += NE;
3046
                r += NE;
3047
                }
3048
        ediv( t, eone, t, ldp );
3049
        }
3050
isone:
3051
/* Find the first (leading) digit. */
3052
emovi( t, w );
3053
emovz( w, t );
3054
emovi( y, w );
3055
emovz( w, y );
3056
eiremain( t, y, ldp );
3057
digit = equot[NI-1];
3058
while( (digit == 0) && (ecmp(y,ezero) != 0) )
3059
        {
3060
        eshup1( y );
3061
        emovz( y, u );
3062
        eshup1( u );
3063
        eshup1( u );
3064
        eaddm( u, y );
3065
        eiremain( t, y, ldp );
3066
        digit = equot[NI-1];
3067
        expon -= 1;
3068
        }
3069
s = string;
3070
if( sign )
3071
        *s++ = '-';
3072
else
3073
        *s++ = ' ';
3074
/* Examine number of digits requested by caller. */
3075
if( outformat == 3 )
3076
        ndigs += expon;
3077
/*
3078
else if( ndigs < 0 )
3079
        ndigs = 0;
3080
*/
3081
if( ndigs > NDEC )
3082
        ndigs = NDEC;
3083
if( digit == 10 )
3084
        {
3085
        *s++ = '1';
3086
        *s++ = '.';
3087
        if( ndigs > 0 )
3088
                {
3089
                *s++ = '0';
3090
                ndigs -= 1;
3091
                }
3092
        expon += 1;
3093
        if( ndigs < 0 )
3094
                {
3095
                ss = s;
3096
                goto doexp;
3097
                }
3098
        }
3099
else
3100
        {
3101
        *s++ = (char )digit + '0';
3102
        *s++ = '.';
3103
        }
3104
/* Generate digits after the decimal point. */
3105
for( k=0; k<=ndigs; k++ )
3106
        {
3107
/* multiply current number by 10, without normalizing */
3108
        eshup1( y );
3109
        emovz( y, u );
3110
        eshup1( u );
3111
        eshup1( u );
3112
        eaddm( u, y );
3113
        eiremain( t, y, ldp );
3114
        *s++ = (char )equot[NI-1] + '0';
3115
        }
3116
digit = equot[NI-1];
3117
--s;
3118
ss = s;
3119
/* round off the ASCII string */
3120
if( digit > 4 )
3121
        {
3122
/* Test for critical rounding case in ASCII output. */
3123
        if( digit == 5 )
3124
                {
3125
                emovo( y, t, ldp );
3126
                if( ecmp(t,ezero) != 0 )
3127
                        goto roun;      /* round to nearest */
3128
                if( (*(s-1) & 1) == 0 )
3129
                        goto doexp;     /* round to even */
3130
                }
3131
/* Round up and propagate carry-outs */
3132
roun:
3133
        --s;
3134
        k = *s & 0x7f;
3135
/* Carry out to most significant digit? */
3136
        if( ndigs < 0 )
3137
                {
3138
                /* This will print like "1E-6". */
3139
                *s = '1';
3140
                expon += 1;
3141
                goto doexp;
3142
                }
3143
        else if( k == '.' )
3144
                {
3145
                --s;
3146
                k = *s;
3147
                k += 1;
3148
                *s = (char )k;
3149
/* Most significant digit carries to 10? */
3150
                if( k > '9' )
3151
                        {
3152
                        expon += 1;
3153
                        *s = '1';
3154
                        }
3155
                goto doexp;
3156
                }
3157
/* Round up and carry out from less significant digits */
3158
        k += 1;
3159
        *s = (char )k;
3160
        if( k > '9' )
3161
                {
3162
                *s = '0';
3163
                goto roun;
3164
                }
3165
        }
3166
doexp:
3167
#ifdef __GO32__
3168
if( expon >= 0 )
3169
        sprintf( ss, "e+%02d", expon );
3170
else
3171
        sprintf( ss, "e-%02d", -expon );
3172
#else
3173
        sprintf( ss, "E%d", expon );
3174
#endif
3175
bxit:
3176
ldp->rndprc = rndsav;
3177
ldp->outexpon =  expon;
3178
}
3179
 
3180
 
3181
 
3182
 
3183
/*
3184
;                                                               ASCTOQ
3185
;               ASCTOQ.MAC              LATEST REV: 11 JAN 84
3186
;                                       SLM, 3 JAN 78
3187
;
3188
;       Convert ASCII string to quadruple precision floating point
3189
;
3190
;               Numeric input is free field decimal number
3191
;               with max of 15 digits with or without
3192
;               decimal point entered as ASCII from teletype.
3193
;       Entering E after the number followed by a second
3194
;       number causes the second number to be interpreted
3195
;       as a power of 10 to be multiplied by the first number
3196
;       (i.e., "scientific" notation).
3197
;
3198
;       Usage:
3199
;               asctoq( string, q );
3200
*/
3201
 
3202
long double _strtold (char *s, char **se)
3203
{
3204
  long double x;
3205
  LDPARMS rnd;
3206
  LDPARMS *ldp = &rnd;
3207
  int lenldstr;
3208
 
3209
  rnd.rlast = -1;
3210
  rnd.rndprc = NBITS;
3211
 
3212
  lenldstr = asctoeg( s, (unsigned short *)&x, LDBL_MANT_DIG, ldp );
3213
  if (se)
3214
    *se = s + lenldstr;
3215
  return x;
3216
}
3217
 
3218
#define REASONABLE_LEN 200
3219
 
3220
static int
3221
asctoeg(char *ss, short unsigned int *y, int oprec, LDPARMS *ldp)
3222
{
3223
unsigned short yy[NI], xt[NI], tt[NI];
3224
int esign, decflg, sgnflg, nexp, exp, prec, lost;
3225
int k, trail, c, rndsav;
3226
long lexp;
3227
unsigned short nsign, *p;
3228
char *sp, *s, *lstr;
3229
int lenldstr;
3230
int mflag = 0;
3231
char tmpstr[REASONABLE_LEN];
3232
 
3233
/* Copy the input string. */
3234
c = strlen (ss) + 2;
3235
if (c <= REASONABLE_LEN)
3236
  lstr = tmpstr;
3237
else
3238
  {
3239
    lstr = (char *) calloc (c, 1);
3240
    mflag = 1;
3241
  }
3242
s = ss;
3243
lenldstr = 0;
3244
while( *s == ' ' ) /* skip leading spaces */
3245
  {
3246
    ++s;
3247
    ++lenldstr;
3248
  }
3249
sp = lstr;
3250
for( k=0; k<c; k++ )
3251
        {
3252
        if( (*sp++ = *s++) == '\0' )
3253
                break;
3254
        }
3255
*sp = '\0';
3256
s = lstr;
3257
 
3258
rndsav = ldp->rndprc;
3259
ldp->rndprc = NBITS; /* Set to full precision */
3260
lost = 0;
3261
nsign = 0;
3262
decflg = 0;
3263
sgnflg = 0;
3264
nexp = 0;
3265
exp = 0;
3266
prec = 0;
3267
ecleaz( yy );
3268
trail = 0;
3269
 
3270
nxtcom:
3271
k = *s - '0';
3272
if( (k >= 0) && (k <= 9) )
3273
        {
3274
/* Ignore leading zeros */
3275
        if( (prec == 0) && (decflg == 0) && (k == 0) )
3276
                goto donchr;
3277
/* Identify and strip trailing zeros after the decimal point. */
3278
        if( (trail == 0) && (decflg != 0) )
3279
                {
3280
                sp = s;
3281
                while( (*sp >= '0') && (*sp <= '9') )
3282
                        ++sp;
3283
/* Check for syntax error */
3284
                c = *sp & 0x7f;
3285
                if( (c != 'e') && (c != 'E') && (c != '\0')
3286
                        && (c != '\n') && (c != '\r') && (c != ' ')
3287
                        && (c != ',') )
3288
                        goto error;
3289
                --sp;
3290
                while( *sp == '0' )
3291
                        *sp-- = 'z';
3292
                trail = 1;
3293
                if( *s == 'z' )
3294
                        goto donchr;
3295
                }
3296
/* If enough digits were given to more than fill up the yy register,
3297
 * continuing until overflow into the high guard word yy[2]
3298
 * guarantees that there will be a roundoff bit at the top
3299
 * of the low guard word after normalization.
3300
 */
3301
        if( yy[2] == 0 )
3302
                {
3303
                if( decflg )
3304
                        nexp += 1; /* count digits after decimal point */
3305
                eshup1( yy );   /* multiply current number by 10 */
3306
                emovz( yy, xt );
3307
                eshup1( xt );
3308
                eshup1( xt );
3309
                eaddm( xt, yy );
3310
                ecleaz( xt );
3311
                xt[NI-2] = (unsigned short )k;
3312
                eaddm( xt, yy );
3313
                }
3314
        else
3315
                {
3316
                /* Mark any lost non-zero digit.  */
3317
                lost |= k;
3318
                /* Count lost digits before the decimal point.  */
3319
                if (decflg == 0)
3320
                        nexp -= 1;
3321
                }
3322
        prec += 1;
3323
        goto donchr;
3324
        }
3325
 
3326
switch( *s )
3327
        {
3328
        case 'z':
3329
                break;
3330
        case 'E':
3331
        case 'e':
3332
                goto expnt;
3333
        case '.':       /* decimal point */
3334
                if( decflg )
3335
                        goto error;
3336
                ++decflg;
3337
                break;
3338
        case '-':
3339
                nsign = 0xffff;
3340
                if( sgnflg )
3341
                        goto error;
3342
                ++sgnflg;
3343
                break;
3344
        case '+':
3345
                if( sgnflg )
3346
                        goto error;
3347
                ++sgnflg;
3348
                break;
3349
        case ',':
3350
        case ' ':
3351
        case '\0':
3352
        case '\n':
3353
        case '\r':
3354
                goto daldone;
3355
        case 'i':
3356
        case 'I':
3357
                goto infinite;
3358
        default:
3359
        error:
3360
#ifdef NANS
3361
                enan( yy, NI*16 );
3362
#else
3363
                mtherr( "asctoe", DOMAIN );
3364
                ecleaz(yy);
3365
#endif
3366
                goto aexit;
3367
        }
3368
donchr:
3369
++s;
3370
goto nxtcom;
3371
 
3372
/* Exponent interpretation */
3373
expnt:
3374
 
3375
esign = 1;
3376
exp = 0;
3377
++s;
3378
/* check for + or - */
3379
if( *s == '-' )
3380
        {
3381
        esign = -1;
3382
        ++s;
3383
        }
3384
if( *s == '+' )
3385
        ++s;
3386
while( (*s >= '0') && (*s <= '9') )
3387
        {
3388
        exp *= 10;
3389
        exp += *s++ - '0';
3390
        if (exp > 4977)
3391
                {
3392
                if (esign < 0)
3393
                        goto zero;
3394
                else
3395
                        goto infinite;
3396
                }
3397
        }
3398
if( esign < 0 )
3399
        exp = -exp;
3400
if( exp > 4932 )
3401
        {
3402
infinite:
3403
        ecleaz(yy);
3404
        yy[E] = 0x7fff;  /* infinity */
3405
        goto aexit;
3406
        }
3407
if( exp < -4977 )
3408
        {
3409
zero:
3410
        ecleaz(yy);
3411
        goto aexit;
3412
        }
3413
 
3414
daldone:
3415
nexp = exp - nexp;
3416
/* Pad trailing zeros to minimize power of 10, per IEEE spec. */
3417
while( (nexp > 0) && (yy[2] == 0) )
3418
        {
3419
        emovz( yy, xt );
3420
        eshup1( xt );
3421
        eshup1( xt );
3422
        eaddm( yy, xt );
3423
        eshup1( xt );
3424
        if( xt[2] != 0 )
3425
                break;
3426
        nexp -= 1;
3427
        emovz( xt, yy );
3428
        }
3429
if( (k = enormlz(yy)) > NBITS )
3430
        {
3431
        ecleaz(yy);
3432
        goto aexit;
3433
        }
3434
lexp = (EXONE - 1 + NBITS) - k;
3435
emdnorm( yy, lost, 0, lexp, 64, ldp );
3436
/* convert to external format */
3437
 
3438
 
3439
/* Multiply by 10**nexp.  If precision is 64 bits,
3440
 * the maximum relative error incurred in forming 10**n
3441
 * for 0 <= n <= 324 is 8.2e-20, at 10**180.
3442
 * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
3443
 * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
3444
 */
3445
lexp = yy[E];
3446
if( nexp == 0 )
3447
        {
3448
        k = 0;
3449
        goto expdon;
3450
        }
3451
esign = 1;
3452
if( nexp < 0 )
3453
        {
3454
        nexp = -nexp;
3455
        esign = -1;
3456
        if( nexp > 4096 )
3457
                { /* Punt.  Can't handle this without 2 divides. */
3458
                emovi( etens[0], tt );
3459
                lexp -= tt[E];
3460
                k = edivm( tt, yy, ldp );
3461
                lexp += EXONE;
3462
                nexp -= 4096;
3463
                }
3464
        }
3465
p = &etens[NTEN][0];
3466
emov( eone, xt );
3467
exp = 1;
3468
do
3469
        {
3470
        if( exp & nexp )
3471
                emul( p, xt, xt, ldp );
3472
        p -= NE;
3473
        exp = exp + exp;
3474
        }
3475
while( exp <= MAXP );
3476
 
3477
emovi( xt, tt );
3478
if( esign < 0 )
3479
        {
3480
        lexp -= tt[E];
3481
        k = edivm( tt, yy, ldp );
3482
        lexp += EXONE;
3483
        }
3484
else
3485
        {
3486
        lexp += tt[E];
3487
        k = emulm( tt, yy, ldp );
3488
        lexp -= EXONE - 1;
3489
        }
3490
 
3491
expdon:
3492
 
3493
/* Round and convert directly to the destination type */
3494
if( oprec == 53 )
3495
        lexp -= EXONE - 0x3ff;
3496
else if( oprec == 24 )
3497
        lexp -= EXONE - 0177;
3498
#ifdef DEC
3499
else if( oprec == 56 )
3500
        lexp -= EXONE - 0201;
3501
#endif
3502
ldp->rndprc = oprec;
3503
emdnorm( yy, k, 0, lexp, 64, ldp );
3504
 
3505
aexit:
3506
 
3507
ldp->rndprc = rndsav;
3508
yy[0] = nsign;
3509
switch( oprec )
3510
        {
3511
#ifdef DEC
3512
        case 56:
3513
                todec( yy, y ); /* see etodec.c */
3514
                break;
3515
#endif
3516
#if LDBL_MANT_DIG == 53
3517
        case 53:
3518
                toe53( yy, y );
3519
                break;
3520
#elif LDBL_MANT_DIG == 24
3521
        case 24:
3522
                toe24( yy, y );
3523
                break;
3524
#elif LDBL_MANT_DIG == 64
3525
        case 64:
3526
                toe64( yy, y );
3527
                break;
3528
#elif LDBL_MANT_DIG == 113
3529
        case 113:
3530
                toe113( yy, y );
3531
                break;
3532
#else
3533
        case NBITS:
3534
                emovo( yy, y, ldp );
3535
                break;
3536
#endif
3537
        }
3538
lenldstr += s - lstr;
3539
if (mflag)
3540
  free (lstr);
3541
return lenldstr;
3542
}
3543
 
3544
 
3545
 
3546
/* y = largest integer not greater than x
3547
 * (truncated toward minus infinity)
3548
 *
3549
 * unsigned short x[NE], y[NE]
3550
 * LDPARMS *ldp
3551
 *
3552
 * efloor( x, y, ldp );
3553
 */
3554
static unsigned short bmask[] = {
3555
0xffff,
3556
0xfffe,
3557
0xfffc,
3558
0xfff8,
3559
0xfff0,
3560
0xffe0,
3561
0xffc0,
3562
0xff80,
3563
0xff00,
3564
0xfe00,
3565
0xfc00,
3566
0xf800,
3567
0xf000,
3568
0xe000,
3569
0xc000,
3570
0x8000,
3571
0x0000,
3572
};
3573
 
3574
static void efloor(short unsigned int *x, short unsigned int *y, LDPARMS *ldp)
3575
{
3576
register unsigned short *p;
3577
int e, expon, i;
3578
unsigned short f[NE];
3579
 
3580
emov( x, f ); /* leave in external format */
3581
expon = (int )f[NE-1];
3582
e = (expon & 0x7fff) - (EXONE - 1);
3583
if( e <= 0 )
3584
        {
3585
        eclear(y);
3586
        goto isitneg;
3587
        }
3588
/* number of bits to clear out */
3589
e = NBITS - e;
3590
emov( f, y );
3591
if( e <= 0 )
3592
        return;
3593
 
3594
p = &y[0];
3595
while( e >= 16 )
3596
        {
3597
        *p++ = 0;
3598
        e -= 16;
3599
        }
3600
/* clear the remaining bits */
3601
*p &= bmask[e];
3602
/* truncate negatives toward minus infinity */
3603
isitneg:
3604
 
3605
if( (unsigned short )expon & (unsigned short )0x8000 )
3606
        {
3607
        for( i=0; i<NE-1; i++ )
3608
                {
3609
                if( f[i] != y[i] )
3610
                        {
3611
                        esub( eone, y, y, ldp );
3612
                        break;
3613
                        }
3614
                }
3615
        }
3616
}
3617
 
3618
 
3619
 
3620
static void eiremain(short unsigned int *den, short unsigned int *num, LDPARMS *ldp)
3621
{
3622
long ld, ln;
3623
unsigned short j;
3624
 unsigned short *equot = ldp->equot;
3625
 
3626
ld = den[E];
3627
ld -= enormlz( den );
3628
ln = num[E];
3629
ln -= enormlz( num );
3630
ecleaz( equot );
3631
while( ln >= ld )
3632
        {
3633
        if( ecmpm(den,num) <= 0 )
3634
                {
3635
                esubm(den, num);
3636
                j = 1;
3637
                }
3638
        else
3639
                {
3640
                j = 0;
3641
                }
3642
        eshup1(equot);
3643
        equot[NI-1] |= j;
3644
        eshup1(num);
3645
        ln -= 1;
3646
        }
3647
emdnorm( num, 0, 0, ln, 0, ldp );
3648
}
3649
 
3650
/* NaN bit patterns
3651
 */
3652
#ifdef MIEEE
3653
static unsigned short nan113[8] = {
3654
  0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
3655
static unsigned short nan64[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
3656
static unsigned short nan53[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
3657
static unsigned short nan24[2] = {0x7fff, 0xffff};
3658
#else /* !MIEEE */
3659
static unsigned short nan113[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0x7fff};
3660
static unsigned short nan64[6] = {0, 0, 0, 0, 0xc000, 0x7fff};
3661
static unsigned short nan53[4] = {0, 0, 0, 0x7ff8};
3662
static unsigned short nan24[2] = {0, 0x7fc0};
3663
#endif /* !MIEEE */
3664
 
3665
 
3666
static void enan (short unsigned int *nan, int size)
3667
{
3668
int i, n;
3669
unsigned short *p;
3670
 
3671
switch( size )
3672
        {
3673
#ifndef DEC
3674
        case 113:
3675
        n = 8;
3676
        p = nan113;
3677
        break;
3678
 
3679
        case 64:
3680
        n = 6;
3681
        p = nan64;
3682
        break;
3683
 
3684
        case 53:
3685
        n = 4;
3686
        p = nan53;
3687
        break;
3688
 
3689
        case 24:
3690
        n = 2;
3691
        p = nan24;
3692
        break;
3693
 
3694
        case NBITS:
3695
        for( i=0; i<NE-2; i++ )
3696
                *nan++ = 0;
3697
        *nan++ = 0xc000;
3698
        *nan++ = 0x7fff;
3699
        return;
3700
 
3701
        case NI*16:
3702
        *nan++ = 0;
3703
        *nan++ = 0x7fff;
3704
        *nan++ = 0;
3705
        *nan++ = 0xc000;
3706
        for( i=4; i<NI; i++ )
3707
                *nan++ = 0;
3708
        return;
3709
#endif
3710
        default:
3711
        mtherr( "enan", DOMAIN );
3712
        return;
3713
        }
3714
for (i=0; i < n; i++)
3715
        *nan++ = *p++;
3716
}
3717
 
3718
 
3719
 
3720
 

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.