Enhanced Computational Support

I l @ ve RuBoard

Historically, FORTRAN has been the premier language for numerical scientific and engineering computation. C90 brought C computational methods into closer agreement with FORTRAN. For example, the specification of floating-point characteristics used in float.h is based on the model developed by the FORTRAN standardization committee. The C9X continues the work of enhancing C's appropriateness for computational work.

The `fenv.h` Header File

The fenv.h header file provides a means of interacting with the floating-point environment. That is, it allows you to set floating-point control mode values that govern how floating-point calculations take place, and it allows you to determine the value of floating-point status flags that report information about the effects of an arithmetic calculation. An example of a control mode setting is specifying the method used to round numbers . An example of a status flag is a flag that is set if an operation produces floating-point overflow.

The status flags and control modes are meaningful only if the hardware supports them. For instance, you can't change the rounding method if the hardware doesn't have that option.

By default, access to the floating-point environment is off. You use a preprocessor directive to turn support on:

 #pragma STDC FENV_ACCESS ON

Support stays on until the program reaches the end of the block containing the pragma, or, if the pragma is external, to the end of the file or translation unit. Or you can use the following directive to turn off support:

 #pragma STDC FENV_ACCESS OFF

This facility is important for those involved in critical floating-point calculations, but of limited interest to the general user , so this appendix doesn't go into the details.

Additions to the `math.h` Library

The C90 math library, for the most part, declares functions with type double arguments and type double return values, such as the following

 double sin(double); double sqrt(double);

The C9X library provides type float and type long float versions of all these functions. These functions use an f or an l suffix in the name , as follows :

 float sinf(float);              /* float version of sin()     */ long double sinl(long double);  /* long double version of sin()  */

Having function families with different levels of precision allows you to choose the most efficient combination of types and functions needed for a particular purpose.

C9X also adds several functions commonly used in scientific, engineering, and mathematical computations . Table H.6 lists the double versions for several of the new functions. All have type float and type long double counterparts, as described previously. In many cases, the functions return values that could be calculated using existing functions, but the new functions do so faster or more accurately. For instance, log1p(x) represents the same value as log(1 + x) , but log1p(x) uses a different algorithm, one that is more accurate for small values of x . So you would use the log() function for calculations in general, but you would use log1p() for small values of x if high accuracy were critical. The table refers to FLT_RADIX . This constant, defined in float.h , is the base used for exponentiation in the internal floating-point representation. The most common value is 2 .

Table 8.6. New math library functions.

Prototype	Description
`double exp2(double x);`	Returns `2` to the `x` power
`double expm1(double x);`	Returns `e ^x - 1`
`double logp1(double x);`	Returns `log(1 + x)`
`double log2(double x);`	Returns the base 2 logarithm of `x`
`double logb(double x);`	Returns the signed exponent of its argument for the underlying base used to represent floating-point values on the system ( `FLT_RADIX` )
`double scalbn (double x, int n);`	Returns `x` — `FLT_RADIX ⁿ`
`double scalbln (double x, long n);`	Returns `x` — `FLT_RADIX ⁿ`
`double ilogb (double x);`	Same as `logb()` except it returns the exponent as an `int`
`double cbrt(double x);`	Returns the cube root of `x`
`double hypot (double x, double y);`	Returns the square root of the sums of the squares of `x` and `y`
`double erf(double x);`	Returns the error function of `x`
`double erfc(double x);`	Returns the complementary error function of `x`
`double gamma(double x);`	Returns the gamma function of `x`
`double lgamma(double x);`	Returns the natural logarithm of the absolute value of the gamma function of `x`
`double copysign (double x, double y);`	Returns the value of `x` with the sign of `y`
`double fdim(double x, double y);`	Returns the absolute value of `x - y`

In addition to these functions, the math library defines several constants and functions relating to classifying numbers and rounding them. For example, a value can be classified as being infinite, not a number ( NaN ), normal, subnormal, and true zero. (A NaN is a special value indicating that a value is not a number; for example, asin(2.0) returns NaN because asin() is defined only for arguments in the range -1 to 1 . A subnormal number is one whose magnitude is smaller than the smallest value that can be represented to full precision.) There are also specialized comparison functions that behave differently from the standard relational operators when one or more arguments are abnormal values. In short, there is expanded support for detailed control of how floating-point calculations are handled.

