|
|
1.1 root 1: /*
2: * Copyright (c) 1982 Regents of the University of California
3: */
4: #ifndef lint
5: static char sccsid[] = "@(#)natof.c 4.3 2/14/82";
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
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.