ext_comp.c

操作系统源代码

/*
  (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
  See the copyright notice in the ACK home directory, in the file "Copyright".
*/

/* $Id: ext_comp.c,v 1.10 1994/06/24 11:53:36 ceriel Exp $ */

/* extended precision arithmetic for the strtod() and cvt() routines */

/* This may require some more work when long doubles get bigger than 8
   bytes. In this case, these routines may become obsolete. ???
*/

#include	"ext_fmt.h"
#include	<float.h>
#include	<errno.h>
#include	<ctype.h>

static int b64_add(struct mantissa *e1, struct mantissa *e2);
static b64_sft(struct mantissa *e1, int n);

static
mul_ext(struct EXTEND *e1, struct EXTEND *e2, struct EXTEND *e3)
{
	/*	Multiply the extended numbers e1 and e2, and put the
		result in e3.
	*/
	register int	i,j;		/* loop control	*/
	unsigned short	mp[4];
	unsigned short	mc[4];
	unsigned short	result[8];	/* result */

	register unsigned short *pres;

	/* first save the sign (XOR)			*/
	e3->sign = e1->sign ^ e2->sign;

	/* compute new exponent */
	e3->exp = e1->exp + e2->exp + 1;

	/* check for overflow/underflow	??? */

	/* 128 bit multiply of mantissas	*/

	/* assign unknown long formats		*/
	/* to known unsigned word formats	*/
	mp[0] = e1->m1 >> 16;
	mp[1] = (unsigned short) e1->m1;
	mp[2] = e1->m2 >> 16;
	mp[3] = (unsigned short) e1->m2;
	mc[0] = e2->m1 >> 16;
	mc[1] = (unsigned short) e2->m1;
	mc[2] = e2->m2 >> 16;
	mc[3] = (unsigned short) e2->m2;
	for (i = 8; i--;) {
		result[i] = 0;
	}
	/*
	 *	fill registers with their components
	 */
	for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
		unsigned short k = 0;
		unsigned long mpi = mp[i];
		for(j=4;j--;) {
			unsigned long tmp = (unsigned long)pres[j] + k;
			if (mc[j]) tmp += mpi * mc[j];
			pres[j] = tmp;
			k = tmp >> 16;
		}
		pres[-1] = k;
	}

	if (! (result[0] & 0x8000)) {
		e3->exp--;
		for (i = 0; i <= 3; i++) {
			result[i] <<= 1;
			if (result[i+1]&0x8000) result[i] |= 1;
		}
		result[4] <<= 1;
	}	
	/*
	 *	combine the registers to a total
	 */
	e3->m1 = ((unsigned long)(result[0]) << 16) + result[1];
	e3->m2 = ((unsigned long)(result[2]) << 16) + result[3];
	if (result[4] & 0x8000) {
		if (++e3->m2 == 0) {
			if (++e3->m1 == 0) {
				e3->m1 = 0x80000000;
				e3->exp++;
			}
		}
	}
}