Support for Complex Numbers

A complex number is a number with a real part and an imaginary part. The real part is an ordinary real number, such as what's represented by the floating-point types. The imaginary part, naturally, is an imaginary number. An imaginary number, in turn, is a multiple of the square root of -1. In mathematics, complex numbers are often written in the form 4.2 + 2.0i ; i symbolically represents the square root of -1.

C9X supports three complex types:

float complex
double complex
long double complex

A float complex value, for example, would be stored using the same memory layout as a two-element array of float , with the real value stored in the first element and the imaginary value in the second element.

C9X implementations may also support three imaginary types:

float imaginary
double imaginary
long double imaginary

You should include the complex.h header file before using the keywords complex and imaginary ; the effect of failing to do so is undefined.

Arithmetic operations are defined for complex types following the usual rules of mathematics. For example, the value of (a+b*I)*(c+d*I) is (a*c-b*d)+(b*c+a*d)*I .

The complex.h header file defines some macros and several functions that accept complex numbers and return complex numbers. In particular, the macro I represents the square root of -1. It enables you to do the following:

 double complex c1 = 4.2 + 2.0 * I; float imaginary c2= -3.0 * I;

The complex.h header file prototypes several complex functions. Many are complex equivalents of math.h functions, using a c prefix. For example, csin() returns the complex sine of its complex argument. Others relate specifically to the features of complex numbers. For example, creal() returns the real part of an complex number, and cimag() returns the imaginary part as a real number. That is, given that z is type double complex , then the following is true:

 z = creal(z) + cimag(z) * I;

If you are familiar with complex numbers and need to use them, you'll want to peruse the contents of complex.h . If you use C++, you should be aware that the C++ complex header file provides a different way, based on classes, of handling complex numbers than does the C complex.h header file.

The `restrict` Keyword

The restrict keyword enhances computational support by giving the compiler permission to optimize certain kinds of code. It can be applied only to pointers, and it indicates that a pointer is the sole initial means of accessing a data object. Consider the following:

 int ar[10]; restrict int * restar = (int *) malloc(10 * sizeof(int)); int * par = ar;

Here, the pointer restar is the sole initial means of access to the memory allocated by malloc() . Therefore, it can be qualified with the keyword restrict . The pointer par , however, is neither the initial nor the sole means of access to the data in the ar array, so it cannot be qualified as restrict .

Now consider the following rather artificial example, in which n is an int :

 for (n = 0; n ;lt 10; n++) {     par[n] += 5;     restar[n] += 5;     ar[n] *= 2;     par[n] += 3;     restar[n] += 3; }

Knowing that restar is the sole initial means of access to the block of data it points to, the compiler can replace the two statements involving restar with a single statement having the same effect:

 restar[n] += 8;   /* ok replacement */

It would be a computational error, however, to condense the two statements involving par to one:

 par[n] += 8;    / * gives wrong answer */

The reason it gives the wrong answer is that the loop uses ar to change the value of the data between the two times par accesses the same data.

Without the restrict keyword, the compiler has to assume the worse case; namely, that some other identifier might have changed the data in between two uses of a pointer. With the restrict keyword used, the compiler is free to look for calculational shortcuts.

You can use the restrict keyword as a qualifier for function parameters that are pointers. This means the compiler can assume that no other identifiers modify the pointed-to data within the body of the function and that the compiler can try optimizations it might not otherwise use. For example, the C library has two functions for copying bytes from one location to another. Under C9X, they have these prototypes:

 void * memcpy(void * restrict s1, const restrict s2, size_t n); void * memmove(void * s1, const s2, size_t n);

Each copies n bytes from location s2 to location s1 . The memcpy() function requires that there be no overlap between the two locations, but memmove() doesn't have that requirement. Declaring s1 and s2 as restrict means each pointer is a sole means of access, so they can't access the same block of data. This matches the requirement that there be no overlap. The memmove() function, which allows overlap, has to be more careful about copying data so that it doesn't overwrite data before it is used.

I l @ ve RuBoard

Enhanced Computational Support