RIFC Experimental Plotter by Chris M. Thomasson

/* Experimental Fractal Data Storage by:
   Chris M. Thomasson
   Copyright 2017
________________________________________________________*/


#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <complex.h>
#include <tgmath.h>
#include <stdbool.h>


#define CT_PI2 6.283185307179586476925286766559

// Our loading threshold, 32 nits
#define CT_LOAD_MAX 32

// Our epsilon for error
#define CT_FIND_EPS 0.00000001

typedef double ct_real;
typedef double complex ct_complex;


/* Generic Plotter by Chris M. Thomasson
________________________________________________________*/
struct ct_axes
{
    ct_real xmin;
    ct_real xmax;
    ct_real ymin;
    ct_real ymax;
};

struct ct_axes
ct_axes_get(
    ct_complex center,
    ct_real diameter
){
    ct_real radius = diameter / 2;

    struct ct_axes axes = {
        creal(center) - radius,
        creal(center) + radius,
        cimag(center) - radius,
        cimag(center) + radius,
    };

    return axes;
}


struct ct_canvas_color
{
    unsigned char red;
    unsigned char green;
    unsigned char blue;
    unsigned char alpha;
};


struct ct_canvas
{
    unsigned long width;
    unsigned long height;
    struct ct_canvas_color* buf;
};

bool
ct_canvas_create(
    struct ct_canvas* const self,
    unsigned long width,
    unsigned long height
){
    size_t size = width * height * sizeof(*self->buf);

    self->buf = calloc(1, size);

    if (self->buf)
    {
        self->width = width;
        self->height = height;

        return true;
    }

    return false;
}

void
ct_canvas_destroy(
    struct ct_canvas const* const self
){
    free(self->buf);
}

bool
ct_canvas_save_ppm(
     struct ct_canvas const* const self,
     char const* fname
){
    FILE* fout = fopen(fname, "w");

    if (fout)
    {
        char const ppm_head[] =
         "P3\n"
         "# Chris M. Thomasson RIFC Cipher Renderer ver:0.0.0.0 (pre-alpha)";

        fprintf(fout, "%s\n%lu %lu\n%u\n",
             ppm_head,
             self->width, self->height,
             255U);

        size_t size = self->width * self->height;

        for (size_t i = 0; i < size; ++i)
        {
            struct ct_canvas_color c = self->buf[i];
            fprintf(fout, "%u %u %u  ", c.red, c.green, c.blue);
        }

        if (! fclose(fout))
        {
            return true;
        }
    }

    return false;
}


struct ct_plane
{
    struct ct_axes axes;
    struct ct_canvas* canvas;
};

struct ct_pixel
{
    size_t x;
    size_t y;
};


struct ct_pixel
ct_plane_project(
    struct ct_plane const* const self,
    ct_complex c
){
     ct_real awidth = self->axes.xmax - self->axes.xmin;
     ct_real aheight = self->axes.ymax - self->axes.ymin;

     ct_real xstep = awidth / (self->canvas->width - 1.0);
     ct_real ystep = aheight / (self->canvas->height - 1.0);

     struct ct_pixel pixel = {
         (creal(c) - self->axes.xmin) / xstep,
         (self->axes.ymax - cimag(c)) / ystep
     };

     return pixel;
}

bool
ct_plane_plot(
    struct ct_plane* const self,
    ct_complex c,
    struct ct_canvas_color color
){
     struct ct_pixel pixel = ct_plane_project(self, c);

     if (pixel.x < self->canvas->width && pixel.y < self->canvas->height)
     {
         size_t cp = pixel.x + pixel.y * self->canvas->height;

         if (cp < self->canvas->height * self->canvas->width)
         {
             self->canvas->buf[cp] = color;
             return true;
         }
     }

     return false;
}


/* RIFC (original) by Chris M. Thomasson
________________________________________________________*/
ct_complex
ct_root(
    ct_complex const z,
    long p,
    unsigned long r
){
    ct_real radius = pow(fabs(z), 1.0 / p);
    ct_real angle_base = carg(z) / p;
    ct_real angle_step = CT_PI2 / p;
    ct_real angle = angle_base + angle_step * r;
    ct_complex root = cos(angle) * radius + sin(angle) * radius * I;

    return root;
}

