mandelbrot-gnofract-compiled-formula

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


#ifndef FRACT_STDLIB_H_
#define FRACT_STDLIB_H_

#ifdef __cplusplus
extern "C"
{
#endif

void fract_rand(double *re, double *im);

typedef struct s_arena *arena_t;

arena_t arena_create(int page_size, int max_pages);

void arena_clear(arena_t arena);

void arena_delete(arena_t arena);

void *arena_alloc(
        arena_t arena,
        int element_size,
        int n_dimensions,
        int *n_elements);

void array_get_int(
        void *allocation, int n_dimensions, int *indexes,
        int *pRetVal, int *pInBounds);

void array_get_double(
        void *allocation, int n_dimensions, int *indexes,
        double *pRetVal, int *pInBounds);

int array_set_int(
        void *allocation, int n_dimensions, int *indexes, int val);

int array_set_double(
        void *allocation, int n_dimensions, int *indexes, double val);

void *alloc_array1D(arena_t arena, int element_size, int size);

void *alloc_array2D(arena_t arena, int element_size, int xsize, int ysize);

void *alloc_array3D(arena_t arena, int element_size, int xsize, int ysize, int zsize);

void *alloc_array4D(arena_t arena, int element_size, int xsize, int ysize, int zsize, int wsize);

int read_int_array_1D(void *array, int x);

int write_int_array_1D(void *array, int x, int val);

int read_int_array_2D(void *array, int x, int y);

int write_int_array_2D(void *array, int x, int y, int val);

double read_float_array_1D(void *array, int x);

int write_float_array_1D(void *array, int x, double val);

double read_float_array_2D(void *array, int x, int y);

int write_float_array_2D(void *array, int x, int y, double val);

#ifdef __cplusplus
}
#endif

#endif


#ifndef PF_H_
#define PF_H_

/* C signature of point-funcs generated by compiler

   This is essentially an opaque object with methods implemented in C.
   Typical usage:

   //#pixel.re, #pixel.im, #z.re, #z.im
   double pparams[] = { 1.5, 0.0, 0.0, 0.0};
   double initparams[] = {5.0, 2.0};
   int nItersDone=0;
   int nFate=0;
   double dist=0.0;
   int solid=0;
   pf_obj *pf = pf_new();
   pf->vtbl->init(pf,0.001,initparams,2);

   pf->vtbl->calc(
        pf,
        pparams,
        100,
        0,0,0,
        &nItersDone, &nFate, &dist, &solid);

   pf->vtbl->kill(pf);
*/

// maximum number of params which can be passed to init
#define PF_MAXPARAMS 200

// number of positional params used
#define N_PARAMS 11

/* the 'fate' of a point. This can be either
   Unknown (255) - not yet calculated
   N - reached an attractor numbered N (up to 30)
   N | FATE_INSIDE - did not escape
   N | FATE_SOLID - color with solid color
   N | FATE_DIRECT - color with DCA
*/

typedef unsigned char fate_t;

#define FATE_UNKNOWN 255
#define FATE_SOLID 0x80
#define FATE_DIRECT 0x40
#define FATE_INSIDE 0x20

typedef enum {
    INT = 0,
    FLOAT = 1,
    GRADIENT = 2,
    PARAM_IMAGE = 3
} e_paramtype;

struct s_param {
    e_paramtype t;
    int intval;
    double doubleval;
    void *gradient;
    void *image;
};

struct s_pf_data;

struct s_pf_vtable {
    /* fill in params with the default values for this formula */
    void (*get_defaults)(
            struct s_pf_data *p,
            double *pos_params,
            struct s_param *params,
            int nparams
    );

    /* fill in fields in pf_data with appropriate stuff */
    void (*init)(
            struct s_pf_data *p,
            double *pos_params,
            struct s_param *params,
            int nparams
    );

    /* calculate one point.
       perform up to nIters iterations,
       return:
       number of iters performed in pnIters
       outcome in pFate: 0 = escaped, 1 = trapped.
       More fates may be generated in future
       dist : index into color table from 0.0 to 1.0
    */
    void (*calc)(
            struct s_pf_data *p,
            // in params
            const double *params, int nIters, int warp_param,
            // tolerance params
            int min_period_iter, double period_tolerance,
            // only used for debugging
            int x, int y, int aa,
            // out params
            int *pnIters, int *pFate, double *pDist, int *pSolid,
            int *pDirectColorFlag, double *pColors
    );

    /* deallocate data in p */
    void (*kill)(
            struct s_pf_data *p
    );
};

struct s_pf_data {
    struct s_pf_vtable *vtbl;
    void *arena;
};

typedef struct s_pf_vtable pf_vtable;
typedef struct s_pf_data pf_obj;


#ifdef __cplusplus
extern "C" {
#endif

/* create a new pf_obj.*/
extern pf_obj *pf_new(void);

#ifdef __cplusplus
}
#endif

#endif /* PF_H_ */


