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

Subversion Repositories openrisc_me

[/] [openrisc/] [trunk/] [gnu-src/] [newlib-1.18.0/] [newlib/] [libc/] [machine/] [powerpc/] [simdldtoa.c] - Blame information for rev 207

Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

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