ext_comp.c

minux的源代码,一个非常小的操作系统

	{ 0,	1395,	0x940919BB,	0xD4620B6D },
	{ 0,	1488,	0x957A4AE1,	0xEBF7F3D4 },
	{ 0,	1581,	0x96EF14C6,	0x454AA840 },
	{ 0,	1674,	0x98678061,	0x27ECE4F5 },
	{ 0,	1767,	0x99E396C1,	0x3A3ACFF2 }
};

static struct EXTEND r_ten_powers[] = { /* representation of 10 ** -i */
	{ 0,	0,	0x80000000,	0 },
	{ 0,	-4,	0xCCCCCCCC,	0xCCCCCCCD },
	{ 0,	-7,	0xA3D70A3D,	0x70A3D70A },
	{ 0,	-10,	0x83126E97,	0x8D4FDF3B },
	{ 0,	-14,	0xD1B71758,	0xE219652C },
	{ 0,	-17,	0xA7C5AC47,	0x1B478423 },
	{ 0,	-20,	0x8637BD05,	0xAF6C69B6 },
	{ 0,	-24,	0xD6BF94D5,	0xE57A42BC },
	{ 0,	-27,	0xABCC7711,	0x8461CEFD },
	{ 0,	-30,	0x89705F41,	0x36B4A597 },
	{ 0,	-34,	0xDBE6FECE,	0xBDEDD5BF },
	{ 0,	-37,	0xAFEBFF0B,	0xCB24AAFF },
	{ 0,	-40,	0x8CBCCC09,	0x6F5088CC },
	{ 0,	-44,	0xE12E1342,	0x4BB40E13 },
	{ 0,	-47,	0xB424DC35,	0x095CD80F },
	{ 0,	-50,	0x901D7CF7,	0x3AB0ACD9 },
	{ 0,	-54,	0xE69594BE,	0xC44DE15B },
	{ 0,	-57,	0xB877AA32,	0x36A4B449 },
	{ 0,	-60,	0x9392EE8E,	0x921D5D07 },
	{ 0,	-64,	0xEC1E4A7D,	0xB69561A5 },
	{ 0,	-67,	0xBCE50864,	0x92111AEB },
	{ 0,	-70,	0x971DA050,	0x74DA7BEF },
	{ 0,	-74,	0xF1C90080,	0xBAF72CB1 },
	{ 0,	-77,	0xC16D9A00,	0x95928A27 },
	{ 0,	-80,	0x9ABE14CD,	0x44753B53 },
	{ 0,	-84,	0xF79687AE,	0xD3EEC551 },
	{ 0,	-87,	0xC6120625,	0x76589DDB },
	{ 0,	-90,	0x9E74D1B7,	0x91E07E48 }
};

static struct EXTEND r_big_ten_powers[] = { /* representation of 10 ** -(28*i) */
	{ 0,	0,	0x80000000,	0 },
	{ 0,	-94,	0xFD87B5F2,	0x8300CA0E },
	{ 0,	-187,	0xFB158592,	0xBE068D2F },
	{ 0,	-280,	0xF8A95FCF,	0x88747D94 },
	{ 0,	-373,	0xF64335BC,	0xF065D37D },
	{ 0,	-466,	0xF3E2F893,	0xDEC3F126 },
	{ 0,	-559,	0xF18899B1,	0xBC3F8CA2 },
	{ 0,	-652,	0xEF340A98,	0x172AACE5 },
	{ 0,	-745,	0xECE53CEC,	0x4A314EBE },
	{ 0,	-838,	0xEA9C2277,	0x23EE8BCB },
	{ 0,	-931,	0xE858AD24,	0x8F5C22CA },
	{ 0,	-1024,	0xE61ACF03,	0x3D1A45DF },
	{ 0,	-1117,	0xE3E27A44,	0x4D8D98B8 },
	{ 0,	-1210,	0xE1AFA13A,	0xFBD14D6E },
	{ 0,	-1303,	0xDF82365C,	0x497B5454 },
	{ 0,	-1396,	0xDD5A2C3E,	0xAB3097CC },
	{ 0,	-1489,	0xDB377599,	0xB6074245 },
	{ 0,	-1582,	0xD91A0545,	0xCDB51186 },
	{ 0,	-1675,	0xD701CE3B,	0xD387BF48 },
	{ 0,	-1768,	0xD4EEC394,	0xD6258BF8 }
};

