xref: /original-bsd/old/as.vax/natof.c (revision 81f6297c)
1 /*
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.
5  */
6 
7 #ifndef lint
8 static char sccsid[] = "@(#)natof.c	5.1 (Berkeley) 04/30/85";
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
326