bignum2.c

<B>Digital的Unix操作系统VAX 4.2源码</B>

/*
 *	Copyright (c) 1982 Regents of the University of California
 */
#ifndef lint
static char sccsid[] = "@(#)bignum2.c 4.3 2/14/82";
#endif not lint

#include <stdio.h>
#include "as.h"
Bignum	Znumber;			/* zero reference */
#define	MINEXP		-32768		/* never generate; reserved for id 0 */

Bignum intconvert(number, convto)
	Bignum	number;
	int	convto;
{
	reg	int	i;
	if (number.num_tag == convto)
		return(number);
	if (ty_nbyte[number.num_tag] > ty_nbyte[convto] && (passno == 2)){
		yywarning("Conversion between %s and %s looses significance",
			ty_string[number.num_tag],
			ty_string[convto]);
	}
	for (i = ty_nbyte[convto]; i < ty_nbyte[TYPO]; i++)
		number.num_uchar[i] = 0;
	return(number);
}

#define	CONV(src, dst)	(((src) << TYPLG) + (dst))

Bignum floatconvert(number, convto)
	Bignum	number;
	int	convto;
{
	reg	u_int	*bp;		/* r11 */
		int	loss = 0;
		int	gain = 0;
		int	mixs = 0;
		Ovf	ovf;

	if (number.num_tag == convto)
		return(number);
	bp = &number.num_uint[0];
#ifdef lint
	*bp = *bp;
#endif lint

	switch(CONV(number.num_tag, convto)){

	case CONV(TYPF, TYPD):	asm("cvtfd (r11), (r11)"); break;
	case CONV(TYPF, TYPG):	mixs++; break;
	case CONV(TYPF, TYPH):	mixs++; break;

	case CONV(TYPD, TYPF):	asm("cvtdf (r11), (r11)"); break;
	case CONV(TYPD, TYPG):	mixs++; break;
	case CONV(TYPD, TYPH):	mixs++; break;

	case CONV(TYPG, TYPF):	mixs++; break;
	case CONV(TYPG, TYPD):	mixs++; break;
	case CONV(TYPG, TYPH):	mixs++; break;

	case CONV(TYPH, TYPF):	mixs++; break;
	case CONV(TYPH, TYPD):	mixs++; break;
	case CONV(TYPH, TYPG):	mixs++; break;
	default:	panic("Bad floating point conversion?");
	}
	if ((gain || mixs || loss) && (passno == 2)){
		yywarning("Converting from %s to %s: %s ",
			ty_string[number.num_tag],
			ty_string[convto],
			gain ? "gains significance" :
			(loss ? "looses significance" : "mixs exponent formats")
		);
	}
	if (mixs){
		number = bignumconvert(number, convto, &ovf);
		if (ovf && passno == 2){
			yywarning("Floating conversion over/underflowed\n");
		}
	} else {
		number.num_tag = convto;
	}
	return(number);
}

/* 
 *	Convert a big number between various representations
 */
Bignum bignumconvert(number, toconv, ovfp)
	Bignum	number;
	int	toconv;
	Ovf	*ovfp;
{
	int	tag;

	*ovfp = 0;
	tag = number.num_tag;
	if (tag == toconv)
		return(number);
	if (tag == TYPUNPACKED){
		return(bignumpack(number, toconv, ovfp));
	}
	if (toconv == TYPUNPACKED){
		return(bignumunpack(number, ovfp));
	}
	return(bignumpack(bignumunpack(number, ovfp), toconv, ovfp));
}

