This discontinuity can be handled gracefully in C99. The evaluation of the expression exp(1/0.0) returns Inf and exp(1/(-0.0)) gives 0.0, which corresponds to mathematical expressions e1/0+ and e1/0- or limx0+e1/x and lim x0-e1/x, respectively. In addition, the evaluation of expressions exp(1.0/Inf) and exp(1.0/(-Inf)) get the value of 1.0. Likewise, the function finite(x) recommended by IEEE 754 is equivalent to the expression -Inf<x&&x<Inf, where x can be a float/double variable or expression. If x is a float, -Inf<x&&x<Inf is equivalent to -FLT_MAXx&&xFLT_MAX; if x is a double, -Inf<x&&x<Inf is equivalent to -DBL_MAXx&&xDBL_MAX;. The mathematical statement "if -<value , then y becomes " can be programmed like this:

if(-Inf < value && value <= Inf) y = Inf;

However, computers can only evaluate expressions step by step. Although the metanumbers are limits of the floating-point numbers, they cannot replace mathematical analysis. For example, the natural number e equal to 2.718281828... is defined as the limit value of the expression in Example 1(a). But the value of the expression pow(1.0 + 1.0/Inf, Inf) is NaN. The evaluation of this expression is carried out like Example 1(b). If the value FLT_MAX is used here instead of Inf, the result is obtained by Example 1(c). Because the metanumber NaN is unordered, a program involving relational operations should be handled cautiously. For example, the expression x>y is not equivalent to !(x<= y) if either x or y is an NaN. Likewise, Example 2(a) is different from the code in Example 2(b). The second if statement should be written as if(x0.0||isnan(x)) to have the same functionality for these two code fragments.

Example 1: Evaluating expressions step by step.

(a)
if(x > 0.0)  function1();
else function2();

(b)
if(x <= 0.0) function2();
else function1();

Example 2: Handling relational operations.

The metanumbers 0.0, -0.0, Inf, -Inf, and NaN are useful for applications in engineering. For example, the discontinuity at the origin can be expressed using signed zeros. The infinity of mechanical advantage at a toggling position for a four-bar linkage can be written as Inf. If no solution exists for output link corresponding to a given input link position of a four-bar linkage, the solution can be represented symbolically as NaN.

Computing in the Entire Complex Domain

Complex numbers zC={(x,y)|x,yR} can be defined as ordered pairs z=(x,y) with specific addition and multiplication rules. The real numbers x and y are called the "real" and "imaginary" parts of z. If you identify the pair of (x,0.0) as the real numbers, the real number R is a subset of C; that is, R={(x,y)|xR, y=0.0} and RC. If a real number is considered either as x or (x,0.0) and let i denote the pure imaginary number (0,1) with i*i=-1, complex numbers can be mathematically represented as z=x+iy.

Complex numbers have wide applications in engineering and science. Due to their importance in scientific programming, numerically oriented languages and packages usually provide complex number support in one way or another. Fortran, for example, has long provided complex data types. Computations involving complex numbers can be accomplished using a data structure in C90. However, programming with such a structure is clumsy because a corresponding function has to be invoked for each operation. With C++, classes for complex numbers can be created. By using operator and function overloadings, built-in operators and functions can be extended to also apply to the complex data type. Therefore, the complex data type can be achieved. Such complex classes have been added in the C++ Standard. Without overloading features of classes, C99 supports complex numbers with built-in data types. There are differences for complex numbers in C/C++. However, a typical C or C++ program with complex numbers can run in Ch without any modification. A complex number z=(x,y)=rei can be constructed by x+I*y, complex(x,y), polar(r,theta), float_complex(x,y), or double_complex(x,y) in Ch. Complex numbers are supported in generic polymorphic mathematical functions in Ch. Textbooks and software often avoid issues related to complex infinity and, other than Ch, no other general-purpose programming language and implementation has provisions for consistent handling of complex infinity. Ch is the only computing environment in existence that can perform numerical computing with consistent numerical results under the IEEE floating-point arithmetic in the entire real domain and complex domain using an extended complex plane for a Riemman sphere. There are two models for implementation of complex numbers under the framework of IEEE floating-point arithmetic. One is to follow the conventional mathematics, while the other is invented in C99's informative Annex G with imaginary types. In this model, the metanumbers (+-Inf,x) (x,+-Inf), (+-Inf,+-0.0) (+-0.0, +-Inf), (+-Inf,NaN), and (NaN,+-Inf) can represent a complex infinity. The sign of zeros is honored. Not-a-number is no longer not a number — it is essentially treated as any number. In short, nonconventional numerical results are invented.

There are major problems in this model:

