Smoking Clover. This code is ported from http://sourceforge.net/projects/smokingclover/

clover.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <GLUT/glut.h> 

#define KEY_ESC 27
#define WINDOW_POS_X 100
#define WINDOW_POS_Y 100
#define DEF_SCREEN_WIDTH 640
#define DEF_SCREEN_HEIGHT 360
#define RADIUS 5432
#define N_COLORS 256
#define UPDATE_PERIOD 16 

#define MIN(x,y) ((x) > (y) ? (y) : (x))
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define  CEIL(a, b) ((a)+(b)-1)/(b)
#define	FLOOR(a, b) CEIL((a)-(b), (b))

struct colors {
    double r, g, b;
};

struct cstate {
    double cur;
    double vel;
};

static double rnd(double lo, double hi);
static void clear_colors(void);
static void rotate_colors(int steps);
static void next_color(int steps);
static int get_count(int x, int y);
static void put_count(int x, int y, int val);
static void plot(int x, int y);
static void redraw_lines (void);
static void redraw_dup(void);
static void clipLine(int x0, int y0, int xn, int yn, int xe, int ye, int xf, int yf);
static void line(int fun, int x0, int y0, int dx, int dy, int xe, int ye, int xf, int yf);

static int curr_width = 0;
static int curr_height = 0;
static GLubyte *screen_buf = NULL;
static int fullscreen = 0;
static clock_t c_old = 0;

static int *count_buf = NULL;
static struct colors colors[N_COLORS];
static struct colors colors_tmp[N_COLORS];

static void
keyboard(unsigned char c, int x, int y)
{
    switch (c) {
    case KEY_ESC:
        exit(0);
        break;
    case 'f':
        if (fullscreen) {
            glutReshapeWindow(DEF_SCREEN_WIDTH, DEF_SCREEN_HEIGHT);
            glutPositionWindow(WINDOW_POS_X, WINDOW_POS_Y);
            fullscreen = 0;
        } else {
            glutFullScreen();
            fullscreen = 1;
        }
        break;
    }
}

void
idle()
{
    glutPostRedisplay();
}

static void
reshape(int w, int h)
{
    curr_width = w;
    curr_height = h;

    if (screen_buf) {
        free(screen_buf);
    }
    screen_buf = (GLubyte *)malloc(curr_width * curr_height * 4);
    if (!screen_buf) {
        abort();
    }

    if (count_buf) {
        free(count_buf);
    }
    count_buf = (int *)calloc(curr_width * curr_height, sizeof(int));
    if (!count_buf) {
        abort();
    }

    redraw_lines();
    redraw_dup();
    clear_colors();
}

static void
draw_pixel(int x, int y, GLubyte r, GLubyte g, GLubyte b)
{
    int base = 4 * ((curr_height - 1 - y) * curr_width + x);
    screen_buf[base + 0] = r;
    screen_buf[base + 1] = g;
    screen_buf[base + 2] = b;
    screen_buf[base + 3] = 255;
}

static void
display()
{
    clock_t c_new;
    int x, y;

    glClear(GL_COLOR_BUFFER_BIT);

    for (y = 0; y < curr_height; y++) {
        for (x = 0; x < curr_width; x++) {
            struct colors c = colors[get_count(x, y) % 256];
            draw_pixel(x, y, 255 * c.r, 255 * c.g, 255 * c.b);
        }
    }

    c_new = clock();
    if (UPDATE_PERIOD < (1000 * (c_new - c_old) / CLOCKS_PER_SEC)) {
        c_old = c_new;
        next_color(1);
    }

    glWindowPos2i(0, 0);
    glDrawPixels(curr_width, curr_height, GL_RGBA, GL_UNSIGNED_BYTE, screen_buf);
    glutSwapBuffers();
}

int
main(int argc, char** argv)
{
    srand(0);

    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH );
    glutInitWindowSize(DEF_SCREEN_WIDTH, DEF_SCREEN_HEIGHT);
    glutInitWindowPosition(WINDOW_POS_X, WINDOW_POS_Y);
    glutCreateWindow(argv[0]);
    glutKeyboardFunc(keyboard);
    glutIdleFunc(idle);
    glutReshapeFunc(reshape);
    glutDisplayFunc(display);

    glClearColor(0.0, 0.0, 0.0, 0.0);
    glShadeModel(GL_FLAT);

    glutMainLoop();
    return 0;
}