Bignum bignumunpack(Packed, ovfp)
	Bignum	Packed;
		Ovf	*ovfp;
{
		Bignum	Mantissa;
		Bignum	Enumber;
	reg	int	i;
		int	j;
		int	k;
	reg	struct	ty_bigdesc	*p;
	reg	chptr	packed;
	reg	chptr	mantissa;
	reg	chptr	enumber;
		u_short	exponent;
		int	sign;
		int	mask;

	p = &ty_bigdesc[Packed.num_tag];

	*ovfp = 0;
	Mantissa = Znumber;
	sign = 0;
	exponent = 0;
	mantissa = CH_FIELD(Mantissa);
	enumber = CH_FIELD(Enumber);
	packed = CH_FIELD(Packed);

	if (isclear(packed)){
		Mantissa.num_tag = TYPUNPACKED;
		Mantissa.num_exponent = MINEXP;
		return(Mantissa);
	}
	/*
	 *	map the packed number into the mantissa, using
	 *	the unpacking map
	 */
	mapnumber(mantissa, packed, 16, p->b_upmmap);
	/*
	 *	perform the mantissa shifting.
	 *	This may appear to overflow; all that is lost
	 *	is low order bits of the exponent.
	 */
	(void)numshift(p->b_mlshift, mantissa, mantissa);
	/*
	 *	handle sign and exponent
	 */
	switch(Packed.num_tag){
	case TYPB:
	case TYPW:
	case TYPL:
	case TYPO:
	case TYPQ:
		sign = 0;
		exponent = p->b_eexcess;
		if (mantissa[HOC] & SIGNBIT){
			sign = -1;
			*ovfp |= numnegate(mantissa, mantissa);
		}
		/*
		 *	Normalize the packed by left shifting,
		 *	adjusting the exponent as we go.
		 *	Do a binary weighted left shift for some speed.
		 */
		k = 0;
		for (j = 4; j >= 0; --j){
			i = 1 << j;	/* 16, 8, 4, 2, 1 */
			while(1){
				if (k >= p->b_msigbits)
					break;
				mask = ONES(i) << (CH_BITS - i);
				if (mantissa[HOC] & mask)
					break;
				(void)numshift(i, mantissa, mantissa);
				k += i;
				exponent -= i;
			}
		}
		assert(mantissa[HOC] & SIGNBIT, "integer <<ing");
		/*
		 *	now, kick the most significant bit off the top
		 */
		(void)numshift(1, mantissa, mantissa);
		break;
	default:
		/*
		 *	map the exponent into the local area.
		 */
		Enumber = Znumber;
		mapnumber(enumber, packed, 2, p->b_upemap);
		/*
		 *	Extract the exponent, and get rid
		 *	of the sign bit
		 */
		exponent = Enumber.num_ushort[0] & ONES(15);
		/*
		 *	shift the exponent, and get rid of high order
		 *	trash
		 */
		exponent >>= p->b_ershift;
		exponent &= ONES(p->b_esigbits);
		/*
		 *	un excess the exponent
		 */
		exponent -= p->b_eexcess;
		/*
		 *	extract and extend the sign bit
		 */
		sign = (Enumber.num_ushort[0] & ~ONES(15)) ? -1 : 0;
	}
	/*
	 *	Assemble the pieces, and return the number
	 */
	Mantissa.num_tag = TYPUNPACKED;
	Mantissa.num_sign = sign;
	Mantissa.num_exponent = exponent;
	return(Mantissa);
}

