11
\$\begingroup\$

Due to software constraints, I cannot use the standard libraries, <math.h>, <algorithm>, templates, inline, or Boost. I am also using standard C (ISO C99) such that array is not a reserved keyword like it is in Visual Studio.

Is this the "best" possible implementation of a min function? This needs to be robust and fast. Is there a more efficient C++ implementation?

double min(double* array, int size){
    // returns the minimum value of array
    static double val = array[0];
    for (int i = 1; i < size; ++i){
        val = val <= array[i] ? val : array[i];
    }
    return val;
}
\$\endgroup\$
  • \$\begingroup\$ Can you elaborate a bit on the software constraints you have to work within? I can understand not using boost but not allowing math.h, algorithms, and templates either? That's ruling out some very useful tools offered by the language as well as reinventing the proverbial wheel as shown by your sample. \$\endgroup\$ – greatwolf Oct 3 '11 at 18:34
  • \$\begingroup\$ @VictorT.,the short answer is due to certification issues and extreme V&V. Trust me, I wish I didn't have these constraints, but really appreciate the resources SE provides \$\endgroup\$ – Elpezmuerto Oct 3 '11 at 19:23
  • 4
    \$\begingroup\$ Is there a guarantee that size > 0? Otherwise you just stepped past the end of your array. \$\endgroup\$ – Martin York Oct 3 '11 at 20:22
  • \$\begingroup\$ @Tux: If this must be optimized for efficiency then the onus will be on the caller to provide only valid data, not the other way around. \$\endgroup\$ – Ed S. Oct 3 '11 at 20:34
  • \$\begingroup\$ @Ed S.: Quote: This needs to be robust \$\endgroup\$ – Martin York Oct 3 '11 at 20:42
19
\$\begingroup\$

It's not quite the best. There's some things IMHO that is hindering its performance and usefulness.

There's no point in declaring val as a static variable. In fact, you've killed any chance of it being usable in multi-threaded programs.

The body of the loop is performing an assignment in every single iteration when it doesn't need to. If you want it to be its best, you should only be doing so when it is required.

Your overall structure is fine, assuming we can expect well-formed inputs where array is nonempty and size is positive, I'd just change it so it's more like this:

double min(const double *arr, size_t length) {
    // returns the minimum value of array
    size_t i;
    double minimum = arr[0];
    for (i = 1; i < length; ++i) {
        if (minimum > arr[i]) {
            minimum = arr[i];
        }
    }
    return minimum;
}
\$\endgroup\$
  • \$\begingroup\$ I'm pretty sure any half-decent optimising compiler will remove the redundant assignment in the original code's ternary operator. \$\endgroup\$ – Rafe Oct 3 '11 at 23:26
  • \$\begingroup\$ Shouldn't it be i++ in stead of ++i? Since you're pre-incrementing and starting at 1, in the first iteration i would already be 2 at the if statement. \$\endgroup\$ – Decent Dabbler Oct 4 '11 at 6:05
  • \$\begingroup\$ @fireeyedboy: It doesn't really matter which one is used, in the context of a for loop, they are exactly the same, increment i. The incrementing expression is evaluated between the execution of the body and the condition check and has no other effect on the two. \$\endgroup\$ – Jeff Mercado Oct 4 '11 at 6:12
  • 1
    \$\begingroup\$ @S.L.Barth: It's the C side in me moving the declaration to the beginning of the scope. I tend to think in (old) C more than I do C++ and the relaxed restrictions. No reason at all really. \$\endgroup\$ – Jeff Mercado Oct 4 '11 at 7:53
  • 1
    \$\begingroup\$ @Elpezmuerto: By declaring the type of the array const double *, we're stating that the contents of the array could not be changed. i.e., The array is immutable. We can still read the contents and store it on a non-const double as long as we don't modify the contents of the array. \$\endgroup\$ – Jeff Mercado Oct 4 '11 at 19:24
6
\$\begingroup\$

