# define INCLUDE_CTYPE # include "dgd.h" # include "xfloat.h" typedef struct { unsigned short sign; /* 0: positive, 0x8000: negative */ unsigned short exp; /* bias: 32767 */ unsigned short high; /* 0 / 1 / 14 bits */ Uint low; /* 0 / 29 bits / 0 / 0 */ } flt; # define NBITS 44 # define BIAS 0x7fff /* constants */ static flt half = { 0x0000, 0x7ffe, 0x4000, 0x00000000L }; /* .5 */ static flt one = { 0x0000, 0x7fff, 0x4000, 0x00000000L }; /* 1 */ /* * NAME: f_edom() * DESCRIPTION: Domain error */ static void f_edom() { error("Math argument"); } /* * NAME: f_erange() * DESCRIPTION: Out of range */ static void f_erange() { error("Result too large"); } static void f_sub P((flt*, flt*)); /* * NAME: f_add() * DESCRIPTION: a = a + b. The result is normalized, but not guaranteed to * be in range. */ static void f_add(a, b) register flt *a, *b; { register unsigned short h, n; register Uint l; flt tmp; if (b->exp == 0) { /* b is 0 */ return; } if (a->exp == 0) { /* a is 0 */ *a = *b; return; } if (a->sign != b->sign) { a->sign ^= 0x8000; f_sub(a, b); a->sign ^= 0x8000; return; } if (a->exp < b->exp) { /* b is the largest; exchange a and b */ tmp = *a; *a = *b; b = &tmp; } n = a->exp - b->exp; if (n <= NBITS) { h = a->high; l = a->low; /* * perform addition */ if (n < 31) { h += (Uint) b->high >> n; l += (((Uint) b->high << (31 - n)) & 0x7fffffffL) | (b->low >> n); } else { l += b->high >> (n - 31); } if ((Int) l < 0) { /* carry */ l &= 0x7fffffffL; h++; } if ((short) h < 0) { /* too large */ l |= (Uint) h << 31; l >>= 1; h >>= 1; a->exp++; } /* * rounding off */ if ((Int) (l += 2) < 0 && (short) ++h < 0) { h >>= 1; a->exp++; } l &= 0x7ffffffcL; a->high = h; a->low = l; } } /* * NAME: f_sub() * DESCRIPTION: a = a - b. The result is normalized, but not guaranteed to be * in range. */ static void f_sub(a, b) register flt *a, *b; { register unsigned short h, n; register Uint l; flt tmp; if (b->exp == 0) { /* b is 0 */ return; } if (a->exp == 0) { *a = *b; a->sign ^= 0x8000; return; } if (a->sign != b->sign) { a->sign ^= 0x8000; f_add(a, b); a->sign ^= 0x8000; return; } if (a->exp <= b->exp && (a->exp < b->exp || (a->high <= b->high && (a->high < b->high || a->low < b->low)))) { /* b is the largest; exchange a and b */ tmp = *a; *a = *b; b = &tmp; a->sign ^= 0x8000; } n = a->exp - b->exp; if (n <= NBITS + 1) { h = a->high; l = a->low; /* * perform subtraction */ if (n < 31) { h -= (Uint) b->high >> n; l -= (((Uint) b->high << (31 - n)) & 0x7fffffffL) | (b->low >> n); if ((b->low >> n) << n != b->low) { --l; } } else { n -= 31; l -= b->high >> n; if (b->low != 0 || (b->high >> n) << n != b->high) { --l; } } if ((Int) l < 0) { /* borrow */ l &= 0x7fffffffL; --h; } /* * normalize */ if (h == 0) { if (l == 0) { a->exp = 0; return; } n = 0; if ((l & 0xffff0000L) == 0) { l <<= 15; n += 15; } h = l >> 16; l <<= 15; l &= 0x7fffffffL; a->exp -= n + 15; } if (h < 0x4000) { n = 0; if ((h & 0xff00) == 0) { h <<= 7; n += 7; } while (h < 0x4000) { h <<= 1; n++; } h |= l >> (31 - n); l <<= n; l &= 0x7fffffffL; a->exp -= n; } /* * rounding off */ if ((Int) (l += 2) < 0 && (short) ++h < 0) { h >>= 1; a->exp++; } l &= 0x7ffffffcL; a->high = h; a->low = l; } } /* * NAME: f_mult() * DESCRIPTION: a = a * b. The result is normalized, but may be out of range. */ static void f_mult(a, b) register flt *a, *b; { register Uint m, l; register unsigned short al, am, ah, bl, bm, bh; if (a->exp == 0) { /* a is 0 */ return; } if (b->exp == 0) { /* b is 0 */ a->exp = 0; return; } al = ((unsigned short) a->low) >> 1; bl = ((unsigned short) b->low) >> 1; am = a->low >> 16; bm = b->low >> 16; ah = a->high; bh = b->high; m = (Uint) al * bl; m >>= 15; m += (Uint) al * bm + (Uint) am * bl; m >>= 15; m += (Uint) al * bh + (Uint) am * bm + (Uint) ah * bl; m >>= 13; l = m & 0x03; m >>= 2; m += (Uint) am * bh + (Uint) ah * bm; l |= (m & 0x7fff) << 2; m >>= 15; m += (Uint) ah * bh; l |= m << 17; ah = m >> 14; a->sign ^= b->sign; a->exp += b->exp - BIAS; if ((short) ah < 0) { ah >>= 1; l >>= 1; a->exp++; } l &= 0x7fffffffL; /* * rounding off */ if ((Int) (l += 2) < 0 && (short) ++ah < 0) { ah >>= 1; a->exp++; } l &= 0x7ffffffcL; a->high = ah; a->low = l; } /* * NAME: f_div() * DESCRIPTION: c = a / b. b must be non-zero. The result is normalized, * but may be out of range. */ static void f_div(a, b) register flt *a, *b; { unsigned short n[3]; register Uint numh, numl, divl, high, low, q; register unsigned short divh, i; if (b->exp == 0) { error("Division by zero"); } if (a->exp == 0) { /* a is 0 */ return; } numh = ((Uint) a->high << 16) | (a->low >> 15); numl = a->low << 17; divh = (b->high << 1) | (b->low >> 30); divl = b->low << 2; n[0] = 0; n[1] = 0; i = 3; do { /* estimate the high word of the quotient */ q = numh / divh; /* highlow = num * q */ low = (unsigned short) (high = q * (unsigned short) divl); high >>= 16; high += q * (divl >> 16); low |= high << 16; high >>= 16; high += q * divh; /* the estimated quotient may be 2 off; correct it if needed */ if (high >= numh && (high > numh || low > numl)) { high -= divh; if (low < divl) { --high; } low -= divl; --q; if (high >= numh && (high > numh || low > numl)) { high -= divh; if (low < divl) { --high; } low -= divl; --q; } } n[--i] = q; if (i == 0) { break; } /* subtract highlow */ numh -= high; if (numl < low) { --numh; } numl -= low; numh <<= 16; numh |= numl >> 16; numl <<= 16; } while (numh != 0 || numl != 0); a->sign ^= b->sign; a->exp -= b->exp - BIAS + 1; high = n[2]; low = ((Uint) n[1] << 15) | (n[0] >> 1); if ((short) high < 0) { low |= high << 31; low >>= 1; high >>= 1; a->exp++; } /* * rounding off */ if ((Int) (low += 2) < 0 && (short) ++high < 0) { high >>= 1; a->exp++; } low &= 0x7ffffffcL; a->high = high; a->low = low; } /* * NAME: f_trunc() * DESCRIPTION: truncate a flt */ static void f_trunc(a) register flt *a; { static unsigned short maskh[] = { 0x4000, 0x6000, 0x7000, 0x7800, 0x7c00, 0x7e00, 0x7f00, 0x7f80, 0x7fc0, 0x7fe0, 0x7ff0, 0x7ff8, 0x7ffc, 0x7ffe, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff }; static Uint maskl[] = { 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L, 0x40000000L, 0x60000000L, 0x70000000L, 0x78000000L, 0x7c000000L, 0x7e000000L, 0x7f000000L, 0x7f800000L, 0x7fc00000L, 0x7fe00000L, 0x7ff00000L, 0x7ff80000L, 0x7ffc0000L, 0x7ffe0000L, 0x7fff0000L, 0x7fff8000L, 0x7fffc000L, 0x7fffe000L, 0x7ffff000L, 0x7ffff800L, 0x7ffffc00L, 0x7ffffe00L, 0x7fffff00L, 0x7fffff80L, 0x7fffffc0L, 0x7fffffe0L, 0x7ffffff0L, 0x7ffffff8L }; if (a->exp < BIAS) { a->exp = 0; } else if (a->exp < BIAS + NBITS) { a->high &= maskh[a->exp - BIAS]; a->low &= maskl[a->exp - BIAS]; } } /* * NAME: f_37bits() * DESCRIPTION: round a flt to 37 binary digits of precision */ static void f_37bits(a) register flt *a; { if ((Int) (a->low += 0x100) < 0) { a->low = 0; if ((short) ++(a->high) < 0) { a->high >>= 1; a->exp++; } } else { a->low &= 0xfffffe00L; } } /* * NAME: f_round() * DESCRIPTION: round off a flt */ static void f_round(a) register flt *a; { static flt half = { 0, 0x7ffe, 0x4000, 0x00000000L }; half.sign = a->sign; f_add(a, &half); f_trunc(a); } /* * NAME: f_ldexp() * DESCRIPTION: add an integer value to the exponent of a flt */ static void f_ldexp(a, exp) register flt *a; register short exp; { if (a->exp == 0) { return; } exp += a->exp; if (exp <= 0) { a->exp = 0; return; } a->exp = exp; } /* * NAME: f_cmp() * DESCRIPTION: compate two flts */ static int f_cmp(a, b) register flt *a, *b; { if (a->exp == 0) { if (b->exp == 0) { return 0; } return (b->sign == 0) ? -1 : 1; } else if (b->exp == 0) { if (a->exp == 0) { return 0; } return (a->sign == 0) ? 1 : -1; } else if (a->sign != b->sign) { return (a->sign == 0) ? 1 : -1; } else { if (a->exp == b->exp && a->high == b->high && a->low == b->low) { return 0; } if (a->exp <= b->exp && (a->exp < b->exp || (a->high <= b->high && (a->high < b->high || a->low < b->low)))) { return (a->sign == 0) ? -1 : 1; } return (a->sign == 0) ? 1 : -1; } } /* * NAME: f_itof() * DESCRIPTION: convert an integer to a flt */ static void f_itof(i, a) Int i; register flt *a; { register Uint n; register unsigned short shift; /* deal with zero and sign */ if (i == 0) { a->exp = 0; return; } else if (i < 0) { a->sign = 0x8000; n = -i; } else { a->sign = 0; n = i; } shift = 0; while ((n & 0xff000000L) == 0) { n <<= 8; shift += 8; } while ((Int) n >= 0) { n <<= 1; shift++; } a->exp = BIAS + 31 - shift; a->high = n >> 17; a->low = (n << 14) & 0x7fffffff; } /* * NAME: f_ftoi() * DESCRIPTION: convert a flt to an integer, disregarding the fractional part */ static Int f_ftoi(a) register flt *a; { Uint i; if (a->exp < BIAS) { return 0; } if (a->exp > BIAS + 30 && (a->sign == 0 || a->exp != BIAS + 31 || a->high != 0x4000 || a->low != 0)) { f_erange(); } i = (((Uint) a->high << 17) | (a->low >> 14)) >> (BIAS + 31 - a->exp); return (a->sign == 0) ? i : -i; } /* * NAME: f_ftoxf() * DESCRIPTION: convert flt to xfloat */ static void f_ftoxf(a, f) register flt *a; register xfloat *f; { register unsigned short exp; register unsigned short high; register Uint low; exp = a->exp; if (exp == 0) { /* zero */ f->high = 0; f->low = 0; return; } high = a->high; low = a->low; /* mantissa */ if ((Int) (low += 0x100) < 0) { low = 0; if ((short) ++high < 0) { high >>= 1; exp++; } } /* exponent */ if (exp > BIAS + 1023) { f_erange(); } if (exp < BIAS - 1022) { /* underflow */ f->high = 0; f->low = 0; return; } f->high = a->sign | ((exp - BIAS + 1023) << 4) | ((high >> 10) & 0x000f); f->low = ((Uint) high << 22) | (low >> 9); } /* * NAME: f_xftof() * DESCRIPTION: convert xfloat to flt */ static void f_xftof(f, a) register xfloat *f; register flt *a; { register unsigned short exp; exp = (f->high >> 4) & 0x07ff; if (exp == 0) { /* zero */ a->exp = 0; return; } a->exp = exp + BIAS - 1023; a->sign = f->high & 0x8000; a->high = 0x4000 | ((f->high & 0x0f) << 10) | (f->low >> 22); a->low = (f->low << 9) & 0x7fffffffL; } static flt tens[] = { { 0, 0x8002, 0x5000, 0x00000000L }, /* 10 ** 1 */ { 0, 0x8005, 0x6400, 0x00000000L }, /* 10 ** 2 */ { 0, 0x800C, 0x4E20, 0x00000000L }, /* 10 ** 4 */ { 0, 0x8019, 0x5F5E, 0x08000000L }, /* 10 ** 8 */ { 0, 0x8034, 0x470D, 0x726FC100L }, /* 10 ** 16 */ { 0, 0x8069, 0x4EE2, 0x6B6A0ADCL }, /* 10 ** 32 */ { 0, 0x80D3, 0x613C, 0x07D27FF4L }, /* 10 ** 64 */ { 0, 0x81A8, 0x49DD, 0x11F2603CL }, /* 10 ** 128 */ { 0, 0x8351, 0x553F, 0x3AFEE780L }, /* 10 ** 256 */ { 0, 0x86A3, 0x718C, 0x682BA984L }, /* 10 ** 512 */ }; static flt tenths[] = { { 0, 0x7FFB, 0x6666, 0x33333334L }, /* 10 ** -1 */ { 0, 0x7FF8, 0x51EB, 0x428F5C28L }, /* 10 ** -2 */ { 0, 0x7FF1, 0x68DB, 0x45D63888L }, /* 10 ** -4 */ { 0, 0x7FE4, 0x55E6, 0x1DC46118L }, /* 10 ** -8 */ { 0, 0x7FC9, 0x734A, 0x652FB114L }, /* 10 ** -16 */ { 0, 0x7F94, 0x67D8, 0x47AB5150L }, /* 10 ** -32 */ { 0, 0x7F2A, 0x543F, 0x7A89E950L }, /* 10 ** -64 */ { 0, 0x7E55, 0x6EE8, 0x119F1930L }, /* 10 ** -128 */ { 0, 0x7CAC, 0x6018, 0x50C958E0L }, /* 10 ** -256 */ { 0, 0x795A, 0x4824, 0x7B8CB6C8L }, /* 10 ** -512 */ }; /* * NAME: float->atof() * DESCRIPTION: Convert a string to a float. The string must be in the * proper format. Return TRUE if the operation was successful, * FALSE otherwise. */ bool flt_atof(s, f) char **s; xfloat *f; { flt a, b, c, *t; register unsigned short e, h; register char *p; p = *s; /* sign */ if (*p == '-') { a.sign = b.sign = 0x8000; p++; } else { a.sign = b.sign = 0; } a.exp = 0; b.low = 0; /* digits before . */ while (isdigit(*p)) { f_mult(&a, &tens[0]); h = (*p++ - '0') << 12; if (h != 0) { e = BIAS + 3; while ((short) h >= 0) { h <<= 1; --e; } b.exp = e; b.high = h >> 1; f_add(&a, &b); } if (a.exp > 0xffff - 10) { return FALSE; } } /* digits after . */ if (*p == '.') { c = tenths[0]; while (isdigit(*++p)) { if (c.exp > 10) { h = (*p - '0') << 12; if (h != 0) { e = BIAS + 3; while ((short) h >= 0) { h <<= 1; --e; } b.exp = e; b.high = h >> 1; b.low = 0; f_mult(&b, &c); f_add(&a, &b); } f_mult(&c, &tenths[0]); } } } /* exponent */ if (*p == 'e' || *p == 'E') { /* sign of exponent */ if (*++p == '-') { t = tenths; p++; } else { t = tens; if (*p == '+') { p++; } } /* get exponent */ e = 0; do { e *= 10; e += *p++ - '0'; if (e >= 1024) { return FALSE; } } while (isdigit(*p)); /* adjust number */ while (e != 0) { if ((e & 1) != 0) { f_mult(&a, t); if (a.exp < 0x1000 || a.exp > 0xf000) { break; } } e >>= 1; t++; } } if (a.exp >= BIAS + 1023 && (a.exp > BIAS + 1023 || (a.high == 0x7fff && a.low >= 0x7fffff00L))) { return FALSE; } f_ftoxf(&a, f); *s = p; return TRUE; } /* * NAME: float->ftoa() * DESCRIPTION: convert a float to a string */ void flt_ftoa(f, buffer) xfloat *f; char *buffer; { static flt tenmillion = { 0, 0x8016, 0x4c4b, 0x20000000L }; register unsigned short i; register short e; register Uint n; register char *p; register flt *t, *t2; char digits[9]; flt a; f_xftof(f, &a); if (a.exp == 0) { strcpy(buffer, "0"); return; } if (a.sign != 0) { *buffer++ = '-'; a.sign = 0; } /* reduce the float to range 1 .. 9.999999999, and extract exponent */ e = 0; if (a.exp >= BIAS) { /* >= 1 */ for (i = 10, t = &tens[9], t2 = &tenths[9]; i > 0; --i, --t, --t2) { e <<= 1; if (a.exp >= t->exp && (a.exp > t->exp || (a.high >= t->high && (a.high > t->high || a.low >= t->low)))) { e |= 1; f_mult(&a, t2); } } } else { /* < 1 */ for (i = 10, t = &tenths[9], t2 = &tens[9]; i > 0; --i, --t, --t2) { e <<= 1; if (a.exp <= t->exp && (a.exp < t->exp || (a.high <= t->high && (a.high < t->high || a.low <= t->low)))) { e |= 1; f_mult(&a, t2); } } if (a.exp < BIAS) { /* still < 1 */ f_mult(&a, &tens[0]); e++; } e = -e; } f_mult(&a, &tenmillion); f_37bits(&a); /* * obtain digits */ f_add(&a, &half); i = a.exp - BIAS + 1 - 15; n = ((Uint) a.high << i) | (a.low >> (31 - i)); if (n == 100000000L) { p = digits + 7; p[0] = '1'; p[1] = '\0'; i = 1; e++; } else { while (n != 0 && n % 10 == 0) { n /= 10; } p = digits + 8; *p = '\0'; i = 0; do { i++; *--p = '0' + n % 10; n /= 10; } while (n != 0); } if (e >= 8 || (e < -3 && i - e > 8)) { buffer[0] = *p; if (i != 1) { buffer[1] = '.'; memcpy(buffer + 2, p + 1, i - 1); i++; } buffer[i++] = 'e'; if (e >= 0) { buffer[i] = '+'; } else { buffer[i] = '-'; e = -e; } p = digits + 8; do { *--p = '0' + e % 10; e /= 10; } while (e != 0); strcpy(buffer + i + 1, p); } else if (e < 0) { e = 1 - e; memcpy(buffer, "0.000000", e); strcpy(buffer + e, p); } else { while (e >= 0) { *buffer++ = (*p == '\0') ? '0' : *p++; --e; } if (*p != '\0') { *buffer = '.'; strcpy(buffer + 1, p); } else { *buffer = '\0'; } } } /* * NAME: float->itof() * DESCRIPTION: convert an integer to a float */ void flt_itof(i, f) Int i; xfloat *f; { flt a; f_itof(i, &a); f_ftoxf(&a, f); f_xftof(f, &a); } /* * NAME: float->ftoi() * DESCRIPTION: convert a float to an integer */ Int flt_ftoi(f) xfloat *f; { flt a; f_xftof(f, &a); f_round(&a); return f_ftoi(&a); } /* * NAME: float->add() * DESCRIPTION: add two floats */ void flt_add(f1, f2) xfloat *f1, *f2; { flt a, b; f_xftof(f2, &b); f_xftof(f1, &a); f_add(&a, &b); f_ftoxf(&a, f1); } /* * NAME: float->sub() * DESCRIPTION: subtract a float from a float */ void flt_sub(f1, f2) xfloat *f1, *f2; { flt a, b; f_xftof(f2, &b); f_xftof(f1, &a); f_sub(&a, &b); f_ftoxf(&a, f1); } /* * NAME: float->mult() * DESCRIPTION: multiply two floats */ void flt_mult(f1, f2) xfloat *f1, *f2; { flt a, b; f_xftof(f1, &a); f_xftof(f2, &b); f_mult(&a, &b); f_ftoxf(&a, f1); } /* * NAME: float->div() * DESCRIPTION: devide a float by a float */ void flt_div(f1, f2) xfloat *f1, *f2; { flt a, b; f_xftof(f2, &b); f_xftof(f1, &a); f_div(&a, &b); f_ftoxf(&a, f1); } /* * NAME: float->cmp() * DESCRIPTION: compare two xfloats */ int flt_cmp(f1, f2) register xfloat *f1, *f2; { if ((short) (f1->high ^ f2->high) < 0) { return ((short) f1->high < 0) ? -1 : 1; } if (f1->high == f2->high && f1->low == f2->low) { return 0; } if (f1->high <= f2->high && (f1->high < f2->high || f1->low < f2->low)) { return ((short) f1->high < 0) ? 1 : -1; } return ((short) f1->high < 0) ? -1 : 1; } /* * NAME: float->floor() * DESCRIPTION: round a float downwards */ void flt_floor(f) xfloat *f; { flt a, b; f_xftof(f, &a); b = a; f_trunc(&b); if (b.sign != 0 && (a.exp != b.exp || a.high != b.high || a.low != b.low)) { f_sub(&b, &one); } f_ftoxf(&b, f); } /* * NAME: float->ceil() * DESCRIPTION: round a float upwards */ void flt_ceil(f) xfloat *f; { flt a, b; f_xftof(f, &a); b = a; f_trunc(&b); if (b.sign == 0 && (a.exp != b.exp || a.high != b.high || a.low != b.low)) { f_add(&b, &one); } f_ftoxf(&b, f); } /* * NAME: float->fmod() * DESCRIPTION: perform fmod */ void flt_fmod(f1, f2) xfloat *f1, *f2; { flt a, b, c; f_xftof(f2, &b); if (b.exp == 0) { f_edom(); } f_xftof(f1, &a); if (a.exp == 0) { return; } c.sign = a.sign; c.high = b.high; c.low = b.low; while (a.exp >= b.exp && (a.exp > b.exp || (a.high >= b.high && (a.high > b.high || a.low >= b.low)))) { c.exp = a.exp; if (a.high <= c.high && (a.high < c.high || a.low < c.low)) { c.exp--; } f_sub(&a, &c); } f_ftoxf(&a, f1); } /* * NAME: float->frexp() * DESCRIPTION: split a float into a fraction and an exponent */ Int flt_frexp(f) register xfloat *f; { short e; if (f->high == 0) { return 0; } e = ((f->high & 0x7ff0) >> 4) - 1022; f->high = (f->high & 0x800f) | (1022 << 4); return e; } /* * NAME: float->ldexp() * DESCRIPTION: make a float from a fraction and an exponent */ void flt_ldexp(f, exp) register xfloat *f; register Int exp; { if (f->high == 0) { return; } exp += (f->high & 0x7ff0) >> 4; if (exp <= 0) { f->high = 0; f->low = 0; return; } if (exp > 1023 + 1023) { f_erange(); } f->high = (f->high & 0x800f) | (exp << 4); } /* * NAME: float->modf() * DESCRIPTION: split float into fraction and integer part */ void flt_modf(f1, f2) xfloat *f1, *f2; { flt a, b; f_xftof(f1, &a); if (a.exp < BIAS) { b.exp = 0; } else { b = a; f_trunc(&b); f_sub(&a, &b); } f_ftoxf(&a, f1); f_ftoxf(&b, f2); }