From 90101c0dec0b3cb2fd06bfca7506472267f2fddf Mon Sep 17 00:00:00 2001 From: ceriel Date: Wed, 26 Jul 1989 11:08:24 +0000 Subject: [PATCH] Added ext_comp.c --- lang/cem/libcc/gen/LIST | 1 + lang/cem/libcc/gen/ecvt.c | 117 ++---- lang/cem/libcc/gen/ext_comp.c | 683 ++++++++++++++++++++++++++++++++++ lang/cem/libcc/gen/strtod.c | 106 +----- 4 files changed, 731 insertions(+), 176 deletions(-) create mode 100644 lang/cem/libcc/gen/ext_comp.c diff --git a/lang/cem/libcc/gen/LIST b/lang/cem/libcc/gen/LIST index 66d70b179..b9251b22c 100644 --- a/lang/cem/libcc/gen/LIST +++ b/lang/cem/libcc/gen/LIST @@ -17,6 +17,7 @@ ffc.c ffs.c gcvt.c ecvt.c +ext_comp.c getlogin.c index.c l3.c diff --git a/lang/cem/libcc/gen/ecvt.c b/lang/cem/libcc/gen/ecvt.c index c23a29117..16cb4d1b4 100644 --- a/lang/cem/libcc/gen/ecvt.c +++ b/lang/cem/libcc/gen/ecvt.c @@ -1,14 +1,39 @@ /* $Header$ */ + #ifndef NOFLOAT -static char *cvt(); -#define NDIGITS 128 +struct mantissa { + unsigned long h_32; + unsigned long l_32; +}; + +struct EXTEND { + short sign; + short exp; + struct mantissa mantissa; +#define m1 mantissa.h_32 +#define m2 mantissa.l_32 +}; + +extern char *_ext_str_cvt(); + +static char * +cvt(value, ndigit, decpt, sign, ecvtflag) + double value; + int ndigit, *decpt, *sign; +{ + struct EXTEND e; + + _dbl_ext_cvt(value, &e); + return _ext_str_cvt(&e, ndigit, decpt, sign, ecvtflag); +} char * ecvt(value, ndigit, decpt, sign) double value; int ndigit, *decpt, *sign; { + return cvt(value, ndigit, decpt, sign, 1); } @@ -20,90 +45,4 @@ fcvt(value, ndigit, decpt, sign) return cvt(value, ndigit, decpt, sign, 0); } -static struct powers_of_10 { - double pval; - double rpval; - int exp; -} p10[] = { - 1.0e32, 1.0e-32, 32, - 1.0e16, 1.0e-16, 16, - 1.0e8, 1.0e-8, 8, - 1.0e4, 1.0e-4, 4, - 1.0e2, 1.0e-2, 2, - 1.0e1, 1.0e-1, 1, - 1.0e0, 1.0e0, 0 -}; - -static char * -cvt(value, ndigit, decpt, sign, ecvtflag) - double value; - int ndigit, *decpt, *sign; -{ - static char buf[NDIGITS+1]; - register char *p = buf; - register char *pe; - - if (ndigit < 0) ndigit = 0; - if (ndigit > NDIGITS) ndigit = NDIGITS; - pe = &buf[ndigit]; - buf[0] = '\0'; - - *sign = 0; - if (value < 0) { - *sign = 1; - value = -value; - } - - *decpt = 0; - if (value != 0.0) { - register struct powers_of_10 *pp = &p10[0]; - - if (value >= 10.0) do { - while (value >= pp->pval) { - value *= pp->rpval; - *decpt += pp->exp; - } - } while ((++pp)->exp > 0); - - pp = &p10[0]; - if (value < 1.0) do { - while (value * pp->pval < 10.0) { - value *= pp->pval; - *decpt -= pp->exp; - } - } while ((++pp)->exp > 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]; - } - while (p <= pe) { - *p++ = (int)value + '0'; - value = 10.0 * (value - (int)value); - } - 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; -} -#endif +#endif NOFLOAT diff --git a/lang/cem/libcc/gen/ext_comp.c b/lang/cem/libcc/gen/ext_comp.c new file mode 100644 index 000000000..e990d541a --- /dev/null +++ b/lang/cem/libcc/gen/ext_comp.c @@ -0,0 +1,683 @@ +/* + (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +/* 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. ??? +*/ + +static int b64_add(); +static int b64_sft(); + +#include + +struct mantissa { + unsigned long h_32; + unsigned long l_32; +}; + +struct EXTEND { + short sign; + short exp; + struct mantissa mantissa; +#define m1 mantissa.h_32 +#define m2 mantissa.l_32 +}; + +static +mul_ext(e1,e2,e3) +struct EXTEND *e1,*e2,*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(e1,e2,e3) +struct EXTEND *e1,*e2,*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(e1, e2) + struct EXTEND *e1, *e2; +{ + int sign = e1->sign ? -1 : 1; + + if (e1->sign > e2->sign) return -1; + if (e1->sign < e2->sign) return 1; + if (e1->exp < e2->exp) return -sign; + if (e1->exp > e2->exp) return sign; + if (e1->m1 < e2->m1) return -sign; + if (e1->m1 > e2->m1) return sign; + if (e1->m2 < e2->m2) return -sign; + if (e1->m2 > e2->m2) return sign; + return 0; +} + +static +b64_sft(e1,n) + 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(e1,e2) + /* + * pointers to 64 bit 'registers' + */ + struct mantissa *e1,*e2; +{ + 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(asign = 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 = s - 1; + + if (c == 'E' || c == 'e') { + int exp1 = 0; + int sign = 1; + + switch(*s) { + case '-': + sign = -1; + case '+': + s++; + } + if (c = *s, isdigit(c)) { + do { + exp1 = 10 * exp1 + (c - '0'); + } while (c = *++s, isdigit(c)); + if (ss) *ss = s; + } + exp += sign * exp1; + } + 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); +} + +extern double ldexp(), frexp(), modf(); + +#define NDIGITS 128 + +char * +_ext_str_cvt(e, ndigit, decpt, sign, ecvtflag) + struct EXTEND *e; + int ndigit, *decpt, *sign; +{ + /* Like cvt(), but for extended precision */ + + static char buf[NDIGITS+1]; + register char *p = buf; + register char *pe; + + 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++; + pp--; + mul_ext(e,&r_big_ten_powers[pp-big_ten_powers],e); + *decpt += (pp - big_ten_powers)* + (sizeof(ten_powers)/sizeof(ten_powers[0])); + pp = &ten_powers[1]; + while(pp<&ten_powers[(sizeof(ten_powers)/sizeof(ten_powers[0]))] && + cmp_ext(e, pp) >= 0) pp++; + pp--; + mul_ext(e, &r_ten_powers[pp-ten_powers], e); + *decpt += pp - ten_powers; + + if (cmp_ext(e, &ten_powers[0]) < 0) { + pp = &r_big_ten_powers[1]; + while(cmp_ext(e,pp) < 0) pp++; + pp--; + mul_ext(e,&big_ten_powers[pp-r_big_ten_powers],e); + *decpt -= (pp - r_big_ten_powers)* + (sizeof(ten_powers)/sizeof(ten_powers[0])); + /* here, value >= 10 ** -28 */ + mul_ext(e, &ten_powers[1], e); + (*decpt)--; + pp = &r_ten_powers[0]; + while(cmp_ext(e, pp) < 0) pp++; + mul_ext(e, &ten_powers[pp-r_ten_powers], e); + *decpt -= pp - r_ten_powers; + } + (*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]; + } + while (p <= pe) { + 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'; + mul_ext(e, &ten_powers[1], e); + } + 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(value, e) + double value; + struct EXTEND *e; +{ + /* Convert double to extended + */ + int exponent; + register int i; + + value = frexp(value, &exponent); + e->sign = value < 0.0; + if (e->sign) value = -value; + e->exp = exponent - 1; + e->m1 = 0; + e->m2 = 0; + for (i = 64; i > 0 && value != 0; i--) { + double ipart; + + b64_sft(&(e->mantissa),-1); + value = modf(2.0*value, &ipart); + if (ipart) { + e->m2 |= 1; + } + } + if (i > 0) b64_sft(&(e->mantissa),-i); +} + +double +_ext_dbl_cvt(e) + struct EXTEND *e; +{ + /* Convert extended to double + */ + double f = ldexp(ldexp((double)e->m1, 32) + (double)e->m2, e->exp-63); + if (e->sign) f = -f; + return f; +} diff --git a/lang/cem/libcc/gen/strtod.c b/lang/cem/libcc/gen/strtod.c index 7982f1f20..2f6f61442 100644 --- a/lang/cem/libcc/gen/strtod.c +++ b/lang/cem/libcc/gen/strtod.c @@ -1,97 +1,29 @@ /* $Header$ */ #ifndef NOFLOAT -#include -extern double ldexp(); +struct mantissa { + unsigned long h_32; + unsigned long l_32; +}; -#define MAX (0x7fffffffL/10) +struct EXTEND { + short sign; + short exp; + struct mantissa mantissa; +#define m1 mantissa.h_32 +#define m2 mantissa.l_32 +}; + + +extern double _ext_dbl_cvt(); double strtod(p, pp) - register char *p; - char **pp; + char *p, **pp; { - register int c; - int exp = 0, sign = 1, expsign = 0; - double fl; - long lowl = 0, highl = 0, pos = 1; - int dotseen = 0; - int digit_seen = 0; - - if (pp) *pp = p; - while (isspace(*p)) p++; - c = *p; - - switch (c) { - case '-': - sign = -1; - case '+': - p++; - } - - while (isdigit(c = *p++) || (c == '.' && ! dotseen++)) { - if (c == '.') continue; - digit_seen = 1; - if (highl < MAX) { - highl = (highl << 3) + (highl << 1) + (c - '0'); - } - else if (pos < MAX) { - pos = (pos << 3) + (pos << 1); - lowl = (lowl << 3) + (lowl << 1) + (c - '0'); - } - else exp++; - if (dotseen) exp--; - } - if (! digit_seen) return 0.0; - fl = highl; - if (pos > 1) { - fl = pos * fl + lowl; - } - - if (pp) *pp = p-1; - - if (c == 'E' || c == 'e') { - int exp1 = 0; - int sign = 1; - - switch (*p) { - case '-': - sign = -1; - case '+': - p++; - } - if (isdigit(c = *p)) { - do { - exp1 = 10 * exp1 + c - '0'; - } while (isdigit(c = *++p)); - if (pp) *pp = p; - } - exp += sign * exp1; - } - - if (fl == 0.0) return 0.0; - - if (exp < 0) { - expsign = 1; - exp = -exp; - } - - if (exp != 0) { - int oldexp = exp; - double exp5 = 5.0; - double correction = 1.0; - - while (exp) { - if (exp % 2) correction *= exp5; - exp /= 2; - if (exp != 0) exp5 *= exp5; - } - if (expsign) fl = fl / correction; - else fl = fl * correction; - - fl = ldexp(fl, expsign ? -oldexp : oldexp); - } + struct EXTEND e; - return sign * fl; + _str_ext_cvt(p, pp, &e); + return _ext_dbl_cvt(&e); } -#endif +#endif NOFLOAT -- 2.34.1