/* cvt-double.c - converting floating point to and from strings
*
****************************************************************
* Copyright (C) 1998, 2000 Thomas Lord
*
* See the file "COPYING" for further information about
* the copyright and warranty status of this work.
*/
#include "hackerlab/bugs/panic.h"
#include "hackerlab/char/char-class.h"
#include "hackerlab/char/str.h"
#include "hackerlab/mem/mem.h"
#include "hackerlab/fmt/cvt.h"
#include "hackerlab/fmt/cvt-double.h"
/****************************************************************
* Serious Bugs (fix NOW!):
*
* Todo:
*
*/
/************************************************************************
*(h0 "Converting Floating Point to and from Strings"
* :include ("fmt/cvt-double.h"))
*
* This functionality is not complete yet.
*/
/************************************************************************
* Converting Floating Point To and From Strings (Internals)
*
* Conversion to the formats used by `printfmt' (similar to `printf')
* is in three stages:
*
* 1. Convert from double-precision floating point to
* binary-encoded base 100 semi-normalized floating point.
* 2. Convert from base 100 to ASCII-encoded decimal scientific
* notation.
* 3. Convert from scientific notation to the requested output
* format.
*
* Not supported (yet) is conversion to minimal-length formats.
* Instead, the precision of the output must be explicitly supplied by
* the caller. More decimal digits may be produced by the conversion
* than are warranted by the underlying floating-point precision (this
* is a bug that will be fixed in the future.)
*
* The initial conversion to base 100 works by adding (in base 100)
* powers of two.
*
* A double precision binary floating point number:
*
* b0 b1 b2 b3 ... x 2^K for k in [-1075...970]
* with bits b0, b1, ... in {0,1}
* is:
*
* b0 * 2^K + b1 * 2^(K-1) + ....
*
* To convert a binary floating point number to decimal, the values
* 2^K are computed using a multi-precision binary encoded base 100
* representation using a combination of table look-ups and easy
* power-of-two, single-digit (base 100) multiplications. Addition
* using that representation is trivial and so is conversion to
* decimal.
*
* The base-100 powers-of-two table is generated by
* "=pow2-base100-tables.scm" and its size can be changed staticly
* without modifying the code. The larger the table, the faster
* conversion runs. The table can be quite small and still produce
* good results. The algorithm also uses tables to handle
* single-digit overflow during addition and multiplication:
*
* quotient100_table[x] == x / 100 // for x in [0 .. 64 * 99]
* mod100_table[x] == x % 100 // for x in [0 .. 64 * 99]
*
* Table lookup in that range seems to be non-trivially faster than
* the divide instruction on the author's machine. Presumably
* table-lookup isn't slow on any machine.
*
* This technique is referred to in "How to Print Floating-Point
* Numbers Accurately" (Guy L. Steele and Jon L. White, `Proceedings
* of the ACM SIGPLAN '90 Conference on Programming Language Design
* and Implementation') as "gross overkill" but (assuming this
* algorithm can be trivially extended to support minimal-length
* formats) "gross overkill" is competitive on modern machines,
* especially in terms of simplicity and implementation size.
*
* This code is not especially optimized. Many uses of `mem_move'
* could be eliminated, probably resulting in a large speed-up. The
* addition and multiplication loops can probably be slightly
* improved.
*
* Conversion from decimal to binary floating point hasn't been
* implemented yet but will probably also work by brute-force: a
* precise conversion form decimal to multi-precision binary encoded
* base 100 floating point-(or some other base 10^n) followed by a
* rounding-conversion to base2 using multi-precision divides by
* powers of two.
*/
#include "hackerlab/fmt/cvt-double-constants.h"
struct base100_float
{
int base2_expt;
int base100_expt;
int n_digits;
t_uchar * digits;
};
#define pow2_enough_base100_digits 1024 /* larger than needed */
#define pow2_enough_base10_digits (2 * pow2_enough_base100_digits) /* larger than needed */
struct base100_float_temp
{
struct base100_float n;
t_uchar digit_space[pow2_enough_base100_digits];
};
#include "hackerlab/fmt/cvt-double-tables.h"
#if 0
#define MOD100(N) ((N)%100)
#define QUOT100(N) ((N)/100)
#else
#define MOD100(N) mod100_table[(N)]
#define QUOT100(N) quotient100_table[(N)]
#endif
static void
weird_pow2_base100 (struct base100_float * answer, int n, int n_in_answer)
{
int biased_n;
int nearest_in_table;
int pow2_from_table;
int pow2_from_answer;
if ((n < pow2_base100_min_expt) || (n > 1023))
panic ("pow2_base100 exponent out of range");
biased_n = n - pow2_base100_min_expt;
nearest_in_table = biased_n / powers_table_frequency;
pow2_from_table = biased_n % powers_table_frequency;
if ( (n_in_answer == pow2_base100_non_expt)
|| ((n - n_in_answer) > pow2_from_table))
{
answer->base2_expt = pow2_table[nearest_in_table].base2_expt;
answer->base100_expt = pow2_table[nearest_in_table].base100_expt;
answer->n_digits = pow2_table[nearest_in_table].n_digits;
mem_move (answer->digits, pow2_table[nearest_in_table].digits, answer->n_digits);
n_in_answer = nearest_in_table;
pow2_from_answer = pow2_from_table;
}
else
pow2_from_answer = n - n_in_answer;
while (pow2_from_answer)
{
int mult_digits;
t_uchar * mult;
t_uchar * prod;
t_uchar carry;
int log2_factor_from_answer;
if (pow2_from_answer <= 6)
{
log2_factor_from_answer = pow2_from_answer;
pow2_from_answer = 0;
}
else
{
log2_factor_from_answer = 6;
pow2_from_answer -= 6;
}
prod = answer->digits + pow2_enough_base100_digits - 1;
mult = answer->digits + answer->n_digits - 1;
mult_digits = answer->n_digits;
carry = 0;
while (mult_digits)
{
int digit_result;
digit_result = carry + (*mult << log2_factor_from_answer);
*prod = MOD100 (digit_result);
carry = QUOT100 (digit_result);
--prod;
--mult;
--mult_digits;
}
if (!carry)
++prod;
else
*prod = carry;
answer->n_digits = (answer->digits + pow2_enough_base100_digits) - prod;
mem_move (answer->digits, prod, answer->n_digits);
answer->base2_expt = n;
while (0 == answer->digits[answer->n_digits - 1])
{
++answer->base100_expt;
--answer->n_digits;
}
}
}
static void
weird_base100_add (struct base100_float * answer, struct base100_float * k)
{
int carry;
t_uchar * sum;
t_uchar * m;
t_uchar * n;
t_uchar * m_stop;
t_uchar * n_stop;
int m_digits;
int n_digits;
int digit;
m = answer->digits + answer->n_digits - 1;
m_stop = answer->digits;
m_digits = answer->n_digits;
n = k->digits + k->n_digits - 1;
n_stop = k->digits;
n_digits = k->n_digits;
sum = answer->digits + pow2_enough_base100_digits - 1;
carry = 0;
if (answer->base100_expt < k->base100_expt)
{
int x;
x = k->base100_expt - answer->base100_expt;
while ((m >= m_stop) && x)
{
*sum = *m;
--m;
--sum;
--x;
}
if (x)
while (x--)
{
*sum = 0;
--sum;
}
}
while ((m >= m_stop) && (n >= n_stop))
{
digit = *m + *n + carry;
*sum = MOD100 (digit);
carry = QUOT100 (digit);
--m;
--n;
--sum;
}
while (m >= m_stop)
{
digit = *m + carry;
*sum = MOD100 (digit);
carry = QUOT100 (digit);
--m;
--sum;
}
while (n >= n_stop)
{
digit = *n + carry;
*sum = MOD100 (digit);
carry = QUOT100 (digit);
--n;
--sum;
}
if (!carry)
++sum;
else
*sum = carry;
/* answer->base100_expt = k->base100_expt; */
answer->n_digits = (answer->digits + pow2_enough_base100_digits - sum);
mem_move (answer->digits, sum, answer->n_digits);
}
static void
double_to_base100 (struct base100_float * answer, double aDouble)
{
int raw_exponent;
int base2_exponent;
int implicit_bit;
unsigned long long fraction;
struct base100_float_temp temp;
union
{
double dbl;
unsigned long long lng;
} double_to_long;
if (aDouble == 0.0)
{
answer->base2_expt = 0;
answer->base100_expt = 0;
answer->n_digits = 1;
answer->digits[0] = 0;
return;
}
temp.n.digits = temp.digit_space;
double_to_long.dbl = aDouble;
raw_exponent = ((*(unsigned long long *)&double_to_long.lng) >> 52) & 0x7ff;
base2_exponent = (!raw_exponent ? 0 : raw_exponent - 1023) - 52;
implicit_bit = (raw_exponent ? 1 : 0);
fraction = ((*(unsigned long long *)&double_to_long.lng) & 0xfffffffffffffLL);
if (implicit_bit)
fraction |= 0x10000000000000LL;
if (!fraction)
panic ("unexpected nan in double_to_base100");
while (!(fraction & 1))
{
++base2_exponent;
fraction >>= 1;
}
weird_pow2_base100 (&temp.n, base2_exponent, pow2_base100_non_expt);
answer->base2_expt = temp.n.base2_expt;
answer->base100_expt = temp.n.base100_expt;
answer->n_digits = temp.n.n_digits;
mem_move (answer->digits, temp.n.digits, temp.n.n_digits);
while (fraction)
{
++base2_exponent;
fraction >>= 1;
if (fraction & 1)
{
weird_pow2_base100 (&temp.n, base2_exponent, temp.n.base2_expt);
weird_base100_add (answer, &temp.n);
}
}
}
/* double_to_e_decimal
*
* Fill buffer with:
*
* (-)?[0-9]\.[0-9]*[eE][+-](0[1-9]|[1-9][0-9]+)
*
* Giving as many digits as necessary to precisely represent
* `k', presuming that `k' is infinite precision. (A smaller
* number of digits may be enough to represent the actual bit
* pattern of (the finite precision value) `k'.)
*/
static int
double_to_e_decimal (t_uchar * buffer,
int * base10_expt,
double k,
char e)
{
struct base100_float_temp answer;
int d;
int b;
t_uchar * digs;
int expt10;
int shift;
answer.n.digits = answer.digit_space;
b = 0;
if (k < 0)
{
buffer[b++] = '-';
k = -k;
}
double_to_base100 (&answer.n, k);
digs = base100_digits_to_ascii_table[answer.n.digits[0]];
if (digs[0] != '0')
{
buffer[b++] = digs[0];
buffer[b++] = '.';
buffer[b++] = digs[1];
shift = 0;
}
else
{
buffer[b++] = digs[1];
buffer[b++] = '.';
shift = 1;
}
d = 1;
while (d < answer.n.n_digits)
{
buffer[b++] = base100_digits_to_ascii_table[answer.n.digits[d]][0];
buffer[b++] = base100_digits_to_ascii_table[answer.n.digits[d]][1];
++d;
}
expt10 = -shift + -1 + 2 * (answer.n.n_digits + answer.n.base100_expt);
*base10_expt = expt10;
buffer[b++] = e;
if (expt10 > 0)
buffer[b++] = '+';
cvt_long_to_decimal (buffer + b, (long)expt10);
{
int l;
l = str_length (buffer + b);
b += l;
if (l == 2)
{
buffer[b] = buffer[b - 1];
buffer[b - 1] = '0';
++b;
}
}
buffer[b] = 0;
return b;
}
static int
fix_e_decimal_precision (t_uchar * output, t_uchar * input, int precision)
{
int x;
int exponent_adjust;
int exponent_position_in_input = -1;
int exponent_position_in_output = -1;
/* From the FreeBSD man-page:
*
* eE The double argument is rounded and converted in the style
* [-]d.ddde+-dd where there is one digit before the decimal-point
* character and the number of digits after it is equal to the pre-
* cision; if the precision is missing, it is taken as 6; if the
* precision is zero, no decimal-point character appears. An E con-
* version uses the letter E (rather than e) to introduce the expo-
* nent. The exponent always contains at least two digits; if the
* value is zero, the exponent is 00.
*
*
* Thanks to `double_to_e_decimal' we already have:
*
* (-)?[0-9]\.[0-9]*[eE][+-](0[1-9]|[1-9][0-9]+)
*
* Where the mantissa contains all of the digits from conversion to decimal
* and the exponent contains as few digits as possible.
*
*/
exponent_adjust = 0;
if (precision == 0)
{
if ((input[2] >= '5') && (input[2] <= '9'))
{
if (input[0] == 9)
{
/* e.g.: 9.9e10 --> 1e11
*/
exponent_adjust = 1;
output[0] = 1;
}
else
{
/* e.g.: 8.9e10 --> 9e10
*/
output[0] = input[0] + 1;
}
}
else
{
output[0] = input[0];
}
for (x = 2; char_to_lower (input[x]) != 'e'; ++x)
;
output[1] = input[x];
exponent_position_in_input = x + 1;
exponent_position_in_output = 2;
/* `exponent_position_in_input' is the position of a decimal number.
*
* `exponent_position_in_output' is where to store that number,
* properly formatted.
*
* `exponent_adjust' is either 0 or 1 and should be added to the
* exponent from the input before formatting the output.
*
* The output contains a single digit decimal number and 'e' or 'E'.
*/
}
else
{
output[0] = input[0];
output[1] = '.';
x = 0;
while ( (x < precision)
&& char_is_digit (input[x + 2]))
{
output[x + 2] = input[x + 2];
++x;
}
if (x < precision)
{
while (x < precision)
output[x++] = '0';
}
else
{
if ( (input[x + 2] >= '5')
&& (input[x + 2] <= '9'))
{
int y;
int carry;
re_round:
carry = 1;
y = x;
while (carry && (y > 1))
{
int old_digit;
int new_carry;
int new_digit;
old_digit = output[y + 2] - '0';
new_carry = ((old_digit + carry) >= 10);
new_digit = (new_carry ? carry : old_digit + 1);
output[y + 2] = new_digit;
carry = new_carry;
--y;
}
if (carry)
{
if (output[0] != '9')
output[0]++;
else
{
mem_move (output + 3, output + 2, x - 2);
output[2] = '0';
output[1] = '1';
exponent_adjust = 1;
}
if ( (output[x + 2] >= '5')
&& (output[x + 2] <= '9'))
goto re_round;
}
}
}
/* `exponent_position_in_input' is the position of a decimal number.
*
* `exponent_position_in_output' is where to store that number,
* properly formatted.
*
* `exponent_adjust' is either 0 or 1 and should be added to the
* exponent from the input before formatting the output.
*
* The output contains a fraction with a single digit to the
* left of the decimal number and followed by 'e' or 'E'.
*/
}
/* `exponent_position_in_input' is the position of a decimal number.
*
* `exponent_position_in_output' is where to store that number,
* properly formatted.
*
* `exponent_adjust' is either 0 or 1 and should be added to the
* exponent from the input before formatting the output.
*/
{
int raw_exponent;
int real_exponent;
int abs_exponent;
int errn;
invariant (exponent_position_in_input != -1);
invariant (exponent_position_in_output != -1);
if (cvt_decimal_to_int (&errn, &raw_exponent,
input + exponent_position_in_input,
str_length (input + exponent_position_in_input)))
panic ("absurd exponent in `fix_e_decimal_precision'");
real_exponent = raw_exponent + exponent_adjust;
if (real_exponent >= 0)
{
output[exponent_position_in_output] = '+';
++exponent_position_in_output;
abs_exponent = real_exponent;
}
else
{
output[exponent_position_in_output] = '-';
++exponent_position_in_output;
abs_exponent = -real_exponent;
}
if (abs_exponent < 10)
{
output[exponent_position_in_output] = '0';
++exponent_position_in_output;
}
cvt_ulong_to_decimal (output + exponent_position_in_output,
abs_exponent);
return str_length (output);
}
}
static int
e_decimal_to_fp_decimal (t_uchar * fpd,
t_uchar * ed,
int base10_expt,
int precision)
{
int x;
int e;
int fix_ed;
x = 0;
e = 0;
fix_ed = 0;
if (base10_expt == 0)
{
fpd[x++] = ed[e++];
if (!precision)
{
done:
if (fix_ed)
{
ed[0] = ed[1];
ed[1] = '.';
}
while ((x > 1) && (fpd[x-2] != '.') && (fpd[x - 1] == '0'))
--x;
fpd[x] = 0;
return x;
}
fpd[x++] = ed[e++]; /* '.' */
}
else if (base10_expt > 0)
{
fpd[x++] = ed[e++];
e++; /* '.' */
while (base10_expt && (char_to_lower (ed[e]) != 'e'))
{
fpd[x++] = ed[e++];
--base10_expt;
}
while (base10_expt)
{
fpd[x++] = '0';
--base10_expt;
}
if (!precision)
goto done;
fpd[x++] = '.';
}
else
{
fpd[x++] = '0';
fpd[x++] = '.';
while (base10_expt < -1)
{
fpd[x++] = '0';
++base10_expt;
}
ed[1] = ed[0];
fix_ed = 1;
e = 1;
}
while (precision && (char_to_lower (ed[e]) != 'e'))
{
fpd[x++] = ed[e++];
--precision;
}
{
int z;
z = x - 1;
if ((ed[e] >= '5') && (ed[e] <= '9'))
{
while (z)
{
if (fpd[z] != '9')
{
++fpd[z];
goto done;
}
else
{
fpd[z] = '0';
--z;
}
}
mem_move (fpd, fpd+1, x - 1);
x += 1;
}
goto done;
}
}
static int
e_decimal_to_g_decimal (t_uchar * answer,
t_uchar * ed,
int base10_expt,
int precision)
{
#if 1
return -1;
#else
if ((base10_expt < -4) || (base10_expt >= precision))
{
str_cpy (answer, ed);
return fix_e_decimal_precision (answer, precision);
}
else
return e_decimal_to_fp_decimal (answer, ed, base10_expt, precision);
#endif
}
/*c
* int cvt_double_to_decimal (t_uchar * buffer,
* int * base10_expt,
* double k,
* int precision,
* char format);
*
* Convert `k' to a decimal number.
*
* [interface not stable]
*/
int
cvt_double_to_decimal (t_uchar * buffer,
int * base10_expt,
double k,
int precision,
char format)
{
t_uchar temp[pow2_enough_base10_digits];
int e;
int l;
int be;
if (!base10_expt)
base10_expt = &be;
if (precision < 0)
precision = 6;
if (format == 'E')
e = 'E';
else
e = 'e';
double_to_e_decimal (temp, base10_expt, k, e);
switch (format)
{
case 'e':
case 'E':
l = fix_e_decimal_precision (buffer, temp, precision);
str_cpy (buffer, temp);
return l;
case 'f':
return e_decimal_to_fp_decimal (buffer, temp, *base10_expt, precision);
case 'g':
return e_decimal_to_g_decimal (buffer, temp, *base10_expt, precision);
default:
while (1)
panic ("unrecognized floating-point format in cvt_double_to_decimal");
}
}
#ifdef TEST_CVT_DOUBLE
/* gcc -DTEST_CVT_DOUBLE -g -I.. -o ,test-cvt-double cvt-double.c ../=build/vu/libvu.a
*/
#include <stdlib.h>
#include <unistd.h>
#include <fcntl.h>
#include <sys/param.h>
#include <sys/types.h>
#include <sys/socket.h>
#include <sys/un.h>
#include <netinet/in.h>
#include <arpa/inet.h>
#include <netdb.h>
#include <errno.h>
#include "hackerlab/cmd/opt.h"
#include "hackerlab/char/char-class.h"
#include "hackerlab/arrays/ar.h"
#include "hackerlab/char/str.h"
#include "hackerlab/rx-posix/regex.h"
#include "hackerlab/bitsets/bitset.h"
#include "hackerlab/vu/printfmt.h"
#include "hackerlab/vu/vu.h"
#include "hackerlab/vu/vfdbuf.h"
#include "hackerlab/vu/url-fd.h"
#include "hackerlab/vu-network/url-socket.h"
#include "hackerlab/mem/must-malloc.h"
#include "hackerlab/hash/hashtree.h"
#include "hackerlab/fs/file-names.h"
#include "hackerlab/vu/safe.h"
#define OPTS(OP, OP2) \
OP (opt_help_msg, "h", "help", 0, \
"Display a help message and exit.") \
OP (opt_version, "V", "version", 0, \
"Display a release identifier string") \
OP2 (opt_version, 0, 0, 0, "and exit.") \
OP (opt_power, "p", "power=n", 1, \
"Display the nth power of 2, base 100.") \
OP (opt_float100, "c", "base100=n", 1, \
"Display N (a double) in base 100.") \
OP (opt_sci, "s", "sci=p,n", 1, \
"Display N (a double) in base 10 scientific") \
OP2 (opt_sci, 0, 0, 1, "notation with precision P.")
enum options
{
OPTS (OPT_ENUM, OPT_IGN)
};
struct opt_desc opts[] =
{
OPTS (OPT_DESC, OPT_DESC)
{-1, 0, 0, 0, 0}
};
static void
usage (int fd, t_uchar * prog, int suggest_help)
{
int errn;
if (0 > printfmt (&errn,
fd,
"usage: %s [options] [input-file]\n",
prog))
panic ("unable to write usage message");
if ( suggest_help
&& (0 > printfmt (&errn,
fd,
"try %s --help\n",
prog)))
panic ("unable to write usage message");
}
static void
help (int fd, char * prog)
{
int errn;
usage (fd, prog, 0);
vu_write (&errn, fd, "\n", 1);
opt_help (fd, opts);
}
static void
version (int fd)
{
int errn;
if (0 > printfmt (&errn, fd, "cvt-double test scaffolding\n"))
panic ("unable to write release identifier");
}
static void
print_base100 (int fd, struct base100_float * n)
{
int x;
/* safe_printfmt (fd, "(00)."); */
for (x = 0; x < n->n_digits; ++x)
safe_printfmt (fd, "(%s)", base100_digits_to_ascii_table[n->digits[x]]);
safe_printfmt (fd, " x 100^%d", n->base100_expt);
}
int
main (int argc, char * argv[])
{
int errn;
int o;
t_uchar * opt_string;
t_uchar * opt_arg;
t_uchar opt_string_space[3];
struct base100_float_temp temp_power;
int last_power_of_2;
last_power_of_2 = pow2_base100_non_expt;
if (!argc)
argv0 = "(argc == 0\?\?\?)";
else
argv0 = file_name_tail (lim_use_must_malloc, argv[0]);
while (1)
{
o = opt_low (&opt_string, &opt_arg, 0, opts, &argc, argv, opt_string_space);
if (o == opt_none)
break;
switch (o)
{
default:
if (argc < 1)
panic ("internal error parsing arguments");
printfmt (&errn, 2, "unhandled option `%s'\n", argv[1]);
usage_error:
usage (2, argv[0], 1);
exit (1);
missing_arg:
printfmt (&errn, 2, "missing argument for `%s'\n", opt_string);
goto usage_error;
bogus_arg:
printfmt (&errn, 2, "ill-formed argument for `%s' (`%s')\n", opt_string, opt_arg);
goto usage_error;
case opt_help_msg:
help (1, argv[0]);
exit (0);
case opt_version:
version (1);
exit (0);
case opt_power:
{
int power_of_2;
if (!opt_arg)
goto missing_arg;
if (cvt_decimal_to_int (&errn, &power_of_2, opt_arg, str_length (opt_arg)))
goto bogus_arg;
temp_power.n.digits = temp_power.digit_space;
weird_pow2_base100 (&temp_power.n, power_of_2, last_power_of_2);
last_power_of_2 = power_of_2;
print_base100 (1, &temp_power.n);
break;
}
case opt_float100:
{
double n;
struct base100_float_temp temp;
if (!opt_arg)
goto missing_arg;
n = strtod (opt_arg, 0);
temp.n.digits = temp.digit_space;
double_to_base100 (&temp.n, n);
print_base100 (1, &temp.n);
break;
}
case opt_sci:
{
t_uchar * comma;
int precision;
double n;
t_uchar cvt[1024];
int base10_expt;
if (!opt_arg)
goto missing_arg;
comma = str_chr_index (opt_arg, ',');
if (!comma)
goto bogus_arg;
if (cvt_decimal_to_uint (&errn, &precision, opt_arg, comma - opt_arg))
goto bogus_arg;
n = strtod (comma + 1, 0);
if (precision > 512)
{
panic ("unsupported precision");
}
cvt_double_to_decimal (cvt, &base10_expt, n, precision, 'e');
printfmt (&errn, 1, "%s\n", cvt);
break;
}
}
}
return 0;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1