I wanted to restart the big-float project I had paused the moment I suddenly realized that the implementation of the main four rounding modes (which I postponed so long, that I got all of the IEEE-754 functions implemented in the meantime) needs three guard bits. And it needs them from the beginning. Which needs a complete refactoring. Yippy.
But that might also be a good chance to clean the whole thing up, I pondered--the function that parses the input is quite a mess, for example. Parsing is normally done with a small program: a lexer: put some regular expressions in and get your stuff out in the right order, very simple. But it has the disadvantage to need a third party program--or an "app" as the kids say today--and that is a dependency to be avoided -- a standard compliant C-compiler must suffice. Which exact C-standard the term "standard compliant" talks about is open for negotiations but it should not be older than, say, 15 years.
Writing a finite state machine manually is probably not the brightest idea I had in my life but the actual amount of such a task is not always easy to approximate and there are worse things to do on a weekend (although I really need to repair the heating before the start of winter).
The format of the input string is best described formally, to avoid any ambiguities, so what follows is a description of the acceptable input in EBNF (ISO 14977).
(* space (0x20), used only for the thousands separators *)
space = ? US-ASCII character 32 ?;
(* sign *)
sign = '+'|'-';
(* thousands separator, must be between two digits *)
tsep = '_' | space;
(* A leading zero is treated as a prefix, hence the special treatment for it *)
zero = '0';
(* Integers. Integers must contain at least one digit of the respective base *)
(* binary digit *)
bindig = '1';
binnum = (bindig | zero) | (bindig | zero), tsep, ( bindig | zero );
(* binary integer *)
binint = binnum, {binnum};
(* octal digit *);
octdig = bindig|'2'|'3'|'4'|'5'|'6'|'7';
octnum = (octdig | zero) | (octdig | zero), tsep, ( octdig | zero );
(* octal integer *)
octint = octnum, {octnum};
(* decimal digit *)
decdig = octdig|'8'|'9';
decnum = (decdig | zero) | (decdig | zero), tsep, ( decdig | zero );
(* decimal integer (no leading zero) *)
decint = decdig, tsep, decnum, {decnum} | decdig, {decnum};
(* hexadecimal digit *)
hexdig = decdig |'A'|'B'|'C'|'D'|'E'|'F'|'a'|'b'|'c'|'e'|'f';
hexnum = (hexdig | zero) | (hexdig | zero), tsep, ( hexdig | zero );
(* hexadecimal integer *)
hexint = hexnum, {hexnum};
(* prefix for the bases *)
prefix = '0';
(* Implementing only these four bases has been deemed sufficient *)
binbase = prefix, ('B'|'b');
octbase = prefix;
hexbase = prefix, ('X'|'x');
(* decimal point, no locale specific variations allowed *)
decpoint = '.';
(* exponents (no thousand delimiters) *)
expodig = (decdig | zero),{(decdig | zero)};
expobin = ('P'|'p'), [sign], expodig;
expodec = ('E'|'e'), [sign], expodig;
(* "xyz", "xyz.", "xyz.zyx", ".zyx" *)
stubrealbin = binint | (binint, decpoint, [binint]) | (decpoint, binint);
stubrealoct = octint | (octint, decpoint, [octint]) | (decpoint, octint);
stubrealdec = decint | (decint, decpoint, [decint]) | (decpoint, decint);
stubrealhex = hexint | (hexint, decpoint, [hexint]) | (decpoint, hexint);
(* No thousands separators after the sign *)
fullrealbin = [sign], binbase, stubrealbin, [expobin];
fullrealoct = [sign], octbase, stubrealoct, [expobin];
fullrealdec = [sign], stubrealdec, [expodec];
fullrealhex = [sign], hexbase, stubrealhex, [expobin];
(* "Inf" and "NaN" are case insensitive in praxi; quiet NaN only, no payload *)
real = ([sign], "Inf") | "NaN" | fullrealbin | fullrealoct | fullrealdec | fullrealhex;
For those you who are more familiar with regular expressions (even if it's only a "sort of" one as the one listed below):
msep = [_ ]
dectd = [0-9] msep [0-9]+ | [0-9]+
decint = [1-9][0-9]* | [1-9] msep dectd | 0
decreal = decint | decint "." dectd? | "." dectd
decexpo = [eE][+-]? [0-9]+
decimal = [+-]? decreal decexpo?
hexint = [0-9a-fA-F] msep [0-9a-fA-F]+ | [0-9a-fA-F]+
hexreal = hexint | hexint "." hexint? | "." hexint
hexexpo = [pP][+-]? [0-9]+
hexadecimal = [+-]? hexreal hexexpo?
octint = [0-7] msep [0-7]+ | [0-7]+
octreal = octint | octint "." octint? | "." octint
octexpo = [pP][+-]? [0-9]+
octal = [+-]? octreal octexpo?
binint = [01] msep [01]+ | [01]+
binreal = binint | binint "." binint? | "." binint
binexpo = [pP][+-]? [0-9]+
binary = [+-]? binreal binexpo?
real = [+-]?"inf" | "nan" | decimal | hexadecimal | octal | binary
Should work in Flex as written when the variables get their correct treatment (e.g.: "{variable2bexpanded}") but I have not tested it.
Compile with (The author's GCC is version 4.9):
gcc -O3 -g3 -W -Wall -Wextra -std=c11 -DDEBUG fsm_parse_real.c -o fsm_parse_real -lm
Or, if you don't have/want to use libmath:
gcc -O3 -g3 -W -Wall -Wextra -std=c11 -DDEBUG -DWITHOUT_MATHLIB fsm_parse_real.c -o fsm_parse_real
The results are different from libmath's powl()
but the actual error is the same, under one ULP. Not much, but still not correctly rounded, the error needs to be smaller than half an ULP to be correctly rounded.
You need to the necessary options to define BIG_ENDIAN
if you run it on a big endian machine and want to print the bits of the result. The rest is endian-agnostic.
The code itself contains more information -- some of it is also hidden in the comments.
The task of the code is to parse a string containing a floating point number at the beginning or fail if the string does not contain a floating point number at the beginning.
The acceptable format of the floating point number is described in the EBNF and the (sort of) regular expression above. All three, EBNF, regexp and code must not differ. In case of discrepancies, EBNF wins.
Note: The code for the actual computation of a numeric result -- the replacement for the big-float -- shall count as not much more than "just slapped on for debugging", the space restrictions of Stack Exchange do not allow me to add the big-integer code necessary to get rid of the dependency of the "long double" which is quite a disappointment. Otherwise it would be a nice extension to strtod
.
Paraphrasing David L. "Dave" Jones of EEV-blog fame:
Don't run it, take it apart!
#ifndef __STDC_IEC_559__
#error "This program needs IEEE-754/IEC-60559 compliant floating point arithmetic."
#endif
// The older, 16-bit versions of MSVS supported it, IIRC, and it is not fully clear why they
// dropped it later (probably a business decision).
// Se e.g.: https://software.intel.com/en-us/forums/intel-c-compiler/topic/277079
//
// The size of a "long double" may vary (it's called "extended precision" in the standard
// whereas all other data-types have a bit length attached) but if the size is not *more*
// than 8 (eight) (assuming an IEEE-754 binary64 for the "double" and 8-bit bytes) it makes
// no sense to use it, we need the extra bits, the results are way off with a binary64 alone.
#ifdef _MSC_VER
#error "Data-type 'long double' not fully supported by MS-Visual Studio"
// For other compilers (but needs C11 to work)
#elif __STDC_VERSION__ >= 201112L
_Static_assert(sizeof(double) < sizeof(long double),
"\"long double\" is not larger than \"double\"");
#else
#include <float.h>
#if LDBL_MANT_DIG == DBL_MANT_DIG
#error "\"long double\" seems not to be larger than \"double\""
#endif
// test above may fail under some circumstances, so at least utter a warning
#warning "Please check manually if \"long double\" is larger than \"double\""
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// only some macros are used from math.h if WITHOUT_MATHLIB is defined
#include <math.h>
// difficult to define INFINITY portably,
// using predefined macros from LibC's math.h instead
#ifndef INFINITY
// Yes, that's correct, a float.
// See ISO/IEC 9899:2011 sec. 7.12 par. 4
#define INFINITY HUGE_VALF
#endif
#ifdef WITHOUT_MATHLIB
static long double local_powl(long double base, int exponent)
{
long double power = 1.0L;
// Not needed/used here
/*
if (base == 0.0L) {
if (exponent > 0) {
return 0.0L;
} else if (exponent == 0) {
return 1.0L;
} else {
// +/- inf
return 1.0L / base;
}
}
if (exponent < 0) {
exponent = -exponent;
base = 1.0L / base;
}
*/
while (exponent) {
if (exponent % 2 == 1) {
power *= base;
}
exponent >>= 1;
base *= base;
}
return power;
}
#define POWER(x,y) local_powl((x),(y))
#else
#define POWER(x,y) powl((x),(long double)(y))
#endif
// strncasecmp() is not in the current (ISO/IEC 9899:2011) C-standard
#if !( defined _BSD_SOURCE \
|| defined _DEFAULT_SOURCE \
|| _POSIX_C_SOURCE >= 200112L \
|| defined YESIHAVESTRNCASECMPSOSHUTUPANDGETOFFMYLAWN)
#include <ctype.h>
static int strncasecmp(const char *s1, const char *s2, size_t n)
{
char c1 = 0;
char c2 = 0;
while (n--) {
c1 = tolower(*s1);
c2 = tolower(*s2);
if (c1 != c2) {
break;
}
if (c1 == '\0') {
break;
}
s1++;
s2++;
}
return (int) (c1 - c2);
}
#endif
// TODO: check if it compiles at all in a recent MSVC version
#ifdef _MSC_VER
#define strncasecmp(x,y,z) _strnicmp((x),(y),(z))
#endif
// Checks for over/underflow needed for computing the exponent
#include <limits.h>
#include <errno.h>
// Full checks not necessary in add_int, both inputs are positive
static inline int add_int(int a, int b)
{
if (((b > 0) && (a > (INT_MAX - b))) || ((b < 0) && (a < (INT_MIN - b)))) {
errno = ERANGE;
return INT_MAX;
} else {
return a + b;
}
}
static inline int sub_int(int a, int b)
{
if (((b > 0) && (a < (INT_MIN + b))) || ((b < 0) && (a > (INT_MAX + b)))) {
errno = ERANGE;
return INT_MIN;
} else {
return a - b;
}
}
// We do not know if the system will have a data-type
// that is larger than 'int', so something like
// long long prod = (long long)a * (long long)b
// may or may not work.
// Algorithm shamelessly stolen from cert.org
static inline int mul_int(int a, int b)
{
if (a > 0) {
if (b > 0) {
if (a > (INT_MAX / b)) {
errno = ERANGE;
return INT_MAX;
}
} else {
if (b < (INT_MIN / a)) {
errno = ERANGE;
return INT_MAX;
}
}
} else {
if (b > 0) {
if (a < (INT_MIN / b)) {
errno = ERANGE;
return INT_MAX;
}
} else {
if ((a != 0) && (b < (INT_MAX / a))) {
errno = ERANGE;
return INT_MAX;
}
}
}
return a * b;
}
#define FSM_OK 1
#define FSM_ERROR 0
#define FSM_EXPRANGE -1
enum fsm_input {
ZERO, // '0'
BINDIG, // '1'
OCTDIG, // '2'|'3'|'4'|'5'|'6'|'7'
DECDIG, // '8'|'9'
HEXDIG, // 'A' |'C'|'D' |'F'
PREHEX, // 'x'
PREBIN, // 'b'
EXPDEC, // 'e'
EXPBIN, // 'p'
DECPNT, // '.'
SGNCHR, // '+' | '-'
MILLSEP, // '_' | ' '
EOS, // '\0'
OTHER
};
enum fsm_states {
START, SIGN, PREFIX, TDBINB, TDOCTB,
TDDECB, TDHEXB, TDBINF, TDOCTF, TDDECF,
TDHEXF, HEXSTART, BINSTART, HEXFRCST, BINFRCST,
OCTFRCST, DECFRCST, DECPSTRT, DECPOINT, DECFRAC,
BINBASE, OCTBASE, DECBASE, HEXBASE, BINFRAC,
OCTFRAC, HEXFRAC, EXPOMARK, EXPOSIGN, EXPONENT,
ERROR, END
};
#ifdef DEBUG
static const char *st2str[32] = {
"START", "SIGN", "PREFIX", "TDBINB", "TDOCTB",
"TDDECB", "TDHEXB", "TDBINF", "TDOCTF", "TDDECF",
"TDHEXF", "HEXSTART", "BINSTART", "HEXFRCST", "BINFRCST",
"OCTFRCST", "DECFRCST", "DECPSTRT", "DECPOINT", "DECFRAC",
"BINBASE", "OCTBASE", "DECBASE", "HEXBASE", "BINFRAC",
"OCTFRAC", "HEXFRAC", "EXPOMARK", "EXPOSIGN", "EXPONENT",
"ERROR", "END"
};
static const char *type2str[14] = {
"ZERO", "BINDIG", "OCTDIG", "DECDIG", "HEXDIG",
"PREHEX", "PREBIN", "EXPDEC", "EXPBIN", "DECPNT",
"SGNCHR", "MILLSEP", "EOS", "OTHER"
};
#endif
// transition-table (transitions of state END not included)
static int fsm_table[31][13] = {
// ZERO, BINDIG, OCTDIG, DECDIG, HEXDIG, PREHEX, PREBIN,
// EXPDEC, EXPBIN, DECPNT, SGNCHR, MILLSEP, EOS
// START
{PREFIX, DECBASE, DECBASE, DECBASE, ERROR, ERROR, ERROR,
ERROR, ERROR, DECPSTRT, SIGN, ERROR, ERROR },
// SIGN
{PREFIX, DECBASE, DECBASE, DECBASE, HEXBASE, ERROR, HEXBASE,
HEXBASE, ERROR, DECPOINT, ERROR, ERROR, ERROR},
//PREFIX
{OCTBASE, OCTBASE, OCTBASE, ERROR, ERROR, HEXSTART, BINSTART,
END, END, DECPOINT, ERROR, ERROR, END},
//TDBINB
{BINBASE, BINBASE, ERROR, ERROR, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//TDOCTB
{OCTBASE, OCTBASE, OCTBASE, ERROR, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//TDDECB
{DECBASE, DECBASE, DECBASE, DECBASE, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//TDHEXB
{HEXBASE, HEXBASE, HEXBASE, HEXBASE, HEXBASE, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//TDBINF
{BINFRAC, BINFRAC, ERROR, ERROR, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//TDOCTF
{OCTFRAC, OCTFRAC, OCTFRAC, ERROR, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//TDDECF
{DECFRAC, DECFRAC, DECFRAC, DECFRAC, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//TDHEXF
{HEXFRAC, HEXFRAC, HEXFRAC, HEXFRAC, HEXFRAC, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//HEXSTART
{HEXBASE, HEXBASE, HEXBASE, HEXBASE, HEXBASE, ERROR, HEXBASE,
HEXBASE, ERROR, HEXFRCST, ERROR, ERROR, END},
//BINSTART
{BINBASE, BINBASE, ERROR, ERROR, ERROR, ERROR, ERROR,
ERROR, ERROR, BINFRCST, ERROR, ERROR, END},
//HEXFRCST
{HEXFRAC, HEXFRAC, HEXFRAC, HEXFRAC, HEXFRAC, ERROR, HEXFRAC,
HEXFRAC, ERROR, ERROR, ERROR, ERROR, END},
//BINFRCST
{BINFRAC, BINFRAC, ERROR, ERROR, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, END},
//OCTFRCST
{OCTFRAC, OCTFRAC, OCTFRAC, ERROR, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, END},
//DECFRCST
{DECFRAC, DECFRAC, DECFRAC, DECFRAC, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, END},
//DECPSTRT
{DECFRAC, DECFRAC, DECFRAC, DECFRAC, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, END},
//DECPOINT
{DECFRAC, DECFRAC, DECFRAC, DECFRAC, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, END},
//DECFRAC
{DECFRAC, DECFRAC, DECFRAC, DECFRAC, ERROR, ERROR, ERROR,
EXPOMARK, ERROR, ERROR, ERROR, TDDECF, END},
//BINBASE
{BINBASE, BINBASE, ERROR, ERROR, ERROR, ERROR, ERROR,
ERROR, EXPOMARK, BINFRCST, ERROR, TDBINB, END},
//OCTBASE
{OCTBASE, OCTBASE, OCTBASE, ERROR, ERROR, ERROR, ERROR,
ERROR, EXPOMARK, OCTFRCST, ERROR, TDOCTB, END},
//DECBASE
{DECBASE, DECBASE, DECBASE, DECBASE, ERROR, ERROR, ERROR,
EXPOMARK, ERROR, DECFRCST, ERROR, TDDECB, END},
//HEXBASE
{HEXBASE, HEXBASE, HEXBASE, HEXBASE, HEXBASE, ERROR, HEXBASE,
HEXBASE, EXPOMARK, HEXFRCST, ERROR, TDHEXB, END},
//BINFRAC
{BINFRAC, BINFRAC, ERROR, ERROR, ERROR, ERROR, ERROR,
ERROR, EXPOMARK, ERROR, ERROR, BINFRAC, END},
//OCTFRAC
{OCTFRAC, OCTFRAC, OCTFRAC, ERROR, ERROR, ERROR, ERROR,
ERROR, EXPOMARK, ERROR, ERROR, TDOCTF, END},
//HEXFRAC
{HEXFRAC, HEXFRAC, HEXFRAC, HEXFRAC, HEXFRAC, ERROR, HEXFRAC,
HEXFRAC, EXPOMARK, ERROR, ERROR, TDHEXF, END},
//EXPOMARK
{EXPONENT, EXPONENT, EXPONENT, EXPONENT, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, EXPOSIGN, ERROR, ERROR},
//EXPOSIGN
{EXPONENT, EXPONENT, EXPONENT, EXPONENT, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR},
//EXPONENT
{EXPONENT, EXPONENT, EXPONENT, EXPONENT, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, END},
//ERROR
{ERROR, ERROR, ERROR, ERROR, ERROR, ERROR, ERROR,
ERROR, ERROR, ERROR, ERROR, ERROR, ERROR}
};
#include <ctype.h>
static int check_type(char c)
{
int type;
switch (tolower(c)) {
case '0':
type = ZERO;
break;
case '1':
type = BINDIG;
break;
case '2':
case '3':
case '4':
case '5':
case '6':
case '7':
type = OCTDIG;
break;
case '8':
case '9':
type = DECDIG;
break;
case 'a':
case 'c':
case 'd':
case 'f':
type = HEXDIG;
break;
case 'x':
type = PREHEX;
break;
case 'b':
type = PREBIN;
break;
case 'e':
type = EXPDEC;
break;
case 'p':
type = EXPBIN;
break;
case '.':
type = DECPNT;
break;
case '+':
case '-':
type = SGNCHR;
break;
case ' ':
case '_':
type = MILLSEP;
break;
case '\0':
type = EOS;
break;
default:
type = OTHER;
break;
};
return type;
}
// Some variables to hold the relevant parts of a real number
// Thread safety is not an issue, parsing gets always done linearly.
// And if it is: just put all of it in one large struct and
// push that around
static int main_sign = 1;
static int expo_sign = 1;
static long double integral_part = 0.0L;
static long double fractional_part = 0.0L;
// Could also be done with one variable holding the number of digits and
// another one holding the type.
static int hex_frac = 0;
static int dec_frac = 0;
static int oct_frac = 0;
static int bin_frac = 0;
static int dec_int = 0;
static int exponent_part = 0;
#define HAS_EXP 1
#define NO_EXP 0
static int has_exponent = NO_EXP;
// Check for 7-bit ASCII (spot-checks only)
// Will not recognize close-but-no-cigar encodings, but almost all of the other encodings
// that are still in "widespread" use are the EBCDIC encodings
// No check for '9' because non-consecutive digits are forbidden by the C-standard
#if ('0' != 0x30) || ('a' != 0x61) || ('z' != 0x7a) || ('A' != 0x41) || ('Z' != 0x5a)
#error "Non-ASCII input encoding found, please change map below accordingly."
#endif
// assuming ASCII. Full map for more bases in the future.
// Far future. Very far future. Very, very far future.
// Close to the heat death of the universe, even.
static const char digit_map[] = {
-1, -1, -1, -1, -1, -1, -1, -1, // 0x00-0x07
-1, -1, -1, -1, -1, -1, -1, -1, // 0x08-0x0f
-1, -1, -1, -1, -1, -1, -1, -1, // 0x10-0x17
-1, -1, -1, -1, -1, -1, -1, -1, // 0x18-0x1f
-1, -1, -1, -1, -1, -1, -1, -1, // 0x20-0x27
-1, -1, -1, 1, -1, -1, -1, -1, // 0x28-0x2f '+' = 0x2b, '-' = 0x2d
0, 1, 2, 3, 4, 5, 6, 7, // 0x30-0x37 '0' - '7'
8, 9, -1, -1, -1, -1, -1, -1, // 0x38-0x3F '8', '9'
-1, 10, 11, 12, 13, 14, 15, -1, // 0x40-0x47 'A' - 'F'
-1, -1, -1, -1, -1, -1, -1, -1, // 0x48-0x4f
-1, -1, -1, -1, -1, -1, -1, -1, // 0x50-0x57
-1, -1, -1, -1, -1, -1, -1, -1, // 0x58-0x5f
-1, 10, 11, 12, 13, 14, 15, -1, // 0x60-0x67 'a' - 'f'
-1, -1, -1, -1, -1, -1, -1, -1, // 0x68-0x6f
-1, -1, -1, -1, -1, -1, -1, -1, // 0x70-0x77
-1, -1, -1, -1, -1, -1, -1, -1 // 0x78-0x7f
};
#define GETDIG(x) ( (long double)digit_map[(size_t)(x)] )
static int fsm(char input, int *state)
{
int res = FSM_OK;
input = tolower(input);
#ifdef DEBUG
printf("INPUT: 0x%2x, map: (%d), STATE: %-10s\n", input,
digit_map[(size_t)input], st2str[*state]);
#endif
switch (*state) {
case SIGN:
main_sign = (int)digit_map[(size_t)input];
break;
case DECBASE:
// no use of fmal() here and below as this is a) for testing only
// and b) will use big-floats in production anyways
integral_part *= 10.0L;
integral_part += GETDIG(input);
// keep number of integer digits (for decimals only)
dec_int++;
break;
case BINBASE:
integral_part *= 2.0L;
integral_part += GETDIG(input);
break;
case OCTBASE:
integral_part *= 8.0L;
integral_part += GETDIG(input);
break;
case HEXBASE:
integral_part *= 16.0L;
integral_part += GETDIG(input);
break;
case BINFRAC:
fractional_part *= 2.0L;
fractional_part += GETDIG(input);
bin_frac++;
break;
case OCTFRAC:
fractional_part *= 8.0L;
fractional_part += GETDIG(input);
oct_frac++;
break;
case DECFRAC:
// put all in one basket to avoid one division
//fractional_part *= 10.0L;
//fractional_part += GETDIG(input);
integral_part *= 10.0L;
integral_part += GETDIG(input);
dec_frac++;
break;
case HEXFRAC:
fractional_part *= 16.0L;
fractional_part += GETDIG(input);
hex_frac++;
break;
case EXPOMARK:
has_exponent = HAS_EXP;
break;
case EXPOSIGN:
expo_sign = (int)digit_map[(size_t)input];
break;
case EXPONENT:
// The variable "exponent_part" is a small native integer, so check for overflow
errno = 0;
exponent_part = mul_int(exponent_part, 10);
if (exponent_part == INT_MAX || errno == ERANGE) {
*state = ERROR;
res = FSM_EXPRANGE;
break;
}
exponent_part = add_int(exponent_part, (int)digit_map[(size_t)input]);
if (exponent_part == INT_MAX && errno == ERANGE) {
*state = ERROR;
res = FSM_EXPRANGE;
}
break;
case END:
// not used
break;
case ERROR:
// "Run in circles, scream and shout!"
// Infantry Journal, Vol. 35, p. 396, United States Infantry Association, 1929
fprintf(stderr, "ERROR state reached\n");
*state = ERROR;
res = FSM_ERROR;
break;
default:
// all the other "intermediate" states
break;
}
return res;
}
// a small helper to take some burden from the parser
static char *trim_both(char *s)
{
char *end, *p;
p = s;
while (isspace(*p)) {
p++;
}
if (*p == '\0') {
return p;
}
end = p + strlen(p) - 1;
while (end > p && isspace(*end)) {
end--;
}
*(end + 1) = '\0';
return p;
}
int str2dbl(char *s, double *d)
{
int cur_state = START;
int type;
int res;
long double ld = 0.0L;
s = trim_both(s);
if (*s == '\0') {
// empty input, would an error be better?
*d = 0.0;
return FSM_OK;
}
// This is not fully IEEE-754 conforming, because "inffoobar" would also parse as
// infinity but GlibC's (2.19 on the author's machine) strtod() does the same
// and the C-standard (current: 9899:2011) is also OK with it.
if (!strncasecmp(s, "-inf", 4)) {
*d = -INFINITY;
return FSM_OK;
}
if (!strncasecmp(s, "+inf", 4) || !strncasecmp(s, "inf", 3)) {
*d = INFINITY;
return FSM_OK;
}
// quiet NAN only, at least for now
if (!strncasecmp(s, "-nan", 4) || !strncasecmp(s, "+nan", 4)
|| !strncasecmp(s, "nan", 3)) {
// Macro NAN might not be defined. See ISO/IEC 9899:2011 sec. 7.2 par. 5
// nan() is in the standard since C99. See ISO/IEC 9899:2011 7.12.11.2
*d = nan("0");
// nan() returns zero if quiet NaNs are not supported. Does that mean
// that signaling ones are?
if(*d == 0.0){
#ifdef DEBUG
fprintf(stderr, "Quiet NaN not supported\n");
#endif
return FSM_ERROR;
}
return FSM_OK;
}
while (cur_state != END) {
type = check_type(*s);
#ifdef DEBUG
printf("CHAR \"%c\" (0x%x), TYPE: %s\n", *s, *s, type2str[type]);
#endif
if (type == OTHER) {
#ifdef DEBUG
fprintf(stderr, "OTHER: %c\n", *s);
#endif
return FSM_ERROR;
}
cur_state = fsm_table[cur_state][type];
res = fsm(*s, &cur_state);
if (res != FSM_OK) {
// exponent over/underflow.
// Still undecided if it better be an error
if (res == FSM_EXPRANGE) {
if (expo_sign < 0) {
*d = -0.0;
} else {
*d = main_sign * INFINITY;
}
}
return res;
}
s++;
}
#ifdef DEBUG
printf("exp; %d, intp: %.20Lg, fracp: %.20Lg\n", exponent_part, integral_part,
fractional_part);
#endif
// The integrity of an integer shall not be questioned but build upon.
ld += integral_part;
/*
* Use of "long double" does not always help, e.g.: 3.571e266 and
* 3.08984926168550152811e-32 are wrong with local_powl();
* 123.123123123123123123e123 and 123.1231231231231231e123 fail
* with libmath's powl() whereas 123.12312312312312312e123 and
* 123.123123123123123123123e123 work. Interestingly, the latter four work
* with local_powl().
* (First two examples from Rick Regan's article at
* http://www.exploringbinary.com/decimal-to-floating-point-needs-arbitrary-precision/)
*
* But it's still less than one ulp--good enough for testing the parser
*/
if (dec_int != 0 || dec_frac != 0) {
exponent_part *= expo_sign;
exponent_part = sub_int(exponent_part, dec_frac);
if (exponent_part == INT_MIN && errno == ERANGE) {
return FSM_EXPRANGE;
}
// NOTE2SELF: 10^x = (2*5)^x = 2^x * 5^x
if (exponent_part < 0) {
// use of absolute value to keep local_powl() simple
ld /= POWER(10.0L, abs(exponent_part));
} else if (exponent_part > 0) {
ld *= POWER(10.0L, abs(exponent_part));
}
} else {
// All numbers here are exact (powers of two) and the actual implementation
// with bigfloats will allow for easy manipulation of the exponent to avoid
// doing actual divisions here and below.
if (hex_frac > 0) {
fractional_part /= POWER(16.0L, hex_frac);
} else if (oct_frac > 0) {
fractional_part /= POWER(8.0L, oct_frac);
} else if (bin_frac > 0) {
fractional_part /= POWER(2.0L, bin_frac);
}
#ifdef DEBUG
printf("fracp: %.20Lg\n", fractional_part);
#endif
ld += fractional_part;
if (has_exponent != NO_EXP) {
if (expo_sign >= 0) {
ld *= POWER(2.0L, exponent_part);
} else {
ld /= POWER(2.0L, exponent_part);
}
}
}
ld *= main_sign;
*d = (double) ld;
// reset global variables
main_sign = 1;
expo_sign = 1;
integral_part = 0.0L;
fractional_part = 0.0L;
hex_frac = 0;
dec_frac = 0;
oct_frac = 0;
bin_frac = 0;
dec_int = 0;
exponent_part = 0;
has_exponent = NO_EXP;
// NOTE2SELF:
// mpf_normalize(&d,c->radix)
return FSM_OK;
}
// for printing the bits of the double
#include <stdint.h>
#ifdef BIG_ENDIAN
typedef union
{
double v;
struct
{
uint32_t m;
uint32_t l;
} p;
} dguts;
// Will fail, of course, if not little endian
#else
typedef union
{
double v;
struct
{
uint32_t l;
uint32_t m;
} p;
} dguts;
#endif
#define DOUBLETOINTS(high, low, d) \
do { \
dguts dg; \
dg.v = (d); \
(high) = dg.p.m; \
(low) = dg.p.l; \
} while (0)
int main(int argc, char **argv)
{
int res;
uint32_t high, low;
char *input, *endptr;
char *fsm_error2str[25] = {
"error",
"ok",
"exponent out of range"
};
double out, libc;
if (argc != 2) {
fprintf(stderr, "Usage: %s float\n", argv[0]);
exit(EXIT_FAILURE);
}
input = malloc(strlen(argv[1]) + 1);
if (input == NULL) {
fprintf(stderr, "Malloc failed to allocate %zu bytes\n",
strlen(argv[1]) + 1);
exit(EXIT_FAILURE);
}
strcpy(input, argv[1]);
// no actual error-check, just printing
errno = 0;
libc = strtod(input, &endptr);
printf("\nerrno after strtod(): %s\n", strerror(errno));
printf("\nINPUT:\n\tstring: \"%s\"\n\tstrtod: %g\tendptr: \"%s\"\n\n", input,
libc, endptr);
out = 0.0;
res = str2dbl(input, &out);
printf("fsm returned %d = \"%s\" and the result (if any) is:", res,
fsm_error2str[res]);
printf(" \n\tinp: %s\n\town: %.20g\n\tlibc: %.20g\n", trim_both(input), out,
libc);
printf("fsm result and libc result are%s equal to the bits\n",
out == libc ? "" : " NOT");
// prints in big-endian form. If you want to check with wolframalpha.com you
// need to click on the "big-endian" button in the upper right corner of the
// "IEEE double-precision number" box
DOUBLETOINTS(high, low, out);
printf("fsm-bits = %08x%08x\n",high, low);
DOUBLETOINTS(high, low, libc);
printf("libc-bits = %08x%08x\n",high, low);
free(input);
exit(EXIT_SUCCESS);
}
long double
related stuff) which has support for a lot from C99 but is basically somewhere in between c90 and c99. The C-compiler in MS-VS is more or less just an add-on, MS-VS is a C++ IDE, but it would be nice if it compiles my stuff, too. \$\endgroup\$flex
to write it for me, based on a formal description similar to what you've given. Did you consider that? If so, what is your objective in rejecting such an approach? \$\endgroup\$