libucontext_fp_test.c

/* standalone test case for https://github.com/kaniini/libucontext/issues/2

 # compile libucontext:
 git clone --depth 1 https://github.com/kaniini/libucontext.git
 cd libucontext
 make FREESTANDING=yes libucontext.a

 # and then compile this test case:
 cc -Wall -Wextra -Wshadow -Os -o libucontext_fp_test libucontext_fp_test.c -I${HOME}/Downloads/libucontext/include/ ${HOME}/Downloads/libucontext/libucontext.a
 */

#ifdef USE_SYSTEM_UCONTEXT
#define _XOPEN_SOURCE
#include <ucontext.h>
#else
#include <libucontext/libucontext.h>
#define ucontext_t libucontext_ucontext_t
#define getcontext libucontext_getcontext
#define swapcontext libucontext_swapcontext
#define makecontext libucontext_makecontext
#define setcontext libucontext_setcontext
#endif

#include <stddef.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>

/* begin very limited coroutine impl using ucontext */

struct channel {
    /* points to the ucontext_t that should be swapcontext'd to on the next parent-to-child or child-to-parent context switch */
    ucontext_t * inactive;

    void (* func)(struct channel *, void *);
    void * arg;

    void * block;
};

static void springboard(struct channel * channel) {
    /* run the main body of the child coroutine */
    channel->func(channel, channel->arg);

    /* reached when coroutine function finishes. the parent may be watching for this condition */
    channel->func = NULL;

    fprintf(stderr, "%s: child body function returned, calling setcontext...\n", __func__);

    /* explicitly jump out of this forever to wherever the parent is waiting (can't use uc_link due to order of operations) */
    setcontext(channel->inactive);
    __builtin_unreachable();
}

static void swapcontext_symmetric(struct channel * channel) {
    ucontext_t buf, * tmp = channel->inactive;
    channel->inactive = &buf;
    swapcontext(&buf, tmp);
}

#define STACK_ALIGNMENT 64
#define BLOCKSIZE 524288

struct channel * coroutine_create(void (* func)(struct channel *, void *), void * arg) {
    /* allocate space for child stack and struct channel */
    unsigned char * block = NULL;
    if (posix_memalign((void *)&block, STACK_ALIGNMENT, BLOCKSIZE)) abort();

    /* offset of initial top of stack within the block, and also the struct channel, respecting alignment reqs */
    const size_t stack_size = (BLOCKSIZE - sizeof(struct channel)) & ~(STACK_ALIGNMENT - 1);

    struct channel * channel = (void *)(block + stack_size);
    *channel = (struct channel) { .func = func, .arg = arg, .block = block };

    /* start the child coroutine */
    ucontext_t dest;
    fprintf(stderr, "%s: calling getcontext...\n", __func__);
    getcontext(&dest);
    dest.uc_stack.ss_sp = block;
    dest.uc_stack.ss_size = stack_size;
    dest.uc_stack.ss_flags = 0;
    fprintf(stderr, "%s: calling makecontext...\n", __func__);
    makecontext(&dest, (void (*)())springboard, 1, channel);

    channel->inactive = &dest;
    fprintf(stderr, "%s: calling swapcontext_symmetric...\n", __func__);
    swapcontext_symmetric(channel);
    fprintf(stderr, "%s: returned from swapcontext_symmetric\n", __func__);
    /* control flow returns here via the first coroutine_switch in the child */

    return channel;
}

void coroutine_switch(struct channel * channel) {
    if (channel->func) swapcontext_symmetric(channel);
}

void close_and_join(struct channel * channel) {
    while (channel->func) swapcontext_symmetric(channel);
    free(channel->block);
}

/* end implementation of coroutines, begin test case code */

#include <string.h>
#include <math.h>
#include <complex.h>

