Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

simdldtoa.c @ 444

Last change on this file since 444 was 444, checked in by satin@…, 6 years ago
add newlib,libalmos-mkh, restructure shared_syscalls.h and mini-libc
File size: 68.5 KB

Line
1
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__ */

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/libs/newlib/src/newlib/libc/machine/powerpc/simdldtoa.c @ 444

Download in other formats: