Type-generic math for everyone.

May 01, 2003
URL:http://www.drdobbs.com/floating-point-math-in-the-intersection/184401653

In our March column, we presented the new complex arithmetic feature of C99 and introduced a small include file that we wrote, math-complex.h, that permits common source using complex arithmetic to be written in the intersection of C99 and C++. Perhaps the most amazing thing about math-complex.h is how only a few lines of code in that file can resolve the differences between C99's built-in complex data type and C++'s class library-based complex.

However, the usefulness of math-complex.h (see Listing 1) extends beyond complex arithmetic. As you might guess from the header's name, C programmers can think of math-complex.h as a shorthand for including both of the two standard headers <math.h> and <complex.h>, while C++ programmers can view it as a shorthand for including both <cmath> and <complex>. In addition, math-complex.h also resolves most of the differences in floating-point math between C and C++ and thus permits both C99 and C++ programmers to write common source in the intersection of the two languages. The March column discussed this intersection with respect to complex arithmetic. This column discusses the intersection with respect to real number arithmetic.

<tgmath.h>

Before discussing the features of the intersection, it is useful to look at a new C99 feature that is used by math-complex.h in order to provide to C programmers one of the useful features in the C++ floating-point library.

When the C++ committee created the <cmath> header from the C language header <math.h>, it added float and long double overloads for the math functions operating on doubles. This allowed the math functions to return the same type of result as their argument. For example, if you call sqrt() with a float, double, or long double argument, a float, double, or long double result respectively will be calculated and returned.

Having only one name for a library function independent of the argument types is a useful feature that has been used by FORTRAN programmers for over two decades. It has two primary advantages:

1. Convenience: programmers only need remember one name for any math function.

2. Flexibility: occasionally, numerical programmers must change the type of their floating-point data. For example, when porting to a new system, float variables might need to become double since float on the new system does not have the required precision. Also, sometimes the easiest way to discover whether float has enough precision for a calculation is to convert the program to double to see if it produces the same results. If the library functions have the same name regardless of their argument types, then there are fewer places in the program that need to be edited if a data type is changed.

C99 did not add a general-purpose overloading feature that allows programmers to write their own overloaded functions. It did however add a new header file, <tgmath.h>, that defines macros that act like the overloaded math functions in C++. These "type-generic" macros inspect the types of their arguments and call the appropriate library function to handle arguments of those types. For example, the type generic macro sqrt() when passed an argument of type float, double, long double, float complex, double complex, or long double complex calls the functions sqrtf, sqrt, sqrtl, csqrtf, csqrt, and csqrtl, respectively. While C99 requires implementations to implement the macros, it does not tell the implementations how to do so. Most implementations will probably invent special-purpose, private extensions to support <tgmath.h>.

When compiled by C99, math-complex.h includes <tgmath.h> so that both C and C++ programmers can use overloaded functions when using the math library. (When discussing the intersection of C99 and C++, we will use the terms "function" and "overloaded function" uniformly rather than dwell on the C99 implementation detail that the "functions" are really macros that achieve their ends without a general-purpose overloading mechanism in the language.)

What Features Are Provided in the Intersection?

Both C and C++ provide the usual arithmetic operators: plus, minus, times, divide, unary-plus, and unary-minus, on the float, double, and long double types. Both provide the corresponding assignment-operators (+=, -=, *=, and /=), as well as ordinary assignment (=). Equality (==), inequality (!=), less-than (<), less-or-equal (<=), greater-than (>), and greater-or-equal (>=) are provided in each language.

Since the three floating types are built-in (in C90, C99, and C++), all these operators accept one operand of one floating type mixed with an operand of any floating or integral type.

Both C and C++ provide the floating-point trigonometric and hyperbolic functions (i.e., sin, cos, tan, sinh, cosh, and tanh) and the corresponding inverse trig functions (asin, acos, atan, and atan2), but C++ does not (yet) provide the inverse hyperbolic functions (asinh, acosh, and atanh). The most common math functions are provided by C90, C99, and C++ (i.e., abs, ceil, exp, exp2, fabs, floor, fmod, frexp, ldexp, log, log10, pow, and sqrt). The type of the result is the same floating-point type as the operand.

What Does C++ provide, above the Intersection Features?

We turn now from our discussion of the C/C++ intersection to consider what is left out when you use the intersection. There is one C++ feature not provided in C99, namely the same-name overloads of modf:

float modf(float, float *);
double modf(double, double *);
long double modf(long double,
  long double *);

inline float modff
  (float xf, float *pxf) {
    return modf(xf, pxf);
  }
inline long double modfl
(long double xl, long double *pxl) {
    return modf(xl, pxl);
  }

What Does C99 Provide, above the Intersection Features?

There are many additional functions in the C99 library, above the library C90 and C++. Here we will list all the double versions; each returns double, unless indicated otherwise; x and y are double, n is int, and ln is long int. The C99 "type-generic math" header also provides float and long double versions with the same name.

• expm1(x) returns e^x- 1 (where e^x means "e to the power x").
• ilogb(x) returns (as int) the exponent of x as a signed int value.
• log1p(x) returns the base-e logarithm of (x+1).
• log2(x) returns the base-2 logarithm of x.
• logb(x) returns the exponent of x as a signed floating-point value.
• scalbn(x, n) and scalbln(x, ln) return x times FLT_RADIX^n.
• cbrt(x) returns the real cube root of x.
• hypot(x, y) returns the square root of the sum of the squares of x and y.
• erf(x) returns the error function of x.
• erfc(x) returns the complementary error function of x (i.e., one minus erf(x)).
• lgamma(x) returns the natural logarithm of the absolute value of gamma of x.
• tgamma(x) returns gamma of x.
• nearbyint(x) returns the rounded integer value of x, in the current rounding mode, without raising the "inexact" floating-point exception.
• rint(x) returns the rounded integer value of x, in the current rounding mode, possibly raising the "inexact" floating-point exception.
• lrint(x) returns (as long int) the rounded integer value of x, in the current rounding mode. (A long long int version is available in C99, but not in Standard C++.)
• round(x) returns the rounded integer value of x, rounding halfway cases away from zero regardless of the current rounding mode.
• lround(x) returns (as long int) the rounded integer value of x, rounding halfway cases away from zero regardless of the current rounding mode. (A long long int version is available in C99, but not in Standard C++.)
• trunc(x) returns the value of x truncated (towards zero) to the nearest integer.
• remainder(x, y) returns the remainder x-n*y, where n is the integer nearest the exact value of x/y (if |x-n*y| equals 1/2, then n is even) as required by IEC 60559.
• remquo(x, y, int *quo) returns remainder(x, y) and stores into *quo a value whose sign is the sign of x/y and whose magnitude is congruent modulo 2*n to the magnitude of the integral quotient of x/y, where n is an implementation-defined integer greater than or equal to 3.
• copysign(x, y) returns a value with the magnitude of x and the sign of y.
• nextafter(x, y) returns the next representable value after x in the direction of y.
• nexttoward(x, long double y) returns the next representable value after x in the direction of y.
• fdim(x, y) returns x-y if x>y, otherwise returns +0.
• fmax(x, y) returns the maximum of x and y.
• fmin(x, y) returns the minimum of x and y.

Cost/Benefit of Working in the Intersection

We have avoided discussion of the special features and behaviors specifically provided by C99 for the boundary cases of IEEE (or "IEC 559") floating-point arithmetic, including NaNs, infinities, signed zeroes, rounding modes, floating-point exceptions, and the functions and macros that manipulate these features. We also avoided discussion of fused-multiply-add and long long integers. If any of these are important to your work, then C99 is clearly your choice, and there would be little value to working in the intersection. Otherwise, we recommend the use of the "intersection" style for your arithmetic applications. As we mentioned in March, in some environments, the C99 implementation will produce faster code execution. On the other hand, there will be environments with no C99 implementation. And a C++ implementation with "small-object" optimization can be as fast as C99. In any event, the application that is written in the "intersection" style can be compiled as C++ or compiled as C99, whichever works best for each environment.

Example Program

Listing 1 shows the math-complex.h header. Listing 2 is a short example program that illustrates that the result type of an overloaded function call depends on the argument type used in the call. Listing 3 shows sample output from running the program.

The CUJ website (<www.cuj.com/code>) contains an expanded version of the test program plus an updated version of the complex arithmetic test from the March column.

About the Authors

Randy Meyers is a consultant providing training and mentoring in C, C++, and Java. He is the current chair of J11, the ANSI C committee, and previously was a member of J16 (ANSI C++) and the ISO Java Study Group. He worked on compilers for Digital Equipment Corporation for 16 years and was Project Architect for DEC C and C++. He can be reached at [email protected].

