Min / Max function of 1D array in C / C++

Question

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

Answer 1

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

Answer 2

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

Answer 3

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

;-)