12
\$\begingroup\$

Here's my first program in C, which is not my first language.

I have been reading Modern C and I attempted this challenge on page 26. The prior pages are pretty much the only exposure I have had to C, for reference. I have a rough understanding of pointers but I have read no formal material on the subject.

#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>

void sort(int num_elements, double array[num_elements]) {
    // the input is trivially sorted
    if (num_elements == 1) {
        return;
    }

    // (from Googling) integer division yields an integer result
    int num_elements_A = num_elements / 2;
    int num_elements_B = num_elements - num_elements_A;

    // copy elements from given array into two smaller sub-arrays
    double A[num_elements_A];
    double B[num_elements_B];
    for (int i = 0; i < num_elements_A; ++i) {
        A[i] = array[i];
    }
    for (int i = 0; i < num_elements_B; ++i) {
        B[i] = array[i+num_elements_A];
    }

    // recursively sort the two small sub-arrays
    sort(num_elements_A, A);
    sort(num_elements_B, B);

    // merge the two sorted sub-arrays
    int current_A_elem = 0;
    int current_B_elem = 0;
    for (int i = 0; i < num_elements; ++i) {
        if ( current_A_elem == num_elements_A ) {
            // we have exhausted subarray A
            array[i] = B[current_B_elem];
            ++current_B_elem;
        } else if ( current_B_elem == num_elements_B ) {
            // we have exhausted subarray B
            array[i] = A[current_A_elem];
            ++current_A_elem;
        } else if ( A[current_A_elem] < B[current_B_elem] ) {
            array[i] = A[current_A_elem];
            ++current_A_elem;
        } else {
            array[i] = B[current_B_elem];
            ++current_B_elem;
        }
    }
}

// returns true if the input is sorted
bool self_test(int num_elements, double array[num_elements]) {
    for (int i = 0; i < num_elements-1; ++i) {
        if (array[i] > array[i+1] ) {
            return false;
        }
    }
    return true;
}

int main(int argc, char* argv[argc+1]) {
    printf("you want me to sort an array of length %d?\n", argc-1);
    
    if (argc < 2) {
        printf("you must provide at least one double to be sorted!\n");
        return EXIT_FAILURE;
    }

    // allocate some memory to store an appropriate number of double values
    double A[argc-1];

    for (int i = 1; i < argc; ++i) {              // process args
        double const a = strtod(argv[i], 0);  // arg -> double
        A[i-1] = a;
    }

    // merge sort!
    sort(argc-1, A);

    printf("here is your sorted array:\n");
    for (int i = 0; i < argc - 1; ++i) {
        printf("%f\n", A[i]);
    }

    if (self_test(argc-1, A)) {
        printf("self-test passed\n");
    } else {
        printf("self-test failed\n");
    }

    return EXIT_SUCCESS;
}
\$\endgroup\$
3
  • \$\begingroup\$ I would replace for (int i = 0; i < num_elements_A; ++i) { A[i] = array[i]; } with memcpy(A, array, num_elements_A*sizeof(double)); and for (int i = 0; i < num_elements_B; ++i) { B[i] = array[i+num_elements_A]; } with memcpy(B, array+num_elements_A, num_elements_B*sizeof(double)); \$\endgroup\$ Feb 23, 2023 at 0:07
  • \$\begingroup\$ @user16217248 Rather than num_elements_A*sizeof(double), sizeof A is simpler. \$\endgroup\$ Feb 23, 2023 at 3:01
  • \$\begingroup\$ Just to help understand why memcpy is faster than loop : this is generally compiled in single instruction of the processor. \$\endgroup\$
    – Ptit Xav
    Mar 5, 2023 at 11:59

4 Answers 4

10
\$\begingroup\$

Handle size 0

// if (num_elements == 1) {
if (num_elements <= 1) {
    return;
}

Sorting double should handle not-a-number

Research isgreater(), isgreaterequal(), ... to well compare FP when NANs are involved.

IMO, all NANs should get sorted to the end of the array.

Minor: Use memcpy() when able

//for (int i = 0; i < num_elements_A; ++i) {
//    A[i] = array[i];
//}
memcpy(A, array, sizeof A);

Consider size_t to handle arrays of all sizes

