C99 and Numeric Computing

Application Examples

Assume that the second order polynomial equation in Example 4(a) can be solved by the formula in Example 4(b).

Example 4: Using metanumbers and complex numbers.

According to this formula, two solutions of x1=2 and x2=3 for Example 4(c) can be obtained by Listing One.

Listing One

#include <stdio.h>
#include <math.h>
int main() {
    double a = 1, b = -5, c = 6, x1, x2;
    x1 = (-b +sqrt(b*b-4*a*c))/(2*a);
    x2 = (-b -sqrt(b*b-4*a*c))/(2*a);
    printf("x1 = %f\n", x1);
    printf("x2 = %f\n", x2);
}

Executing Listing One yields the output:

x1=3.000000

x2=2.000000

For the equation x²-4x+13=0, two complex solutions of x₁=2+i3 and x₂=2-i3 cannot be solved in the real domain. These complex numbers cannot be represented in double data type.

Listing Two

#include <stdio.h>
#include <math.h>
int main() {
    double a = 1, b = -4, c = 13, x1, x2;
    x1 = (-b +sqrt(b*b-4*a*c))/(2*a);
    x2 = (-b -sqrt(b*b-4*a*c))/(2*a);
    printf("x1 = %f\n", x1);
    printf("x2 = %f\n", x2);
}

In the domain of real numbers, the values of these solutions are Not-a-Number. Therefore, executing Listing Two gives this output:

x1=NaN

x2=NaN

However, using complex numbers, x²-4x+13=0 can be solved by the C99-conforming program in Listing Three.

Listing Three

#include <stdio.h>
#include <math.h>
#include <complex.h>
int main() {
    double complex a = 1, b = -4, c = 13, x1, x2;
    x1 = (-b +csqrt(b*b-4*a*c))/(2*a);
    x2 = (-b -csqrt(b*b-4*a*c))/(2*a);
    printf("x1 = complex(%f,%f)\n", creal(x1), cimag(x1));
    printf("x2 = complex(%f,%f)\n", creal(x2), cimag(x2));
}

Executing Listing Three gives this output:

x1=complex(2.000000,3.000000)

x2=complex(2.000000,-3.000000)

The same problem can be solved by the C++-conforming program in Listing Four with the same output as Listing Three. Note that complex type is declared with the type specifier double_complex in C++ instead of double complex in C99. Real and imaginary parts of a complex number are obtained by functions real() and imag() in C++, instead of creal() and cimag() in C99.

Listing Four

#include <stdio.h>
#include <math.h>
#include <complex.h>
int main() {
    double_complex a = 1, b = -4, c = 13, x1, x2;
    x1 = (-b +sqrt(b*b-4*a*c))/(2*a);
    x2 = (-b -sqrt(b*b-4*a*c))/(2*a);
    printf("x1 = complex(%f,%f)\n", real(x1), imag(x1));
    printf("x2 = complex(%f,%f)\n", real(x2), imag(x2));
}

Listing Five

#include <stdio.h>
#include <math.h>
#ifndef INFINITY
#define INFINITY (1.0/0.0)
#endif
int main() {
    double a = 1, b = INFINITY, c = 13, x1, x2;
    x1 = (-b +sqrt(b*b-4*a*c))/(2*a);
    x2 = (-b -sqrt(b*b-4*a*c))/(2*a);
    printf("x1 = %f\n", x1);
    printf("x2 = %f\n", x2);
}

If the coefficient b in ax²+bx+c=0 is , which can be represented as INFINITY in C99, the solution can be calculated by Listing Five with this output:

x1=NaN

x2=-Inf

Listing Six

#include <stdio.h>
#include <math.h>
#include <complex.h>
#ifndef INFINITY
#define INFINITY (1.0/0.0)
#endif
int main() {
    double complex a = 1, b = INFINITY, c = 13, x1, x2;
    x1 = (-b +csqrt(b*b-4*a*c))/(2*a);
    x2 = (-b -csqrt(b*b-4*a*c))/(2*a);
    printf("x1 = complex(%f,%f)\n", creal(x1), cimag(x1));
    printf("x2 = complex(%f,%f)\n", creal(x2), cimag(x2));
}

Note that the C99 macro INFINITY for is not supported in gcc; it is obtained conveniently in the program by dividing 1.0 by 0. If complex numbers are used in the calculation by Listing Six, which is C99 conforming, the output becomes:

x1=complex(NaN,NaN)

x2=complex(NaN,NaN)

Conclusion

C99 makes the power of IEEE 754 floating-point arithmetic easily available to you for applications in the real domain. Still, a few special values of mathematical functions related to Not-a-Number — improperly defined in C99 — should be redefined. For example, the results of hypot(Inf,NaN)==Inf, hypot(-Inf,NaN)==Inf, pow(1,NaN)==1.0, pow(NaN,0.0)==1.0, pow(NaN,-0.0)==1.0 in C99 should be redefined as NaN.

Complex numbers are a generalization of real numbers. C99 supports complex numbers as built-in data types, whereas C++ uses classes. There are differences for complex numbers in C and C++. But a typical C/C++ conforming program using complex numbers can run in Ch without any modification. Still, I recommend that the familiar complex-construction function polar() be added to resolve the compatibility of C/C++.

Ch conforms to the C99 standard related to the IEEE 754 floating-point arithmetic and complex numbers. Following the conventional mathematics, Ch handles numerical numbers consistently in the entire real and complex domains. Ch treats floating-point real numbers with signed zeros, signed infinities (Inf and -Inf), and a Not-a-Number (NaN); and complex numbers with unsigned zeros, a unique complex infinity (ComplexInf), and a unique complex Not-a-Number (ComplexNaN) in a consistent manner. In addition, generic mathematical functions in Ch have optional arguments to handle different branches of complex functions. Commonly used high-level functions with complex numbers for numerical analysis have been implemented in Ch. Furthermore, Ch has computational arrays for numerical computing of linear algebra and matrices. For example, linear equation b=A*x can be written verbatim in Ch. Users do not need to worry about the underlying optimization with fast and accurate numerical algorithms. Numerical computing with high-level numerical analysis functions in the entire real and complex domains under the framework of C99 and IEEE 754 Standards are very useful for solving problems in engineering and science.

Harry is an associate professor and director of the Integration Engineering Lab at the University of California Davis. He can be contacted at [email protected].