Implementing the modulus operator for floating-point numbers

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

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;
}

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

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

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.

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

Where AbsMod represents the modulus of \(|x|\) and \(|y|\):

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

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);
}