typedef struct {
    pf_obj parent;
    struct s_param p[PF_MAXPARAMS];
    double pos_params[N_PARAMS];
    int f__bailout;
    void *t__a__gradient;
    double t__a_fbailout;

    double z_re;
    double z_im;

    double t__f0;
    double t__f1;
    double t__f2;
    double t__f3;
    double t__f4;
    double t__f5;
    double t__f6;
    double t__f7;
    double t__f8;
    double t__f9;
    int t__f10;
    double t__a_cf0_offset;
    double t__a_cf0_density;
    double t__a_cf0bailout;
    double cf0ed;


    double t__cf00;
    double t__cf01;
    double t__cf02;
    double t__cf03;
    double t__cf04;
    double t__cf05;
    double t__cf06;
    double t__cf07;
    double t__cf08;
    double t__cf09;
    double t__a_cf1_offset;
    double t__a_cf1_density;

    double t__cf10;
    double t__cf11;
} pf_real;

static void pf_init(
        struct s_pf_data *p_stub,
        double *pos_params,
        struct s_param *params,
        int nparams) {
    pf_real *pfo = (pf_real *) p_stub;
    int i;

    if (nparams > PF_MAXPARAMS) {
        nparams = PF_MAXPARAMS;
    }
    for (i = 0; i < nparams; ++i) {
        pfo->p[i] = params[i];
        /* printf("param %d = %.17g\n",i,params[i]); */
    }
    for (i = 0; i < N_PARAMS; ++i) {
        pfo->pos_params[i] = pos_params[i];
    }
}

static void pf_get_defaults(
        // "object" pointer
        struct s_pf_data *t__p_stub,
        // in params
        double *pos_params,
        // out params
        struct s_param *params,
        // in param
        int nparams
) {
    pf_real *t__pfo = (pf_real *) t__p_stub;
    /*
    t__start_fdefault: ;

t__a_fbailout = 4.00000000000000000;
goto t__end_fdefault;
t__end_fdefault: ;

    t__pfo->p[3].doubleval = t__a_cf0_density;
t__pfo->p[2].doubleval = t__a_cf0_offset;
t__pfo->p[4].doubleval = t__a_cf0bailout;
t__pfo->p[6].doubleval = t__a_cf1_density;
t__pfo->p[5].doubleval = t__a_cf1_offset;
t__pfo->p[0].gradient = t__a__gradient;
t__pfo->p[1].doubleval = t__a_fbailout;
    */
}