// void sort(int num_elements, double array[num_elements]) {
void sort(size_t num_elements, double array[num_elements]) {

Note: This is a minor issue as long as temporary space is local to the function. For me, once n == 100,000, I ran out of stack space.

Make cleaner faster code that executes 2 loops

Reduce index testing. Something like:

size_t left = 0;
size_t right = 0;
size_t i = 0;
while (1) {
  if (fcmpd(A[left], B[right]) < 0) {
    array[i++] = A[left++];
    if (left >= nmemb_left) {
      memcpy(&array[i], &B[right], sizeof array[0] * (nmemb_right - right));
      break;
    }
  } else {
    array[i++] = B[right++];
    if (right >= nmemb_right) {
      memcpy(&array[i], &A[left], sizeof array[0] * (nmemb_left - left));
      break;
    }
  }
}

Advanced: Sorting double issues: +0, -0, NAN

If code wants all -0.0 < +0.0, more tests needed. Research signbit().

If code wants to sort NANs, consider memcmp().

General

A generic sort would pass in the compare function rather than embed that test within the sort function.

\$\endgroup\$
2
  • \$\begingroup\$ It seems to me that num_elements < 1 is an impossible state. What is the reasoning for adding this logic? My guess: in more complex examples, it may not be obvious that num_elements >= 1 always, so making the assertion num_elements > 1 at the beginning of the subroutine is helpful to the reader. \$\endgroup\$ Feb 23, 2023 at 15:55
  • 2
    \$\begingroup\$ @user14464173 As used in this main(), sort(0, array) is not going to happen. In general use, sort(0, array) is possible. Plan for the future. If staying with int num_elements, if (num_elements <= 1) also prevents UB when num_elements < 0 - something code analyzers may whine about too otherwise. \$\endgroup\$ Feb 23, 2023 at 16:35
8
\$\begingroup\$

int main(int argc, char* argv[argc+1])

You don't need to put a size in the argv array. The compiler ignores it.

You should #include <stdio.h> for the printf definition.

double const a = strtod(argv[i], 0); // arg -> double

strtod can fail if you give it non-numeric data. Use the 2nd argument to detect failures and alert the user. See this SO question

You can collapse the if-else branches in your merge to

for (int i = 0; i < num_elements; ++i) {
    if (( current_A_elem < num_elements_A ) &&
        ((current_B_elem == num_elements_B) ||
         (A[current_A_elem] < B[current_B_elem]))) {
        array[i] = A[current_A_elem++];
    }
    else {
        array[i] = B[current_B_elem++];
    }
}

although I don't know if that's clearer. Another way is to do 3 loops:

while ( current_A_elem < num_elements_A  && current_B_elem < num_elements_B) {
    if (A[current_A_elem] < B[current_B_elem]) {
        array[i++] = A[current_A_elem++];
    }
    else {
        array[i++] = B[current_B_elem++];
    }
}
// only one of the following 2 while loops will execute.
while ( current_A_elem < num_elements_A) {
    array[i++] = A[current_A_elem++];
}
while ( current_B_elem < num_elements_B) {
    array[i++] = B[current_B_elem++];
}

double A[num_elements_A];

As @user16217248 mentions, you could use memcpy to populate the working arrays. However, declaring these arrays on the stack has the potential of a StackOverflow. The memory needed grows O(N^2): 10 items needs storage space for 100; 1000 items need space for 1,000,000 etc.

It might be a better design to pass in an auxiliary array that can get reused.

Change the function signature to void sort(int n, double array[], double work_array[]) . Inside the function, instead of copying to local storage, copy array to work_array and then do the recursive part like:

double* work_A = &(work_array[0]);
double* work_B = &(work_array[num_elements_A]);
// note that the roles of `array` and `work_array` are swapped.
sort(num_elements_A, work_A, array);
sort(num_elements_B, work_B, &(array[num_elements_A]));

Afterwords, use work_A and work_B where you had A and B.

\$\endgroup\$
3
  • 5
    \$\begingroup\$ Code analyzers can use int main(int argc, char* argv[argc+1]) to assess proper argv[] access. \$\endgroup\$ Feb 23, 2023 at 2:27
  • 1
    \$\begingroup\$ I do #include <stdio.h>! I'm going to leave the size of argv in, because of the suggestion by @chux-ReinstateMonica and because my book does this too. Thank you for the link to the explanation of strtod, I have added rudimentary input validation as shown. I like your implementation with 3 loops! Finally, as for the suggestion about array copying, I'm going to leave that issue until I have a better understanding of pointers. \$\endgroup\$ Feb 23, 2023 at 15:50
  • \$\begingroup\$ collapse ... given i < num_elements, current_A_elem < num_elements_A && A[current_A_elem] < B[current_B_elem] || num_elements_B <= current_B_elem should be enough. \$\endgroup\$
    – greybeard
    Mar 5, 2023 at 11:20
5
\$\begingroup\$

To well investigate a function like OP's sort(), it is useful to incorporate a test harness.

Interestingly a test harness is often more code than the code-under-test, yet does not need the robustness of production code. It is OK for it to have all sorts of assert()s, slow tests, etc.

Below is some sample code using a modified version of OP's sort.

It allocates space rather than use local function space.

#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>

int fcmpd(double a, double b) {
  if (isgreater(a, b))
    return 1;
  if (isless(a, b))
    return -1;
  if (isnan(a)) {
    if (isnan(b)) {
      return memcmp(&a, &b, sizeof a);
    }
    return 1;
  }
  if (isnan(b)) {
    return -1;
  }
  // same value, but sort out -0.0 from +0.0
  return !signbit(a) - !signbit(b);
}

int fcmpdqs(const void *ap, const void *bp) {
  double a = *(const double*) ap;
  double b = *(const double*) bp;
  return fcmpd(a, b);
}

double *dt = 0;

void sort_double_helper(size_t nmemb, double array[nmemb], double *temp) {
  if (nmemb > 1) {
    size_t nmemb_left = nmemb / 2;
    size_t nmemb_right = nmemb - nmemb_left;

    // Copy elements from given array into 2 smaller (but adjacent) sub-arrays.
    memcpy(temp, array, sizeof temp[0] * nmemb);
    double *A = temp;
    temp += nmemb_left;
    double *B = temp;
    temp += nmemb_right;

    if (temp > dt) {
      dt = temp;
    }

    // recursively sort the two small sub-arrays
    sort_double_helper(nmemb_left, A, temp);
    sort_double_helper(nmemb_right, B, temp);

    // merge the two sorted sub-arrays
    size_t i_left = 0;
    size_t i_right = 0;
    size_t i = 0;
    while (1) {
      if (fcmpd(A[i_left], B[i_right]) < 0) {
        array[i++] = A[i_left++];
        if (i_left >= nmemb_left) {
          memcpy(&array[i], &B[i_right],
              sizeof array[0] * (nmemb_right - i_right));
          break;
        }
      } else {
        array[i++] = B[i_right++];
        if (i_right >= nmemb_right) {
          memcpy(&array[i], &A[i_left],
              sizeof array[0] * (nmemb_left - i_left));
          break;
        }
      }
    }
  }
}

void sort_double(size_t nmemb, double array[nmemb]) {
  size_t temp_n = 2 * nmemb + 256; // WAG (needs further analysis)
  double *temp = malloc(sizeof temp[0] * temp_n);
  dt = temp;
  assert(temp);
  sort_double_helper(nmemb, array, temp);
  // printf("TEMP %zu %td %g\n", nmemb, dt - temp, 1.0 * (dt - temp) - 2 * nmemb);
  free(temp);
}

///////////
#include <stdint.h>

double rand_double(void) {
  union {
    double d;
    unsigned char uc[sizeof(double)];
  } u;
  for (unsigned i = 0; i < sizeof u.uc; i++) {
    u.uc[i] = (unsigned char) rand();
  }
  return u.d;
}

void print_double_array(size_t nmemb, const double a[nmemb]) {
  const char *sep = "";
  for (size_t i = 0; i < nmemb; i++) {
    if (i <= 7 || i + 7 >= nmemb) {
      printf("%s%- 13.5g", sep, a[i]);
      sep = ", ";
    } else {
      ;
    }
  }
  puts("");
}

unsigned long long checkcode_double_array(size_t nmemb, const double a[nmemb]) {
  unsigned long long cc = 0;
  while (nmemb > 0) {
    union {
      double d;
      uint64_t u64;
    } u = {.d = a[--nmemb]};
    cc ^= u.u64;
  }
  return cc;
}

int sort_test(size_t n) {
  int error = 0;
  double *a = malloc((sizeof a[0] * n) | 1);
  assert(a);
  long double sum = 0.0;
  for (size_t i = 0; i < n; i++) {
    a[i] = rand_double();
    sum += a[i];
  }
  double *q = malloc((sizeof q[0] * n) | 1);
  assert(q);
  memcpy(q, a, sizeof q[0] * n);
  qsort(q, n, sizeof q[0], fcmpdqs);
  unsigned long long cc0 = checkcode_double_array(n, a);
  print_double_array(n, a);
  sort_double(n, a);
  print_double_array(n, a);
  unsigned long long cc1 = checkcode_double_array(n, a);
  // Simplistic checks
  if (cc0 != cc1) {
    printf("Mis-match! %016llX  %016llX\n", cc0, cc1);
    error = 1;
  }
  for (size_t i = 1; i < n; i++) {
    if (isgreater(a[i - 1], a[i])) {
      printf("%zu %g %g\n", i, a[i - 1], a[i]);
      error = 2;
      break;
    }
  }
  // Complete check
  if (memcmp(a, q, sizeof a[0] * n)) {
    print_double_array(n, q);
    printf("qsort() mis-match!");
    error = 3;
  }
  printf("Result for test size: %zu, error code: %d\n", n, error);
  fflush(stdout);
  free(a);
  free(q);
  return error;
}

int main(void) {
  // srand((unsigned)time(0));
  sort_test(0);
  for (size_t i = 1; i < INT_MAX / 10; i = i * 10) {
    sort_test(i);
  }
}

Output

Result for test size: 0, error code: 0
-1.5574e+118 
-1.5574e+118 
Result for test size: 1, error code: 0
-3.0612e+125 ,  8.0501e+175 , -1.0608e+132 ,  3.3927e-209 , -2.3682e-244 ,  7.9119e+301 , -1.4051e-53  , -1.379e-14   , -2.1587e+125 , -6.4137e+13  
-1.0608e+132 , -3.0612e+125 , -2.1587e+125 , -6.4137e+13  , -1.379e-14   , -1.4051e-53  , -2.3682e-244 ,  3.3927e-209 ,  8.0501e+175 ,  7.9119e+301 
Result for test size: 10, error code: 0
 4.9009e-109 , -1.1004e-267 , -4.9987e-84  ,  1.2384e-196 , -4.7458e+240 ,  7.877e+259  , -1.1082e+184 ,  3.0209e+145 ,  2.7375e+85  ,  1.1423e+195 ,  3.563e+245  ,  5.7227e+244 ,  3.5208e+156 , -6.4709e-167 ,  1.2248e+260 
-4.0149e+257 , -4.7458e+240 , -5.6195e+235 , -4.1516e+224 , -7.0454e+202 , -3.7949e+202 , -3.5251e+186 , -1.1082e+184 ,  5.7227e+244 ,  3.563e+245  ,  5.1718e+246 ,  1.3249e+258 ,  7.877e+259  ,  1.2248e+260 ,  1.3056e+301 
Result for test size: 100, error code: 0
-1.4878e-238 , -6.3365e+224 , -1.2271e+167 , -3.8502e+248 , -4.2007e+26  ,  3.2845e+250 , -3.0233e+132 , -1.061e+150  ,  1.6955e+57  , -5.7906e-248 ,  7.7691e+236 ,  4.512e-67   , -1.8815e-69  ,  1.9407e-21  ,  7.0348e+282 
-4.3422e+307 , -1.643e+307  , -1.4944e+306 , -1.4327e+305 , -4.4285e+302 , -7.8393e+301 , -6.982e+301  , -6.8221e+301 ,  1.4194e+296 ,  3.5929e+296 ,  1.7932e+297 ,  4.2531e+297 ,  1.3644e+299 ,  8.9723e+302 ,  8.7216e+307 
Result for test size: 1000, error code: 0
-6.0499e-245 , -9.7001e+171 ,  8.5045e-135 ,  2.0394e+170 ,  2.3656e+114 ,  3.6892e+135 ,  7.5027e+247 , -6.1797e-267 , -7.203e-130  ,  2.2239e+293 ,  7.363e-282  ,  9.7721e-87  ,  3.9824e-59  , -9.2624e-226 ,  8.4479e-217 
-1.388e+308  , -8.6533e+307 , -7.3786e+307 , -7.3449e+307 , -5.8173e+307 , -4.9323e+307 , -1.6201e+307 , -1.2287e+307 , -nan         ,  nan         ,  nan         ,  nan         ,  nan         ,  nan         ,  nan         
Result for test size: 10000, error code: 0
-9.5033e+149 ,  1.8574e-124 , -1.3543e+34  ,  3.112e-197  , -4.7614e-219 ,  4.0548e-85  , -8.8576e-166 , -2.2164e-176 , -1.0475e-76  , -2.906e-93   , -1.2964e+166 ,  3.056e-55   ,  1.6323e-74  ,  2.801e-142  ,  3.1944e-251 
-1.7492e+308 , -1.7355e+308 , -1.7194e+308 , -1.7152e+308 , -1.6971e+308 , -1.658e+308  , -1.617e+308  , -1.5684e+308 ,  nan         , -nan         ,  nan         , -nan         , -nan         ,  nan         , -nan         
Result for test size: 100000, error code: 0
 8.146e+87   , -7.0071e+129 , -9.1492e+68  ,  1.2535e+50  , -2.0746e-34  , -2.9128e-303 ,  3.9061e+302 ,  5.97e-127   , -2.5632e-192 ,  3.0825e-251 , -3.2951e-41  , -3.339e+164  ,  1.9077e-19  ,  1.1938e+250 ,  2.9245e-159 
-1.7973e+308 , -1.797e+308  , -1.7712e+308 , -1.7701e+308 , -1.7674e+308 , -1.7592e+308 , -1.7531e+308 , -1.7507e+308 , -nan         , -nan         ,  nan         ,  nan         , -nan         ,  nan         , -nan         
Result for test size: 1000000, error code: 0
-1.7589e+258 , -5.3367e+118 ,  2.6179e-155 ,  2.7818e-28  , -2.0408e-257 ,  2.5443e-193 ,  1.5116e+219 , -3.481e+103  ,  2.3648e+297 ,  5.3499e-173 ,  6.6864e-42  , -7.7392e-62  ,  5.2705e+212 , -1.6984e-199 , -1.705e-244  
-1.7976e+308 , -1.7974e+308 , -1.7969e+308 , -1.7964e+308 , -1.7961e+308 , -1.796e+308  , -1.7954e+308 , -1.7951e+308 ,  nan         , -nan         , -nan         , -nan         ,  nan         ,  nan         , -nan         
Result for test size: 10000000, error code: 0

\$\endgroup\$
4
  • \$\begingroup\$ Why post this as a separate answer, instead of as an edit to your existing one? \$\endgroup\$ Feb 24, 2023 at 8:50
  • \$\begingroup\$ @CodyGray Other answer was stable, this one volatile and focused on whole other aspect: testing and memory allocation. As stated "test harness is often more code than the code-under-test", it seemed prudent to post this code-heavy answer by itself to allow OP & others to quickly use it to test their alternatives. \$\endgroup\$ Feb 24, 2023 at 17:11
  • \$\begingroup\$ Could you expand a little bit on what the results imply about the performance of sort()? \$\endgroup\$ Feb 24, 2023 at 22:32
  • \$\begingroup\$ @user14464173 The 2 lines of numbers are values (up to the first 8 and last 7) of the array before and then after sorting. The error code 0 is good, else there was a detected problem. No problem even with sizes up to 10,000,000. \$\endgroup\$ Feb 24, 2023 at 22:36
3
\$\begingroup\$

This is a Great First Program in C

I don’t think I’ve seen a better one. However, I’ll suggest some improvements.

Write Contracts for your Functions

These can be short and sweet; names like sort are quite self-explanatory. You already have a good one for the only other function in the program. I find that it helps me to write these before I implement them.

Refactor Into Helper Functions

If you’re setting out to write a recursive solution, you might as well go all the way and replace all the loops with recursion! But, seriously, this is one of the many times that the most efficient recursive implementation has a different interface than the one you want to make public. (I’ll have more to say about this later.)

If you’re taking my advice to write static single assignments, it’s also handy to have to refactor your computations that require a lot of local temporary variables, especially mutable ones, into phi functions (or “phony functions”). This simplifies the code that calls the computation by moving all the temporary local state into the phi function. The only data the caller needs to keep around is the result, frozen as a const so it will not be modified again. The calling function becomes simpler and safer.

A lot of programmers also swear by refactoring their functions into smaller ones with no more than one level of recursion. If you give the building blocks good names, the code will often come out close to a high-level description of what it does, and be easier to understand. That’s not always a good idea, but I do think it’s good coding style in general to keep local variables in scope only for as long as they’re needed. This reduces the number of things to keep track of at once.

Array Indices are Never, Ever int!

On many systems, int is a signed 32-bit type that will overflow at a little more than two billion elements—which can easily happen if you’re crunching big data. On some older or embedded systems, int will wrap around after a mere 32,767 elements!

Worse, signed integer overflow is formally undefined behavior, which many compilers take as a license to break your program and introduce security bugs.

An array index or offset should have type size_t. If you want to use a signed quantity instead, so that the difference of two pointers or indices will compare properly as negative values, use ptrdiff_t. They’re both in stddef.h.

Use Static Single Assignments Where Possible

Or, in other words, declare everything as const that you can. This is a program where it’s not always possible, since you’re sorting the array in place.

Writing static single assignments eliminate whole categories of bugs, such as losing track of what the current value of a variable is, using a variable before it is set or on a code path where it didn’t get set, or modifying it from the wrong place.

Static single assignments also help the optimizer out, especially when it’s trying to convince itself that it doesn’t need to recalculate variables inside a loop, or at all.

So, for example, num_elements_A and num_elements_B should be set once and never modified. If you did, it was an mistake, and you want the compiler to stop you. So you can declare them both const.

Your Merge Loop is a Great Start

I mean the one that begins

for (int i = 0; i < num_elements; ++i) {

I especially approve of how you use the invariant that one element is written to the result array on each iteration, and how you check that you always fill the buffer without overflowing it.

There are still some ways it can be improved. One, I just mentioned, is that i should be a size_t, not int. chux also had a good suggestion to check for Not-a-Number (which make conditionals behave oddly).

You should consider checking for logic errors in the program with assertions. Here, for example, whenever the loop completes, you know that both current_A_element and current_B_element should both be at the end of their arrays. If they aren’t, your algorithm has gone wrong, and the closer to the source of the error you So, you could assert that both these things are true.

If you don’t think this is necessary, they’re never used outside the loop, so you could instead declare them as:

for ( size_t i = 0, left = 0, right = 0; i < num_elements; ++i ) {

This gives them the same lifetime as i: they stick around until the loop terminates, then are destroyed.

Now, because I’m weird, I went ahead and wrote a tail-recursive function to replace the loop. I personally like those because they turn the local state of the loop body into static single assignments that all get updated together, at the same time, when function recurses. It’s impossible to update any of it in any other place, and impossible not to update any of it. But, like I said, that’s weird. That’s how things are done in functional languages, and C isn’t really designed for it.

Besides, you already know how to do it with a for loop. What you did is much more normal for C code.

Test that the Output is Correct, not Merely that it is Sorted

For example, the test array[i] > array[i+1] will always fail if one of the two elements is a NaN value, so a NaN in the middle of the output would never be detected.

You could in theory work around this (first rewriting your check for ordering to fail if it encounters a non-NaN value after a NaN, then also checking whether the array of results is a permutation of the original), but that’s very complicated and error-prone.

What test frameworks actually do is have you specify the expected result for each test case, and compare that to the result from the library.

Beware of Variable-Length Arrays

In your function, you write

double A[num_elements_A];
double B[num_elements_B];

Since num_elements_A and num_elements_B are not compile-time constants, these are variable-length arrays. Those are very poorly supported! They’ve formally been in the standard since ’99, but that doesn’t mean much. Some compilers don’t allow them at all, and others say they’re deprecated.

If you need a dynamic buffer, your best option is calloc.

However, you should actually

Allocate Buffers Once and Switch Between Them

Right now, you make two calls to the heap per iteration, and then copy the entire contents of the array to the two subarrays. These are not freed until the function returns, so you end up keeping lg n copies of the data in memory at once. These allocations also won’t have good locality of reference, so you’ll have a lot of cache misses.

It would be far more efficient if you allocated a second buffer, at the beginning of the function, and then ping-pong between the original buffer and the second buffer with each step of the algorithm.

That is, if the input data is currently in its original location, you want to partition it into two subarrays and store the two sorted subarrays in the secondary buffer.

Then, when you partition the two subarrays into four subarrays, in order to sort them, you can overwrite the original input data with the four sorted subarrays. This is fine, because you no longer need the original data to compute the final output, only the two sorted subarrays, which you will compute from the four sorted subarrays once you have those.

To compute the four sorted subarrays, you want to partition them into eight subarrays and sort those. You can overwrite the secondary buffer with the eight sorted subarrays first; after you have used those to compute the four sorted subarrays in the primary buffer, you will then use the four sorted subarrays to write two sorted subarrays to the secondary buffer, which you will finally use to write the sorted array back to its original location. When that is done, you free the secondary buffer.

In keeping with my advice to refactor into helper functions, I turned this into two helper functions, both of which receive pointers to both buffers. Both divide the input array, call the other version on both subarrays, and merge the results into the appropriate location.

Doing things this way also eliminates the extra copy of the entire array on each iteration. You just write each intermediate result into the correct place to begin with.

To implement this cleanly, it helps immensely to use the right data structures.

Use Array Slices

A struct containing a pointer to the start of a subarray, and the number of elements it contains, simplifies your life in a number of ways. The main one is that you can now write a helper function that returns a slice. This lets you compose your helper functions together, passing the output of one function to the next.

You also always keep track of what size your buffers are supposed to be, allowing you to prevent and test for buffer overruns.

Putting it All Together

What I came up with arguably has too many fiddly little helpers, but each one is short and simple, and the resulting code is very efficient. You can check it out on Godbolt.

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#if __clang__ || __INTEL_LLVM_COMPILER 
#  define MUSTTAIL __attribute((musttail))
#else
#  define MUSTTAIL /**/
#endif

/* Represents a slice of an array of double, whose base pointer is p and
 * whose length is n.
 */
typedef struct slice_double_t {
  double* p;
  size_t n;
} slice_double_t;


/* Tail-recursive helper function to merge the two sorted arrays at left_ptr
 * and right_ptr to the array out_ptr.  Checks that it wrote exactly as many
 * elements as expected.
 */
static void merge_double( const double* const left_ptr,
                          const double* const left_end,
                          const double* const right_ptr,
                          const double* const right_end,
                          double* const out_ptr,
                          const double* const out_end )
{
  if ( left_ptr == left_end && right_ptr == right_end ) {
    assert(out_ptr == out_end); // Always check for buffer overruns!
    return;
  } else if (left_ptr == left_end) { // Only the left subarray is empty.
    assert( right_end - right_ptr == out_end - out_ptr );
    memcpy( out_ptr, right_ptr, sizeof(*out_ptr)*(size_t)(out_end - out_ptr) );
    return;
  } else if (right_ptr == right_end) { // Only the right subarray is empty.
    assert( left_end - left_ptr == out_end - out_ptr );
    memcpy( out_ptr, left_ptr, sizeof(*out_ptr)*(size_t)(out_end - out_ptr) );
    return;
  } else if (*left_ptr <= *right_ptr) { // Never true for NaN!
    *out_ptr = *left_ptr;
    MUSTTAIL return merge_double( left_ptr+1, left_end,
                                  right_ptr, right_end,
                                  out_ptr+1, out_end );
  } else if (*left_ptr > *right_ptr) { // Also never true for NaN.
    *out_ptr = *right_ptr;
    MUSTTAIL return merge_double( left_ptr, left_end,
                                  right_ptr+1, right_end,
                                  out_ptr+1, out_end );
  } else if (isnan(*right_ptr)) {
    *out_ptr = *left_ptr;
    MUSTTAIL return merge_double( left_ptr+1, left_end,
                                  right_ptr, right_end,
                                  out_ptr+1, out_end );
  } else {
    *out_ptr = *right_ptr;
    MUSTTAIL return merge_double( left_ptr, left_end,
                                  right_ptr+1, right_end,
                                  out_ptr+1, out_end );
  }
}


/* Wrapper for merge_double.
 */
static slice_double_t merge_slices_double( const slice_double_t left,
                                           const slice_double_t right,
                                           const slice_double_t output )
{
  assert(left.n + right.n == output.n);
  merge_double( left.p,
                left.p + left.n,
                right.p,
                right.p + right.n,
                output.p,
                output.p + output.n );
  return output;
}

// Forward declaration for mutually-recursive call.
static slice_double_t merge_sort_in_place_slice_double( const slice_double_t inputs,
                                                        const slice_double_t scratchpad );


/* Sorts the input slice, storing the result in the output slice, which must
 * be the same length.  May overwrite the inputs.  Returns the sorted slice.
 */
static slice_double_t merge_sort_to_slice_double( const slice_double_t inputs,
                                                  const slice_double_t outputs )
{
  assert (inputs.n == outputs.n);

  if (inputs.n == 0) {
    return outputs;
  }

  if (inputs.n == 1) {
    outputs.p[0] = inputs.p[0];
    return outputs;
  }

  // Divide:
  const size_t pivot = inputs.n / 2U;
  const slice_double_t left_in = { inputs.p, pivot };
  const slice_double_t right_in = { inputs.p + pivot, inputs.n - pivot };
  const slice_double_t left_out = { outputs.p, pivot };
  const slice_double_t right_out = { outputs.p + pivot,
                                     outputs.n - pivot };

  // Conquer:
  return merge_slices_double( merge_sort_in_place_slice_double( left_in, left_out ),
                              merge_sort_in_place_slice_double( right_in, right_out ),
                              outputs );
}


/* Sorts the input slice in place.  Requires another array slice the same size
 * as the inputs, to hold intermediate results.  Returns the sorted slice.
 */
static slice_double_t merge_sort_in_place_slice_double( const slice_double_t inputs,
                                                        const slice_double_t scratchpad )
{
  assert(inputs.n == scratchpad.n);

  if (inputs.n <= 1) {
    return inputs;
  }
  // When this line is reached, sorting is non-trivial.

  // Divide:
  const size_t pivot = inputs.n / 2U;
  const slice_double_t left_in = { inputs.p, pivot };
  const slice_double_t right_in = { inputs.p + pivot, inputs.n - pivot };
  const slice_double_t left_scratch = { scratchpad.p, pivot };
  const slice_double_t right_scratch = { scratchpad.p + pivot,
                                         scratchpad.n - pivot };

  // Conquer:
  return merge_slices_double( merge_sort_to_slice_double( left_in, left_scratch ),
                              merge_sort_to_slice_double( right_in, right_scratch ),
                              inputs );
}


/* Merge-sorts the input array slice in place.  Returns the sorted slice.
 */
slice_double_t merge_sort_slice_double(const slice_double_t inputs)
{
  if (inputs.n <= 1) { // Input this small is ordered, by definition.
    return inputs;
  }

  // Could also get the scratchpad from alloca().
  const slice_double_t scratchpad = { calloc( inputs.n, sizeof(*inputs.p) ),
                                      inputs.n };
  assert(scratchpad.p); // Check for allocation failure.
  (void)merge_sort_in_place_slice_double( inputs, scratchpad );
  free(scratchpad.p);
  return inputs;
}


/* Merge-sorts the input array of length n in place, by calling
 * merge_sort_double_slice.
 */
void merge_sort_double( const size_t n, double data[n] )
{
  const slice_double_t inputs = { data, n };
  (void)merge_sort_slice_double(inputs);
}

There’s one piece of boilerplate here that isn’t standard (and maybe I should’ve kept simpler). I wrote a tail-recursive version of merge_double, and C is not really designed to do that. For it to work, I #define MUSTTAIL to expand to a non-standard extension on Clang or ICX, and other compilers ignore it. (Therefore, for this code to work on GCC, you need to give it one of the flags -O2, -O3, -Os or -foptimize-sibling-calls.)

If you’d rather stick to a for loop, which you could still move into a helper function, you don’t need to worry about any of that.

You should also write some kind of test harness for the code:

// Test Harness:
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

/* Returns true if and only if the input arrays (of the same length) are
 * identical.  Checks for NAN.
 */
bool compare_doubles( const double* const p,
                      const double* const q,
                      const size_t n )
{
  for ( size_t i = 0; i < n; ++i ) {
    
    if ( (p[i] != q[i]) && (!isnan(p[i]) || !isnan(q[i])) ) {
      return false;
    }
  }

  return true;
}


/* Prints an array of doubles to the file in the format [1, 2, 3].
 */
void fprint_doubles( FILE* const file,
                     const double* const inputs,
                     const size_t n )
{
  fputs( "[", file );

  if (n > 0) {
    fprintf( file, "%f", inputs[0] );
  }

  for ( size_t i = 1; i < n; ++i ) {
    fprintf( file, ", %f", inputs[i] );
  }

  fputs( "]", file );
}


/* Either prints that the test succeeded to stdout and returns, or prints
 * information about what failed to stderr and aborts.
 */
void expect_test( const double* const expected,
                  const double* const got,
                  const size_t n,
                  const unsigned test_no )
{
  if (!compare_doubles( expected, got, n )) {
    fflush(stdout);
    fprintf( stderr, "Test %u FAILED!\nExpected: ", test_no );
    fprint_doubles( stderr, expected, n );
    fputs( "\nGot: ", stderr );
    fprint_doubles( stderr, got, n );
    exit(EXIT_FAILURE);
  }

  printf( "Test %u passed.\n", test_no );
}


int main(void)
{
  unsigned test_no = 1;

  {
    static double test_array1[15] = { 11, 12, 13, 14, 15,
                                       6,  7,  8,  9, 10,
                                       1,  2,  3,  4,  5 };
    static const size_t test_array1_size = sizeof(test_array1) / sizeof(*test_array1);
    static const double expected_array1[15] = {  1,  2,  3,  4,  5,
                                                 6,  7,  8,  9, 10,
                                                11, 12, 13, 14, 15 };

    merge_sort_double( test_array1_size, test_array1 );
    expect_test( expected_array1, test_array1, test_array1_size, test_no++ );
  }

  {
    static double test_array2[16] = { 13, 14, 15, 16,
                                       9, 10, 11, 12,
                                       5,  6,  7,  8,
                                       1,  2,  3,  4 };
    static const size_t test_array2_size = sizeof(test_array2)/sizeof(*test_array2);
    static const double expected_array2[16] = {  1,  2,  3,  4,
                                                 5,  6,  7,  8,
                                                 9, 10, 11, 12,
                                                13, 14, 15, 16 };

    merge_sort_double( test_array2_size, test_array2 );
    expect_test( expected_array2, test_array2, test_array2_size, test_no++ );
  }

  {
    static double test_array3[6] = { 2, NAN, 1, NAN, 3, NAN };
    static const size_t test_array3_size = sizeof(test_array3)/sizeof(*test_array3);
    static const double expected_array3[6] = { 1, 2, 3, NAN, NAN, NAN };
  
    merge_sort_double( test_array3_size, test_array3 );
    expect_test( expected_array3, test_array3, test_array3_size, test_no++ );
  }

  return EXIT_SUCCESS;
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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