2

Question

What is the computational cost of the remainder function, is there a specific instruction to calculate it in a cheap way in a specific case?

Description

I need to transform a mathematical variable x into the range of I=[-0.5; 0.5) from R=[-2; 2). While x is not element of I, then x is shifted towards I by repeatedly adding or subtracting 1 to the value of x. x is represented with double x in my code. I need the fastest way of this transformation for this I and R values but wider R ranges can be also interesting.

Ideas and speed comparison

The function I was suggested to use was the naive implementation following the description:

void shift_to_I(double& x) // version 1
{
    while (x < -0.5)
            x += 1;

    while (x >= 0.5)
            x -= 1;
}

Not only for speed considerations but also for code quality I was thinking of using remainder from <cmath> introduced in c++11. With remainder the code shortens to

void shift_to_I(double& x) // version 2
{
    x = remainder(x,1);
}

I had to realize though that it was slower than the original function on my architecture (Intel i7 whatver with VC++). I believed there was a dedicated instruction for this purpose, but either the compiler doesn't know it or it doesn't exist. For wider R interval (on my architecture it is around [-25; 25)) the second version will be faster but I need a code that is fast for narrow intervals too. clang and gcc specific solutions are also welcomed.

DanielTuzes
  • 2,494
  • 24
  • 40
  • note that there is also `std::fmod`, i didnt understand the difference, but it is not the same as `std::remainder` – 463035818_is_not_an_ai Nov 11 '19 at 14:44
  • Note that special functions line `remainder` and `fmod` have to deal properly with corner cases like `inf` and `nan`. – Evg Nov 11 '19 at 14:49
  • `x -= int(x);` would reduce to `(-1, 1)` interval without a loop; it would take another operation to get to (-0.5, 0.5). No idea how well this will perform; you may want to benchmark. This should perform the same independent of the width of `R` (as long as it's within `int` range). – Igor Tandetnik Nov 11 '19 at 14:59
  • 2
    The Intel x86 and x64 CPUs (FPUs) have the `fprem` and `fprem1` [instructions](https://www.felixcloutier.com/x86/fprem1) (single op, gets remainder). But whether or not a particular compiler uses these for `std::remainder()` (or `std::fmod()`) depends on their implementers' decisions. – Adrian Mole Nov 11 '19 at 15:01
  • @formerlyknownas_463035818 The only difference seem to be that _"In contrast to std::fmod(), the returned value is not guaranteed to have the same sign as x."_ (where x is the first parameter) from (https://en.cppreference.com/w/cpp/numeric/math/remainder) – Lukas-T Nov 11 '19 at 15:10
  • Why not use both solutions? Check the range, if it is small, use the naive approach, if it is large, use the other. A good way to tune this is to use google benchmark to test the performance of the code. – NathanOliver Nov 11 '19 at 15:15
  • "I had to realize though that it was slower than the original function" - Did you test that with an optimized build or only with a debug build? – Jesper Juhl Nov 11 '19 at 15:37
  • @JesperJuhl with Release mode. I also paid attention not to optimize out variables (collecting sum e.g.) and to measure on 1e8 random numbers uniformly distributed in *R*. – DanielTuzes Nov 11 '19 at 15:58
  • @Evg I believe due to IEEE specs, the implementer must take care of those silly cases except one turns on `-ffast-math`. In this case, do you think the implementer removes those code parts checking for `nan`s and `inf`s? But clearly, `fmod` is not the right function here. – DanielTuzes Nov 11 '19 at 16:09
  • With `-ffast-math` GCC issues just a single assembly instruction, no special checks for `inf` and `nan` are performed. But Clang still calls `__ieee754_remainder`. – Evg Nov 11 '19 at 16:17
  • 1
    @Evg While `-ffast-math` is nice and fast and all, the *unpleasant* part of it is that it breaks standard conformance. Personally I avoid it. I'll rather have standard guaranteed behaviour than slightly better performance with sometimes surprising results. But of course, *sometimes* you know what you are doing and know that the `fast-math` shortcuts don't affect you and the extra perf is super important, in that case, sure; use it. – Jesper Juhl Nov 11 '19 at 17:00

1 Answers1

3

This question is compiler and implementation dependent.

For example, on my machine with GCC 8.3:

  1. Without -ffast-math, std::remainder translates into a call to this function:

    double __remainder(double x, double y)
{
    if (((__builtin_expect (y == 0.0, 0) && ! isnan(x)) || (__builtin_expect(isinf(x), 0) && ! isnan(y))) && _LIB_VERSION != _IEEE_)
        return __kernel_standard(x, y, 28);
    return __ieee754_remainder(x, y);
}

    with __ieee754_remainder looking like this:

    double __ieee754_remainder(double x, double y)
{
    double z, d, xx;
    int4 kx, ky, n, nn, n1, m1, l;
    mynumber u, t, w = {{0, 0}}, v = {{0, 0}}, ww = {{0, 0}}, r;
    u.x = x;
    t.x = y;
    kx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign  for x*/
    t.i[HIGH_HALF] &= 0x7fffffff;     /*no sign for y */
    ky = t.i[HIGH_HALF];
    /*------ |x| < 2^1023  and   2^-970 < |y| < 2^1024 ------------------*/
    if (kx < 0x7fe00000 && ky < 0x7ff00000 && ky >= 0x03500000)
    {
        SET_RESTORE_ROUND_NOEX(FE_TONEAREST);
        if (kx + 0x00100000 < ky)
            return x;
        if ((kx - 0x01500000) < ky)
        {
            z = x / t.x;
            v.i[HIGH_HALF] = t.i[HIGH_HALF];
            d = (z + big.x) - big.x;
            xx = (x - d * v.x) - d * (t.x - v.x);
            if (d - z != 0.5 && d - z != -0.5)
                return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x);
            else
            {
                if (fabs(xx) > 0.5 * t.x)
                    return (z > d) ? xx - t.x : xx + t.x;
                else
                    return xx;
            }
        } /*    (kx<(ky+0x01500000))         */
        else
        {
            r.x = 1.0 / t.x;
            n = t.i[HIGH_HALF];
            nn = (n & 0x7ff00000) + 0x01400000;
            w.i[HIGH_HALF] = n;
            ww.x = t.x - w.x;
            l = (kx - nn) & 0xfff00000;
            n1 = ww.i[HIGH_HALF];
            m1 = r.i[HIGH_HALF];
            while (l > 0)
            {
                r.i[HIGH_HALF] = m1 - l;
                z = u.x * r.x;
                w.i[HIGH_HALF] = n + l;
                ww.i[HIGH_HALF] = (n1) ? n1 + l : n1;
                d = (z + big.x) - big.x;
                u.x = (u.x - d * w.x) - d * ww.x;
                l = (u.i[HIGH_HALF] & 0x7ff00000) - nn;
            }
            r.i[HIGH_HALF] = m1;
            w.i[HIGH_HALF] = n;
            ww.i[HIGH_HALF] = n1;
            z = u.x * r.x;
            d = (z + big.x) - big.x;
            u.x = (u.x - d * w.x) - d * ww.x;
            if (fabs(u.x) < 0.5 * t.x)
                return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x);
                else if (fabs(u.x) > 0.5 * t.x)
                return (d > z) ? u.x + t.x : u.x - t.x;
            else
            {
                z = u.x / t.x;
                d = (z + big.x) - big.x;
                return ((u.x - d * w.x) - d * ww.x);
            }
        }
    } /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
    else
    {
        if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0))
        {
            y = fabs(y) * t128.x;
            z = __ieee754_remainder(x, y) * t128.x;
            z = __ieee754_remainder(z, y) * tm128.x;
            return z;
        }
        else
        {
            if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 &&
                (ky > 0 || t.i[LOW_HALF] != 0))
            {
                y = fabs(y);
                z = 2.0 * __ieee754_remainder(0.5 * x, y);
                d = fabs(z);
                if (d <= fabs(d - y))
                    return z;
                else if (d == y)
                    return 0.0 * x;
                else
                    return (z > 0) ? z - y : z + y;
            }
            else /* if x is too big */
            {
                if (ky == 0 && t.i[LOW_HALF] == 0) /* y = 0 */
                    return (x * y) / (x * y);
                else if (kx >= 0x7ff00000    /* x not finite */
                         || (ky > 0x7ff00000 /* y is NaN */
                                || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0)))
                    return (x * y) / (x * y);
                else
                    return x;
            }
        }
    }
}

    Pretty far from a single machine instruction.

  2. With -ffast-math, single fprem1 assembly instruction is used.

Evg
  • 25,259
  • 5
  • 41
  • 83
  • With a little benchmark using `-ffast-math`, I got that in case *R* = [-2; 2), `remainder` function is faster than the loop-based one by 10%. For *R* = [-1.5; 1.5), the `remainder` is slower and the loop-based is faster, again, around 10%. – DanielTuzes Nov 11 '19 at 16:16
  • @DanielTuzes, and without `-ffast-math`? – Evg Nov 11 '19 at 16:18