xref: /original-bsd/old/as.vax/bignum2.c (revision fbed46ce)
1 /*
2  *	Copyright (c) 1982 Regents of the University of California
3  */
4 #ifndef lint
5 static char sccsid[] = "@(#)bignum2.c 4.3 02/14/82";
6 #endif not lint
7 
8 #include <stdio.h>
9 #include "as.h"
10 Bignum	Znumber;			/* zero reference */
11 #define	MINEXP		-32768		/* never generate; reserved for id 0 */
12 
13 Bignum intconvert(number, convto)
14 	Bignum	number;
15 	int	convto;
16 {
17 	reg	int	i;
18 	if (number.num_tag == convto)
19 		return(number);
20 	if (ty_nbyte[number.num_tag] > ty_nbyte[convto] && (passno == 2)){
21 		yywarning("Conversion between %s and %s looses significance",
22 			ty_string[number.num_tag],
23 			ty_string[convto]);
24 	}
25 	for (i = ty_nbyte[convto]; i < ty_nbyte[TYPO]; i++)
26 		number.num_uchar[i] = 0;
27 	return(number);
28 }
29 
30 #define	CONV(src, dst)	(((src) << TYPLG) + (dst))
31 
32 Bignum floatconvert(number, convto)
33 	Bignum	number;
34 	int	convto;
35 {
36 	reg	u_int	*bp;		/* r11 */
37 		int	loss = 0;
38 		int	gain = 0;
39 		int	mixs = 0;
40 		Ovf	ovf;
41 
42 	if (number.num_tag == convto)
43 		return(number);
44 	bp = &number.num_uint[0];
45 #ifdef lint
46 	*bp = *bp;
47 #endif lint
48 
49 	switch(CONV(number.num_tag, convto)){
50 
51 	case CONV(TYPF, TYPD):	asm("cvtfd (r11), (r11)"); break;
52 	case CONV(TYPF, TYPG):	mixs++; break;
53 	case CONV(TYPF, TYPH):	mixs++; break;
54 
55 	case CONV(TYPD, TYPF):	asm("cvtdf (r11), (r11)"); break;
56 	case CONV(TYPD, TYPG):	mixs++; break;
57 	case CONV(TYPD, TYPH):	mixs++; break;
58 
59 	case CONV(TYPG, TYPF):	mixs++; break;
60 	case CONV(TYPG, TYPD):	mixs++; break;
61 	case CONV(TYPG, TYPH):	mixs++; break;
62 
63 	case CONV(TYPH, TYPF):	mixs++; break;
64 	case CONV(TYPH, TYPD):	mixs++; break;
65 	case CONV(TYPH, TYPG):	mixs++; break;
66 	default:	panic("Bad floating point conversion?");
67 	}
68 	if ((gain || mixs || loss) && (passno == 2)){
69 		yywarning("Converting from %s to %s: %s ",
70 			ty_string[number.num_tag],
71 			ty_string[convto],
72 			gain ? "gains significance" :
73 			(loss ? "looses significance" : "mixs exponent formats")
74 		);
75 	}
76 	if (mixs){
77 		number = bignumconvert(number, convto, &ovf);
78 		if (ovf && passno == 2){
79 			yywarning("Floating conversion over/underflowed\n");
80 		}
81 	} else {
82 		number.num_tag = convto;
83 	}
84 	return(number);
85 }
86 
87 /*
88  *	Convert a big number between various representations
89  */
90 Bignum bignumconvert(number, toconv, ovfp)
91 	Bignum	number;
92 	int	toconv;
93 	Ovf	*ovfp;
94 {
95 	int	tag;
96 
97 	*ovfp = 0;
98 	tag = number.num_tag;
99 	if (tag == toconv)
100 		return(number);
101 	if (tag == TYPUNPACKED){
102 		return(bignumpack(number, toconv, ovfp));
103 	}
104 	if (toconv == TYPUNPACKED){
105 		return(bignumunpack(number, ovfp));
106 	}
107 	return(bignumpack(bignumunpack(number, ovfp), toconv, ovfp));
108 }
109 
110 Bignum bignumunpack(Packed, ovfp)
111 	Bignum	Packed;
112 		Ovf	*ovfp;
113 {
114 		Bignum	Mantissa;
115 		Bignum	Enumber;
116 	reg	int	i;
117 		int	j;
118 		int	k;
119 	reg	struct	ty_bigdesc	*p;
120 	reg	chptr	packed;
121 	reg	chptr	mantissa;
122 	reg	chptr	enumber;
123 		u_short	exponent;
124 		int	sign;
125 		int	mask;
126 
127 	p = &ty_bigdesc[Packed.num_tag];
128 
129 	*ovfp = 0;
130 	Mantissa = Znumber;
131 	sign = 0;
132 	exponent = 0;
133 	mantissa = CH_FIELD(Mantissa);
134 	enumber = CH_FIELD(Enumber);
135 	packed = CH_FIELD(Packed);
136 
137 	if (isclear(packed)){
138 		Mantissa.num_tag = TYPUNPACKED;
139 		Mantissa.num_exponent = MINEXP;
140 		return(Mantissa);
141 	}
142 	/*
143 	 *	map the packed number into the mantissa, using
144 	 *	the unpacking map
145 	 */
146 	mapnumber(mantissa, packed, 16, p->b_upmmap);
147 	/*
148 	 *	perform the mantissa shifting.
149 	 *	This may appear to overflow; all that is lost
150 	 *	is low order bits of the exponent.
151 	 */
152 	(void)numshift(p->b_mlshift, mantissa, mantissa);
153 	/*
154 	 *	handle sign and exponent
155 	 */
156 	switch(Packed.num_tag){
157 	case TYPB:
158 	case TYPW:
159 	case TYPL:
160 	case TYPO:
161 	case TYPQ:
162 		sign = 0;
163 		exponent = p->b_eexcess;
164 		if (mantissa[HOC] & SIGNBIT){
165 			sign = -1;
166 			*ovfp |= numnegate(mantissa, mantissa);
167 		}
168 		/*
169 		 *	Normalize the packed by left shifting,
170 		 *	adjusting the exponent as we go.
171 		 *	Do a binary weighted left shift for some speed.
172 		 */
173 		k = 0;
174 		for (j = 4; j >= 0; --j){
175 			i = 1 << j;	/* 16, 8, 4, 2, 1 */
176 			while(1){
177 				if (k >= p->b_msigbits)
178 					break;
179 				mask = ONES(i) << (CH_BITS - i);
180 				if (mantissa[HOC] & mask)
181 					break;
182 				(void)numshift(i, mantissa, mantissa);
183 				k += i;
184 				exponent -= i;
185 			}
186 		}
187 		assert(mantissa[HOC] & SIGNBIT, "integer <<ing");
188 		/*
189 		 *	now, kick the most significant bit off the top
190 		 */
191 		(void)numshift(1, mantissa, mantissa);
192 		break;
193 	default:
194 		/*
195 		 *	map the exponent into the local area.
196 		 */
197 		Enumber = Znumber;
198 		mapnumber(enumber, packed, 2, p->b_upemap);
199 		/*
200 		 *	Extract the exponent, and get rid
201 		 *	of the sign bit
202 		 */
203 		exponent = Enumber.num_ushort[0] & ONES(15);
204 		/*
205 		 *	shift the exponent, and get rid of high order
206 		 *	trash
207 		 */
208 		exponent >>= p->b_ershift;
209 		exponent &= ONES(p->b_esigbits);
210 		/*
211 		 *	un excess the exponent
212 		 */
213 		exponent -= p->b_eexcess;
214 		/*
215 		 *	extract and extend the sign bit
216 		 */
217 		sign = (Enumber.num_ushort[0] & ~ONES(15)) ? -1 : 0;
218 	}
219 	/*
220 	 *	Assemble the pieces, and return the number
221 	 */
222 	Mantissa.num_tag = TYPUNPACKED;
223 	Mantissa.num_sign = sign;
224 	Mantissa.num_exponent = exponent;
225 	return(Mantissa);
226 }
227 
228 Bignum bignumpack(Unpacked, toconv, ovfp)
229 	Bignum	Unpacked;
230 	int	toconv;
231 	Ovf	*ovfp;
232 {
233 	Bignum	Back;
234 	Bignum	Enumber;
235 	Bignum	Temp;
236 
237 		short	exponent;
238 		char	sign;
239 	reg	struct	ty_bigdesc	*p;
240 	reg	chptr	back;
241 	reg	chptr	enumber;
242 	reg	chptr	temp;
243 	reg	chptr	unpacked;
244 
245 		int	i,j;
246 
247 	if (Unpacked.num_tag != TYPUNPACKED)
248 		panic("Argument to bignumpack is not unpacked");
249 
250 	*ovfp = 0;
251 	Back = Znumber;
252 	Temp = Znumber;
253 	Back.num_tag = toconv;
254 
255 	back = CH_FIELD(Back);
256 	temp = CH_FIELD(Temp);
257 	enumber = CH_FIELD(Enumber);
258 	unpacked = CH_FIELD(Unpacked);
259 	p = &ty_bigdesc[toconv];
260 
261 	exponent = Unpacked.num_exponent;
262 	sign = Unpacked.num_sign;
263 	if (exponent == MINEXP)
264 		return(Back);	/* identically zero */
265 
266 	switch(toconv){
267 	case TYPB:
268 	case TYPW:
269 	case TYPL:
270 	case TYPQ:
271 	case TYPO:
272 		/*
273 		 *	Put back in the assumed high order fraction
274 		 *	bit that is always a 1.
275 		 */
276 		(void)numshift(-1, temp, unpacked);
277 		temp[HOC] |= SIGNBIT;
278 		if (exponent > p->b_eexcess){
279 			/*
280 			 *	Construct the largest positive integer
281 			 */
282 			(void)numclear(temp);
283 			(void)num1comp(temp, temp);
284 			temp[HOC] &= ~SIGNBIT;
285 			sign = sign;
286 			*ovfp |= OVF_OVERFLOW;
287 		} else
288 		if (exponent <= 0){
289 			/*
290 			 *	chop the temp; underflow to integer 0
291 			 */
292 			(void)numclear(temp);
293 			sign = 0;
294 			*ovfp |= OVF_UNDERFLOW;
295 		} else {
296 			/*
297 			 *	denormalize the temp.
298 			 *	This will again chop, by shifting
299 			 *	bits off the right end into oblivion.
300 			 */
301 			for (j = 4; j >= 0; --j){
302 				i = 1 << j;	/* 16, 8, 4, 2, 1 */
303 				while(exponent + i <= p->b_eexcess){
304 					numshift(-i, temp, temp);
305 					exponent += i;
306 				}
307 			}
308 		}
309 		/*
310 		 *	negate the temp if the sign is set
311 		 */
312 		if (sign)
313 			*ovfp |= numnegate(temp, temp);
314 		/*
315 		 *	Stuff the temp number into the return area
316 		 */
317 		mapnumber(back, temp, 16, p->b_pmmap);
318 		return(Back);
319 	default:
320 		/*
321 		 *	Shift the mantissa to the right, filling in zeroes on
322 		 *	the left.  This aligns the least significant bit
323 		 *	on the bottom of a byte, something that upround
324 		 *	will use.
325 		 *	Put the result into a temporary.
326 		 *	Even though the shift may be zero, there
327 		 *	is still a copy involved here.
328 		 */
329 		(void)numshift(-(p->b_mlshift), temp, unpacked);
330 		/*
331 		 *	Perform the rounding by adding in 0.5 ulp's
332 		 */
333 		exponent = upround(&Temp, p, exponent);
334 		/*
335 		 *	Do a range check on the exponent, in preparation
336 		 *	to stuffing it in.
337 		 */
338 		if ((short)(exponent + p->b_eexcess) == 0){
339 			/*
340 			 *	Sorry, no gradual underflow on the
341 			 *	VAX.  Chop this beasty totally to zero
342 			 */
343 			goto zeroret;
344 		} else
345 		if ((short)(exponent + p->b_eexcess) < 0){
346 			/*
347 			 *	True underflow will happen;
348 			 *	Chop everything to positive zero
349 			 */
350 		  zeroret:
351 			(void)numclear(temp);
352 			exponent = 0;
353 			sign = 0;	/* avoid reserved operand! */
354 			*ovfp |= OVF_UNDERFLOW;
355 		} else
356 		if ((unsigned)(exponent + p->b_eexcess)
357 		    >= (unsigned)(1 << p->b_esigbits)){
358 			/*
359 			 *	Construct the largest magnitude possible
360 			 *	floating point unpacked: 0.{1}111111111
361 			 */
362 			(void)numclear(temp);
363 			(void)num1comp(temp, temp);
364 			exponent = ONES(p->b_esigbits);
365 			sign = sign;
366 			*ovfp |= OVF_OVERFLOW;
367 		} else {
368 			/*
369 			 *	The exponent will fit.
370 			 *	Bias it up, and the common code will stuff it.
371 			 */
372 			exponent += p->b_eexcess;
373 		}
374 		exponent <<= p->b_ershift;
375 		/*
376 		 *	mask out trash for the sign, and put in the sign.
377 		 */
378 		exponent &= ONES(15);
379 		if (sign)
380 			exponent |= ~ONES(15);
381 		Enumber.num_ushort[0] = exponent;
382 		/*
383 		 *	Map the unpacked exponent into the value going back
384 		 */
385 		mapnumber(back, enumber, 2, p->b_pemap);
386 		/*
387 		 *	Stuff the unpacked mantissa into the return area
388 		 */
389 		mapnumber(back, temp, 16, p->b_pmmap);
390 		return(Back);
391 	}
392 	/*NOTREACHED*/
393 }
394 
395 mapnumber(chp1, chp2, nbytes, themap)
396 		chptr	chp1, chp2;
397 		int	nbytes;
398 		char	*themap;
399 {
400 	reg	int	i;
401 	reg	u_char	*p1, *p2;
402 
403 	p1 = (u_char *)chp1;
404 	p2 = (u_char *)chp2;
405 	for (i = 0; i < nbytes; i++){
406 		switch(themap[i]){
407 		case NOTAKE:
408 			break;
409 		default:
410 			p1[themap[i]] |= p2[i];
411 			break;
412 		}
413 	}
414 }
415 
416 #define	UPSHIFT	2
417 /*
418  *	round in 1/2 ulp in the number, possibly modifying
419  *	the binary exponent if there was total carry out.
420  *	Return the modified exponent
421  */
422 int upround(numberp, p, exponent)
423 	reg	Bignum	*numberp;
424 	reg	struct	ty_bigdesc	*p;
425 		int	exponent;
426 {
427 	reg	u_char	*bytep;
428 		int	nbytes;
429 		int	byteindex;
430 		int	hofractionbit;
431 		int	ovffractionbit;
432 	reg	int	ovfbitindex;
433 	reg	chptr	number;
434 	static	Bignum	ulp;
435 
436 	/*
437 	 *	Find out the byte index of the byte containing the ulp
438 	 */
439 	number = CH_FIELD(numberp[0]);
440 	bytep = numberp->num_uchar;
441 
442 	nbytes = (p->b_msigbits - 1) + p->b_mlshift;
443 	assert((nbytes % 8) == 0, "mantissa sig bits");
444 	nbytes /= 8;
445 	byteindex = 15 - nbytes;
446 	assert(byteindex >= 0, "ulp in outer space");
447 	/*
448 	 *	Shift the number to the right by two places,
449 	 *	so that we can do full arithmetic without overflowing
450 	 *	to the left.
451 	 */
452 	numshift(-UPSHIFT, number, number);
453 	/*
454 	 *	Construct the missing high order fraction bit
455 	 */
456 	ovfbitindex = 8 - (p->b_mlshift + UPSHIFT);
457 	assert(ovfbitindex >= 0, "Shifted byte 15 into byte 14");
458 	hofractionbit = (0x01 << ovfbitindex);
459 	ovffractionbit = (0x02 << ovfbitindex);
460 	bytep[15] |= hofractionbit;
461 	/*
462 	 *	construct the unit in the last place, and it
463 	 *	to the fraction
464 	 */
465 	ulp.num_uchar[byteindex] |= (0x80 >> UPSHIFT);
466 	numaddv(number, number, CH_FIELD(ulp));
467 	ulp.num_uchar[byteindex] &= ~(0x80 >> UPSHIFT);
468 	/*
469 	 *	Check if there was an overflow,
470 	 *	and adjust by shifting.
471 	 *	Also, bring the number back into canonical
472 	 *	unpacked form by left shifting by two to undeo
473 	 *	what we did before.
474 	 */
475 	if (bytep[15] & ovffractionbit){
476 		exponent += 1;
477 		numshift(UPSHIFT - 1, number, number);
478 	} else {
479 		numshift(UPSHIFT, number, number);
480 	}
481 	/*
482 	 *	Clear off trash in the unused bits of the high
483 	 *	order byte of the number
484 	 */
485 	bytep[15] &= ONES(8 - p->b_mlshift);
486 	return(exponent);
487 }
488 #ifdef DEBUG
489 bignumprint(number)
490 	Bignum	number;
491 {
492 	printf("Bignum: %s (exp: %d, sign: %d) 0x%08x%08x%08x%08x",
493 		ty_string[number.num_tag],
494 		number.num_exponent,
495 		number.num_sign,
496 		number.num_uint[3],
497 		number.num_uint[2],
498 		number.num_uint[1],
499 		number.num_uint[0]);
500 	switch(number.num_tag){
501 	case TYPB:
502 	case TYPW:
503 	case TYPL:
504 	case TYPQ:
505 	case TYPO:
506 	case TYPUNPACKED:
507 		break;
508 	case TYPF:
509 		printf(" == %10.8e", number.num_num.numFf_float.Ff_value);
510 		break;
511 	case TYPD:
512 		printf(" == %20.17e", number.num_num.numFd_float.Fd_value);
513 		break;
514 	case TYPG:
515 	case TYPH:
516 		break;
517 	}
518 }
519 
520 numprintovf(ovf)
521 	Ovf	ovf;
522 {
523 	int	i;
524 	static struct	ovftab{
525 		Ovf	which;
526 		char	*print;
527 	} ovftab[] = {
528 		OVF_POSOVF,	"posovf",
529 		OVF_MAXINT,	"maxint",
530 		OVF_ADDV,	"addv",
531 		OVF_LSHIFT,	"lshift",
532 		OVF_F,		"F float",
533 		OVF_D,		"D float",
534 		OVF_G,		"G float",
535 		OVF_H,		"H float",
536 		OVF_OVERFLOW,	"cvt overflow",
537 		OVF_UNDERFLOW,	"cvt underflow",
538 		0,		0
539 	};
540 	for(i = 0; ovftab[i].which; i++){
541 		if (ovf & ovftab[i].which)
542 			printf("Overflow(%s) ", ovftab[i].print);
543 	}
544 }
545 #endif DEBUG
546