[unix-history] / usr / src / old / as.vax / natof.c

/*
 *	Copyright (c) 1982 Regents of the University of California
 */
#ifndef lint
static char sccsid[] = "@(#)natof.c 4.3 %G%";
#endif not lint

#include <stdio.h>
#include <ctype.h>
#include <errno.h>

#include "as.h"

Bignum bigatof(str, radix)
	reg	char	*str;		/* r11 */
		int	radix;		/* TYPF ... TYPH */
{
		int	msign;
		int	esign;
		int	decpt;
	reg	chptr	temp;		/* r10 */
	reg	u_int	quotient;	/* r9 */	/* must be here */
	reg	u_int	remainder;	/* r8 */	/* must be here */
	reg	chptr	acc;
	reg	int	dividend;	/* for doing division */
	reg	u_int	i;	
		short	*sptr;		/* for doing division */
		int	ch;
		int	dexponent;	/* decimal exponent */
		int	bexponent;	/* binary exponent */
		Bignum	Acc;
		Bignum	Temp;
	static	Bignum	znumber;
		Ovf	ovf;
		u_int	j;
	extern	int	errno;
		u_int	ediv();

#ifdef lint
	quotient = 0;
	remainder = 0;
#endif lint
	msign = 0;
	esign = 0;
	decpt = 0;
	dexponent = 0;
	Acc = znumber;
	Acc.num_tag = radix;
	acc = CH_FIELD(Acc);
	temp = CH_FIELD(Temp);

	do{
		ch = *str++;
	} while(isspace(ch));

	switch(ch){
	case '-':
		msign = -1;
		/* FALLTHROUGH */
	case '+':
		ch = *str++;
		break;
	}
  dofract:
	for(; isdigit(ch); ch = *str++){
		assert(((acc[HOC] & SIGNBIT) == 0), "Negative HOC");
		if (acc[HOC] < MAXINT_10){
			ovf = numshift(3, temp, acc);
			ovf |= numshift(1, acc, acc);
			ovf |= numaddv(acc, temp, acc);
			ovf |= numaddd(acc, acc, ch - '0');
			assert(ovf == 0, "Overflow building mantissa");
		} else {
			/*
			 *	Then, the number is too large anyway
			 */
			dexponent++;
		}
		if (decpt)
			dexponent--;
	}
	switch(ch){
	case '.':
		if (decpt == 0){
			decpt++;
			ch = *str++;
			goto dofract;
		}
		break;
	/*
	 *	only 'e' and 'E' are recognized by atof()
	 */
	case 'e':
	case 'E':
	/*
	 *	we include the remainder for compatability with as formats
	 *	in as, the radix actual paramater agrees with the character
	 *	we expect; consequently, no checking is done.
	 */
	case 'd':
	case 'D':
	case 'g':
	case 'G':
	case 'h':
	case 'H':
		j = 0;
		ch = *str++;
		esign = 0;
		switch(ch){
		case '-':
			esign = 1;
			/* FALLTHROUGH */
		case '+':
			ch = *str++;
		}
		for(; isdigit(ch); ch = *str++){
			if (j < MAXINT_10){
				j *= 10;
				j += ch - '0';
			} else {
				/*
				 *	outrageously large exponent
				 */
				 /*VOID*/
			}
		}
		if (esign)
			dexponent -= j;
		else
			dexponent += j;
		/*
		 *	There should be a range check on dexponent here
		 */
	}
	/*
	 *	The number has now been reduced to a mantissa
	 *	and an exponent.
	 *	The mantissa is an n bit number (to the precision
	 *	of the extended words) in the acc.
	 *	The exponent is a signed power of 10 in dexponent.
	 *	msign is on if the resulting number will eventually
	 *	be negative.
	 *
	 *	We now must convert the number to standard format floating
	 *	number, which will be done by accumulating
	 *	a binary exponent in bexponent, as we gradually
	 *	drive dexponent towards zero, one count at a time.
	 */
	if (isclear(acc)){
		return(Acc);
	}
	bexponent = 0;

	/*
	 *	Scale the number down.
	 *	We must divide acc by 10 as many times as needed.
	 */
	for (; dexponent < 0; dexponent++){
		/*
		 *	Align the number so that the most significant
		 *	bits are aligned in the most significant
		 *	bits of the accumulator, adjusting the
		 *	binary exponent as we shift.
		 *	The goal is to get the high order bit (NOT the
		 *	sign bit) set.
		 */
		assert(((acc[HOC] & SIGNBIT) == 0), "Negative HOC");
		ovf = 0;

		for (j = 5; j >= 1; --j){
			i = 1 << (j - 1); 		/* 16, 8, 4, 2, 1 */
			quotient = ONES(i);
			quotient <<= (CH_BITS - 1) - i;
			while((acc[HOC] & quotient) == 0){
				ovf |= numshift((int)i, acc, acc);
				bexponent -= i;
			}
		}
		/*
		 *	Add 2 to the accumulator to effect rounding,
		 *	and get set up to divide by 5.
		 */
		ovf = numaddd(acc, acc, 2);
		assert(ovf == 0, "Carry out of left rounding up by 2");
		/*
		 *	Divide the high order chunks by 5;
		 *	The last chunk will be divided by 10,
		 *	(to see what the remainder is, also to effect rounding)
		 *	and then multipiled by 2 to effect division by 5.
		 */
		remainder = 0;
#if DEBUGNATOF
		printf("Dividing: ");
		bignumprint(Acc);
		printf("\n");
#endif DEBUGNATOF
		sptr = (short *)acc;
		for (i = (CH_N * 2 - 1); i >= 1; --i){
			/*
			 *	Divide (remainder:16).(acc[i]:16)
			 *	by 5, putting the quotient back
			 *	into acc[i]:16, and save the remainder
			 *	for the next iteration.
			 */
			dividend = (remainder << 16) | (sptr[i] & ONES(16));
			assert(dividend >= 0, "dividend < 0");
			quotient = dividend / 5;
			remainder = dividend - (quotient * 5);
			sptr[i] = quotient;
			remainder = remainder;
		}
		/*
		 *	Divide the lowest order chunk by 10,
		 *	saving the remainder to decide how to round.
		 *	Then, multiply by 2, making it look as
		 *	if we divided by 10.
		 *	This multiply fills in a 0 on the least sig bit.
		 */
		dividend = (remainder << 16) | (sptr[0] & ONES(16));
		assert(dividend >= 0, "dividend < 0");
		quotient = dividend / 10;
		remainder = dividend - (quotient * 10);
		sptr[0] = quotient + quotient;

		if (remainder >= 5)
			ovf = numaddd(acc, acc, 1);
		/*
		 *	Now, divide by 2, effecting division by 10,
		 *	merely by adjusting the binary exponent.
		 */
		bexponent--;
	}
	/*
	 *	Scale the number up by multiplying by 10 as 
	 *	many times as necessary
	 */
	for (; dexponent > 0; dexponent--){
		/*
		 *	Compare high word to (2**31)/5,
		 *	and scale accordingly
		 */
		while ( ((unsigned)acc[HOC]) > MAXINT_5){
			(void)numshift(-1, acc, acc);
			bexponent++;
		}
		/*
		 *	multiply the mantissa by 5,
		 *	and scale the binary exponent by 2
		 */
		ovf = numshift(2, temp, acc);
		ovf |= numaddv(acc, acc, temp);
		assert(ovf == 0, "Scaling * 10 of manitissa");
		bexponent++;
	}
	/*
	 *	We now have:
	 *	a CH_N chunk length binary integer, right
	 *		justified (in native format).
	 *	a binary exponent.
	 *
	 *	Now, we treat this large integer as an octa word
	 *	number, and unpack it into standard unpacked
	 *	format.  That unpacking will give us yet
	 *	another binary exponent, which we adjust with
	 *	the accumulated binary exponent.
	 */
	Acc.num_tag = TYPO;
#if DEBUGNATOF
	printf("Octal number: ");
	bignumprint(Acc);
	printf("\n");
#endif DEBUGNATOF
	Acc = bignumunpack(Acc, &ovf);

	if (ovf)
		errno = ERANGE;
#if DEBUGNATOF
	printf("Unpacked octal number: ");
	bignumprint(Acc);
	printf("bexponent == %d\n", bexponent);
#endif DEBUGNATOF
	Acc.num_exponent += bexponent;
	assert(Acc.num_sign == 0, "unpacked integer is < 0");
	Acc.num_sign = msign;
	/*
	 *	We now pack the number back into a radix format number.
	 *	This checks for overflow, underflow,
	 *	and rounds by 1/2 ulp.
	 */
	ovf = 0;
	Acc = bignumpack(Acc, radix, &ovf);
	if (ovf)
		errno = ERANGE;
#if DEBUGNATOF
	printf("packed number: ");
	bignumprint(Acc);
	printf("\n");
#endif DEBUGNATOF
	return(Acc);
}
#if 0
/*
 *	Unfortunately, one can't use the ediv instruction to do
 *	division on numbers with > 64 bits.
 *	This is because ediv returns signed quantities;
 *	if the quotient is an unsigned number > 2^31,
 *	ediv sets integer overflow.
 */
