diff options
Diffstat (limited to 'src/util/xldio.c')
| -rw-r--r-- | src/util/xldio.c | 603 |
1 files changed, 603 insertions, 0 deletions
diff --git a/src/util/xldio.c b/src/util/xldio.c new file mode 100644 index 0000000..ab40058 --- /dev/null +++ b/src/util/xldio.c @@ -0,0 +1,603 @@ +/* This file contain long double I/O routines for Windows (because Windows API + don't support long double at all. + + They don't work on other architectures. So be curefull. */ + + +/* This source comes from the DJGPP runtime library. It has been hacked + to work with XaoS */ +/* Copyright (C) 1998 DJ Delorie, see COPYING.DJ for details */ +/* Copyright (C) 1997 DJ Delorie, see COPYING.DJ for details */ +/* Copyright (C) 1996 DJ Delorie, see COPYING.DJ for details */ +/* Copyright (C) 1994 DJ Delorie, see COPYING.DJ for details */ +#include <config.h> +#include <stdlib.h> +#include <ctype.h> +#include <xio.h> +/*#include <libc/unconst.h>*/ +#include <sys/types.h> +#include <ctype.h> +#include <string.h> +#include <stdio.h> +#include <math.h> +#include <xldio.h> +#ifdef USE_XLDIO + +static long double powten[] = { + 1e1L, 1e2L, 1e4L, 1e8L, 1e16L, 1e32L, 1e64L, 1e128L, 1e256L, + 1e512L, 1e1024L, 1e2048L, 1e4096L +}; + +long double x_strtold(CONST char *s, CONST char **sret) +{ + long double r; /* result */ + int e, ne; /* exponent */ + int sign; /* +- 1.0 */ + int esign; + int flags = 0; + int l2powm1; + + r = 0.0L; + sign = 1; + e = ne = 0; + esign = 1; + + while (*s && isspace(*s)) + s++; + + if (*s == '+') + s++; + else if (*s == '-') { + sign = -1; + s++; + } + + while ((*s >= '0') && (*s <= '9')) { + flags |= 1; + r *= 10.0L; + r += *s - '0'; + s++; + } + + if (*s == '.') { + s++; + while ((*s >= '0') && (*s <= '9')) { + flags |= 2; + r *= 10.0L; + r += *s - '0'; + s++; + ne++; + } + } + if (flags == 0) { + if (sret) + *sret = /*unconst(s, char *) */ s; + return 0.0L; + } + + if ((*s == 'e') || (*s == 'E')) { + s++; + if (*s == '+') + s++; + else if (*s == '-') { + s++; + esign = -1; + } + while ((*s >= '0') && (*s <= '9')) { + e *= 10; + e += *s - '0'; + s++; + } + } + if (esign < 0) { + esign = -esign; + e = -e; + } + e = e - ne; + if (e < -4096) { + /* possibly subnormal number, 10^e would overflow */ + r *= 1.0e-2048L; + e += 2048; + } + if (e < 0) { + e = -e; + esign = -esign; + } + if (e >= 8192) + e = 8191; + if (e) { + long double d = 1.0L; + l2powm1 = 0; + while (e) { + if (e & 1) + d *= powten[l2powm1]; + e >>= 1; + l2powm1++; + } + if (esign > 0) + r *= d; + else + r /= d; + } + if (sret) + *sret = /*unconst(s, char *) */ s; + return r * sign; +} + +#if 0 +main() +{ + printf("%E", (float) x_strtold("1.4E15", NULL)); +} +#endif + + +#define MAXEXPLD 4952 /* this includes subnormal numbers */ +static int is_nan = 0; +static char decimal = '.'; +static long double pten[] = { + 1e1L, 1e2L, 1e4L, 1e8L, 1e16L, 1e32L, 1e64L, 1e128L, 1e256L, + 1e512L, 1e1024L, 1e2048L, 1e4096L +}; + +static long double ptenneg[] = { + 1e-1L, 1e-2L, 1e-4L, 1e-8L, 1e-16L, 1e-32L, 1e-64L, 1e-128L, 1e-256L, + 1e-512L, 1e-1024L, 1e-2048L, 1e-4096L +}; + +static inline char tochar(int n) +{ + if (n >= 9) + return '9'; + if (n <= 0) + return '0'; + return n + '0'; +} + +static inline int todigit(char c) +{ + if (c <= '0') + return 0; + if (c >= '9') + return 9; + return c - '0'; +} + +#define LONGINT 0x01 /* long integer */ +#define LONGDBL 0x02 /* long double */ +#define SHORTINT 0x04 /* short integer */ +#define ALT 0x08 /* alternate form */ +#define LADJUST 0x10 /* left adjustment */ +#define ZEROPAD 0x20 /* zero (as opposed to blank) pad */ +#define HEXPREFIX 0x40 /* add 0x or 0X prefix */ + +#define MAXP 4096 +#define NP 12 +#define P (4294967296.0L * 4294967296.0L * 2.0L) /* 2^65 */ +static long double INVPREC = P; +static long double PREC = 1.0L / P; +#undef P +/* + * Defining FAST_LDOUBLE_CONVERSION results in a little bit faster + * version, which might be less accurate (about 1 bit) for long + * double. For 'normal' double it doesn't matter. + */ +/* #define FAST_LDOUBLE_CONVERSION */ +#if 1 +#define modfl mymodfl +inline long double m_floor(long double x) +{ + register long double __value; + volatile unsigned short int __cw, __cwtmp; + + asm volatile ("fnstcw %0":"=m" (__cw)); + __cwtmp = (__cw & 0xf3ff) | 0x0400; /* rounding down */ + asm volatile ("fldcw %0"::"m" (__cwtmp)); + asm volatile ("frndint":"=t" (__value):"0"(x)); + asm volatile ("fldcw %0"::"m" (__cw)); + + return __value; + +} + +static inline long double mymodfl(long double x, long double *pint) +{ + /*int p=(int) x; */ + long double p = m_floor(x); + long double frac = x - p; + if (x < 0) + p = p + 1, frac = frac - 1; + *pint = p; + return frac; +} +#endif +static char *exponentl(char *p, int expv, unsigned char fmtch) +{ + char *t; + char expbuf[MAXEXPLD]; + + *p++ = fmtch; + if (expv < 0) { + expv = -expv; + *p++ = '-'; + } else + *p++ = '+'; + t = expbuf + MAXEXPLD; + if (expv > 9) { + do { + *--t = tochar(expv % 10); + } + while ((expv /= 10) > 9); + *--t = tochar(expv); + for (; t < expbuf + MAXEXPLD; *p++ = *t++); + } else { + *p++ = '0'; + *p++ = tochar(expv); + } + return p; +} + +static int isspeciall(long double d, char *bufp) +{ + struct IEEExp { + unsigned manl:32; + unsigned manh:32; + unsigned exp:15; + unsigned sign:1; + } *ip = (struct IEEExp *) &d; + + is_nan = 0; /* don't assume the static is 0 (emacs) */ + if (ip->exp != 0x7fff) + return (0); + if ((ip->manh & 0x7fffffff) || ip->manl) { + strcpy(bufp, "NaN"); + is_nan = ip->sign ? -1 : 1; /* kludge: we don't need the sign, it's not nice + but it should work */ + } else + (void) strcpy(bufp, "Inf"); + return (3); +} + +static char *my_roundl(long double fract, int *expv, char *start, + char *end, char ch, char *signp) +{ + long double tmp; + + if (fract) { + if (fract == 0.5L) { + char *e = end; + if (*e == '.') + e--; + if (*e == '0' || *e == '2' || *e == '4' || *e == '6' + || *e == '8') { + tmp = 3.0; + goto start; + } + } + (void) modfl(fract * 10.0L, &tmp); + } else + tmp = todigit(ch); + start: + if (tmp > 4) + for (;; --end) { + if (*end == decimal) + --end; + if (++*end <= '9') + break; + *end = '0'; + if (end == start) { + if (expv) { /* e/E; increment exponent */ + *end = '1'; + ++*expv; + } else { /* f; add extra digit */ + *--end = '1'; + --start; + } + break; + } + } + /* ``"%.3f", (double)-0.0004'' gives you a negative 0. */ + else if (*signp == '-') + for (;; --end) { + if (*end == decimal) + --end; + if (*end != '0') + break; + if (end == start) + *signp = 0; + } + return start; +} + + +static int +cvtl(long double number, int prec, int flags, char *signp, + unsigned char fmtch, char *startp, char *endp) +{ + char *p, *t; + long double fract = 0; + int dotrim, expcnt, gformat; + int doextradps = 0; /* Do extra decimal places if the precision needs it */ + int doingzero = 0; /* We're displaying 0.0 */ + long double integer, tmp; + + if ((expcnt = isspeciall(number, startp))) + return (expcnt); + + dotrim = expcnt = gformat = 0; + /* fract = modfl(number, &integer); */ + integer = number; + + /* get an extra slot for rounding. */ + t = ++startp; + + p = endp - 1; + if (integer) { + int i, lp = NP, pt = MAXP; +#ifndef FAST_LDOUBLE_CONVERSION + long double oint = integer, dd = 1.0L; +#endif + if (integer > INVPREC) { + integer *= PREC; + while (lp >= 0) { + if (integer >= pten[lp]) { + expcnt += pt; + integer *= ptenneg[lp]; +#ifndef FAST_LDOUBLE_CONVERSION + dd *= pten[lp]; +#endif + } + pt >>= 1; + lp--; + } +#ifndef FAST_LDOUBLE_CONVERSION + integer = oint / dd; +#else + integer *= INVPREC; +#endif + } + /* + * Do we really need this ? + */ + for (i = 0; i < expcnt; i++) + *p-- = '0'; + } + number = integer; + fract = modfl(number, &integer); + /* If integer is zero then we need to look at where the sig figs are */ + if (integer < 1) { + /* If fract is zero the zero before the decimal point is a sig fig */ + if (fract == 0.0) + doingzero = 1; + /* If fract is non-zero all sig figs are in fractional part */ + else + doextradps = 1; + } + /* + * get integer portion of number; put into the end of the buffer; the + * .01 is added for modf(356.0 / 10, &integer) returning .59999999... + */ + for (; integer; ++expcnt) { + tmp = modfl(integer * 0.1L, &integer); + *p-- = tochar((int) ((tmp + .01L) * 10)); + } + switch (fmtch) { + case 'f': + /* reverse integer into beginning of buffer */ + if (expcnt) + for (; ++p < endp; *t++ = *p); + else + *t++ = '0'; + /* + * if precision required or alternate flag set, add in a + * decimal point. + */ + if (prec || flags & ALT) + *t++ = decimal; + /* if requires more precision and some fraction left */ + if (fract) { + if (prec) + do { + fract = modfl(fract * 10.0L, &tmp); + *t++ = tochar((int) tmp); + } + while (--prec && fract); + if (fract) + startp = my_roundl(fract, (int *) NULL, startp, + t - 1, (char) 0, signp); + } + for (; prec--; *t++ = '0'); + break; + case 'e': + case 'E': + eformat: + if (expcnt) { + *t++ = *++p; + if (prec || flags & ALT) + *t++ = decimal; + /* if requires more precision and some integer left */ + for (; prec && ++p < endp; --prec) + *t++ = *p; + /* + * if done precision and more of the integer component, + * round using it; adjust fract so we don't re-round + * later. + */ + if (!prec && ++p < endp) { + fract = 0; + startp = my_roundl((long double) 0.0L, &expcnt, + startp, t - 1, *p, signp); + } + /* adjust expcnt for digit in front of decimal */ + --expcnt; + } + /* until first fractional digit, decrement exponent */ + else if (fract) { + int lp = NP, pt = MAXP; +#ifndef FAST_LDOUBLE_CONVERSION + long double ofract = fract, dd = 1.0L; +#endif + expcnt = -1; + if (fract < PREC) { + fract *= INVPREC; + while (lp >= 0) { + if (fract <= ptenneg[lp]) { + expcnt -= pt; + fract *= pten[lp]; +#ifndef FAST_LDOUBLE_CONVERSION + dd *= pten[lp]; +#endif + } + pt >>= 1; + lp--; + } +#ifndef FAST_LDOUBLE_CONVERSION + fract = ofract * dd; +#else + fract *= PREC; +#endif + } + /* adjust expcnt for digit in front of decimal */ + for ( /* expcnt = -1 */ ;; --expcnt) { + fract = modfl(fract * 10.0L, &tmp); + if (tmp) + break; + } + *t++ = tochar((int) tmp); + if (prec || flags & ALT) + *t++ = decimal; + } else { + *t++ = '0'; + if (prec || flags & ALT) + *t++ = decimal; + } + /* if requires more precision and some fraction left */ + if (fract) { + if (prec) + do { + fract = modfl(fract * 10.0L, &tmp); + *t++ = tochar((int) tmp); + } + while (--prec && fract); + if (fract) + startp = + my_roundl(fract, &expcnt, startp, t - 1, (char) 0, + signp); + } + /* if requires more precision */ + for (; prec--; *t++ = '0'); + + /* unless alternate flag, trim any g/G format trailing 0's */ + if (gformat && !(flags & ALT)) { + while (t > startp && *--t == '0'); + if (*t == decimal) + --t; + ++t; + } + t = exponentl(t, expcnt, fmtch); + break; + case 'g': + case 'G': + if (prec) { + /* If doing zero and precision is greater than 0 count the + * 0 before the decimal place */ + if (doingzero) + --prec; + } else { + /* a precision of 0 is treated as precision of 1 unless doing zero */ + if (!doingzero) + ++prec; + } + /* + * ``The style used depends on the value converted; style e + * will be used only if the exponent resulting from the + * conversion is less than -4 or greater than the precision.'' + * -- ANSI X3J11 + */ + if (expcnt > prec || (!expcnt && fract && fract < .0001)) { + /* + * g/G format counts "significant digits, not digits of + * precision; for the e/E format, this just causes an + * off-by-one problem, i.e. g/G considers the digit + * before the decimal point significant and e/E doesn't + * count it as precision. + */ + --prec; + fmtch -= 2; /* G->E, g->e */ + gformat = 1; + goto eformat; + } + /* + * reverse integer into beginning of buffer, + * note, decrement precision + */ + if (expcnt) + for (; ++p < endp; *t++ = *p, --prec); + else + *t++ = '0'; + /* + * if precision required or alternate flag set, add in a + * decimal point. If no digits yet, add in leading 0. + */ + if (prec || flags & ALT) { + dotrim = 1; + *t++ = decimal; + } else + dotrim = 0; + /* if requires more precision and some fraction left */ + while (prec && fract) { + fract = modfl(fract * 10.0L, &tmp); + *t++ = tochar((int) tmp); + /* If we're not adding 0s + * or we are but they're sig figs: + * decrement the precision */ + if ((doextradps != 1) || ((int) tmp != 0)) { + doextradps = 0; + prec--; + } + } + if (fract) + startp = + my_roundl(fract, (int *) NULL, startp, t - 1, (char) 0, + signp); + /* alternate format, adds 0's for precision, else trim 0's */ + if (flags & ALT) + for (; prec--; *t++ = '0'); + else if (dotrim) { + while (t > startp && *--t == '0'); + if (*t != decimal) + ++t; + } + } + return t - startp; +} + +#if 0 +main() +{ + static char buf[4096]; + int i; + cvtl(0.00000000000000000005, 4, 0, &softsign, 'G', buf, + buf + sizeof(buf) - 1); + printf("%s\n", buf + 1); + printf("%.30LG\n", (long double) 234236723234234231235324.47239L); +} + +#endif +void x_ldout(long double param, int prec, xio_file stream) +{ + static char buf[4095]; + char softsign = 0; + int l; + if (param < 0) + xio_putc('-', stream), param = -param; + l = cvtl(param, prec, 0, &softsign, 'G', buf, buf + sizeof(buf)); + /*printf("a:%s %i\n",buf+1, prec); */ + buf[l + 2] = 0; + l = strlen(buf + 1); + if (buf[l] == '.') + buf[l] = 0; + /*printf("b:%s %i\n",buf+1, prec); */ + xio_puts(buf + 1, stream); +} +#endif |