double
rnd(double lo, double hi)
{
    static char state[8];
    static int beenhere;
    char *old;
    int seed;
    double r;

    if (beenhere == 0) {
        beenhere = 1;

        seed = 1;
        old = initstate (seed, state, sizeof state);
        setstate (old);
    }

    old = setstate (state);
    r = random ();
    setstate (old);

    return ((r / RAND_MAX) * (hi - lo) + lo);
}

void
clear_colors ()
{
    int i;

    for (i = 0; i < N_COLORS; i++) {
        colors[i].r = 0;
        colors[i].g = 0;
        colors[i].b = 0;
    }
}

void
rotate_colors (int steps)
{
    if (steps >= N_COLORS) {
        return;
    }

    memcpy(&colors_tmp[0], &colors[steps],
            (N_COLORS - steps) * sizeof colors_tmp[0]);
    memcpy(&colors_tmp[N_COLORS - steps], &colors[0],
            steps * sizeof colors_tmp[0]);
    memcpy(colors, colors_tmp, N_COLORS * sizeof *colors);
}

void
next_color (int steps)
{
    static struct cstate cstates[3];
    static int beenhere;
    int i;
    struct cstate *cp;
    int step;
    double maxvel, minvel;

    minvel = .003;
    maxvel = .15;
    if (beenhere == 0) {
        beenhere = 1;
        cstates[0].vel = rnd(minvel, maxvel);
        cstates[1].vel = rnd(minvel, maxvel);
        cstates[2].vel = rnd(minvel, maxvel);
    }

    for (step = 0; step < steps; step++) {
        for (i = 0, cp = cstates; i < 3; i++, cp++) {
            cp->cur += cp->vel;
            if (cp->cur > 1) {
                cp->cur = 1;
                cp->vel = rnd(-maxvel, -minvel);
            }
            if (cp->cur < 0) {
                cp->cur = 0;
                cp->vel = rnd(minvel, maxvel);
            }
        }

        rotate_colors(1);
        colors[N_COLORS - 1].r = cstates[0].cur;
        colors[N_COLORS - 1].g = cstates[1].cur;
        colors[N_COLORS - 1].b = cstates[2].cur;
    }
}

void
redraw_lines(void)
{
    int maxX = curr_width - 1;
    int maxY = curr_height - 1;
    int midX = maxX / 2;
    int midY = maxY / 2;
    int line_x = RADIUS;
    int line_y = 0;
    int line_f = 0;

    for (;;) {
        if (line_f > line_x) {
            line_x--;
            line_f = line_f - line_x - (line_x - 1);
        }

        clipLine(
                midX, midY,
                line_x + midX, line_y + midY,
                0, 0,
                maxX, maxY);

        line_f = line_f + line_y + line_y + 1;
        line_y++;
        if (line_y >= line_x) {
            break;
        }
    }
}

void
redraw_dup(void)
{
    int maxX = curr_width - 1;
    int maxY = curr_height - 1;
    int midX = maxX / 2;
    int midY = maxY / 2;
    int dup_x = midX;
    int nsteps, step;

    for (;;) {
        nsteps = 4;
        for (step = 0; step < nsteps; step++) {
            int yy, y, val, x1, y1, i, o;

            /* set vals on diagonal to 2*v-1 */
            if (dup_x - midX + midY <= maxY)
                put_count(dup_x, dup_x-midX+midY,
                        (get_count(dup_x, dup_x-midX+midY)<<1)-1);

            /* now do a column from horizontal, down to diag */
            yy = MIN(maxY, dup_x - midX + midY);
            for (y = midY; y <= yy; y++) {
                val = get_count(dup_x, y);
                x1 = dup_x;
                y1 = y;
                for (i = 0; i < 4; i++) {
                    if ((y1 < maxY) && (y1 > 0)) {
                        put_count(midX + midX - x1, y1, val);
                        put_count(x1, y1, val);
                    }
                    o = x1;
                    x1 = midX + midY - y1;
                    y1 = midY + o - midX;
                }
            }

            dup_x++;
            if (dup_x >= maxX) {
                break;
            }
        }

        if (dup_x >= maxX) {
            break;
        }
    }
}