static void pf_calc(
        // "object" pointer
        struct s_pf_data *t__p_stub,
        // in params
        const double *t__params, int maxiter, int t__warp_param,
        // periodicity params
        int min_period_iter, double period_tolerance,
        // only used for debugging
        int t__p_x, int t__p_y, int t__p_aa,
        // out params
        int *t__p_pnIters, int *t__p_pFate, double *t__p_pDist, int *t__p_pSolid,
        int *t__p_pDirectColorFlag, double *t__p_pColors
) {
    pf_real *t__pfo = (pf_real *) t__p_stub;

    double pixel_re = t__params[0];
    double pixel_im = t__params[1];
    double t__h_zwpixel_re = t__params[2];
    double t__h_zwpixel_im = t__params[3];

    double t__h_index = 0.0;
    int t__h_solid = 0;
    int t__h_fate = 0;
    int t__h_inside = 0;
    double t__h_color_re = 0.0;
    double t__h_color_i = 0.0;
    double t__h_color_j = 0.0;
    double t__h_color_k = 0.0;

    *t__p_pDirectColorFlag = 0;

    if (t__warp_param != -1) {
        t__pfo->p[t__warp_param].doubleval = t__h_zwpixel_re;
        t__pfo->p[t__warp_param + 1].doubleval = t__h_zwpixel_im;
        t__h_zwpixel_re = t__h_zwpixel_im = 0.0;
    }

    /* variable declarations */
    int f__bailout = 0;
    void *t__a__gradient = t__pfo->p[0].gradient;
    double t__a_fbailout = t__pfo->p[1].doubleval;

    double z_re = 0.00000000000000000;
    double z_im = 0.00000000000000000;

    double t__f0 = 0.00000000000000000;
    double t__f1 = 0.00000000000000000;
    double t__f2 = 0.00000000000000000;
    double t__f3 = 0.00000000000000000;
    double t__f4 = 0.00000000000000000;
    double t__f5 = 0.00000000000000000;
    double t__f6 = 0.00000000000000000;
    double t__f7 = 0.00000000000000000;
    double t__f8 = 0.00000000000000000;
    double t__f9 = 0.00000000000000000;
    int t__f10 = 0;
    double t__a_cf0_offset = t__pfo->p[2].doubleval;
    double t__a_cf0_density = t__pfo->p[3].doubleval;
    double t__a_cf0bailout = t__pfo->p[4].doubleval;
    double cf0ed = 0.00000000000000000;


    double t__cf00 = 0.00000000000000000;
    double t__cf01 = 0.00000000000000000;
    double t__cf02 = 0.00000000000000000;
    double t__cf03 = 0.00000000000000000;
    double t__cf04 = 0.00000000000000000;
    double t__cf05 = 0.00000000000000000;
    double t__cf06 = 0.00000000000000000;
    double t__cf07 = 0.00000000000000000;
    double t__cf08 = 0.00000000000000000;
    double t__cf09 = 0.00000000000000000;
    double t__a_cf1_offset = t__pfo->p[5].doubleval;
    double t__a_cf1_density = t__pfo->p[6].doubleval;

    double t__cf10 = 0.00000000000000000;
    double t__cf11 = 0.00000000000000000;

    double old_z_re;
    double old_z_im;
    int period_iters = 0;
    int save_mask = 9;
    int save_incr = 1;
    int next_save_incr = 4;

    int t__h_numiter = 0;


    t__start_finit:;

    z_re = t__h_zwpixel_re;
    z_im = t__h_zwpixel_im;
    goto t__end_finit;
    t__end_finit:;


    old_z_re = z_re;
    old_z_im = z_im;
    do {
        t__start_floop:;

        t__f0 = z_re * z_re;
        t__f1 = z_im * z_im;
        t__f2 = z_re * z_im;
        t__f3 = t__f0 - t__f1;
        t__f4 = 2.00000000000000000 * t__f2;
        t__f5 = t__f3 + pixel_re;
        t__f6 = t__f4 + pixel_im;
        z_re = t__f5;
        z_im = t__f6;
        goto t__end_floop;
        t__end_floop:;


        t__start_fbailout:;

        t__f7 = z_re * z_re;
        t__f8 = z_im * z_im;
        t__f9 = t__f7 + t__f8;
        t__f10 = t__f9 < t__a_fbailout;
        f__bailout = t__f10;
        goto t__end_fbailout;
        t__end_fbailout:;


        if (!f__bailout) break;

        if (t__h_numiter >= min_period_iter) {
            if ((t__h_numiter & save_mask) == 0) {
                /* save a value */
                old_z_re = z_re;
                old_z_im = z_im;

                if (--save_incr == 0) {
                    /* lengthen period check */
                    save_mask = (save_mask << 1) + 1;
                    save_incr = next_save_incr;
                }
            } else {
                /* compare to an older value */
                if ((fabs(z_re - old_z_re) < period_tolerance)
                    && (fabs(z_im - old_z_im) < period_tolerance)) {
                    period_iters = t__h_numiter;
                    //t__h_numiter = maxiter;
                    t__h_inside = 1;
                    *t__p_pFate = 1;

                    goto loop_done;
                }
            }
        }


        t__h_numiter++;
    } while (t__h_numiter < maxiter);

    /* fate of 0 = escaped, 1 = trapped */
    t__h_inside = (t__h_numiter >= maxiter);
    *t__p_pFate = (t__h_numiter >= maxiter);

    loop_done:


    *t__p_pnIters = t__h_numiter;
    if (t__h_inside == 0) {
        t__start_cf0final:;

        t__cf00 = z_re * z_re;
        t__cf01 = z_im * z_im;
        t__cf02 = t__cf00 + t__cf01;
        t__cf03 = t__cf02 + 0.00000000100000000;
        t__cf04 = t__a_cf0bailout / t__cf03;
        cf0ed = t__cf04;
        t__cf05 = ((double) t__h_numiter);
        t__cf06 = t__cf05 + cf0ed;
        t__cf07 = t__cf06 / 256.00000000000000000;
        t__h_index = t__cf07;
        t__cf08 = t__a_cf0_density * t__h_index;
        t__cf09 = t__cf08 + t__a_cf0_offset;
        t__h_index = t__cf09;
        goto t__end_cf0final;
        t__end_cf0final:;

        ;
    } else {
        t__start_cf1final:;

        t__h_solid = 1;
        t__cf10 = t__a_cf1_density * t__h_index;
        t__cf11 = t__cf10 + t__a_cf1_offset;
        t__h_index = t__cf11;
        goto t__end_cf1final;
        t__end_cf1final:;

        ;
    }
    *t__p_pFate = t__h_fate | (t__h_inside ? FATE_INSIDE : 0);
    *t__p_pDist = t__h_index;
    *t__p_pSolid = t__h_solid;


    arena_clear((arena_t) (t__p_stub->arena));
    return;
}

static void pf_kill(
        struct s_pf_data *p_stub) {
    arena_delete((arena_t) (p_stub->arena));
    free(p_stub);
}

static struct s_pf_vtable vtbl =
        {
                pf_get_defaults,
                pf_init,
                pf_calc,
                pf_kill
        };

pf_obj *pf_new() {
    pf_real *p = (pf_real *) malloc(sizeof(pf_real));
    if (!p) return NULL;
    p->parent.vtbl = &vtbl;
    p->parent.arena = arena_create(100000, 1);
    return (pf_obj *) p;
}