static void fft8_with_intermission(struct channel * peer, float complex y[8], const float complex x[8], const char in_child_for_stderr_only) {
    /* perform four dfts of size 2, two of which are multiplied by a twiddle factor (a -90 degree phase shift) */
    const float complex a0 = x[0] + x[4];
    const float complex a1 = x[0] - x[4];
    const float complex a2 = x[2] + x[6];
    const float complex a3 = CMPLXF( cimagf(x[2]) - cimagf(x[6]), crealf(x[6]) - crealf(x[2]) );
    const float complex a4 = x[1] + x[5];
    const float complex a5 = x[1] - x[5];
    const float complex a6 = x[3] + x[7];
    const float complex a7 = CMPLXF( cimagf(x[3]) - cimagf(x[7]), crealf(x[7]) - crealf(x[3]) );

    /* perform two more dfts of size 2 */
    const float complex c0 = a0 + a2;
    const float complex c1 = a1 + a3;
    const float complex c2 = a0 - a2;
    const float complex c3 = a1 - a3;
    const float complex c4 = a4 + a6;
    const float complex b5 = a5 + a7;
    const float complex b6 = a4 - a6;
    const float complex b7 = a5 - a7;

    /* intermission */
    fprintf(stderr, "%s(%s): calling coroutine_switch...\n", __func__, in_child_for_stderr_only ? "child" : "parent");
    coroutine_switch(peer);
    fprintf(stderr, "%s(%s): returned from coroutine_switch\n", __func__, in_child_for_stderr_only ? "child" : "parent");

    /* apply final twiddle factors */
    const float complex c5 = CMPLXF((cimagf(b5) + crealf(b5)) * (float)M_SQRT1_2,  (cimagf(b5) - crealf(b5)) * (float)M_SQRT1_2);
    const float complex c6 = CMPLXF( cimagf(b6), -crealf(b6) );
    const float complex c7 = CMPLXF((cimagf(b7) - crealf(b7)) * (float)M_SQRT1_2, -(crealf(b7) + cimagf(b7)) * (float)M_SQRT1_2);

    /* intermission */
    fprintf(stderr, "%s(%s): calling coroutine_switch...\n", __func__, in_child_for_stderr_only ? "child" : "parent");
    coroutine_switch(peer);
    fprintf(stderr, "%s(%s): returned from coroutine_switch\n", __func__, in_child_for_stderr_only ? "child" : "parent");

    /* perform four dfts of length two */
    y[0] = c0 + c4;
    y[1] = c1 + c5;
    y[2] = c2 + c6;
    y[3] = c3 + c7;
    y[4] = c0 - c4;
    y[5] = c1 - c5;
    y[6] = c2 - c6;
    y[7] = c3 - c7;
}

int validate(const float complex y[], const float complex ref[], const size_t Y) {
    float error_squared_sum = 0, denominator = 0;
    for (size_t iy = 0; iy < Y; iy++) {
        const float complex error = ref[iy] - y[iy];
        error_squared_sum += crealf(error) * crealf(error) + cimagf(error) * cimagf(error);
        denominator += crealf(ref[iy]) * crealf(ref[iy]) + cimagf(ref[iy]) * cimagf(ref[iy]);
    }
    return error_squared_sum * 1e3f < denominator;
}

void child(struct channel * channel_to_parent, void * arg) {
    (void)arg;

    float complex y[8], x[8] = { 1, I, -1, -I, 1, I, -1, -I };
    fft8_with_intermission(channel_to_parent, y, x, 1);

    for (size_t ix = 0; ix < 8; ix++)
        fprintf(stderr, "%s: y[%zu] = %g %+gi\n", __func__, ix, crealf(y[ix]), cimagf(y[ix]));
    fprintf(stderr, "\n");

    /* communicate result to parent */
    int * child_passed_p = arg;
    *child_passed_p = validate(y, (float complex[8]) { [2] = 8.0f }, 8);

    fprintf(stderr, "%s: child returning...\n", __func__);
}

int main(void) {
    int child_passed = 0;
    struct channel * channel_to_child = coroutine_create(child, &child_passed);
    fprintf(stderr, "%s: returned from coroutine_create\n", __func__);

    float complex y[8], x[8] = { 0.25f, 0.25f, 1.25f, 0.25f, 0.25f, 0.25f, 0.25f, 0.25f };
    fft8_with_intermission(channel_to_child, y, x, 0);

    fprintf(stderr, "%s: calling close_and_join...\n", __func__);
    close_and_join(channel_to_child);
    fprintf(stderr, "%s: returned from close_and_join\n", __func__);

    for (size_t ix = 0; ix < 8; ix++)
        fprintf(stderr, "%s: y[%zu] = %g %+gi\n", __func__, ix, crealf(y[ix]), cimagf(y[ix]));
    fprintf(stderr, "\n");

    const int parent_passed = validate(y, (float complex[8]) { 3.0f, -I, -1, I, 1, -I, -1, I }, 8);

    printf("%s\n", child_passed && parent_passed ? "pass" : "fail");

    return (child_passed && parent_passed) ? EXIT_SUCCESS : EXIT_FAILURE;
}