Get a floating point number from a string with a finite state machine

Question

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);
}

Answer 1

There is a lot in your post. So far only some minor ideas.

1. Fundamentally, code is relying on long double with a greater precision than double. Better to code a solution without this reliance. Else how to code a long double version?

2. str2dbl(char *s, double *d) only uses char *s as a non-const to do trimming. I'd expect code to cope with a const string. The code modifications needed are small.

3. Re-order str2dbl() for clarity by starting with space and sign handling. Suggested layout:

   int str2dbl(char *s, double *d) {
     ...
     while (isspace(*s)) s++;
     int sign = *s;
     if (sign == '+' || sign == '-') s++;
     // At this point sign is either '-' or not.
   
     // Handle NAN
   
     // Handle Infinity
   
     // main converison
   
     // detect is end-of-string is \0
   
     if (sign == '-') result = -result;
   

4. Small simplification with add_int(int a, int b). Using a b >= 0 instead of b > 0 allows a simplification that a compiler may not catch. Similar with sub_int()

   // if (((b > 0) && (a > (INT_MAX - b))) || ((b < 0) && (a < (INT_MIN - b)))) {
   //      v-------------- opposite ------------v
   if (((b >= 0) && (a > (INT_MAX - b))) || ((b < 0) && (a < (INT_MIN - b)))) {
   // or
   if (b >= 0) { if (a > INT_MAX - b) Over(); } 
   else { if (a < INT_MIN - b) Under(); }. 
   

5. Wrong return value on overflow (2 places).

   if (b < (INT_MIN / a)) {
     errno = ERANGE;
     // return INT_MAX;
     return INT_MIN;
   

6. See little value is type matching simple constants.

   // long double power = 1.0L;
   long double power = 1.0;
   

7. Questionable code when exponent < 0 (if that occurs in this code.) Odd mixture of using % and >>. I'd expect either (% and /) or (& and >>). With signed arithmetic, these have slightly different functionalities.

   while (exponent) {
     // if (negative_exponent % 2 --> -1) 
     // if (exponent % 2 == 1) {  
     if (exponent % 2) {
       power *= base;
     }
     // exponent >>= 1;
     exponent /= 2;
     base *= base;
   }
   

8. Casting to unsigned char is better for string is....() and to...() functions. These functions are UB for negative values. Of course, ASCII is always positive, but not much work to prevent.

     const char *s1,
     ...
     // c1 = tolower(*s1);
     c1 = tolower((unsigned char) *s1);
   

9. Questionable code when char is unsigned. Suggest signed char. Else main_sig will the wrong value even is one assumes ASCII-only in the string.

   // static const char digit_map[] = {
   static const signed char digit_map[] = {
     -1, -1, -1, -1, -1, -1, -1, -1, //  0x00-0x07
   

10. strncasecmp() has the same problem. Better to use unsigned char here.

    // char c1 = 0;
    // char c2 = 0;
    unsigned char c1 = 0;
    unsigned char c2 = 0;
    

11. For locale issues, the decimal point may be something else than '.'.

12. Avoid magic numbers.

    // static int fsm_table[31][13] = {
    static int fsm_table[][OTHER] = {
    

13. Concerning // empty input, would an error be better?. Any parsing that lacks a scanned digit (or inf or NaN) should flag an error. E. g.: "", "-", "-+0", "in", " " etc.

Answer 2

I think there are several issues around your code's dependency on long double and its characteristics, and around uses of that type. Inasmuch as you have made it clear that the whole long double bit is a placeholder for an alternative internal representation, however, I'll not devote time or space to those issues. I mention it only to clarify the scope of this review.

In addition to @chux's astute observations, then ...

Character class assignment

You perform character class assignment via a big switch statement in function check_type(). I suggest instead implementing a (bounds-checked) table lookup, which you clearly know how to do. All of the characters you assign to classes other than "other" belong to C's basic character set, so they are guaranteed to have positive char values even on systems that have signed default chars, in case that was a concern. You could even eliminate your calls to tolower() that way.

Also, it might be convenient in this and other respects to declare enum fsm_input so that OTHER corresponds to the value 0.

Globals

As @JS1 remarked first (and since deleted), reliance on file-scope variables makes your code non-thread-safe. There's no particular need for that -- yes, if different functions need to access shared parse state then passing around a pointer to a local struct that contains the needed state variables is a reasonable solution. As a bonus, the compiler may be able to optimize the thread-safe code better.

Half-hearted use of enums

You go to the trouble to declare enum types for input classes and machine states, but you then use only the declared enum constants, not the types themselves (you declare objects expected to hold instances of those types as type int). It would be clearer and would afford compilers more opportunity for type checking and optimization to declare objects as the appropriate one of those enum types wherever that is in fact what they are expected to represent.

Numbers with trailing garbage

Your specification says that the input must start with a floating-point number, and parsing must fail if it does not contain a floating point number at the beginning. It appears, however, that the program will report an error for strings that start with a parseable number, followed by a character from class "OTHER", even though it seems from your description that such inputs should be accepted.

Similarly, but separately, your code rejects any input that causes the machine to transition to the ERROR state, even if the input parsed before that state is reached constitutes a valid representation of a number. This, too, seems at odds with the specification.

Signature of str2dbl()

If you really mean it about parsing a number from the beginning of the input, as opposed to parsing one from the whole input, then you should consider adding a mechanism for telling the caller where the parsed initial substring stops. You might emulate strtod() in this regard, for example.

Modifying the the input

You modify the input string via function trim_both(). This is inherently undesirable, if only because it means your function cannot work on const strings, as @chux already observed. Even if const strings were not a concern, however, modifying the input still needlessly narrows the use cases to which your function can be applied. Returning an indication of the extent or end of the parsed portion, as already discussed, combined with the changes necessary to accept inputs with trailing garbage, would provide a reasonable alternative.

Return values of add_int() and sub_int()

When these functions detect that the sum (difference) would not be representable as an int, the former always returns INT_MAX, and the latter always INT_MIN, regardless of whether the result over- or underflows. That's a bit surprising given that both alternatives are possible for each function, and that you seem to take a different approach in mul_int(). You already signal the error by setting errno, so these return values convey no additional information in the error case. You certainly can and do work with the implementation as it is, but this is the kind of little quirk that breeds bugs later.