C99 and Numeric Computing

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.