I have a .h file that implements my generic by defining LINK_* beforehand. This is an implementation of natural merge sort on doubly-linked lists using log n extra space by greedily merging as it receives the data. I also insert, (ie, insertion sort,) the first element in the merge and then merge from there; empirically, it is much faster on partially ordered lists, which is what I'm going for.I would like some input on the user-interface for the includer. Is it too confusing? What would make it better? Would you use it? There is a lot of pre-processor stuff that's going on; can I simplify a bit?LinkRun contains containing weakly increasing sub-lists; LinkRuns is a global array of the maximum possible LinkRuns, log n = sizeof(size_t) * 8; way too much, in most cases. Should I dynamically allocate, as in Natural merge sort in C? Each time this is included, a new LinkRuns static array is created, but concurrently calling the same instance of this function uses the same array and fails. I would like something that doesn't do that, but also doesn't waste a lot of space, is transparent to the caller, is fast, and is simple.Suggestions for algorithm improvement in general? Anywhere this could be faster?LinkSort.h:/** Generic sort of doubly linked-List. You must #define the required constants before including this file; they are undefined at the end of the this file for convenience when including multiple link types. @param LINK_NAME Name NameLinkSort() will be generated, only 1 word; required; @param LINK_TYPE struct Name because in C, a type does not have to be a single word; required; @param LINK_PREV prev @param LINK_NEXT next within the structure LINK_TYPE, previous and next fields; required. Does not do checks that they are valid; @param LINK_COMPARATOR LINK_NAMEComparator a comparator function taking two pointers to LINK_TYPE and returning an int; required. @author    Neil @version   1.0; 2016-11 @since     1.0; 2016-11 *//* check defines */#ifndef LINK_NAME#error Link generic LINK_NAME undefined.#endif#ifndef LINK_TYPE#error Link generic LINK_TYPE undefined.#endif#ifndef LINK_PREV#error Link generic LINK_PREV undefined.#endif#ifndef LINK_NEXT#error Link generic LINK_NEXT undefined.#endif#ifndef LINK_COMPARATOR#error Link generic LINK_COMPARATOR undefined.#endif/* After this block, the preprocessor replaces T with LINK_TYPE, T_(X) with LINK_NAMEX, PRIVATE_T_(X) with _LINK_NAME_X, and T_NAME with the string version. <p> http://c-faq.com/decl/namespace.html You may use identifiers consisting of an underscore followed by a digit or lower case letter at function, block, or prototype scope. So PRIVATE_T_, which makes use of PCAT_, is intended for prototype scope. <p> http://stackoverflow.com/questions/16522341/pseudo-generics-in-c */#ifdef CAT#undef CAT#endif#ifdef CAT_#undef CAT_#endif#ifdef PCAT#undef PCAT#endif#ifdef PCAT_#undef PCAT_#endif#ifdef T#undef T#endif#ifdef T_#undef T_#endif#ifdef PRIVATE_T_#undef PRIVATE_T_#endif#ifdef T_NAME#undef T_NAME#endif#ifdef QUOTE#undef QUOTE#endif#ifdef QUOTE_#undef QUOTE_#endif#define CAT_(x, y) x ## y#define CAT(x, y) CAT_(x, y)#define PCAT(x, y) PCAT_(x, y)#define PCAT_(x, y) _ ## x ## _ ## y#define QUOTE_(name) #name#define QUOTE(name) QUOTE_(name)#define T_(thing) CAT(LINK_NAME, thing)#define PRIVATE_T_(thing) PCAT(LINK_NAME, thing)#define T_NAME QUOTE(LINK_NAME)typedef LINK_TYPE PRIVATE_T_(Type);#define T PRIVATE_T_(Type)/* typedef a comparison function */typedef int (*T_(LinkCompare))(const T *, const T *);/* global comparison function */static const T_(LinkCompare) PRIVATE_T_(compare) = LINK_COMPARATOR;/* A run is a temporary sequence of values in the array that is weakly increasing. */struct PRIVATE_T_(LinkRun) {    T *head, *tail;    size_t size;};/* Store the maximum capacity for the indexing with size_t. (Overkill, really.) range(runs) = Sum_{k=0}^runs 2^{runs-k} - 1             = 2^{runs+1} - 2 2^bits      = 2 (r^runs - 1) runs        = log(2^{bits-1} + 1) / log 2 runs       <= 2^{bits - 1}, 2^{bits + 1} > 0 fixme: will crash if you call the same function multiple, simultaneous, times; it's okay to call different functions via different includes in the same file. There needs to be some management of resources for multi-threading, perhaps some dynamic allocation? */static struct PRIVATE_T_(LinkRuns) {    struct PRIVATE_T_(LinkRun) run[(sizeof(size_t) << 3) - 1];    size_t run_no;} PRIVATE_T_(runs);/* generated prototypes */void T_(LinkSort)(T **const);/** Inserts the first element from the larger of two sorted runs, then merges the rest. \cite{Peters2002Timsort}, via \cite{McIlroy1993Optimistic}, does long merges by galloping, but we don't have random access to the data. In practice, this is 2% slower on randomly distributed keys when the linked-list size is over 500 000; randomly distributed keys have high insertion times that to well in standard merging. However, it's (potentially much) faster when the keys have structure: observed, [-2%, 500%]. <p> Assumes array contains at least 2 elements and there are at least two runs. */static void PRIVATE_T_(natural_merge)(void) {    struct PRIVATE_T_(LinkRun) *const run_a = PRIVATE_T_(runs).run        + PRIVATE_T_(runs).run_no - 2;    struct PRIVATE_T_(LinkRun) *const run_b = run_a + 1;    T *a = run_a->tail, *b = run_b->head;    T *chosen;    if(run_a->size <= run_b->size) {        T *prev_chosen;        /* run a is smaller: downwards insert b.head followed by upwards         merge */        /* insert the first element of b downwards into a */        for( ; ; ) {            if(PRIVATE_T_(compare)(a, b) <= 0) {                chosen = a;                a = a->LINK_NEXT;                break;            }            if(!a->LINK_PREV) {                run_a->head = run_b->head;                chosen = b;                b = b->LINK_NEXT;                break;            }            a = a->LINK_PREV;        }        /* merge upwards, while the lists are interleaved */        while(chosen->LINK_NEXT) {            prev_chosen = chosen;            if(PRIVATE_T_(compare)(a, b) > 0) {                chosen = b;                b = b->LINK_NEXT;            } else {                chosen = a;                a = a->LINK_NEXT;            }            prev_chosen->LINK_NEXT = chosen;            chosen->LINK_PREV = prev_chosen;        }        /* splice the one list left */        if(!a) {            b->LINK_PREV = chosen;            chosen->LINK_NEXT = b;            run_a->tail = run_b->tail;        } else {            a->LINK_PREV = chosen;            chosen->LINK_NEXT = a;        }    } else {        T *next_chosen;        int is_a_tail = 0;        /* run b is smaller; upwards insert followed by downwards merge */        /* insert the last element of a upwards into b */        for( ; ; ) {            if(PRIVATE_T_(compare)(a, b) <= 0) {                chosen = b;                b = b->LINK_PREV;                break;            }            /* here, a > b */            if(!b->LINK_NEXT) {                is_a_tail = -1;                chosen = a;                a = a->LINK_PREV;                break;            }            b = b->LINK_NEXT;        }        if(!is_a_tail) run_a->tail = run_b->tail;        /* merge downwards, while the lists are interleaved */        while(chosen->LINK_PREV) {            next_chosen = chosen;            if(PRIVATE_T_(compare)(a, b) > 0) {                chosen = a;                a = a->LINK_PREV;            } else {                chosen = b;                b = b->LINK_PREV;            }            next_chosen->LINK_PREV = chosen;            chosen->LINK_NEXT = next_chosen;        }        /* splice the one list left */        if(!a) {            b->LINK_NEXT = chosen;            chosen->LINK_PREV = b;            run_a->head = run_b->head;        } else {            a->LINK_NEXT = chosen;            chosen->LINK_PREV = a;        }    }    run_a->size += run_b->size;    PRIVATE_T_(runs).run_no--;}/** Greedy natural insertion-merge sort on doubly-linked lists. @param phead   A pointer to the head of the list, which is a pointer to                LIST_TYPE; the head of the list will, in general, change,                unless it's the smallest item. */void T_(LinkSort)(T **const phead) {    T *head;    /* new_run is an index into link_runs, a temporary sorting structure;     head is first smallest, tail is last largest */    struct PRIVATE_T_(LinkRun) *new_run;    /* part of the state machine for classifying points wrt their neighbours */    enum { UNSURE, INCREASING, DECREASING } mono;    /* the data that we are sorting */    T *a, *b, *c, *first_iso_a;    /* run_count is different from link_runs.run_no in that it only increases;     only used for calculating the path up the tree */    size_t run_count, rc;    /* the value of the comparison */    int comp;    if(!phead || !(head = *phead)) return; /* ensure we have an 'a' */    /* reset the state machine and output to just 'a' in the first run */    mono = UNSURE;    PRIVATE_T_(runs).run_no = 1;    new_run = PRIVATE_T_(runs).run + 0, run_count = 1;    new_run->size = 1;    a = first_iso_a = new_run->head = new_run->tail = head;    for(b = a->LINK_NEXT; b; a = b, b = c) {        /* b.next can be modified, and we always want the iteration original */        c = b->LINK_NEXT;        comp = PRIVATE_T_(compare)(a, b);        /* state machine that considers runs in both directions -- in practice,         slightly slower than only considering increasing runs on most cases;         however, I would hate to see my code replaced with one line; reverse         order is 15 times faster, but it's not likely */        if(comp < 0) { /* a < b, increasing -- good */            if(mono != DECREASING) { /* if decreasing, inflection */                mono = INCREASING;                new_run->size++;                continue;            }        } else if(comp > 0) { /* decreasing; reverse preserving stability */            if(mono != INCREASING) { /* if increasing, inflection */                mono = DECREASING;                b->LINK_NEXT = first_iso_a;                first_iso_a->LINK_PREV = b;                new_run->head = first_iso_a = b;                new_run->size++;                continue;            }            new_run->tail = a; /* terminating an increasing sequence */        } else { /* a == b */            if(mono == DECREASING) { /* extend */                T *const a_next = a->LINK_NEXT;                b->LINK_NEXT = a_next;                a_next->LINK_PREV = b;                a->LINK_NEXT = b;                b->LINK_PREV = a;            } else { /* weakly increasing */                new_run->tail = b;            }            new_run->size++;            continue;        }        /* head and tail don't necessarily correspond to the first and last */        new_run->head->LINK_PREV = new_run->tail->LINK_NEXT = 0;        /* greedy merge: keeps space to O(log n) instead of O(n) */        for(rc = run_count; !(rc & 1) && PRIVATE_T_(runs).run_no >= 2; rc >>= 1)            PRIVATE_T_(natural_merge)();        /* reset the state machine and output to just 'b' at the next run */        mono = UNSURE;        new_run = PRIVATE_T_(runs).run + PRIVATE_T_(runs).run_no++, run_count++;        new_run->size = 1;        new_run->head = new_run->tail = first_iso_a = b;    }    /* terminating the last increasing sequence */    if(mono == INCREASING) new_run->tail = a;    new_run->tail->LINK_NEXT = new_run->head->LINK_PREV = 0;    /* clean up the rest; when only one run, propagate link_runs[0] to head */    while(PRIVATE_T_(runs).run_no > 1) PRIVATE_T_(natural_merge)();    *phead = PRIVATE_T_(runs).run[0].head;}#undef LINK_NAME#undef LINK_TYPE#undef LINK_PREV#undef LINK_NEXT#undef LINK_COMPARATORTest.c:/* This is a test for timing of linked-list sorts and LinkSort.h, a generic linked-list sorting header. Outputs concurrent.tsv, qsort.tsv, serial.tsv, sort.p (gnuplot script,) and, assuming you have gnuplot installed, sort.eps. @author    Neil @version   1.0; 2016-11 @since     1.0; 2016-11 */#include <stdlib.h> /* EXIT_SUCCESS */#include <stdio.h>  /* fprintf */#include <time.h>   /* clock */#include <math.h>   /* pow *//* constants */static const char *const gnu_name     = sort.p;static const char *const graph_name   = sort.eps;static const double impatient_ms      = 500.0;static const double sample_grow_power = 4.0;static const unsigned replicas        = 5;/* a data structure with two doubly-linked lists */struct Foo {    int x, y;    struct Foo *x_prev, *x_next;    struct Foo *y_prev, *y_next;};/* the first elements of one list */struct FooFirsts {    struct Foo *x_first, *y_first;};/* for qsort_keys */struct Key {    int key;    struct Foo *reference;};/* function type used in TimeData */typedef void (*SortFn)(struct FooFirsts *const);/* prototypes */static int Foo_compare_x(const struct Foo *, const struct Foo *);static int Foo_compare_y(const struct Foo *, const struct Foo *);/* generate FooXLinkSort */#define LINK_NAME FooX#define LINK_TYPE struct Foo#define LINK_PREV x_prev#define LINK_NEXT x_next#define LINK_COMPARATOR &Foo_compare_x#include LinkSort.h/* generate FooYLinkSort */#define LINK_NAME FooY#define LINK_TYPE struct Foo#define LINK_PREV y_prev#define LINK_NEXT y_next#define LINK_COMPARATOR &Foo_compare_y#include LinkSort.h/* scratch space */static struct Foo foo[1000000];static const unsigned foo_capacity = sizeof foo / sizeof *foo;/* random seed */static unsigned seed;/** @implements FooXLinkCompare */static int Foo_compare_x(const struct Foo *a, const struct Foo *b) {    return a->x - b->x;}/** @implements FooYLinkCompare */static int Foo_compare_y(const struct Foo *a, const struct Foo *b) {    return a->y - b->y;}/** @implements int (*)(const void *, const void *) */static int Key_compare(const void *a_key, const void *b_key) {    const struct Key *a = (struct Key *)a_key, *b = (struct Key *)b_key;    return a->key - b->key;}/** @implements SortFn */static void concurrent(struct FooFirsts *const this) {    #pragma omp parallel sections    {        #pragma omp section        FooXLinkSort(&this->x_first);        #pragma omp section        FooYLinkSort(&this->y_first);    }}/** @implements SortFn */static void serial(struct FooFirsts *const this) {    FooXLinkSort(&this->x_first);    FooYLinkSort(&this->y_first);}/** @implements SortFn */static void qsort_keys(struct FooFirsts *const this) {    struct Foo *a;    struct Key *keys = 0, *key;    size_t keys_size = 0, keys_capacity[2] = { 13, 21 }, i;    enum { E_NO_ERR, E_ERRNO } error = E_NO_ERR;    /* try-catch-finally */    do {        /* initial malloc */        if(!(keys = malloc(sizeof *keys * keys_capacity[0])))            { error = E_ERRNO; break; }        /* copy all the keys into keys */        for(a = this->x_first; a; a = a->x_next) {            /* make sure we have space; fibonacci some more if need be */            if(keys_size >= keys_capacity[0]) {                struct Key *new_keys;                keys_capacity[0] ^= keys_capacity[1];                keys_capacity[1] ^= keys_capacity[0];                keys_capacity[0] ^= keys_capacity[1];                keys_capacity[1] += keys_capacity[0];                if(keys_capacity[1] < keys_capacity[0])                    keys_capacity[1] = (size_t)-1;                if(!(new_keys = realloc(keys, sizeof *keys * keys_capacity[0])))                    { error = E_ERRNO; break; }                keys = new_keys;            }            /* copy the element */            key = keys + keys_size++;            key->key = a->x;            key->reference = a;        }        if(a) break; /* <- continued break from inner */        /* qsort the keys */        qsort(keys, keys_size, sizeof *keys, &Key_compare);        /* map the order onto the linked-list */        this->x_first = keys->reference;        for(i = 0; i < keys_size; i++) {            key = keys + i;            key->reference->x_prev =              (i) ? keys[i-1].reference : 0;            key->reference->x_next = (i!=keys_size-1) ? keys[i+1].reference : 0;        }        keys_size = 0;        /* now sort y */        for(a = this->y_first; a; a = a->y_next) {            /* copy the element */            key = keys + keys_size++;            key->key = a->y;            key->reference = a;        }        /* assert(keys_size == keys_size from before) */        /* qsort the keys */        qsort(keys, keys_size, sizeof *keys, &Key_compare);        /* map the order onto the linked-list */        this->y_first = keys->reference;        for(i = 0; i < keys_size; i++) {            key = keys + i;            key->reference->y_prev =              (i) ? keys[i-1].reference : 0;            key->reference->y_next = (i!=keys_size-1) ? keys[i+1].reference : 0;        }    } while(0);    if(error == E_ERRNO) {        perror(sort);        free(keys);        exit(EXIT_FAILURE);    } {        free(keys);    }}/** Sets up a random linked list of size elements and calls sort. @param size    Must be (0, foo_capacity). @return        The amount of time taken by sort. */static double time_sort_ms(const SortFn sort, const unsigned size) {    clock_t begin, end;    unsigned i;    struct Foo *a, *b;    struct FooFirsts firsts;    /* build up the linked list */    firsts.x_first = firsts.y_first = foo;    for(i = 0, a = foo; i < size; ++i, a++) {        a->x_prev =            (!i) ? 0 : foo + i - 1;        a->y_prev =            (!i) ? 0 : foo + i - 1;        a->x_next = (i == size - 1) ? 0 : foo + i + 1;        a->y_next = (i == size - 1) ? 0 : foo + i + 1;#ifdef NOISE        a->x = rand() / (RAND_MAX / 100.0f);        a->y = rand() / (RAND_MAX / 100.0f);#elif defined(NOISE_INCRESE)        a->x = rand() / (RAND_MAX / 100.0f) + i;        a->y = rand() / (RAND_MAX / 100.0f) + i;#else        a->x = size - i;        a->y = size - i;#endif    }    /* time */    begin = clock();    sort(&firsts);    end   = clock();    /* test that they're really in-order; assumes at least 1 datum */    for(a = firsts.x_first, b = a->x_next, i = 1; b; a = b, b = b->x_next, i++) {        if(a->x > b->x) {            fprintf(stderr,                Not x-sorted at [ ... , %d, %d, ... ], seed %u; :0.\n,                a->x, b->x, seed);            exit(EXIT_FAILURE);        }    }    if(i != size) {        fprintf(stderr, Input size %u, but output size %u on x.\n, size, i);        exit(EXIT_FAILURE);    }    for(a = firsts.y_first, b = a->y_next, i = 1; b; a = b, b = b->y_next, i++) {        if(a->y > b->y) {            fprintf(stderr,                Not y-sorted at [ ... , %d, %d, ... ], seed %u; :0.\n,                a->y, b->y, seed);            exit(EXIT_FAILURE);        }    }    if(i != size) {        fprintf(stderr, Input size %u, but output size %u on y.\n, size, i);        exit(EXIT_FAILURE);    }    return (double)(end - begin) / (CLOCKS_PER_SEC / 1000.0);}/** Does an experiment to see how new your computer is. */int main(void) {    /* look-up table */    struct TimeData {        char *const name;        SortFn sort;    } const time_data[] = {        { concurrent, &concurrent },        { serial,     &serial },        { qsort,      &qsort_keys }    }, *td;    const size_t time_data_size = sizeof time_data / sizeof *time_data;    unsigned td_i;    /* experiment */    FILE *fp = 0, *gnu = 0;    char fn[32] = not_a_file;    unsigned r, samples, base;    /* error handling */    enum { E_NO_ERR, E_ERRNO, E_EXTERNAL } error = E_NO_ERR;    /* seed */    srand(seed = (unsigned)clock()); rand();    /* try-catch-finally */    do {        if(!(gnu = fopen(gnu_name, w))) { error = E_ERRNO; break; }        fprintf(gnu, set term postscript eps enhanced color\n            set output \%s\\n            set grid\n            set xlabel \size of linked list\\n            set ylabel \time, t (ms)\\n            set yrange [0:]\n            # set xrange [0:1000] # zooming in\n            # seed %u\n            \n            plot, graph_name, seed); /* . . . */        for(td_i = 0; td = time_data + td_i, td_i < time_data_size; td_i++) {            /* open <sort>.tsv for writing */            if(snprintf(fn, sizeof fn, %s.tsv, td->name) < 0                || !(fp = fopen(fn, w))) { error = E_ERRNO; break; }            fprintf(stderr, Created/overwrote, \%s,\ to store data on %s                 sort.\n, fn, td->name);            /* do several experiments, increasing number of elements until we             hit impatient_ms or foo_capacity is reached */            fprintf(fp, # %s: elements ms\n, td->name);            for(samples = 1, base = 2;                samples < foo_capacity;                samples = (unsigned)pow((double)base++, sample_grow_power)) {                int n = 0;                double dt_ms = 0.0, mean_ms = 0.0, delta_ms, ssdm = 0.0;                fprintf(stderr, %s: sorting %u elements distributed uniformly                         at random.\n, td->name, samples);                /* \cite{Welford1962Note} */                for(r = 0; r < replicas; r++) {                    fprintf(stderr, #);                    dt_ms = time_sort_ms(td->sort, samples);                    n++;                    delta_ms = dt_ms - mean_ms;                    mean_ms += delta_ms / n;                    ssdm += delta_ms * (dt_ms - mean_ms);                }                fprintf(stderr,  done.\n);                fprintf(fp, %u\t%f\t%f\n, samples, mean_ms, sqrt(ssdm / (n - 1)));                /* loop until the process is taking too long */                if(mean_ms >= impatient_ms) break;            }            /* close the file */            if(fclose(fp)) { fp = 0; error = E_ERRNO; break; }            fp = 0;            /* write the graph */            fprintf(gnu,                %s\t\%s\ using 1:2:3 with errorlines lw 3 title \%s\,                td_i ? , \\\n :  , fn, td->name);        }    } while(0);    if(error == E_ERRNO) {        perror(fn);    } {        if(fp  && fclose(fp)) perror(fn);        if(gnu) {            fprintf(gnu, \n);            if(fclose(gnu)) perror(gnu_name);        }    }    if(error) return EXIT_FAILURE;    /* it doesn't matter at this point if we fail, but let's try gnuplot? and     opening it? */    error = E_NO_ERR;    do {        int s;        fprintf(stderr, Running gnuplot to get a graph of, \%s,\             (http://www.gnuplot.info/.)\n, gnu_name);        if(snprintf(fn, sizeof fn, gnuplot %s, gnu_name) < 0            || (s = system(fn)) == -1) { error = E_ERRNO; break; }        else if(s != EXIT_SUCCESS) { error = E_EXTERNAL; break; }        fprintf(stderr, Output, \%s;\ on a Mac, this opens automatically.\n,            graph_name);        if(snprintf(fn, sizeof fn, open %s, graph_name) < 0            || (s = system(fn)) == -1) { error = E_ERRNO; break; }        else if(s != EXIT_SUCCESS) { error = E_EXTERNAL; break; }    } while(0);    if(error == E_ERRNO) {        perror(fn);    } else if(error == E_EXTERNAL) {        fprintf(stderr, System call, \%s,\ returned an error; the machine             probably lacks that functionality.\n, fn);    }    return EXIT_SUCCESS;}These are the results of running this on a 2007 laptop which is 'simulated' multi-core. It is faster to sort noisy data by copying the keys to a temporary array, qsort, and copy the order back,2007 #define NOISEOn the other hand, if the data has structure, natural merge sort seems better,2007 #define NOISE_INCREASE2007 #define REVERSEOn a 2015 laptop, the concurrency is good,2015 #define NOISE
Generic natural merge sort a linked-list in C
c;sorting;linked list;generics;mergesort