/*
 * (xe, ye) and (xf, yf) are the corners of a rectangle to clip a line to.
 * (x0, y0) and (xn, yn) are the endpoints of the line to clip.
 * The function argument that's being computed is the semi-quadrant;
 *  dx and dy are used to determine whether we're above or below the diagonal,
 *  since (x0, y0) is always the midpoint of the pattern.
 * (The LispM has the origin at lower left, instead of upper left, so
 * the numbers don't correspond to the normal Cartesian plane quadrants.)
 *
 * This routine is very general, but the calling code only builds lines in the
 * first semi-quadrant and then copies them everywhere else.
 */

void
clipLine(int x0, int y0, int xn, int yn, int xe, int ye, int xf, int yf)
{
    int	dx, dy;

    dx = abs(xn - x0);
    dy = abs(yn - y0);

    if (xn > x0) {			/* moving right */
        if (yn >= y0) {			/* moving up */
            if (dx > dy)		/* below diagonal */
                line(0, x0, y0, dx, dy, xe, ye, xf, yf);
            else
                line(1, y0, x0, dy, dx, ye, xe, yf, xf);
        } else {
            if (dx > dy)
                line(7, x0, -y0, dx, dy, xe, -yf, xf, -ye);
            else
                line(6, -y0, x0, dy, dx, -yf, xe, -ye, xf);
        }
    } else {
        if (yn >= y0) {
            if (dx > dy)
                line(3, -x0, y0, dx, dy, -xf, ye, -xe, yf);
            else
                line(2, y0, -x0, dy, dx, ye, -xf, yf, -xe);
        } else {
            if (dx > dy)
                line(4, -x0, -y0, dx, dy, -xf, -yf, -xe, -ye);
            else
                line(5, -y0, -x0, dy, dx, -yf, -xf, -ye, -xe);
        }
    }
}

/*
 * Clip symmetric segment (x0, y0) thru (xn, yn) to the rectangle
 * (xe, ye) < (xf, yf).
 *
 * The original says:
 *
 * "This routine incorrectly assumes that the subsegment starts prior to the
 * midpoint of the supersegment.  The 'divide for nearest integer' (i.e.,
 * divide for remainder of minimum magnitude), which is simulated by the FLOOR
 * and CEIL of num and (dx <<1), always rounds up on the half integer case, but
 * should round down (for symmetry) if startup is in 2nd half. It would be
 * nice to have these other flavors of divide.'
 */

void
line(int fun, int x0, int y0, int dx, int dy, int xe, int ye, int xf, int yf)
{
    int	x, num, lx;
    int	xx, y, x00, f;
    int	x11;

    x = MAX(x0,  MAX(xe,
                (dy == 0) ? xe :
                x0 + CEIL(dx * (((ye - y0)<<1) - 1), (dy << 1))));
    num = dx + 2 * dy * (x - x0);
    lx = MIN(xf, (dy == 0) ? xf :
            x0 + CEIL(dx * (((yf - y0)<<1) - 1), (dy << 1)));
    xx = MIN(lx, x0 + (dx>>1));
    y = y0 + FLOOR(num, (dx<<1));
    f = (FLOOR(num, (dx<<1)) - dx) >> 1;

    for (x00 = x; x00 < xx; x00++, f += dy) {
        if (f + f > dx) {
            f -= dx;
            y++;
        }
        switch(fun) {
        case 0:	plot(x00, y);	break;
        case 1:	plot(y, x00);	break;
        case 2:	plot(-y, x00);	break;
        case 3:	plot(-x00, y);	break;
        case 4:	plot(-x00, -y);	break;
        case 5:	plot(-y, -x00);	break;
        case 6:	plot(y, -x00);	break;
        case 7:	plot(x00, -y);	break;
        }
    }

    for (x11 = x00; x11 < lx; x11++, f += dy) {
        if (f + f > dx) {
            f -= dx;
            y++;
        }
        switch(fun) {
        case 0:	plot(x11, y);	break;
        case 1:	plot(y, x11);	break;
        case 2:	plot(-y, x11);	break;
        case 3:	plot(-x11, y);	break;
        case 4:	plot(-x11, -y);	break;
        case 5:	plot(-y, -x11);	break;
        case 6:	plot(y, -x11);	break;
        case 7:	plot(x11, -y);	break;
        }
    }
}

int
get_count(int x, int y) {
    return count_buf[curr_width * y + x];
}

void
put_count(int x, int y, int val) {
    count_buf[curr_width * y + x] = val;
}

void
plot(int x, int y) {
    count_buf[curr_width * y + x]++;
}