You can use SSE2/SSE3 minps / pminsd or relevant instruction set for your processor/architectoure since it is supported directly in GCC / MASM / TASM (In case MASM or TASM is not supported such SSE2/SSE3 instruction set there are also the .inc files to create macros simulating instruction sets on the web for MASM), create .OBJ file by your favorite linker then link it as usual and use in you favorite IDE. You will get from 4x to 16x performance boost compared to the traditional "classic" algorithm. It depend on data size (old compilers treats double not in IEEE format, bout like float in several configurations, on 16x systems, particularly, double means 32 bit data structure, not 64 bit data structure, in modern languages it is correlated to "double" and "long double" data structures, respectively)

The idea is simple: if you have k elements, [k=4n+p, 4>p=>0], complete it with n-p elements or just load last 4 doubles resetting to 0 last p elements, so you can fast evaluate n candidates. evaluate candidates n times comparing to the accumulator, you will get a minimum.

If your processor supports SSE5 or is a brand new, most likely you also will be using one of the HD instructions, which really handy, because it can find maximum (not minimum yet) in array of double values.

Sample of using SSE to calculate peak values of a float array:

#include <xmmintrin.h>

double min(double* array, int size) { 
    // returns the minimum value of array 
    int i; 
    double val = array[0]; 
    for (i = 1; i < size; ++i) { 
        if (val > array[i]) { 
            val = array[i]; 
        } 
    } 
    return val; 
} 

#define ARRAY_SIZE 16

float m_fArray[ARRAY_SIZE];

void x86_sse_find_peaks(float *buf, unsigned nframes, float *min, float *max)
{
    __m128 current_max, current_min, work;

    // Load max and min values into all four slots of the XMM registers
    current_min = _mm_set1_ps(*min);
    current_max = _mm_set1_ps(*max);

    // Work input until "buf" reaches 16 byte alignment
    while ( ((unsigned long)buf) % 16 != 0 && nframes > 0) {

        // Load the next float into the work buffer
        work = _mm_set1_ps(*buf);

        current_min = _mm_min_ps(current_min, work);
        current_max = _mm_max_ps(current_max, work);

        buf++;
        nframes--;
    }

    // use 64 byte prefetch for quadruple quads
    while (nframes >= 16) {
        //__builtin_prefetch(buf+64,0,0); // for GCC 4.3.2+

        work = _mm_load_ps(buf);
        current_min = _mm_min_ps(current_min, work);
        current_max = _mm_max_ps(current_max, work);
        buf+=4;
        work = _mm_load_ps(buf);
        current_min = _mm_min_ps(current_min, work);
        current_max = _mm_max_ps(current_max, work);
        buf+=4;
        work = _mm_load_ps(buf);
        current_min = _mm_min_ps(current_min, work);
        current_max = _mm_max_ps(current_max, work);
        buf+=4;
        work = _mm_load_ps(buf);
        current_min = _mm_min_ps(current_min, work);
        current_max = _mm_max_ps(current_max, work);
        buf+=4;
        nframes-=16;
    }

    // work through aligned buffers
    while (nframes >= 4) {

        work = _mm_load_ps(buf);

        current_min = _mm_min_ps(current_min, work);
        current_max = _mm_max_ps(current_max, work);

        buf+=4;
        nframes-=4;
    }

    // work through the rest < 4 samples
    while ( nframes > 0) {

        // Load the next float into the work buffer
        work = _mm_set1_ps(*buf);

        current_min = _mm_min_ps(current_min, work);
        current_max = _mm_max_ps(current_max, work);

        buf++;
        nframes--;
    }

    // Find min & max value in current_max through shuffle tricks

    work = current_min;
    work = _mm_shuffle_ps(work, work, _MM_SHUFFLE(2, 3, 0, 1));
    work = _mm_min_ps (work, current_min);
    current_min = work;
    work = _mm_shuffle_ps(work, work, _MM_SHUFFLE(1, 0, 3, 2));
    work = _mm_min_ps (work, current_min);

    _mm_store_ss(min, work);

    work = current_max;
    work = _mm_shuffle_ps(work, work, _MM_SHUFFLE(2, 3, 0, 1));
    work = _mm_max_ps (work, current_max);
    current_max = work;
    work = _mm_shuffle_ps(work, work, _MM_SHUFFLE(1, 0, 3, 2));
    work = _mm_max_ps (work, current_max);

    _mm_store_ss(max, work);
}

