Where the Operators Run Out

C's arithmetic operators stop at %. There is no operator for a square root, a logarithm, or a power, and no operator for the remainder of two double values. Those live in the math library, and reaching them takes one #include and, on most systems, one extra word on the compiler command line.

The second half of this lesson is the uncomfortable half. Real arithmetic on a computer is not the arithmetic of the algebra textbook: some results are approximations, some have no defined value at all, and both kinds go wrong quietly. Knowing exactly where the guarantees stop is the difference between a program you can trust and one that happens to work on the numbers you tried.

The Math Library

Declaring the functions takes one header:

#include <math.h>

The C89 set covers what exam questions ask for. Every one of these takes double arguments and returns a double:

  • sqrt(x) square root, pow(x, y) for x to the power y
  • fabs(x) absolute value of a floating point number
  • floor(x) rounds down to a whole number, ceil(x) rounds up
  • fmod(x, y) remainder of x / y, the % that floating point values are denied
  • log(x) natural logarithm, log10(x) base 10, exp(x) for e to the x
  • sin(x), cos(x), tan(x) and their inverses, all working in radians

A quadratic's discriminant and roots, the classic exam program, is sqrt plus arithmetic:

#include <stdio.h>
#include <math.h>

int main(void)
{
    double a = 0.0;
    double b = 0.0;
    double c = 0.0;
    double discriminant = 0.0;

    printf("Enter a, b and c: ");
    scanf("%lf %lf %lf", &a, &b, &c);
    discriminant = b * b - 4.0 * a * c;

    printf("root1 = %.2f\n", (-b + sqrt(discriminant)) / (2.0 * a));
    printf("root2 = %.2f\n", (-b - sqrt(discriminant)) / (2.0 * a));
    return 0;
}

Given 2, 4 and -16 that prints roots of 2.00 and -4.00. It is shown rather than run here for a reason, which is the next section. (A negative discriminant has no real square root, and handling that case needs the if of the next chapter.)

The Linker Step Everyone Forgets

On Unix systems the math functions live in a library separate from the rest of the C runtime, so the compiler needs to be told to link it. Compiling the program above without that produces no complaint about your code at all, only this:

