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

/*
 * Copyright (c) 1982 Regents of the University of California.
 * All rights reserved.  The Berkeley software License Agreement
 * specifies the terms and conditions for redistribution.
 */

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