Dr. Thomas Plum has authored four books on C, and co-authored Efficient C (with Jim Brodie) and C++ Programming Guidelines (with Daniel Saks). He has been an officer of the United States and international C and C++ standards committees. His company Plum Hall Inc. provides test suites for C, C++, Java, and C#. His address is [email protected].

Listing 1: math-complex.h

// Compatibility file for C99 and C++ complex from May 2003
// C/C++ Users Journal.
//
// This header can be included by either C99 or ANSI C++
// programs to allow both real and complex floating point
// arithmetic to be written in a common subset of C99 & C++.
// Note that overloads for both the real and complex math
// functions are available after this header has been
// included.

#ifndef MATH_COMPLEX_H_INCLUDED
#define MATH_COMPLEX_H_INCLUDED

#ifdef __cplusplus

#include <cmath>
#include <complex>

using namespace std;

typedef complex<float> float_complex;
typedef complex<double> double_complex;
typedef complex<long double> long_double_complex;

// Since C99 does not have type-generics for modf of float
// and long double and C++ does not have the suffixed names,
// provide suffixed names for C++.
// Note: Some C++ implementations do not yet support the
// overloaded modf function called by these inlines.
inline float modff(float xf, float *pxf)
    {return modf(xf, pxf);}
inline long double modfl(long double xl, long double *pxl)
    {return modf(xl, pxl);}

#else

// Note that <tgmath.h> includes <math.h> and <complex.h>
#include <tgmath.h>

typedef float complex float_complex;
typedef double complex double_complex;
typedef long double complex long_double_complex;

#define float_complex(r,i) ((float)(r) + ((float)(i))*I)
#define double_complex(r,i) ((double)(r) + ((double)(i))*I)
#define long_double_complex(r,i) ((long double)(r) + ((long double)(i))*I)

#define real(x) creal(x)
#define imag(x) cimag(x)
#define abs(x) fabs(x)
#define arg(x) carg(x)

#endif  // #ifdef __cplusplus

#endif  // #ifndef MATH_COMPLEX_H_INCLUDED

Listing 2: Test program

// Test program for math-complex.h. This file can be

// compiled by either C99 or C++.

// An extended version of this test is available on

// the CUJ website at www.cuj.com/code.

#include <stdio.h>

#include "math-complex.h"

static const long double pi =

3.141592653589793238462643383279502884197L;

int main()

{

// Establish the size of floating point type on

// this platform. Note on some platforms, two or

// more of these datatypes may have the same size.

printf("sizeof(float) = %d\n", (int) sizeof(float));

printf("sizeof(double) = %d\n", (int) sizeof(double));

printf("sizeof(long double) = %d\n", (int) sizeof(long double));

// The overloaded library function sin() should return

// the same type as its argument. So, take the sizeof

// function calls passing a float, double, and long

// double to see if the function result has sizeof

// a float, double, and long double, respectively.

printf("sizeof sin(0.0F) = %d\n", (int) sizeof sin(0.0F));

printf("sizeof sin(0.0) = %d\n", (int) sizeof sin(0.0));

printf("sizeof sin(0.0L) = %d\n", (int) sizeof sin(0.0L));

// arc sine of 1 is pi/2. Try this for the three real

// floating point types and see if the result is more

// precise as the float point type gets more precise.

printf("pi/2 = %.20Lf\n", pi/2);

printf("asin(1.0F) = %.20f error=%Lg\n",

asin(1.0F), asin(1.0F) - pi/2);

printf("asin(1.0) = %.20f error=%Lg\n",

asin(1.0), asin(1.0) - pi/2);

printf("asin(1.0L) = %.20Lf error=%Lg\n",

asin(1.0L), asin(1.0L) - pi/2);

return 0;

}

-- End of Listing --

Listing 3: Sample output

sizeof(float) = 4

sizeof(double) = 8

sizeof(long double) = 12

sizeof sin(0.0F) = 4

sizeof sin(0.0) = 8

sizeof sin(0.0L) = 12

pi/2 = 1.57079632679489661928

asin(1.0F) = 1.57079637050628662109 error=4.37114e-08

asin(1.0) = 1.57079632679489655799 error=-6.12574e-17

asin(1.0L) = 1.57079632679489661928 error=0

-- End of Listing --

Terms of Service | Privacy Statement | Copyright © 2024 UBM Tech, All rights reserved.