Fastest way to clamp a real (fixed/floating point) value?

Question

Hi, Is there a more efficient way to clamp real numbers than using if statements or ternary operators? I want to do this both for doubles and for a 32-bit fixpoint implementation (16.16). I'm not asking for code that can handle both cases; they will be handled in separate functions.

Obviously, I can do something like:

double clampedA;
double a = calculate();
clampedA = a > MY_MAX ? MY_MAX : a;
clampedA = a < MY_MIN ? MY_MIN : a;

or

double a = calculate();
double clampedA = a;
if(clampedA > MY_MAX)
    clampedA = MY_MAX;
else if(clampedA < MY_MIN)
    clampedA = MY_MIN;

The fixpoint version would use functions/macros for comparisons.

This is done in a performance-critical part of the code, so I'm looking for an as efficient way to do it as possible (which I suspect would involve bit-manipulation)

EDIT: It has to be standard/portable C, platform-specific functionality is not of any interest here. Also, MY_MIN and MY_MAX are the same type as the value I want clamped (doubles in the examples above).

Answer 1

If I understand properly, you want to limit a value "a" to a range between MY_MIN and MY_MAX. The type of "a" is a double. You did not specify the type of MY_MIN or MY_MAX.

The simple expression:

clampedA = (a > MY_MAX)? MY_MAX : (a < MY_MIN)? MY_MIN : a;

should do the trick.

I think there may be a small optimization to be made if MY_MAX and MY_MIN happen to be integers:

int b = (int)a;
clampedA = (b > MY_MAX)? (double)MY_MAX : (b < MY_MIN)? (double)MY_MIN : a;

By changing to integer comparisons, it is possible you might get a slight speed advantage.

Answer 2

I think you could use SSE3 or some similar technology for this, but do not know exactly which commands/how... You can take a look at: Saturation arithmetic

Answer 3

Realistically, no decent compiler will make a difference between an if() statement and a ?: expression. The code is simple enough that they'll be able to spot the possible paths. That said, your two examples are not identical. The equivalent code using ?: would be

a = (a > MAX) ? MAX : ((A < MIN) ? MIN : a);

as that avoid the A < MIN test when a > MAX. Now that could make a difference, as the compiler otherwise would have to spot the relation between the two tests.

If clamping is rare, you can test the need to clamp with a single test:

if (abs(a - (MAX+MIN)/2) > ((MAX-MIN)/2)) ...

E.g. with MIN=6 and MAX=10, this will first shift a down by 8, then check if it lies between -2 and +2. Whether this saves anything depends a lot on the relative cost of branching.

Answer 4

The bits of IEEE 754 floating point are ordered in a way that if you compare the bits interpreted as an integer you get the same results as if you would compare them as floats directly. So if you find or know a way to clamp integers you can use it for (IEEE 754) floats as well. Sorry, I don't know a faster way.

