cosmopolitan/libc/stdio/dtoa.c
2023-04-27 05:44:32 -07:00

534 lines
15 KiB
C

/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
╞══════════════════════════════════════════════════════════════════════════════╡
│ Copyright (C) 1997, 1999, 2001 Lucent Technologies │
│ All Rights Reserved │
│ │
│ Permission to use, copy, modify, and distribute this software and │
│ its documentation for any purpose and without fee is hereby │
│ granted, provided that the above copyright notice appear in all │
│ copies and that both that the copyright notice and this │
│ permission notice and warranty disclaimer appear in supporting │
│ documentation, and that the name of Lucent or any of its entities │
│ not be used in advertising or publicity pertaining to │
│ distribution of the software without specific, written prior │
│ permission. │
│ │
│ LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, │
│ INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. │
│ IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY │
│ SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES │
│ WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER │
│ IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, │
│ ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF │
│ THIS SOFTWARE. │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/assert.h"
#include "libc/fmt/fmt.h"
#include "libc/fmt/fmt.internal.h"
#include "libc/fmt/internal.h"
#include "libc/intrin/bsr.h"
#include "libc/macros.internal.h"
#include "libc/str/str.h"
#include "third_party/gdtoa/gdtoa.h"
/**
* @fileoverview Floating-Point Formatting
*
* This implements most of ANSI C's printf floating-point directives.
* omitting L, with %.0g and %.0G giving the shortest decimal string
* that rounds to the number being converted, and with negative
* precisions allowed for %f.
*/
struct FPBits {
uint32_t bits[4];
const FPI *fpi;
int sign;
int ex; // exponent
int kind;
};
union U {
double d;
uint64_t q;
long double ld;
unsigned int ui[4];
unsigned short us[5];
};
static const FPI kFpiDbl = {
.nbits = 53,
.emin = 1 - 1023 - 53 + 1,
.emax = 2046 - 1023 - 53 + 1,
.rounding = FPI_Round_near,
.sudden_underflow = 0,
};
static const FPI kFpiLdbl = {
.nbits = 64,
.emin = 1 - 16383 - 64 + 1,
.emax = 32766 - 16383 - 64 + 1,
.rounding = FPI_Round_near,
.sudden_underflow = 0,
};
static const char kSpecialFloats[2][2][4] = {
{"INF", "inf"},
{"NAN", "nan"},
};
static void dfpbits(union U *u, struct FPBits *b) {
int ex, i;
uint32_t *bits;
b->fpi = &kFpiDbl;
b->sign = u->ui[1] & 0x80000000L;
bits = b->bits;
bits[1] = u->ui[1] & 0xfffff;
bits[0] = u->ui[0];
if ((ex = (u->ui[1] & 0x7ff00000L) >> 20) != 0) {
if (ex == 0x7ff) {
// Infinity or NaN
i = bits[0] | bits[1] ? STRTOG_NaN : STRTOG_Infinite;
} else {
i = STRTOG_Normal;
bits[1] |= 0x100000;
}
} else if (bits[0] | bits[1]) {
i = STRTOG_Denormal;
ex = 1;
} else {
i = STRTOG_Zero;
}
b->kind = i;
b->ex = ex - (0x3ff + 52);
}
static void xfpbits(union U *u, struct FPBits *b) {
uint32_t *bits;
int ex, i;
b->fpi = &kFpiLdbl;
b->sign = u->us[4] & 0x8000;
bits = b->bits;
bits[1] = ((unsigned)u->us[3] << 16) | u->us[2];
bits[0] = ((unsigned)u->us[1] << 16) | u->us[0];
if ((ex = u->us[4] & 0x7fff) != 0) {
i = STRTOG_Normal;
if (ex == 0x7fff) // Infinity or NaN
i = bits[0] | bits[1] ? STRTOG_NaN : STRTOG_Infinite;
} else if (bits[0] | bits[1]) {
i = STRTOG_Denormal;
ex = 1;
} else {
i = STRTOG_Zero;
}
b->kind = i;
b->ex = ex - (0x3fff + 63);
}
// returns number of hex digits minus 1, or 0 for zero
static int fpiprec(struct FPBits *b) {
FPI *fpi;
int i, j, k, m;
uint32_t *bits;
if (b->kind == STRTOG_Zero) return (b->ex = 0);
fpi = b->fpi;
bits = b->bits;
for (k = (fpi->nbits - 1) >> 2; k > 0; --k) {
if ((bits[k >> 3] >> 4 * (k & 7)) & 0xf) {
m = k >> 3;
for (i = 0; i <= m; ++i)
if (bits[i]) {
if (i > 0) {
k -= 8 * i;
b->ex += 32 * i;
for (j = i; j <= m; ++j) {
bits[j - i] = bits[j];
}
}
break;
}
for (i = 0; i < 28 && !((bits[0] >> i) & 0xf); i += 4) donothing;
if (i) {
b->ex += i;
m = k >> 3;
k -= (i >> 2);
for (j = 0;; ++j) {
bits[j] >>= i;
if (j == m) break;
bits[j] |= bits[j + 1] << (32 - i);
}
}
break;
}
}
return k;
}
// round to prec hex digits after the "."
// prec1 = incoming precision (after ".")
static int bround(struct FPBits *b, int prec, int prec1) {
uint32_t *bits, t;
int i, inc, j, k, m, n;
m = prec1 - prec;
bits = b->bits;
inc = 0;
k = m - 1;
if ((t = bits[k >> 3] >> (j = (k & 7) * 4)) & 8) {
if (t & 7) goto inc1;
if (j && bits[k >> 3] << (32 - j)) goto inc1;
while (k >= 8) {
k -= 8;
if (bits[k >> 3]) {
inc1:
inc = 1;
goto haveinc;
}
}
}
haveinc:
b->ex += m * 4;
i = m >> 3;
k = prec1 >> 3;
j = i;
if ((n = 4 * (m & 7)))
for (;; ++j) {
bits[j - i] = bits[j] >> n;
if (j == k) break;
bits[j - i] |= bits[j + 1] << (32 - n);
}
else
for (;; ++j) {
bits[j - i] = bits[j];
if (j == k) break;
}
k = prec >> 3;
if (inc) {
for (j = 0; !(++bits[j] & 0xffffffff); ++j) donothing;
if (j > k) {
onebit:
bits[0] = 1;
b->ex += 4 * prec;
return 1;
}
if ((j = prec & 7) < 7 && bits[k] >> (j + 1) * 4) goto onebit;
}
for (i = 0; !(bits[i >> 3] & (0xf << 4 * (i & 7))); ++i) donothing;
if (i) {
b->ex += 4 * i;
prec -= i;
j = i >> 3;
i &= 7;
i *= 4;
for (m = j;; ++m) {
bits[m - j] = bits[m] >> i;
if (m == k) break;
bits[m - j] |= bits[m + 1] << (32 - i);
}
}
return prec;
}
int __fmt_dtoa(int (*out)(const char *, void *, size_t), void *arg, int d,
int flags, int prec, int sign, int width, bool longdouble,
char qchar, unsigned char signbit, const char *alphabet,
va_list va) {
double x;
union U u;
struct FPBits fpb;
char *s, *q, *se, *s0, special[8];
int c, k, i1, ui, bw, bex, sgn, prec1, decpt;
x = 0;
switch (d) {
case 'F':
case 'f':
if (!(flags & FLAGS_PRECISION)) prec = 6;
if (!longdouble) {
x = va_arg(va, double);
s = s0 = dtoa(x, 3, prec, &decpt, &fpb.sign, &se);
if (decpt == 9999) {
if (s && s[0] == 'N') {
fpb.kind = STRTOG_NaN;
} else {
fpb.kind = STRTOG_Infinite;
}
}
} else {
u.ld = va_arg(va, long double);
xfpbits(&u, &fpb);
s = s0 =
gdtoa(fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, 3, prec, &decpt, &se);
}
if (decpt == 9999) {
Format9999:
if (s0) freedtoa(s0);
bzero(special, sizeof(special));
s = q = special;
if (fpb.sign) {
*q++ = '-';
} else if (flags & FLAGS_PLUS) {
*q++ = '+';
} else if (flags & FLAGS_SPACE) {
*q++ = ' ';
}
memcpy(q, kSpecialFloats[fpb.kind == STRTOG_NaN][d >= 'a'], 4);
flags &= ~(FLAGS_PRECISION | FLAGS_PLUS | FLAGS_HASH | FLAGS_SPACE);
prec = 0;
return __fmt_stoa(out, arg, s, flags, prec, width, signbit, qchar);
}
FormatReal:
if (fpb.sign /* && (x || sign) */) sign = '-';
if (prec > 0) width -= prec;
if (width > 0) {
if (sign) --width;
if (decpt <= 0) {
--width;
if (prec > 0) --width;
} else {
if (s == se) decpt = 1;
width -= decpt;
if (prec > 0 || (flags & FLAGS_HASH)) --width;
}
}
if (width > 0 && !(flags & FLAGS_LEFT)) {
if ((flags & FLAGS_ZEROPAD)) {
if (sign) __FMT_PUT(sign);
sign = 0;
do __FMT_PUT('0');
while (--width > 0);
} else
do __FMT_PUT(' ');
while (--width > 0);
}
if (sign) __FMT_PUT(sign);
if (decpt <= 0) {
__FMT_PUT('0');
if (prec > 0 || (flags & FLAGS_HASH)) __FMT_PUT('.');
while (decpt < 0) {
__FMT_PUT('0');
prec--;
decpt++;
}
} else {
do {
if ((c = *s)) {
s++;
} else {
c = '0';
}
__FMT_PUT(c);
} while (--decpt > 0);
if (prec > 0 || (flags & FLAGS_HASH)) __FMT_PUT('.');
}
while (--prec >= 0) {
if ((c = *s)) {
s++;
} else {
c = '0';
}
__FMT_PUT(c);
}
while (--width >= 0) __FMT_PUT(' ');
if (s0) freedtoa(s0);
break;
case 'G':
case 'g':
if (!(flags & FLAGS_PRECISION)) prec = 6;
if (prec < 1) prec = 1;
if (!longdouble) {
x = va_arg(va, double);
s = s0 = dtoa(x, 2, prec, &decpt, &fpb.sign, &se);
if (decpt == 9999) {
if (s && s[0] == 'N') {
fpb.kind = STRTOG_NaN;
} else {
fpb.kind = STRTOG_Infinite;
}
}
} else {
u.ld = va_arg(va, long double);
xfpbits(&u, &fpb);
s = s0 = gdtoa(fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, prec ? 2 : 0, prec,
&decpt, &se);
}
if (decpt == 9999) goto Format9999;
c = se - s;
prec1 = prec;
if (!prec) {
prec = c;
prec1 = c + (s[1] || (flags & FLAGS_HASH) ? 5 : 4);
// %.0g gives 10 rather than 1e1
}
if (decpt > -4 && decpt <= prec1) {
if ((flags & FLAGS_HASH))
prec -= decpt;
else
prec = c - decpt;
if (prec < 0) prec = 0;
goto FormatReal;
}
d -= 2;
if (!(flags & FLAGS_HASH) && prec > c) prec = c;
--prec;
goto FormatExpo;
case 'e':
case 'E':
if (!(flags & FLAGS_PRECISION)) prec = 6;
if (prec < 0) prec = 0;
if (!longdouble) {
x = va_arg(va, double);
s = s0 = dtoa(x, 2, prec + 1, &decpt, &fpb.sign, &se);
if (decpt == 9999) {
if (s && s[0] == 'N') {
fpb.kind = STRTOG_NaN;
} else {
fpb.kind = STRTOG_Infinite;
}
}
} else {
u.ld = va_arg(va, long double);
xfpbits(&u, &fpb);
s = s0 = gdtoa(fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, prec ? 2 : 0, prec,
&decpt, &se);
}
if (decpt == 9999) goto Format9999;
FormatExpo:
if (fpb.sign /* && (x || sign) */) sign = '-';
if ((width -= prec + 5) > 0) {
if (sign) --width;
if (prec || (flags & FLAGS_HASH)) --width;
}
if ((c = --decpt) < 0) c = -c;
while (c >= 100) {
--width;
c /= 10;
}
if (width > 0 && !(flags & FLAGS_LEFT)) {
if ((flags & FLAGS_ZEROPAD)) {
if (sign) __FMT_PUT(sign);
sign = 0;
do __FMT_PUT('0');
while (--width > 0);
} else
do __FMT_PUT(' ');
while (--width > 0);
}
if (sign) __FMT_PUT(sign);
__FMT_PUT(*s++);
if (prec || (flags & FLAGS_HASH)) __FMT_PUT('.');
while (--prec >= 0) {
if ((c = *s))
s++;
else
c = '0';
__FMT_PUT(c);
}
__FMT_PUT(d);
if (decpt < 0) {
__FMT_PUT('-');
decpt = -decpt;
} else
__FMT_PUT('+');
for (c = 2, k = 10; 10 * k <= decpt; c++, k *= 10) donothing;
for (;;) {
i1 = decpt / k;
__FMT_PUT(i1 + '0');
if (--c <= 0) break;
decpt -= i1 * k;
decpt *= 10;
}
while (--width >= 0) __FMT_PUT(' ');
freedtoa(s0);
break;
case 'A':
alphabet = "0123456789ABCDEFPX";
goto FormatBinary;
case 'a':
alphabet = "0123456789abcdefpx";
FormatBinary:
if (longdouble) {
u.ld = va_arg(va, long double);
xfpbits(&u, &fpb);
} else {
u.d = va_arg(va, double);
dfpbits(&u, &fpb);
}
if (fpb.kind == STRTOG_Infinite || fpb.kind == STRTOG_NaN) {
s0 = 0;
goto Format9999;
}
prec1 = fpiprec(&fpb);
if ((flags & FLAGS_PRECISION) && prec < prec1) {
prec1 = bround(&fpb, prec, prec1);
}
bw = 1;
bex = fpb.ex + 4 * prec1;
if (bex) {
if ((i1 = bex) < 0) i1 = -i1;
while (i1 >= 10) {
++bw;
i1 /= 10;
}
}
if (fpb.sign /* && (sign || fpb.kind != STRTOG_Zero) */) {
sign = '-';
}
if ((width -= bw + 5) > 0) {
if (sign) --width;
if (prec1 || (flags & FLAGS_HASH)) --width;
}
if ((width -= prec1) > 0 && !(flags & FLAGS_LEFT) &&
!(flags & FLAGS_ZEROPAD)) {
do __FMT_PUT(' ');
while (--width > 0);
}
if (sign) __FMT_PUT(sign);
__FMT_PUT('0');
__FMT_PUT(alphabet[17]);
if ((flags & FLAGS_ZEROPAD) && width > 0 && !(flags & FLAGS_LEFT)) {
do __FMT_PUT('0');
while (--width > 0);
}
i1 = prec1 & 7;
k = prec1 >> 3;
__FMT_PUT(alphabet[(fpb.bits[k] >> 4 * i1) & 0xf]);
if (prec1 > 0 || (flags & FLAGS_HASH)) __FMT_PUT('.');
if (prec1 > 0) {
prec -= prec1;
while (prec1 > 0) {
if (--i1 < 0) {
if (--k < 0) break;
i1 = 7;
}
__FMT_PUT(alphabet[(fpb.bits[k] >> 4 * i1) & 0xf]);
--prec1;
}
if ((flags & FLAGS_HASH) && prec > 0) {
do __FMT_PUT(0);
while (--prec > 0);
}
}
__FMT_PUT(alphabet[16]);
if (bex < 0) {
__FMT_PUT('-');
bex = -bex;
} else
__FMT_PUT('+');
for (c = 1; 10 * c <= bex; c *= 10) donothing;
for (;;) {
i1 = bex / c;
__FMT_PUT('0' + i1);
if (!--bw) break;
bex -= i1 * c;
bex *= 10;
}
while (--width >= 0) __FMT_PUT(' ');
break;
default:
unreachable;
}
return 0;
}