#define	TP	(int)(sizeof(ten_powers)/sizeof(ten_powers[0]))
#define BTP	(int)(sizeof(big_ten_powers)/sizeof(big_ten_powers[0]))
#define MAX_EXP	(TP * BTP - 1)

static
add_exponent(struct EXTEND *e, int exp)
{
	int neg = exp < 0;
	int divsz, modsz;
	struct EXTEND x;

	if (neg) exp = -exp;
	divsz = exp / TP;
	modsz = exp % TP;
	if (neg) {
		mul_ext(e, &r_ten_powers[modsz], &x);
		mul_ext(&x, &r_big_ten_powers[divsz], e);
	}
	else {
		mul_ext(e, &ten_powers[modsz], &x);
		mul_ext(&x, &big_ten_powers[divsz], e);
	}
}

_str_ext_cvt(const char *s, char **ss, struct EXTEND *e)
{
	/*	Like strtod, but for extended precision */
	register int	c;
	int		dotseen = 0;
	int		digitseen = 0;
	int		exp = 0;

	if (ss) *ss = (char *)s;
	while (isspace(*s)) s++;

	e->sign = 0;
	e->exp = 0;
	e->m1 = e->m2 = 0;

	c = *s;
	switch(c) {
	case '-':
		e->sign = 1;
	case '+':
		s++;
	}
	while (c = *s++, isdigit(c) || (c == '.' && ! dotseen++)) {
		if (c == '.') continue;
		digitseen = 1;
		if (e->m1 <= (unsigned long)(0xFFFFFFFF)/10) {
			struct mantissa	a1;

			a1 = e->mantissa;
			b64_sft(&(e->mantissa), -3);
			b64_sft(&a1, -1);
			b64_add(&(e->mantissa), &a1);
			a1.h_32 = 0;
			a1.l_32 = c - '0';
			b64_add(&(e->mantissa), &a1);
		}
		else exp++;
		if (dotseen) exp--;
	}
	if (! digitseen) return;

	if (ss) *ss = (char *)s - 1;

	if (c == 'E' || c == 'e') {
		int	exp1 = 0;
		int	sign = 1;
		int	exp_overflow = 0;

		switch(*s) {
		case '-':
			sign = -1;
		case '+':
			s++;
		}
		if (c = *s, isdigit(c)) {
			do {
				int tmp;

				exp1 = 10 * exp1 + (c - '0');
				if ((tmp = sign * exp1 + exp) > MAX_EXP ||
				     tmp < -MAX_EXP) {
					exp_overflow = 1;
				}
			} while (c = *++s, isdigit(c));
			if (ss) *ss = (char *)s;
		}
		exp += sign * exp1;
		if (exp_overflow) {
			exp = sign * MAX_EXP;
			if (e->m1 != 0 || e->m2 != 0) errno = ERANGE;
		}
	}
	if (e->m1 == 0 && e->m2 == 0) return;
	e->exp = 63;
	while (! (e->m1 & 0x80000000)) {
		b64_sft(&(e->mantissa),-1);
		e->exp--;
	}
	add_exponent(e, exp);
}

#include	<math.h>

static
ten_mult(struct EXTEND *e)
{
	struct EXTEND e1 = *e;

	e1.exp++;
	e->exp += 3;
	add_ext(e, &e1, e);
}

#define NDIGITS 128
#define NSIGNIFICANT 19