int _tmain(int argc, _TCHAR* argv[])
{
    float   min = FLT_MAX;
    float   max = FLT_MIN;

    m_fArray[0] = 0;
    m_fArray[1] = 1;
    m_fArray[2] = 2;
    m_fArray[3] = 3;
    m_fArray[4] = 4;
    m_fArray[5] = 3;
    m_fArray[6] = 2;
    m_fArray[7] = 1;
    m_fArray[8] = -1;
    m_fArray[9] = -2;
    m_fArray[10] = -3;
    m_fArray[11] = -4;
    m_fArray[12] = -5;
    m_fArray[13] = -6;
    m_fArray[14] = -7;
    m_fArray[15] = -8;

    x86_sse_find_peaks(m_fArray, ARRAY_SIZE, &min, &max);

    printf("value = %.2f, max = %.2f\n", min, max); // output is: value = -8.00, max = 4.00
    return 0;
}
\$\endgroup\$
  • 1
    \$\begingroup\$ This is great peaks implementation but I am working with double and don't want to downcast. If I recall, SSE was originally designed for floats. I've read that SSE can now handle double precision values, but not sure how to change your above code to handle doubles. \$\endgroup\$ – Elpezmuerto Oct 5 '11 at 14:09
  • \$\begingroup\$ @Elpezmuerto: I will definitely make it a try, would you mind that i will use Visual Studio 11 to compile C++ code, I'm also having a Borland C++, TASM in Windows XP virtual machine under Windows 7? \$\endgroup\$ – Artur Mustafin Oct 6 '11 at 1:53
  • \$\begingroup\$ sure I am happy with the code current state. Anything beyond this is really for research curiosity and the good for all, etc. \$\endgroup\$ – Elpezmuerto Oct 6 '11 at 18:14
3
\$\begingroup\$

Recommended form for robustness:

Ok, you want robust and fast. If this was needed in production code, I would probably write it like this:

double min(double array[], int size)
{
    assert(array != NULL);
    assert(size >= 0);

    if ((array == NULL) || (size <= 0))
       return NAN;

    double val = -DBL_MAX;
    for (int i = 0; i < size; i++)
        if (val < array[i])
            val = array[i];
    return val;
}

Looping: I noticed you optimized your loop to begin with val = array[0] initially and i = 1 instead of i = 0. This may have saved machine cycles 20 years ago, but these days memory accesses are the biggest bottleneck, and you're stuck with the same number of memory accesses into array[] no matter how you write the loop, so I would try to be as non-clever as possible here and just stick to a more traditional loop.

Argument checking: Is it valid to call this function with size of 0? If so, then you also have a bug if you start off with val = array[0] when size is 0, because you might be reading invalid memory depending on where array is pointing. At the very least, it's logically incorrect to examine array[0] when size is 0.

Alternative formulation (not recommended):

Just for fun, if you have a sufficiently large stack and lots of machine cycles to burn, you could write it like this:

double min(double a[], int n)
{
    if ((a == NULL) || (n <= 0))
    {
        return NAN;
    }
    else if (n == 1)
    {
        return array[0];
    }
    else
    {
        double x = a[0], y = min(&a[1], n-1);
        return x < y? x : y;
    }
}

Or, a recursive binary search using only O(log₂ n) stack space:

inline double min2(double x, double y)
{
    return x < y? x : y;
}

double min(double a[], int n)
{
    return (a == NULL) || (n <= 0)? NAN
         : (n == 1)? a[0]
         : min2(min(a, n/2), min(a+n/2, n-n/2));
}

Or, in avoidance of error-checking and in adherence to the principle of maximum terseness:

inline double min2(double x, double y) { return x<y? x:y; }
double min(double a[], int n) { return n==1? a[0] : min2(min(a,n/2),min(a+n/2,n-n/2)); }

;-)

\$\endgroup\$

protected by Jamal Feb 14 '17 at 18:03

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

Not the answer you're looking for? Browse other questions tagged or ask your own question.