static
add_ext(struct EXTEND *e1, struct EXTEND *e2, struct EXTEND *e3)
{
	/*	Add two extended numbers e1 and e2, and put the result
		in e3
	*/
	struct EXTEND ce2;
	int diff;

	if ((e2->m1 | e2->m2) == 0L) {
		*e3 = *e1;
		return;
	}
	if ((e1->m1 | e1->m2) == 0L) {
		*e3 = *e2;
		return;
	}
	ce2 = *e2;
	*e3 = *e1;
	e1 = &ce2;

	/* adjust mantissas to equal power */
	diff = e3->exp - e1->exp;
	if (diff < 0) {
		diff = -diff;
		e3->exp += diff;
		b64_sft(&(e3->mantissa), diff);
	}
	else if (diff > 0) {
		e1->exp += diff;
		b64_sft(&(e1->mantissa), diff);
	}
	if (e1->sign != e3->sign) {
		/* e3 + e1 = e3 - (-e1) */
		if (e1->m1 > e3->m1 ||
                    (e1->m1 == e3->m1 && e1->m2 > e3->m2)) {
                	/*      abs(e1) > abs(e3) */
                	if (e3->m2 > e1->m2) {
                        	e1->m1 -= 1;    /* carry in */
                	}
                	e1->m1 -= e3->m1;
                	e1->m2 -= e3->m2;
                	*e3 = *e1;
        	}
        	else {
                	if (e1->m2 > e3->m2)
                        	e3->m1 -= 1;    /* carry in */
                	e3->m1 -= e1->m1;
                	e3->m2 -= e1->m2;
        	}
	}
	else {
		if (b64_add(&e3->mantissa,&e1->mantissa)) {/* addition carry */
			b64_sft(&e3->mantissa,1);/* shift mantissa one bit RIGHT */
			e3->m1 |= 0x80000000L;	/* set max bit	*/
			e3->exp++;		/* increase the exponent */
		}
	}
	if ((e3->m2 | e3->m1) != 0L) {
		/* normalize */
		if (e3->m1 == 0L) {
			e3->m1 = e3->m2; e3->m2 = 0L; e3->exp -= 32;
		}
		if (!(e3->m1 & 0x80000000)) {
			unsigned long l = 0x40000000;
			int cnt = -1;

			while (! (l & e3->m1)) {
				l >>= 1; cnt--;
			}
			e3->exp += cnt;
			b64_sft(&(e3->mantissa), cnt);
		}
	}
}

static int
cmp_ext(struct EXTEND *e1, struct EXTEND *e2)
{
        struct EXTEND tmp; 
         
        e2->sign = ! e2->sign; 
        add_ext(e1, e2, &tmp);  
        e2->sign = ! e2->sign;
        if (tmp.m1 == 0 && tmp.m2 == 0) return 0; 
        if (tmp.sign) return -1;
        return 1;
}

static
b64_sft(struct mantissa *e1, int n)
{
	if (n > 0) {
		if (n > 63) {
			e1->l_32 = 0;
			e1->h_32 = 0;
			return;
		}
		if (n >= 32) {
			e1->l_32 = e1->h_32;
			e1->h_32 = 0;
			n -= 32;
		}
		if (n > 0) {
			e1->l_32 >>= n;
			if (e1->h_32 != 0) {
				e1->l_32 |= (e1->h_32 << (32 - n));
				e1->h_32 >>= n;
			}
		}
		return;
	}
	n = -n;
	if (n > 0) {
		if (n > 63) {
			e1->l_32 = 0;
			e1->h_32 = 0;
			return;
		}
		if (n >= 32) {
			e1->h_32 = e1->l_32;
			e1->l_32 = 0;
			n -= 32;
		}
		if (n > 0) {
			e1->h_32 <<= n;
			if (e1->l_32 != 0) {
				e1->h_32 |= (e1->l_32 >> (32 - n));
				e1->l_32 <<= n;
			}
		}
	}
}

static int
b64_add(struct mantissa *e1, struct mantissa *e2)
		/*
		 * pointers to 64 bit 'registers'
		 */
{
	register int	overflow;
	int		carry;

			/* add higher pair of 32 bits */
	overflow = ((unsigned long) 0xFFFFFFFF - e1->h_32 < e2->h_32);
	e1->h_32 += e2->h_32;

			/* add lower pair of 32 bits */
	carry = ((unsigned long) 0xFFFFFFFF - e1->l_32 < e2->l_32);
	e1->l_32 += e2->l_32;
	if ((carry) && (++e1->h_32 == 0))
		return(1);		/* had a 64 bit overflow */
	else
		return(overflow);	/* return status from higher add */
}