• First, many results in Annex G are pure invention and mathematically incorrect, according to the conventional complex analysis. For instance, in a well-known example in complex analysis, exp(ComplexInf) is undefined because the complex infinity can be approached from different directions, leading to different numerical results. But Annex G invented exp(Inf,NaN)==(+-0,+-0), exp(Inf,NaN)==(+-Inf,NaN), exp(-Inf,Inf)==(+-0,+-0), exp(Inf,Inf)==(Inf,NaN), and the like. In this model, real and imaginary parts of a complex number are treated as a whole in some cases. In other cases, they are treated as if completely independent; for example, imag(complex(NaN,x)) gives x.
• Second, the model has to introduce a new data type of imaginary in the type system to achieve the aforementioned nonconventional numerical results. In fact, three new data types — float imaginary, double imaginary, and long double imaginary — have been added to the type system. It is expensive to add a type system to C because many operators and type conversion rules have to be overloaded internally, which in turn slows compilation and execution speed, even for programs that don't use imaginary types.

• Third, the new imaginary type is not orthogonal with other types. The set of imaginaries is closed under addition but not multiplication, so it isn't a proper number system under the field operations. They do not form a self-contained subspace. For example, the square root of an imaginary does not return a pure imaginary. Few (if any) real-world applications use only pure-imaginary numbers. Imaginaries are typically embedded in the whole complex plane. Nevertheless, additional functions are needed to handle mathematical functions with arguments of imaginary type.

• Fourth, this model dramatically complicates complex arithmetic operations and complex functions, which leads to inefficiency for numerical computing.

• Fifth, this model also shifts the burden of keeping track of zeros to users. It is almost impossible to get the correct computing results without constantly consulting a manual, especially for complex functions with multiple arguments. If the sign of zero is ignored in complex numbers, you can ignore the sign of zeros. If the sign is honored, you have to worry about the sign of zeros.

And most importantly, inventions such as exp(-Inf,NaN)==(+-0,+-0) in Annex G have no practical use in engineering and science. These inventions, which run counter to conventional complex mathematical analysis, are even harmful; see Example 3.

double x, y, r;
double complex z;
double tolerance = DBL_EPSILON;
 ...
z = cexp(complex(x/y, sqrt(r));
if(cabs(z) < tolerance)
  move_needle_one_inch_for_robot_brain_surgery();
else
  error_handling_for_some_thing_wrong();

Example 3: Harmful Annex G inventions.

If it happens that the numerical value for x is any positive number such as 3.0, y is -0.0, and r is NaN at a singularity point of a robot. This code works fine for implementations such as Ch, based on an alternative complex model. But based on the invention in Annex G, a patient's brain might be damaged. Propagating NaN to a valid, but erroneous, numerical result according to Annex G has a serious consequence for applications in engineering and science.

In short, Annex G is neither normative nor mandated for C99-conforming implementations. Consequently, the C99-conforming implementation of Ch uses a different complex model in the spirit of IEEE 754 floating-point arithmetic — one that follows the conventional complex analysis under a Riemann sphere. In this model, there is only one complex infinity (Inf, Inf)==ComplexInf in the extended complex domain, one Complex-Not-a-Number (NaN, NaN)==ComplexNaN, and zero is unsigned. This complex model is easy to implement and use.

In conventional complex analysis under a Riemann sphere, complex numbers can be represented in the extended complex plane. In Figure 2, there is a one-to-one correspondence between the points on the Riemann sphere and the points on the extended complex plane C. The point p on the surface of the sphere is determined by the intersection of the line through the point z and the north pole N of the sphere. There is only one complex infinity in the extended complex plane. The north pole N corresponds to the point at infinity.

Figure 2: Riemann sphere G and extended complex plane.

Because of the finite representation of floating-point numbers, complex numbers in Ch are programmed in the extended finite complex plane; see Figure 3.

Figure 3: Unit sphere and extended finite complex plane.

Any complex or float complex values inside the ranges of |x|<FLT_MAX and |y|<FLT_MAX are representable in finite floating-point numbers. Variable x is used to represent the real part of a complex number and y the imaginary part; FLT_MAX, a predefined system constant in header file float.h, is the maximum representable finite floating-point number in the float data type. For double complex, the boundary for the extended finite complex plane will be expanded from FLT_MAX to DBL_MAX. Outside the extended finite complex plane, a complex number is treated as a complex infinity represented as ComplexInf or complex(Inf,Inf) in Ch. The origin of the extended finite complex plane is complex(0.0, 0.0). In Ch, an undefined or mathematically indeterminate complex number is denoted as complex(NaN, NaN) or ComplexNaN (short for "Complex-Not-a-Number"). The special complex numbers of ComplexInf and ComplexNaN are referred to as "complex metanumbers." Because of the mathematical infinities of +-, it becomes necessary to distinguish a positive zero 0.0 from a negative zero -0.0 for real numbers. Unlike the real line, along which real numbers can approach the origin through the positive or negative numbers, the origin of the complex plane can be reached in any direction in terms of the limit value of limr0rei, where r is the modulus and the phase of a complex number. Therefore, complex operations and complex functions in Ch do not distinguish 0.0 from -0.0 for real and imaginary parts of complex numbers.

Under this complex model, the conventional numerical results such as ComplexInf+ComplexInf==ComplexNaN, ComplexInf-ComplexInf==ComplexNaN, ComplexInf*ComplexInf==ComplexInf, and exp(ComplexInf)==ComplexNaN can be obtained. The ex- pression (i)*(i) is equivalent to ComplexInf*ComplexInf with the result of ComplexInf.

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].

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