If you have the floats stored in an arrays you can consider to use some CPU extensions like SSE3, as rkj said. You can take a look at liboil it does all the dirty work for you. Keeps your program portable and uses faster cpu instructions if possible. (I'm not sure tho how OS/compiler-independent liboil is).

Answer 5

For the 16.16 representation, the simple ternary is unlikely to be bettered speed-wise.

And for doubles, because you need it standard/portable C, bit-fiddling of any kind will end badly.

Even if a bit-fiddle was possible (which I doubt), you'd be relying on the binary representation of doubles. THIS (and their size) IS IMPLEMENTATION-DEPENDENT.

Possibly you could "guess" this using sizeof(double) and then comparing the layout of various double values against their common binary representations, but I think you're on a hiding to nothing.

The best rule is TELL THE COMPILER WHAT YOU WANT (ie ternary), and let it optimise for you.

EDIT: Humble pie time. I just tested quinmars idea (below), and it works - if you have IEEE-754 floats. This gave a speedup of about 20% on the code below. IObviously non-portable, but I think there may be a standardised way of asking your compiler if it uses IEEE754 float formats with a #IF...?

  double FMIN = 3.13;
  double FMAX = 300.44;

  double FVAL[10] = {-100, 0.23, 1.24, 3.00, 3.5, 30.5, 50 ,100.22 ,200.22, 30000};
  uint64  Lfmin = *(uint64 *)&FMIN;
  uint64  Lfmax = *(uint64 *)&FMAX;

    DWORD start = GetTickCount();

    for (int j=0; j<10000000; ++j)
    {
     uint64 * pfvalue = (uint64 *)&FVAL[0];
     for (int i=0; i<10; ++i)
      *pfvalue++ = (*pfvalue < Lfmin) ? Lfmin : (*pfvalue > Lfmax) ? Lfmax : *pfvalue;
    }

    volatile DWORD hacktime = GetTickCount() - start;

    for (int j=0; j<10000000; ++j)
    {
     double * pfvalue = &FVAL[0];
     for (int i=0; i<10; ++i)
      *pfvalue++ = (*pfvalue < FMIN) ? FMIN : (*pfvalue > FMAX) ? FMAX : *pfvalue;
    }

    volatile DWORD normaltime = GetTickCount() - (start + hacktime);

Answer 6

Here's a possibly faster implementation similar to @Roddy's answer:

typedef int64_t i_t;
typedef double  f_t;

static inline
i_t i_tmin(i_t x, i_t y) {
  return (y + ((x - y) & -(x < y))); // min(x, y)
}

static inline
i_t i_tmax(i_t x, i_t y) {
  return (x - ((x - y) & -(x < y))); // max(x, y)
}

f_t clip_f_t(f_t f, f_t fmin, f_t fmax)
{
#ifndef TERNARY
  assert(sizeof(i_t) == sizeof(f_t));
  //assert(not (fmin < 0 and (f < 0 or is_negative_zero(f))));
  //XXX assume IEEE-754 compliant system (lexicographically ordered floats)
  //XXX break strict-aliasing rules
  const i_t imin = *(i_t*)&fmin;
  const i_t imax = *(i_t*)&fmax;
  const i_t i    = *(i_t*)&f;
  const i_t iclipped = i_tmin(imax, i_tmax(i, imin));

#ifndef INT_TERNARY
  return *(f_t *)&iclipped;
#else /* INT_TERNARY */
  return i < imin ? fmin : (i > imax ? fmax : f); 
#endif /* INT_TERNARY */

#else /* TERNARY */
  return fmin > f ? fmin : (fmax < f ? fmax : f);
#endif /* TERNARY */
}

See Compute the minimum (min) or maximum (max) of two integers without branching and Comparing floating point numbers

The IEEE float and double formats were designed so that the numbers are “lexicographically ordered”, which – in the words of IEEE architect William Kahan means “if two floating-point numbers in the same format are ordered ( say x < y ), then they are ordered the same way when their bits are reinterpreted as Sign-Magnitude integers.”

A test program:

/** gcc -std=c99 -fno-strict-aliasing -O2 -lm -Wall *.c -o clip_double && clip_double */
#include <assert.h> 
#include <iso646.h>  // not, and
#include <math.h>    // isnan()
#include <stdbool.h> // bool
#include <stdint.h>  // int64_t
#include <stdio.h>

static 
bool is_negative_zero(f_t x) 
{
  return x == 0 and 1/x < 0;
}

static inline 
f_t range(f_t low, f_t f, f_t hi) 
{
  return fmax(low, fmin(f, hi));
}

static const f_t END = 0./0.;

#define TOSTR(f, fmin, fmax, ff) ((f) == (fmin) ? "min" :    \
         ((f) == (fmax) ? "max" :  \
          (is_negative_zero(ff) ? "-0.": \
           ((f) == (ff) ? "f" : #f))))

static int test(f_t p[], f_t fmin, f_t fmax, f_t (*fun)(f_t, f_t, f_t)) 
{
  assert(isnan(END));
  int failed_count = 0;
  for ( ; ; ++p) {
    const f_t clipped  = fun(*p, fmin, fmax), expected = range(fmin, *p, fmax);
    if(clipped != expected and not (isnan(clipped) and isnan(expected))) {
      failed_count++;
      fprintf(stderr, "error: got: %s, expected: %s\t(min=%g, max=%g, f=%g)\n", 
          TOSTR(clipped,  fmin, fmax, *p), 
          TOSTR(expected, fmin, fmax, *p), fmin, fmax, *p);
    }
    if (isnan(*p))
      break;
  }
  return failed_count;
}  

int main(void)
{
  int failed_count = 0;
  f_t arr[] = { -0., -1./0., 0., 1./0., 1., -1., 2, 
     2.1, -2.1, -0.1, END};
  f_t minmax[][2] = { -1, 1,  // min, max
            0, 2, };

  for (int i = 0; i < (sizeof(minmax) / sizeof(*minmax)); ++i) 
    failed_count += test(arr, minmax[i][0], minmax[i][1], clip_f_t);      

  return failed_count & 0xFF;
}

In console:

$ gcc -std=c99 -fno-strict-aliasing -O2 -lm *.c -o clip_double && ./clip_double

It prints:

error: got: min, expected: -0.  (min=-1, max=1, f=0)
error: got: f, expected: min    (min=-1, max=1, f=-1.#INF)
error: got: f, expected: min    (min=-1, max=1, f=-2.1)
error: got: min, expected: f    (min=-1, max=1, f=-0.1)

Answer 7

Ternary operator is really the way to go, because most compilers are able to compile them into a native hardware operation that uses a conditional move instead of a branch (and thus avoids the mispredict penalty and pipeline bubbles and so on). Bit-manipulation is likely to cause a load-hit-store.

In particular, PPC and x86 with SSE2 have a hardware op that could be expressed as an intrinsic something like this:

double fsel( double a, double b, double c ) {
  return a >= 0 ? b : c; 
}

The advantage is that it does this inside the pipeline, without causing a branch. In fact, if your compiler uses the intrinsic, you can use it to implement your clamp directly:

inline double clamp ( double a, double min, double max ) 
{
   a = fsel( a - min , a, min );
   return fsel( a - max, max, a );
}

I strongly suggest you avoid bit-manipulation of doubles using integer operations. On most modern CPUs there is no direct means of moving data between double and int registers other than by taking a round trip to the dcache. This will cause a data hazard called a load-hit-store which basically empties out the CPU pipeline until the memory write has completed (usually around 40 cycles or so).

The exception to this is if the double values are already in memory and not in a register: in that case there is no danger of a load-hit-store. However your example indicates you've just calculated the double and returned it from a function which means it's likely to still be in XMM1.

Answer 8

Rather than testing and branching, I normally use this format for clamping:

clampedA = fmin(fmax(a,MY_MIN),MY_MAX);

Although I have never done any pref analysis on the compiled code.

Answer 9

I'm surprised nobody commented on Linasses solution : fmin(fmax(a,MY_MIN),MY_MAX)

It looks the most elegant/portable of all and shouldn't be slower than ternaries/branching ?