|
|
1.1 root 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) 4/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
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.