/* (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. See the copyright notice in the ACK home directory, in the file "Copyright". */ /* $Id$ */ #include #include "flt_misc.h" /* 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= BIGSZ) { flt_mul(&x, neg ? &r_big_10pow[BIGSZ-1] : &big_10pow[BIGSZ-1],&x); if (!status) status = flt_status; divsz -= BIGSZ-1; } flt_mul(&x, (neg ? r_big_10pow : big_10pow) + divsz, e); if (!status) status = flt_status; flt_status = status; } void flt_str2flt(char *s, flt_arith *e) { register int c; int dotseen = 0; int digitseen = 0; int exp = 0; while (isspace(*s)) s++; flt_status = 0; e->flt_sign = 0; e->flt_exp = 0; e->m1 = e->m2 = 0; c = *s; switch(c) { case '-': e->flt_sign = 1; case '+': s++; } while (c = *s++, isdigit(c) || (c == '.' && ! dotseen++)) { if (c == '.') continue; digitseen = 1; if (e->m1 >= 0 && e->m1 <= 0x7FFFFFFF/5) { struct flt_mantissa a1; a1 = e->flt_mantissa; flt_b64_sft(&(e->flt_mantissa), -3); flt_b64_sft(&a1, -1); flt_b64_add(&(e->flt_mantissa), &a1); a1.flt_h_32 = 0; a1.flt_l_32 = c - '0'; flt_b64_add(&(e->flt_mantissa), &a1); } else exp++; if (dotseen) exp--; } if (! digitseen) { flt_status = FLT_NOFLT; return; } 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)); } exp += sign * exp1; } if (e->m1 == 0 && e->m2 == 0) return; e->flt_exp = 63; flt_nrm(e); add_exponent(e, exp); flt_chk(e); } #define NDIG 18 static char *flt_ecvt(flt_arith *e, int *decpt, int *sign) { /* Like ecvt(), but for extended precision */ static char buf[NDIG+1]; register char *p = buf; register char *pe; register int findex = 0; pe = &buf[NDIG]; buf[0] = '\0'; *sign = 0; if (e->flt_sign) { *sign = 1; e->flt_sign = 0; } *decpt = 0; if (e->m1 != 0) { register flt_arith *pp = &big_10pow[1]; findex = 1; while (flt_cmp(e, &big_10pow[BIGSZ-1]) >= 0) { flt_mul(e,&r_big_10pow[BIGSZ-1],e); *decpt += (BIGSZ-1)*SMALLSZ; } while (flt_cmp(e,pp) >= 0) { pp++; findex++; } findex--; flt_mul(e,&r_big_10pow[findex],e); *decpt += findex*SMALLSZ; pp = &s10pow[1]; findex = 1; while (pp < &s10pow[SMALLSZ] && flt_cmp(e, pp) >= 0) { pp++; findex++; } findex--; *decpt += findex; if (flt_cmp(e, &s10pow[0]) < 0) { while (flt_cmp(e, &r_big_10pow[BIGSZ-1]) < 0) { flt_mul(e,&big_10pow[BIGSZ-1],e); *decpt -= (BIGSZ-1)*SMALLSZ; } pp = &r_big_10pow[1]; findex = 1; while(flt_cmp(e,pp) < 0) { pp++; findex++; } findex--; flt_mul(e,&big_10pow[findex],e); *decpt -= findex*SMALLSZ; /* here, value >= 10 ** -28 */ flt_mul(e, &s10pow[1], e); (*decpt)--; pp = &r_10pow[0]; findex = 0; while(flt_cmp(e, pp) < 0) { pp++; findex++; } flt_mul(e, &s10pow[findex], e); *decpt -= findex; findex = 0; } (*decpt)++; /* because now value in [1.0, 10.0) */ } while (p <= pe) { if (findex) { flt_arith tc, oldtc; int count = 0; oldtc.flt_exp = 0; oldtc.flt_sign = 0; oldtc.m1 = 0; oldtc.m2 = 0; tc = s10pow[findex]; while (flt_cmp(e, &tc) >= 0) { oldtc = tc; flt_add(&tc, &s10pow[findex], &tc); count++; } *p++ = count + '0'; oldtc.flt_sign = 1; flt_add(e, &oldtc, e); findex--; continue; } if (e->flt_exp >= 0 && e->m1 != 0) { flt_arith x; x.m2 = 0; x.flt_exp = e->flt_exp; x.flt_sign = 1; x.m1 = (e->m1 >> 1) & 0x7FFFFFFF; x.m1 = x.m1>>(30-e->flt_exp); *p++ = (x.m1) + '0'; if (x.m1) { x.m1 = x.m1 << (31-e->flt_exp); flt_add(e, &x, e); } } else *p++ = '0'; if (e->m1) flt_mul(e, &s10pow[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; } } *pe = '\0'; while (--pe > buf && *pe == '0') *pe = '\0'; } return buf; } void flt_flt2str(flt_arith *e, char *buf, int bufsize) { int sign, dp; register int i; register char *s1; char Xbuf[NDIG+12]; register char *s = Xbuf; flt_arith e1; e1 = *e; flt_status = 0; s1 = flt_ecvt(&e1,&dp,&sign); if (sign) *s++ = '-'; *s++ = *s1++; *s++ = '.'; for (i = NDIG-1; i > 0; i--) { if (*s1) *s++ = *s1++; else { if (i == NDIG-1) *s++ = '0'; break; } } if (e->m1 | e->m2) { --dp ; } if (dp != 0) { *s++ = 'e'; if ( dp<0 ) { *s++ = '-' ; dp= -dp ; } else { *s++ = '+' ; } s1 = &Xbuf[NDIG+12]; *--s1 = '\0'; do { *--s1 = dp % 10 + '0'; dp /= 10; } while (dp != 0); while (*s1) *s++ = *s1++; } *s++ = '\0'; if (s - Xbuf > bufsize) { flt_status = FLT_BTSM; return; } s = Xbuf; s1 = buf; do { *s1++ = *s; } while (*s++); }