Bignum bignumpack(Unpacked, toconv, ovfp)
	Bignum	Unpacked;
	int	toconv;
	Ovf	*ovfp;
{
	Bignum	Back;
	Bignum	Enumber;
	Bignum	Temp;

		short	exponent;
		char	sign;
	reg	struct	ty_bigdesc	*p;
	reg	chptr	back;
	reg	chptr	enumber;
	reg	chptr	temp;
	reg	chptr	unpacked;

		int	i,j;

	if (Unpacked.num_tag != TYPUNPACKED)
		panic("Argument to bignumpack is not unpacked");

	*ovfp = 0;
	Back = Znumber;
	Temp = Znumber;
	Back.num_tag = toconv;

	back = CH_FIELD(Back);
	temp = CH_FIELD(Temp);
	enumber = CH_FIELD(Enumber);
	unpacked = CH_FIELD(Unpacked);
	p = &ty_bigdesc[toconv];

	exponent = Unpacked.num_exponent;
	sign = Unpacked.num_sign;
	if (exponent == MINEXP)
		return(Back);	/* identically zero */

	switch(toconv){
	case TYPB:
	case TYPW:
	case TYPL:
	case TYPQ:
	case TYPO:
		/*
		 *	Put back in the assumed high order fraction
		 *	bit that is always a 1.
		 */
		(void)numshift(-1, temp, unpacked);
		temp[HOC] |= SIGNBIT;
		if (exponent > p->b_eexcess){
			/*
			 *	Construct the largest positive integer
			 */
			(void)numclear(temp);
			(void)num1comp(temp, temp);
			temp[HOC] &= ~SIGNBIT;
			sign = sign;
			*ovfp |= OVF_OVERFLOW;
		} else
		if (exponent <= 0){
			/*
			 *	chop the temp; underflow to integer 0
			 */
			(void)numclear(temp);
			sign = 0;
			*ovfp |= OVF_UNDERFLOW;
		} else {
			/*
			 *	denormalize the temp.
			 *	This will again chop, by shifting
			 *	bits off the right end into oblivion.
			 */
			for (j = 4; j >= 0; --j){
				i = 1 << j;	/* 16, 8, 4, 2, 1 */
				while(exponent + i <= p->b_eexcess){
					numshift(-i, temp, temp);
					exponent += i;
				}
			}
		}
		/*
		 *	negate the temp if the sign is set
		 */
		if (sign)
			*ovfp |= numnegate(temp, temp);
		/*
		 *	Stuff the temp number into the return area
		 */
		mapnumber(back, temp, 16, p->b_pmmap);
		return(Back);
	default:
		/*
		 *	Shift the mantissa to the right, filling in zeroes on
		 *	the left.  This aligns the least significant bit
		 *	on the bottom of a byte, something that upround
		 *	will use.
		 *	Put the result into a temporary.
		 *	Even though the shift may be zero, there
		 *	is still a copy involved here.
		 */
		(void)numshift(-(p->b_mlshift), temp, unpacked);
		/*
		 *	Perform the rounding by adding in 0.5 ulp's
		 */
		exponent = upround(&Temp, p, exponent);
		/*
		 *	Do a range check on the exponent, in preparation
		 *	to stuffing it in.
		 */
		if ((short)(exponent + p->b_eexcess) == 0){
			/*
			 *	Sorry, no gradual underflow on the
			 *	VAX.  Chop this beasty totally to zero
			 */
			goto zeroret;
		} else
		if ((short)(exponent + p->b_eexcess) < 0){
			/*
			 *	True underflow will happen;
			 *	Chop everything to positive zero
			 */
		  zeroret:
			(void)numclear(temp);
			exponent = 0;
			sign = 0;	/* avoid reserved operand! */
			*ovfp |= OVF_UNDERFLOW;
		} else
		if ((unsigned)(exponent + p->b_eexcess)
		    >= (unsigned)(1 << p->b_esigbits)){
			/*
			 *	Construct the largest magnitude possible
			 *	floating point unpacked: 0.{1}111111111
			 */
			(void)numclear(temp);
			(void)num1comp(temp, temp);
			exponent = ONES(p->b_esigbits);
			sign = sign;
			*ovfp |= OVF_OVERFLOW;
		} else {
			/*
			 *	The exponent will fit.
			 *	Bias it up, and the common code will stuff it.
			 */
			exponent += p->b_eexcess;
		}
		exponent <<= p->b_ershift;
		/*
		 *	mask out trash for the sign, and put in the sign.
		 */
		exponent &= ONES(15);
		if (sign)
			exponent |= ~ONES(15);
		Enumber.num_ushort[0] = exponent;
		/*
		 *	Map the unpacked exponent into the value going back
		 */
		mapnumber(back, enumber, 2, p->b_pemap);
		/*
		 *	Stuff the unpacked mantissa into the return area
		 */
		mapnumber(back, temp, 16, p->b_pmmap);
		return(Back);
	}
	/*NOTREACHED*/
}

mapnumber(chp1, chp2, nbytes, themap)
		chptr	chp1, chp2;
		int	nbytes;
		char	*themap;
{
	reg	int	i;
	reg	u_char	*p1, *p2;

	p1 = (u_char *)chp1;
	p2 = (u_char *)chp2;
	for (i = 0; i < nbytes; i++){
		switch(themap[i]){
		case NOTAKE:
			break;
		default:
			p1[themap[i]] |= p2[i];
			break;
		}
	}
}