/* The following tables can be computed with the following bc(1)
   program:

obase=16
scale=0
define t(x){
	auto a, b, c
	a=2;b=1;c=2^32;n=1
	while(a<x) {
		b=a;n+=n;a*=a
	}
	n/=2
	a=b
	while(b<x) {
		a=b;b*=c;n+=32
	}
	n-=32
	b=a
	while(a<x) {
		b=a;a+=a;n+=1
	}
	n-=1
	x*=16^16
	b=x%a
	x/=a
	if(a<=(2*b)) x+=1
	obase=10
	n
	obase=16
	return(x)
}
for (i=1;i<28;i++) {
	t(10^i)
}
0
for (i=1;i<20;i++) {
	t(10^(28*i))
}
0
define r(x){
	auto a, b, c
	a=2;b=1;c=2^32;n=1
	while(a<x) {
		b=a;n+=n;a*=a
	}
	n/=2
	a=b
	while(b<x) {
		a=b;b*=c;n+=32
	}
	n-=32
	b=a
	while(a<x) {
		b=a;a+=a;n+=1
	}
	a=b
	a*=16^16
	b=a%x
	a/=x
	if(x<=(2*b)) a+=1
	obase=10
	-n
	obase=16
	return(a)
}
for (i=1;i<28;i++) {
	r(10^i)
}
0
for (i=1;i<20;i++) {
	r(10^(28*i))
}
0

*/
static struct EXTEND ten_powers[] = {	/* representation of 10 ** i */
	{ 0,	0,	0x80000000,	0 },
	{ 0,	3,	0xA0000000,	0 },
	{ 0,	6,	0xC8000000,	0 },
	{ 0,	9,	0xFA000000,	0 },
	{ 0,	13,	0x9C400000,	0 },
	{ 0,	16,	0xC3500000,	0 },
	{ 0,	19,	0xF4240000,	0 },
	{ 0,	23,	0x98968000,	0 },
	{ 0,	26,	0xBEBC2000,	0 },
	{ 0,	29,	0xEE6B2800,	0 },
	{ 0,	33,	0x9502F900,	0 },
	{ 0,	36,	0xBA43B740,	0 },
	{ 0,	39,	0xE8D4A510,	0 },
	{ 0,	43,	0x9184E72A,	0 },
	{ 0,	46,	0xB5E620F4,	0x80000000 },
	{ 0,	49,	0xE35FA931,	0xA0000000 },
	{ 0,	53,	0x8E1BC9BF,	0x04000000 },
	{ 0,	56,	0xB1A2BC2E,	0xC5000000 },
	{ 0,	59,	0xDE0B6B3A,	0x76400000 },
	{ 0,	63,	0x8AC72304,	0x89E80000 },
	{ 0,	66,	0xAD78EBC5,	0xAC620000 },
	{ 0,	69,	0xD8D726B7,	0x177A8000 },
	{ 0,	73,	0x87867832,	0x6EAC9000 },
	{ 0,	76,	0xA968163F,	0x0A57B400 },
	{ 0,	79,	0xD3C21BCE,	0xCCEDA100 },
	{ 0,	83,	0x84595161,	0x401484A0 },
	{ 0,	86,	0xA56FA5B9,	0x9019A5C8 },
	{ 0,	89,	0xCECB8F27,	0xF4200F3A }
};
static struct EXTEND big_ten_powers[] = {  /* representation of 10 ** (28*i) */
	{ 0,	0,	0x80000000,	0 },
	{ 0,	93,	0x813F3978,	0xF8940984 },
	{ 0,	186,	0x82818F12,	0x81ED44A0 },
	{ 0,	279,	0x83C7088E,	0x1AAB65DB },
	{ 0,	372,	0x850FADC0,	0x9923329E },
	{ 0,	465,	0x865B8692,	0x5B9BC5C2 },
	{ 0,	558,	0x87AA9AFF,	0x79042287 },
	{ 0,	651,	0x88FCF317,	0xF22241E2 },
	{ 0,	744,	0x8A5296FF,	0xE33CC930 },
	{ 0,	837,	0x8BAB8EEF,	0xB6409C1A },
	{ 0,	930,	0x8D07E334,	0x55637EB3 },
	{ 0,	1023,	0x8E679C2F,	0x5E44FF8F },
	{ 0,	1116,	0x8FCAC257,	0x558EE4E6 },
	{ 0,	1209,	0x91315E37,	0xDB165AA9 },
	{ 0,	1302,	0x929B7871,	0xDE7F22B9 },