Implementing the modulus operator for floating-point numbers

Index | Up


Table of Contents

1. Motivation

This article was very remotely inspired by Jacob Smith’s article about color palettes. I started writing some Scheme code for generating color palettes, and I decided to first write some functions for converting between RGB and HSV formats. Then, I realized that Scheme (unlike Emacs Lisp, for example) doesn’t support floating-point inputs in its modulus function.

I decided to write my own fmod function, starting with a simple version that only supported positive values, and eventually adding support for negative values, even though this was not needed when converting between RGB and HSV. I started writing the process in my scratch repository, but I think it has become interesting enough to deserve its own place in this blog.

2. Simple version for positive values

The first function I wrote simply keeps subtracting the divisor to the dividend until the dividend is smaller than the divisor.

;; Calculate remainder of X by Y. Supports floating-point, as long as they are
;; positive.
(define (basic-fmod x y)
  (if (< x y)
      x
      (basic-fmod (- x y) y)))

I also wrote an iterative C version.

/* Calculate remainder of X and Y, as long as they are positive. */
double basic_fmod(double x, double y) {
    while (x >= y)
        x -= y;
    return x;
}

3. Supporting negative values

There are two main approaches for handling negative values: the one used by the fmod function in the math.h C header and the one used by the mod function in Emacs Lisp. Let’s have a look at the outputs of both functions, and how they can be implemented.

3.1. C’s fmod

Let’s have a look at some outputs from the fmod function in the math header.

X Y Result
9.5 2.5 2.00
9.5 -2.5 2.00
-9.5 2.5 -2.00
-9.5 -2.5 -2.00

As you can see, the result is always the modulus of the absolute values, while keeping the original sign of x. We can easily implement this by adding some calls to abs and wrapping our previous function with a sign conditional.

;; Scheme equivalent of C's `fmod'.
(define (c-fmod x y)
  (define (basic-fmod x y)
    (if (< x y)
        x
        (basic-fmod (- x y) y)))

  (let ((result (basic-fmod (abs x) (abs y))))
    (if (negative? x)
        (- result)
        result)))

This is the iterative C version. I decided to use an int instead of a bool for the bNegativeResult variable because I didn’t want to include stdbool.h just for this. I used hungarian notation to emphasize that it acts as a boolean.

#include <math.h> /* fabs() */

/* Equivalent of `fmod' from math.h  */
double my_fmod(double x, double y) {
    int bNegativeResult = x < 0;

    x = fabs(x);
    y = fabs(y);

    while (x >= y)
        x -= y;

    return bNegativeResult ? -x : x;
}

3.2. Emacs’ mod

Now let’s compare the previous results with the ones from Emacs.

X Y Result
9.5 2.5 2.00
9.5 -2.5 -0.50
-9.5 2.5 0.50
-9.5 -2.5 -2.00

These results are better in my opinion, since they satisfy the following formula.

\begin{equation*}
  \lfloor a / b \rfloor \times b + a \bmod b = a
\end{equation*}

Where \lfloor a / b \rfloor indicates the integer division of a and b.

Before explaining my approach, I want to show how Emacs’ mod actually works. The actual C function is called fmod_float and is located in the src/floatfns.c file. Omitting the emacs-specific parts, we get the following function.

#include <math.h> /* fmod() */

double my_emacs_fmod(double x, double y) {
    x = fmod(x, y);

    /* If the "remainder" comes out with the wrong sign, fix it. */
    if (y < 0 ? x > 0 : x < 0)
        x += y;

    return x;
}

I want to note that, although Emacs’ obviously uses the real fmod from math.h, the previous my_fmod function can be used here as well.

As you can see, the only part that differences the Emacs modulus from the C modulus is the conditional in the middle. We could simply implement this behavior in Scheme by adding the missing conditional, but I think it’s better if we adapt our previous function.

3.3. My approach

If we look again at the outputs from Emacs’ mod, we can see that the changes in the output match the following table.

X Y Result
Positive Positive AbsMod(x, y)
Positive Negative y + AbsMod(x, y)
Negative Positive y - AbsMod(x, y)
Negative Negative -AbsMod(x, y)

Where AbsMod represents the modulus of |x| and |y|:

\begin{equation*}
  \text{AbsMod}(x, y) = |x| \bmod |y|
\end{equation*}

The table can also be expressed as a conditional formula, if you are into that.

\begin{equation*}
  x \bmod b =
  \begin{cases}
    \text{AbsMod}(x, y),     &amp; x \geq 0 \land y \geq 0 \\
    y + \text{AbsMod}(x, y), &amp; x \geq 0 \land y &lt; 0 \\
    y - \text{AbsMod}(x, y), &amp; x &lt; 0 \land y \geq 0 \\
    -\text{AbsMod}(x, y),    &amp; x &lt; 0 \land y &lt; 0
  \end{cases}
\end{equation*}

With this, we can make a final fmod version.

;; Calculate remainder of X by Y, supporting floating point and negative values.
(define (fmod x y)
  (define (basic-fmod x y)
    (if (< x y)
        x
        (basic-fmod (- x y) y)))

  (let ((abs-result (fmod-positive (abs x) (abs y))))
    (cond ((and (positive? x) (positive? y))
           abs-result)
          ((and (positive? x) (negative? y))
           (+ y abs-result))
          ((and (negative? x) (positive? y))
           (- y abs-result))
          ((and (negative? x) (negative? y))
           (- abs-result)))))

There are some unnecessary calls to positive? and negative?, but I think it’s clearer this way. This issue does not happen in the following C version.

#include <math.h> /* fabs() */

double emacs_fmod(double x, double y) {
    const double abs_x = fabs(x);
    const double abs_y = fabs(y);

    /* Calculate fmod(fabs(x), fabs(y)) */
    double abs_result = abs_x;
    while (abs_result >= abs_y)
        abs_result -= abs_y;

    return (x >= 0) ? ((y >= 0) ? abs_result : y + abs_result)
                    : ((y >= 0) ? y - abs_result : -abs_result);
}

4. Final note

Most of these functions were made just by observing the output of already existing functions, so I don’t know if this is the most optimal or mathematical approach. If you have any suggestions or improvements, feel free to contribute to this page.