#define	UPSHIFT	2
/*
 *	round in 1/2 ulp in the number, possibly modifying
 *	the binary exponent if there was total carry out.
 *	Return the modified exponent
 */
int upround(numberp, p, exponent)
	reg	Bignum	*numberp;
	reg	struct	ty_bigdesc	*p;
		int	exponent;
{
	reg	u_char	*bytep;
		int	nbytes;
		int	byteindex;
		int	hofractionbit;
		int	ovffractionbit;
	reg	int	ovfbitindex;
	reg	chptr	number;
	static	Bignum	ulp;

	/*
	 *	Find out the byte index of the byte containing the ulp
	 */
	number = CH_FIELD(numberp[0]);
	bytep = numberp->num_uchar;

	nbytes = (p->b_msigbits - 1) + p->b_mlshift;
	assert((nbytes % 8) == 0, "mantissa sig bits");
	nbytes /= 8;
	byteindex = 15 - nbytes;
	assert(byteindex >= 0, "ulp in outer space");
	/*
	 *	Shift the number to the right by two places,
	 *	so that we can do full arithmetic without overflowing
	 *	to the left.
	 */
	numshift(-UPSHIFT, number, number);
	/*
	 *	Construct the missing high order fraction bit
	 */
	ovfbitindex = 8 - (p->b_mlshift + UPSHIFT);
	assert(ovfbitindex >= 0, "Shifted byte 15 into byte 14");
	hofractionbit = (0x01 << ovfbitindex);
	ovffractionbit = (0x02 << ovfbitindex);
	bytep[15] |= hofractionbit;
	/*
	 *	construct the unit in the last place, and it
	 *	to the fraction
	 */
	ulp.num_uchar[byteindex] |= (0x80 >> UPSHIFT);
	numaddv(number, number, CH_FIELD(ulp));
	ulp.num_uchar[byteindex] &= ~(0x80 >> UPSHIFT);
	/*
	 *	Check if there was an overflow,
	 *	and adjust by shifting.
	 *	Also, bring the number back into canonical
	 *	unpacked form by left shifting by two to undeo
	 *	what we did before.
	 */
	if (bytep[15] & ovffractionbit){
		exponent += 1;
		numshift(UPSHIFT - 1, number, number);
	} else {
		numshift(UPSHIFT, number, number);
	}
	/*
	 *	Clear off trash in the unused bits of the high
	 *	order byte of the number
	 */
	bytep[15] &= ONES(8 - p->b_mlshift);
	return(exponent);
}
#ifdef DEBUG
bignumprint(number)
	Bignum	number;
{
	printf("Bignum: %s (exp: %d, sign: %d) 0x%08x%08x%08x%08x",
		ty_string[number.num_tag],
		number.num_exponent,
		number.num_sign,
		number.num_uint[3],
		number.num_uint[2],
		number.num_uint[1],
		number.num_uint[0]);
	switch(number.num_tag){
	case TYPB:
	case TYPW:
	case TYPL:
	case TYPQ:
	case TYPO:
	case TYPUNPACKED:
		break;
	case TYPF:
		printf(" == %10.8e", number.num_num.numFf_float.Ff_value);
		break;
	case TYPD:
		printf(" == %20.17e", number.num_num.numFd_float.Fd_value);
		break;
	case TYPG:
	case TYPH:
		break;
	}
}

numprintovf(ovf)
	Ovf	ovf;
{
	int	i;
	static struct	ovftab{
		Ovf	which;
		char	*print;
	} ovftab[] = {
		OVF_POSOVF,	"posovf",
		OVF_MAXINT,	"maxint",
		OVF_ADDV,	"addv",
		OVF_LSHIFT,	"lshift",
		OVF_F,		"F float",
		OVF_D,		"D float",
		OVF_G,		"G float",
		OVF_H,		"H float",
		OVF_OVERFLOW,	"cvt overflow",
		OVF_UNDERFLOW,	"cvt underflow",
		0,		0
	};
	for(i = 0; ovftab[i].which; i++){
		if (ovf & ovftab[i].which)
			printf("Overflow(%s) ", ovftab[i].print);
	}
}
#endif DEBUG