/* 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 #include #include #include #include #include #include #include #include #include #include #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