char *
_ext_str_cvt(struct EXTEND *e, int ndigit, int *decpt, int *sign, int ecvtflag)
{
	/*	Like cvt(), but for extended precision */

	static char buf[NDIGITS+1];
	struct EXTEND m;
	register char *p = buf;
	register char *pe;
	int findex = 0;

	if (ndigit < 0) ndigit = 0;
	if (ndigit > NDIGITS) ndigit = NDIGITS;
	pe = &buf[ndigit];
	buf[0] = '\0';

	*sign = 0;
	if (e->sign) {
		*sign = 1;
		e->sign = 0;
	}

	*decpt = 0;
	if (e->m1 != 0) {
		register struct EXTEND *pp = &big_ten_powers[1];

		while(cmp_ext(e,pp) >= 0) {
			pp++;
			findex = pp - big_ten_powers;
			if (findex >= BTP) break;
		}
		pp--;
		findex = pp - big_ten_powers;
		mul_ext(e,&r_big_ten_powers[findex],e);
		*decpt += findex * TP;
		pp = &ten_powers[1];
		while(pp < &ten_powers[TP] && cmp_ext(e, pp) >= 0) pp++;
		pp--;
		findex = pp - ten_powers;
		*decpt += findex;

		if (cmp_ext(e, &ten_powers[0]) < 0) {
			pp = &r_big_ten_powers[1];
			while(cmp_ext(e,pp) < 0) pp++;
			pp--;
			findex = pp - r_big_ten_powers;
			mul_ext(e, &big_ten_powers[findex], e);
			*decpt -= findex * TP;
			/* here, value >= 10 ** -28 */
			ten_mult(e);
			(*decpt)--;
			pp = &r_ten_powers[0];
			while(cmp_ext(e, pp) < 0) pp++;
			findex = pp - r_ten_powers;
			mul_ext(e, &ten_powers[findex], e);
			*decpt -= findex;
			findex = 0;
		}
		(*decpt)++;	/* because now value in [1.0, 10.0) */
	}
	if (! ecvtflag) {
		/* for fcvt() we need ndigit digits behind the dot */
		pe += *decpt;
		if (pe > &buf[NDIGITS]) pe = &buf[NDIGITS];
	}
	m.exp = -62;
	m.sign = 0;
	m.m1 = 0xA0000000;
	m.m2 = 0;
	while (p <= pe) {
		struct EXTEND oneminm;

		if (p - pe > NSIGNIFICANT) {
			findex = 0;
			e->m1 = 0;
		}
		if (findex) {
			struct EXTEND tc, oldtc;
			int count = 0;

			oldtc.exp = 0;
			oldtc.sign = 0;
			oldtc.m1 = 0;
			oldtc.m2 = 0;
			tc = ten_powers[findex];
			while (cmp_ext(e, &tc) >= 0) {
				oldtc = tc;
				add_ext(&tc, &ten_powers[findex], &tc);
				count++;
			}
			*p++ = count + '0';
			oldtc.sign = 1;
			add_ext(e, &oldtc, e);
			findex--;
			continue;
		}
		if (e->m1) {
			m.sign = 1;
			add_ext(&ten_powers[0], &m, &oneminm);
			m.sign = 0;
			if (e->exp >= 0) {
				struct EXTEND x;

				x.m2 = 0; x.exp = e->exp;
				x.sign = 1;
				x.m1 = e->m1>>(31-e->exp);
				*p++ = (x.m1) + '0';
				x.m1 = x.m1 << (31-e->exp);
				add_ext(e, &x, e);
			}
			else *p++ = '0';
			/* Check that remainder is still significant */
			if (cmp_ext(&m, e) > 0 || cmp_ext(e, &oneminm) > 0) {
				if (e->m1 && e->exp >= -1) *(p-1) += 1;
				e->m1 = 0;
				continue;
			}
			ten_mult(&m);
			ten_mult(e);
		}
		else *p++ = '0';
	}
	if (pe >= buf) {
		p = pe;
		*p += 5;	/* round of at the end */
		while (*p > '9') {
			*p = '0';
			if (p > buf) ++*--p;
			else {
				*p = '1';
				++*decpt;
				if (! ecvtflag) {
					/* maybe add another digit at the end,
					   because the point was shifted right
					*/
					if (pe > buf) *pe = '0';
					pe++;
				}
			}
		}
		*pe = '\0';
	}
	return buf;
}

_dbl_ext_cvt(double value, struct EXTEND *e)
{
	/*	Convert double to extended
	*/
	int exponent;

	value = frexp(value, &exponent);
	e->sign = value < 0.0;
	if (e->sign) value = -value;
	e->exp = exponent - 1;
	value *= 4294967296.0;
	e->m1 = value;
	value -= e->m1;
	value *= 4294967296.0;
	e->m2 = value;
}

static struct EXTEND max_d;

double
_ext_dbl_cvt(struct EXTEND *e)
{
	/*	Convert extended to double
	*/
	double f;
	int sign = e->sign;

	e->sign = 0;
	if (e->m1 == 0 && e->m2 == 0) {
		return 0.0;
	}
	if (max_d.exp == 0) {
		_dbl_ext_cvt(DBL_MAX, &max_d);
	}
	if (cmp_ext(&max_d, e) < 0) {
		f = HUGE_VAL;
		errno = ERANGE;
	}
	else	f = ldexp((double)e->m1*4294967296.0 + (double)e->m2, e->exp-63);
	if (sign) f = -f;
	if (f == 0.0 && (e->m1 != 0 || e->m2 != 0)) {
		errno = ERANGE;
	}
	return f;
}