/usr/bin/ld: /tmp/cc6j5xuT.o: in function `main':
roots.c:(.text.startup+0x64): undefined reference to `sqrt'
/usr/bin/ld: roots.c:(.text.startup+0x94): undefined reference to `sqrt'
collect2: error: ld returned 1 exit status

Nothing is wrong with the source. The compiler accepted every line, and the linker then failed to find the machine code for the function, once for each of the two calls. The fix is -lm, and it goes at the end of the command:

gcc -std=c89 -pedantic-errors -Wall -Wextra roots.c -o roots -lm

Read "undefined reference" to a name you did not write yourself as a missing library, not a missing declaration. It is the most common first encounter with the difference between compiling and linking.

This course's runnable examples compile without -lm, so the sqrt program above is one to run with your own compiler. fabs, floor and ceil are the exceptions that do run here, because GCC recognizes them and computes them inline rather than calling the library; fabs appears in a runnable box further down. One more reason the program above reads its values instead of hard-coding them: given constants, GCC computes sqrt at compile time and the link succeeds without -lm, which teaches the wrong lesson.

Always Include the Header

Leaving out #include <math.h> is not caught by C89 the way you might hope. The language allows a call to a function it has never seen and assumes it returns int, so the code compiles, with a warning:

warning: implicit declaration of function 'sqrt' [-Wimplicit-function-declaration]
    7 |     r = sqrt(2.0);
      |         ^~~~

Then the program calls a function that really returns double as though it returned int, which is undefined behaviour: the value you get back is not a wrong number so much as no promised number at all. Treat that warning as an error, always.

Round-off: Where Decimal Fractions Go

Type 0.1 into a program and the machine stores the nearest double it can represent, which is not exactly one tenth. Binary fractions represent halves, quarters and eighths exactly; a tenth is a repeating fraction in binary, exactly as a third is in decimal. Print enough digits and the gap shows:

0.1 + 0.2 = 0.30000000000000004
0.3       = 0.29999999999999999
sum == 0.3 gives 0

Both values are slightly wrong, they are wrong by different amounts, and so the comparison that any algebra student would call true comes back 0. The usual %f would have printed 0.300000 for both and hidden the whole affair, which is precisely why this bites in real programs.

Every operation rounds, and rounding accumulates:

Ten tenths comes to 0.99999999999999989. Nothing here is broken; each addition gave the closest double to the true sum, and nine roundings in a row drifted.

An exam note, since textbook versions of this section like to assert that (1.0 / 3.0) * 3.0 differs from 1.0. Run it and it compares equal, in float and in double both: the division rounds one way and the multiplication rounds back. The real guarantee is weaker and stranger than "the answer will be slightly off". It is that no particular answer is promised, so you must not depend on either outcome.

Never Compare Floating Point with ==

The rule that follows from all of this: compare floating point values by asking whether they are close enough, using fabs to measure the distance between them.

The tolerance is a judgement about your problem, not a universal constant: money to the nearest cent tolerates far more slack than a physics simulation. What is never right is ==. The same warning covers < and > at the boundary, where a value that "should" be exactly the limit may sit either side of it.

Division by Zero

Integer division by zero, n / 0 or n % 0, is undefined behaviour. Not an exception, not a guaranteed crash, not infinity: the standard imposes no requirement whatsoever, so a program containing it has no defined meaning at all.

/* WRONG: undefined behaviour whenever count is zero */
average = total / count;

The guard is the short-circuit && from earlier in this chapter, which promises not to evaluate its right side when the left is false:

count != 0 && total / count > threshold

Floating point division by zero is a different question with a different answer. Where the implementation follows IEEE 754, as the platforms you will meet do, 1.0 / 0.0 yields infinity and 0.0 / 0.0 yields a not-a-number value, both of which propagate through later arithmetic instead of ending the program. Defined, but rarely what you wanted.

Overflow and Underflow

Signed integer overflow, a result too large for the type, is also undefined behaviour. INT_MAX + 1 is not "wraps around to INT_MIN" in C, whatever the machine's instruction happens to do; the compiler is entitled to optimize on the assumption that it never occurs. Unsigned arithmetic is the exception: it wraps, and the standard says so, which is why unsigned types are the right ones for bit patterns and hash values.

Floating point has its own two ends. A result too large in magnitude overflows to infinity, and one too small underflows toward zero, losing all its significant digits on the way. Both are defined and both destroy your answer, so keep intermediate values inside the range the type can hold: sqrt(b * b - 4.0 * a * c) can overflow on inputs where the roots themselves are perfectly ordinary.

One last thing worth knowing about this whole family. The sandbox running your code checks memory errors with AddressSanitizer, and arithmetic undefined behaviour is invisible to it. A clean run is not evidence that your integer division is safe.

Key Takeaways

  • math.h declares sqrt, pow, fabs, floor, ceil, fmod, log, exp and the trigonometric functions; in C89 all of them take and return double.
  • fmod(x, y) is the remainder for floating point values, since % accepts integers only.
  • Link with -lm on Unix systems: "undefined reference" to a library function is a linker error, not a compiler error.
  • A missing #include <math.h> compiles with a warning and calls a double function as though it returned int, which is undefined behaviour.
  • 0.1 and 0.2 are not exact in binary, so 0.1 + 0.2 == 0.3 is false and roundings accumulate over repeated operations.
  • Never compare floating point values with ==; test fabs(a - b) < tolerance with a tolerance chosen for the problem.
  • Integer division by zero and signed integer overflow are undefined behaviour; unsigned overflow wraps, and floating point division by zero gives infinity or not-a-number under IEEE 754.