Channels ▼
RSS

C/C++

C99 and Numeric Computing

Source Code Accompanies This Article. Download It Now.


Although C is a powerful systems programming language that can deliver much of the same control over devices as assembly language, it has deficiencies when it comes to scientific and engineering applications that require extensive numerical computing. While some numerical computing deficiencies in the original K&R C were addressed in C90, many limitations still remain.

C99, ratified as the ANSI/ISO C Standard (ISO/IEC IS 9899), is a milestone in C's evolution into a viable programming language for scientific and numerical computing. Among other features, C99 supports IEEE floating-point arithmetic, complex numbers, and variable-length arrays (VLAs) for numerical programming. Complex numbers and VLAs were added mainly based on the prior art of implementation of Ch from SoftIntegration, SCC from Cray Research, gcc from Free Software Foundation, and others.

Although support for C99 is limited, more compilers are adding these new features. For example, Comeau C 4.2.45.2 from Comeau Computing supports VLAs without complex numbers. The C compiler from Hewlett-Packard supports VLAs and complex numbers. GCC 2.95 and later provide limited support of VLAs and complex numbers. The Dinkum C99 Library from Dinkumware is a complete library for C99.

Complex numbers are handled as built-in data types in C. With C++, on the other hand, complex numbers are treated as classes. For example, Forte C++ 6 Update 2 (formerly Sun Visual WorkShop C++) provides some support for complex arithmetic. Still, there is no provision for IEEE floating-point arithmetic for both real numbers and complex numbers in C++.

In this article, I'll examine what's involved with C99 compliance by looking at Ch, a C virtual machine produced by SoftIntegration, a company I founded. Ch has provisions for consistent handling of numerical numbers in the entire real and complex domains with VLAs under the framework of IEEE floating-point arithmetic. As a superset of C interpreter, Ch also supports classes in C++. In particular, I'll focus on the design and implementation of IEEE floating-point arithmetic and complex numbers.

Ch was designed for script computing in the framework of C/C++. C or C++ conforming programs with complex numbers will run in this virtual machine without modification. I've tested all C programs presented here with both Ch and GCC 2.96, and C++ programs with Ch and G++ 2.95.

Computing in the Entire Real Domain

The IEEE 754 standard for binary floating-point arithmetic is significant for consistent floating-point arithmetic with respect to real numbers. IEEE 754 distinguishes +0.0 from -0.0, which introduces additional programming complexity. Another important IEEE 754 feature is the internal representations for the mathematical infinities and invalid values. The mathematical infinity is represented by Inf. A mathematically indeterminate or an undefined value such as division of zero by zero is represented by NaN (short for "Not-a-Number"). Many computers support in-hardware signed zeros, infinity, and NaN. However, information about low-level and limited high-level instruction sets provided by hardware vendors may not be relevant to application programmers, and most features of a final system depend on the software implementation. Even for IEEE machines, there is no provision for propagating via software the sign of zeros, infinity, and NaN in a consistent and useful manner. They must be programmed as if zeros were unsigned, without infinity and NaN. Based on IEEE machines, some vendors provide software support for IEEE 754 through libraries. However, these special values in libraries are not transparent to programmers. Although the application of symbols such as Inf and NaN can be found in some mathematical packages and libraries, vendor handling of these special numbers is often flawed. These are the gray areas in which IEEE 754 is not supported in many hardware and software systems.

To make the power of IEEE 754 easily available to you, C99 introduced the floating-point numbers INFINITY, -INFINITY, NAN, and signed zeros -0.0 and 0.0. In Ch, INFINITY and NAN — which correspond to the built-in metanumbers Inf and NaN — are defined as macros in the header file math.h. For convenience, I use the metanumbers Inf and NaN, which are transparent to programmers. Signed zeros (+0.0 and -0.0) in C99 behave like correctly signed infinitesimal quantities 0+ and 0-, whereas symbols Inf and -Inf correspond to mathematical infinities and -, respectively. IEEE 754 only addresses the arithmetic involving these metanumbers. These metanumbers are extended in C99 to commonly used mathematical functions in the spirit of IEEE 754. Ch includes provisions for consistent handling of metanumbers in I/O, arithmetic, relational and logic operations, and polymorphic mathematical functions. An NaN is propagated consistently through subsequent computations. Many people believe the C99 committee errored in handling some mathematical functions. For example, the values of function calls for hypot(Inf,NaN), hypot(-Inf,NaN), pow(1,NaN), and pow(NaN,+/-0.0) are defined in C99 as Inf, Inf, 1.0, and 1.0, respectively. In Ch, I implement them to return NaN because these functions are mathematically undefined for the arguments with the aforementioned values.

For real numbers, C99 distinguishes -0.0 from 0.0. The metanumbers 0.0, -0.0, Inf, -Inf, and NaN are useful for scientific computing. For example, the function f(x)=e1/x is not continuous at the origin; see Figure 1.


Figure 1: Function f(x)=e1/x.

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.


Related Reading


More Insights






Currently we allow the following HTML tags in comments:

Single tags

These tags can be used alone and don't need an ending tag.

<br> Defines a single line break

<hr> Defines a horizontal line

Matching tags

These require an ending tag - e.g. <i>italic text</i>

<a> Defines an anchor

<b> Defines bold text

<big> Defines big text

<blockquote> Defines a long quotation

<caption> Defines a table caption

<cite> Defines a citation

<code> Defines computer code text

<em> Defines emphasized text

<fieldset> Defines a border around elements in a form

<h1> This is heading 1

<h2> This is heading 2

<h3> This is heading 3

<h4> This is heading 4

<h5> This is heading 5

<h6> This is heading 6

<i> Defines italic text

<p> Defines a paragraph

<pre> Defines preformatted text

<q> Defines a short quotation

<samp> Defines sample computer code text

<small> Defines small text

<span> Defines a section in a document

<s> Defines strikethrough text

<strike> Defines strikethrough text

<strong> Defines strong text

<sub> Defines subscripted text

<sup> Defines superscripted text

<u> Defines underlined text

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task. However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

 
Disqus Tips To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.
 

Video