unsigned int ediv(high, low, divisor, qp, i)
	register	unsigned int	high;		/* r11 */
	register	unsigned int	low;		/* r10 */
	register	unsigned int	divisor;	/* r9 */
			unsigned int	*qp;
{
	register	unsigned int	remainder;	/* r8 */
	register	unsigned int	quotient;	/* r7 */

	asm("ediv r9, r10, r7, r8	# Divide.  q->r7, r->r8 (discarded)");
	*qp = quotient;
	return(remainder);
}
#endif 0
Commit	Line	Data
7a47add7 RH	1	/*
	2	* Copyright (c) 1982 Regents of the University of California
	3	*/
	4	#ifndef lint
f70ab843	5	static char sccsid[] = "@(#)natof.c 4.3 %G%";
7a47add7 RH	6	#endif not lint
	7
	8	#include <stdio.h>
	9	#include <ctype.h>
	10	#include <errno.h>
	11
	12	#include "as.h"
	13
	14	Bignum bigatof(str, radix)
	15	reg char str; / r11 */
	16	int radix; /* TYPF ... TYPH */
	17	{
	18	int msign;
	19	int esign;
	20	int decpt;
	21	reg chptr temp; /* r10 */
	22	reg u_int quotient; /* r9 / / must be here */
	23	reg u_int remainder; /* r8 / / must be here */
	24	reg chptr acc;
	25	reg int dividend; /* for doing division */
	26	reg u_int i;
	27	short sptr; / for doing division */
	28	int ch;
	29	int dexponent; /* decimal exponent */
	30	int bexponent; /* binary exponent */
	31	Bignum Acc;
	32	Bignum Temp;
	33	static Bignum znumber;
	34	Ovf ovf;
	35	u_int j;
	36	extern int errno;
	37	u_int ediv();
	38
	39	#ifdef lint
	40	quotient = 0;
	41	remainder = 0;
	42	#endif lint
	43	msign = 0;
	44	esign = 0;
	45	decpt = 0;
	46	dexponent = 0;
	47	Acc = znumber;
	48	Acc.num_tag = radix;
	49	acc = CH_FIELD(Acc);
	50	temp = CH_FIELD(Temp);
	51
	52	do{
	53	ch = *str++;
	54	} while(isspace(ch));
	55
	56	switch(ch){
	57	case '-':
	58	msign = -1;
	59	/* FALLTHROUGH */
	60	case '+':
	61	ch = *str++;
	62	break;
	63	}
	64	dofract:
	65	for(; isdigit(ch); ch = *str++){
	66	assert(((acc[HOC] & SIGNBIT) == 0), "Negative HOC");
	67	if (acc[HOC] < MAXINT_10){
	68	ovf = numshift(3, temp, acc);
	69	ovf \|= numshift(1, acc, acc);
70	ovf \|= numaddv(acc, temp, acc);
71	ovf \|= numaddd(acc, acc, ch - '0');
72	assert(ovf == 0, "Overflow building mantissa");
73	} else {
74	/*
75	* Then, the number is too large anyway
76	*/
77	dexponent++;
78	}
79	if (decpt)
80	dexponent--;
81	}
82	switch(ch){
83	case '.':
84	if (decpt == 0){
85	decpt++;
86	ch = *str++;
87	goto dofract;
88	}
89	break;
90	/*
91	* only 'e' and 'E' are recognized by atof()
92	*/
93	case 'e':
94	case 'E':
95	/*
96	* we include the remainder for compatability with as formats
97	* in as, the radix actual paramater agrees with the character
98	* we expect; consequently, no checking is done.
99	*/
100	case 'd':
101	case 'D':
102	case 'g':
103	case 'G':
104	case 'h':
105	case 'H':
106	j = 0;
107	ch = *str++;
108	esign = 0;
109	switch(ch){
110	case '-':
111	esign = 1;
112	/* FALLTHROUGH */
113	case '+':
114	ch = *str++;
115	}
116	for(; isdigit(ch); ch = *str++){
117	if (j < MAXINT_10){
118	j *= 10;
119	j += ch - '0';
120	} else {
121	/*
122	* outrageously large exponent
123	*/
124	/VOID/
125	}
126	}
127	if (esign)
128	dexponent -= j;
129	else
130	dexponent += j;
131	/*
132	* There should be a range check on dexponent here
133	*/
134	}
135	/*
136	* The number has now been reduced to a mantissa
137	* and an exponent.
138	* The mantissa is an n bit number (to the precision
139	* of the extended words) in the acc.
140	* The exponent is a signed power of 10 in dexponent.
141	* msign is on if the resulting number will eventually
142	* be negative.
143	*
144	* We now must convert the number to standard format floating
145	* number, which will be done by accumulating
146	* a binary exponent in bexponent, as we gradually
147	* drive dexponent towards zero, one count at a time.
148	*/
149	if (isclear(acc)){
150	return(Acc);
151	}
152	bexponent = 0;
153
154	/*
155	* Scale the number down.
156	* We must divide acc by 10 as many times as needed.
157	*/
158	for (; dexponent < 0; dexponent++){
159	/*
160	* Align the number so that the most significant
161	* bits are aligned in the most significant
162	* bits of the accumulator, adjusting the
163	* binary exponent as we shift.
164	* The goal is to get the high order bit (NOT the
165	* sign bit) set.
166	*/
167	assert(((acc[HOC] & SIGNBIT) == 0), "Negative HOC");
168	ovf = 0;
169
170	for (j = 5; j >= 1; --j){
171	i = 1 << (j - 1); /* 16, 8, 4, 2, 1 */
172	quotient = ONES(i);
173	quotient <<= (CH_BITS - 1) - i;
174	while((acc[HOC] & quotient) == 0){
175	ovf \|= numshift((int)i, acc, acc);
176	bexponent -= i;
177	}
178	}
179	/*
180	* Add 2 to the accumulator to effect rounding,
181	* and get set up to divide by 5.
182	*/
183	ovf = numaddd(acc, acc, 2);
184	assert(ovf == 0, "Carry out of left rounding up by 2");
185	/*
186	* Divide the high order chunks by 5;
187	* The last chunk will be divided by 10,
188	* (to see what the remainder is, also to effect rounding)
189	* and then multipiled by 2 to effect division by 5.
190	*/
191	remainder = 0;
192	#if DEBUGNATOF
193	printf("Dividing: ");
194	bignumprint(Acc);
195	printf("\n");
196	#endif DEBUGNATOF
197	sptr = (short *)acc;
198	for (i = (CH_N * 2 - 1); i >= 1; --i){
199	/*
200	* Divide (remainder:16).(acc[i]:16)
201	* by 5, putting the quotient back
202	* into acc[i]:16, and save the remainder
203	* for the next iteration.
204	*/
205	dividend = (remainder << 16) \| (sptr[i] & ONES(16));
206	assert(dividend >= 0, "dividend < 0");
207	quotient = dividend / 5;
208	remainder = dividend - (quotient * 5);
209	sptr[i] = quotient;
210	remainder = remainder;
211	}
212	/*
213	* Divide the lowest order chunk by 10,
214	* saving the remainder to decide how to round.
215	* Then, multiply by 2, making it look as
216	* if we divided by 10.
217	* This multiply fills in a 0 on the least sig bit.
218	*/
219	dividend = (remainder << 16) \| (sptr[0] & ONES(16));
220	assert(dividend >= 0, "dividend < 0");
221	quotient = dividend / 10;
222	remainder = dividend - (quotient * 10);
223	sptr[0] = quotient + quotient;
224
225	if (remainder >= 5)
226	ovf = numaddd(acc, acc, 1);
227	/*
228	* Now, divide by 2, effecting division by 10,
229	* merely by adjusting the binary exponent.
230	*/
231	bexponent--;
232	}
233	/*
234	* Scale the number up by multiplying by 10 as
235	* many times as necessary
236	*/
237	for (; dexponent > 0; dexponent--){
238	/*
239	* Compare high word to (2**31)/5,
240	* and scale accordingly
241	*/
242	while ( ((unsigned)acc[HOC]) > MAXINT_5){
243	(void)numshift(-1, acc, acc);
244	bexponent++;
245	}
246	/*
247	* multiply the mantissa by 5,
248	* and scale the binary exponent by 2
249	*/
250	ovf = numshift(2, temp, acc);
251	ovf \|= numaddv(acc, acc, temp);
252	assert(ovf == 0, "Scaling * 10 of manitissa");
253	bexponent++;
254	}
255	/*
256	* We now have:
257	* a CH_N chunk length binary integer, right
258	* justified (in native format).
259	* a binary exponent.
260	*
261	* Now, we treat this large integer as an octa word
262	* number, and unpack it into standard unpacked
263	* format. That unpacking will give us yet
264	* another binary exponent, which we adjust with
265	* the accumulated binary exponent.
266	*/
267	Acc.num_tag = TYPO;
268	#if DEBUGNATOF
269	printf("Octal number: ");
270	bignumprint(Acc);
271	printf("\n");
272	#endif DEBUGNATOF
273	Acc = bignumunpack(Acc, &ovf);
274
275	if (ovf)
276	errno = ERANGE;
277	#if DEBUGNATOF
278	printf("Unpacked octal number: ");
279	bignumprint(Acc);
280	printf("bexponent == %d\n", bexponent);
281	#endif DEBUGNATOF
282	Acc.num_exponent += bexponent;
283	assert(Acc.num_sign == 0, "unpacked integer is < 0");
284	Acc.num_sign = msign;
285	/*
286	* We now pack the number back into a radix format number.
287	* This checks for overflow, underflow,
288	* and rounds by 1/2 ulp.
289	*/
290	ovf = 0;
291	Acc = bignumpack(Acc, radix, &ovf);
292	if (ovf)
293	errno = ERANGE;
294	#if DEBUGNATOF
295	printf("packed number: ");
296	bignumprint(Acc);
297	printf("\n");
298	#endif DEBUGNATOF
299	return(Acc);
300	}
301	#if 0
302	/*
303	* Unfortunately, one can't use the ediv instruction to do
304	* division on numbers with > 64 bits.
305	* This is because ediv returns signed quantities;
306	* if the quotient is an unsigned number > 2^31,
307	* ediv sets integer overflow.
308	*/
309	unsigned int ediv(high, low, divisor, qp, i)
310	register unsigned int high; /* r11 */
311	register unsigned int low; /* r10 */
312	register unsigned int divisor; /* r9 */
313	unsigned int *qp;
314	{
315	register unsigned int remainder; /* r8 */
316	register unsigned int quotient; /* r7 */
317
318	asm("ediv r9, r10, r7, r8 # Divide. q->r7, r->r8 (discarded)");
319	*qp = quotient;
320	return(remainder);
321	}
322	#endif 0