Proper design of C code that handles both single- and double-precision floating point?

Question

I am developing a library of special-purpose math functions in C. I need to provide a capability for the library to handle both single-precision and double-precision. The important point here is that the "single" functions should use ONLY "single" arithmetic internally (resp. for the "double" functions).

As an illustration, take a look at LAPACK (Fortran), which provides two versions of each of its function (SINGLE and DOUBLE). Also the C math library (example, expf and exp).

To clarify, I want to support something similar to the following (contrived) example:

float MyFloatFunc(float x) {
    return expf(-2.0f * x)*logf(2.75f*x);
}

double MyDoubleFunc(double x) {
    return exp(-2.0 * x)*log(2.75*x);
}

I've thought about the following approaches:

1. Using macros for the function name. This still requires two separate source codebases:

   #ifdef USE_FLOAT
   #define MYFUNC MyFloatFunc
   #else
   #define MYFUNC MyDoubleFunc
   #endif
   

2. Using macros for the floating point types. This allows me to share the codebase across the two different versions:

   #ifdef USE_FLOAT
   #define NUMBER float
   #else
   #define NUMBER double
   #endif
   

3. Just developing two separate libraries, and forgetting about trying to save headaches.

Does anyone have a recommendation or additional suggestions?

Answer 1

The <tgmath.h> header, standardized in C 1999, provides type-generic calls to the routines in <math.h> and <complex.h>. After you include <tgmath.h>;, the source text sin(x) will call sinl if x is long double, sin if x is double, and sinf if x is float.

You will still need to conditionalize your constants, so that you use 3.1 or 3.1f as appropriate. There are a variety of syntactic techniques for this, depending on your needs and what appears more aesthetic to you. For constants that are exactly represented in float precision, you can simply use the float form. E.g., y = .5f * x will automatically convert .5f to double if x is double. However, sin(.5f) will produce sinf(.5f), which is less accurate than sin(.5).

You might be able to reduce the conditionalization to a single clear definition:

#if defined USE_FLOAT
    typedef float Float;
#else
    typedef double Float;
#endif

Then you can use constants in ways like this:

const Float pi = 3.14159265358979323846233;
Float y = sin(pi*x);
Float z = (Float) 2.71828182844 * x;

That might not be completely satisfactory because there are rare cases where a numeral converted to double and then to float is less accurate than a numeral converted directly to float. So you may be better off with a macro described above, where C(numeral) appends a suffix to the numeral if necessary.

Answer 2

(Partially inspired by Pascal Cuoq's answer) If you want one library with float and double versions of everything, you could use recursive #includes in combination with macros. It doesn't result in the clearest of code, but it does let you use the same code for both versions, and the obfuscation is thin enough it's probably manageable:

mylib.h:

#ifndef MYLIB_H_GUARD
  #ifdef MYLIB_H_PASS2
    #define MYLIB_H_GUARD 1
    #undef C
    #undef FLT
    #define C(X) X
    #define FLT double
  #else
    /* any #include's needed in the header go here */

    #undef C
    #undef FLT
    #define C(X) X##f
    #define FLT float
  #endif

  /* All the dual-version stuff goes here */
  FLT C(MyFunc)(FLT x);

  #ifndef MYLIB_H_PASS2
    /* prepare 2nd pass (for 'double' version) */
    #define MYLIB_H_PASS2 1
    #include "mylib.h"
  #endif
#endif /* guard */

mylib.c:

#ifdef MYLIB_C_PASS2
  #undef C
  #undef FLT
  #define C(X) X
  #define FLT double
#else
  #include "mylib.h"
  /* other #include's */

  #undef C
  #undef FLT
  #define C(X) X##f
  #define FLT float
#endif

/* All the dual-version stuff goes here */
FLT C(MyFunc)(FLT x)
{
  return C(exp)(C(-2.0) * x) * C(log)(C(2.75) * x);
}

#ifndef MYLIB_C_PASS2
  /* prepare 2nd pass (for 'double' version) */
  #define MYLIB_C_PASS2 1
  #include "mylib.c"
#endif

Each file #includes itself one additional time, using different macro definitions on the second pass, to generate two versions of the code that uses the macros.

Some people may object to this approach, though.

Answer 3

The big question for you will be:

• Is it easier to maintain two separate unobfuscated source trees, or one obfuscated one?

If you have the proposed common coding, you will have to write the code in a stilted fashion, being very careful not to write any undecorated constants or non-macro function calls (or function bodies).

If you have separate source code trees, the code will be simpler to maintain in that each tree will look like normal (non-obfuscated) C code, but if there is a bug in YourFunctionA in the 'float' version, will you always remember to make the matching change in the 'double' version.

I think this depends on the complexity and volatility of the functions. My suspicion is that once written and debugged the first time, there will seldom be a need to go back to it. This actually means it doesn't matter much which mechanism you use - both will be workable. If the function bodies are somewhat volatile, or the list of functions is volatile, then the single code base may be easier overall. If everything is very stable, the clarity of the two separate code bases may make that preferable. But it is very subjective.

I'd probably go with a single code base and wall-to-wall macros. But I'm not certain that's best, and the other way has its advantages too.

Answer 4

For polynomial approximations, interpolations, and other inherently approximative math functions, you cannot share code between a double-precision and a single-precision implementation without either wasting time in the single-precision version or being more approximative than necessary in the double-precision one.

Nevertheless, if you go the route of the single codebase, the following should work for constants and standard library functions:

#ifdef USE_FLOAT
#define C(x) x##f
#else
#define C(x) x
#endif

... C(2.0) ... C(sin) ...