// Stores a number v using a base
ct_complex
ct_store(
    struct ct_plane* const plane,
    ct_complex z,
    ct_complex c,
    unsigned long v,
    long p
){
    unsigned long b = abs(p);

    /*
    printf("ct_store:((%lf%+lfi), (%lf%+lfi), %lu, %ld)\n"
           "_______________________________\n",
           creal(z), cimag(z), creal(c), cimag(c), v, p);
*/
    while (v)
    {
        // Gain the base
        unsigned long d = v / b;
        unsigned long r = v - d * b;

        // Compute the proper root for the base
        z = ct_root(z - c, p, r);

        //printf("(%lf%+lfi):%lu\n", creal(z), cimag(z), r);

        struct ct_canvas_color color = {
            v & 0x000000FF,
            (v & 0x0000FF00) >> 8,
            //(v & 0x00FF0000) >> 16,
            255,
            0
        };

        // Plot the point
        ct_plane_plot(plane, z, color);

        v = d;
    }

   // printf("result:(%lf%+lfi)\n", creal(z), cimag(z));
    //printf("_______________________________\n\n");

    return z;
}

// Loads a number at a base from z.
ct_complex
ct_load(
    ct_complex z0,
    ct_complex z,
    ct_complex c,
    unsigned long* pv,
    long p
){
    unsigned long v = 0;
    unsigned long b = abs(p);

    /*
    printf("ct_load:((%lf%+lfi), (%lf%+lfi), %p, %ld)\n"
       "_______________________________\n",
       creal(z), cimag(z), creal(c), cimag(c), (void*)pv, p);
       */

    for (unsigned long i = 0; i < CT_LOAD_MAX; ++i)
    {
        // Check for termination condition
        ct_real d = fabs(z - z0);
        if (d < CT_FIND_EPS) break;

        // Raise z to the power and add c
        ct_complex zn = pow(z, p) + c;

        ct_complex r = 0 + 0 * I;
        unsigned long ii = 0;

        // Search for the stored root
        for (ii = 0; ii < b; ++ii)
        {
            r = ct_root(zn - c, p, ii);
            ct_real d = fabs(z - r);
            if (d < CT_FIND_EPS) break;
        }

        // Recreate the number
        v = b * v + ii;

        //printf("(%lf%+lfi):%lu\n", creal(z), cimag(z), ii);

        // Set the next iteration
        z = zn;
    }

    *pv = v;

    /*
    printf("result:((%lf%+lfi), %lu)\n", creal(z), cimag(z), *pv);
    printf("_______________________________\n\n");
    */

    return z;
}


// Canvas width and height
#define CT_WIDTH 512
#define CT_HEIGHT CT_WIDTH

int main(void)
{
    struct ct_canvas canvas;

    bool status = ct_canvas_create(&canvas, CT_WIDTH, CT_HEIGHT);
    assert(status);

    // Setup our plotting axes and plane
    ct_complex center = 0 + 0 * I;
    ct_real zoom_diameter = 3.5;
    struct ct_axes axes = ct_axes_get(center, zoom_diameter);
    struct ct_plane plane = { axes, &canvas };

    // Our initial conditions
    ct_complex z = -1+.0;
    ct_complex c = -.75+.06*I;

    // Power of -2
    long p = 2;
    unsigned long e = 0;

    // Try to store/load values 0 through vmax-1
    unsigned long vmax = 65536;
    for (unsigned long v = 0; v < vmax; ++v)
    {
        // Store the data!
        ct_complex z0 = z;
        z = ct_store(&plane, z, c, v, p);

        // Load the data and check for errors!
        unsigned long lv = 0;
        ct_load(z0, z, c, &lv, p);

        if (lv != v)
        {
            //printf("ERROR: %lu != %lu       \n", lv, v);
            e += 1;
        }

        // Display progress
        if (! (v % (vmax / 2048)))
        {
            printf("plotted:%lu of %lu, error: %lu\r", v, vmax, e);
        }
    }

    printf("Plotting Complete:%lu of %lu, %lu errors\r", vmax, vmax, e);

    status = ct_canvas_save_ppm(&canvas, "ct_cipher_rifc.ppm");
    assert(status);
    printf("\ncreated: ct_cipher_rifc.ppm\n");
    ct_canvas_destroy(&canvas);

    return 0;
}