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

Subversion Repositories openrisc_me

[/] [openrisc/] [trunk/] [gnu-src/] [newlib-1.17.0/] [newlib/] [libc/] [stdlib/] [ldtoa.c] - Blame information for rev 158

Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

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