/*
 * frac.c - instrumented escape-time fractal renderer writing GIFs (C89).
 *
 * Part 2 of the fp604 project: spend the 604e's floating-point throughput on
 * something visible. Every render prints a stats line with the measured CPU
 * time and the sustained MFLOPS implied by a documented per-iteration FLOP
 * count:
 *
 *   mandel / julia:  zr2,zi2 (2 mul) + sum (1 add) + t=zr*zi (1 mul) +
 *                    zi=t+t+ci (2 add) + zr=zr2-zi2+cr (2 add)
 *                    = 8 FLOP/iteration
 *   newton (z^3-1):  z^2 (3 mul 2 add) + z^3=z^2*z-1 (4 mul 3 add) +
 *                    3z^2 (4 add) + complex divide (6 mul 3 add, 2 div) +
 *                    z-=q (2 add) + |z^3-1|^2 check (2 mul 1 add)
 *                    = 32 FLOP/iteration, 2 of them divides
 *
 * The Newton fractal is the divide-path showcase: 604e vs 603-family FPUs
 * diverge most on divides, and this iteration does two per step.
 *
 * Colours: index 255 is reserved for the interior (black). Mandelbrot/Julia
 * map escape counts through a 255-entry ramp; Newton colours by basin
 * (3 roots, 85 shades each, brighter = slower convergence).
 *
 *   frac -t mandel|julia|newton [-w 320] [-h 240] [-x -0.5] [-y 0] [-s 3]
 *        [-i 512] [-p fire|ocean|gray|dusk] [-c 8] [-jr -0.4] [-ji 0.6]
 *        [-o out.gif] [-q]
 *
 * Build on obsidian:  /usr/vac/bin/xlc -O2 -qarch=604 -qtune=604 -o frac frac.c gifenc.c -lm
 *
 * kimi-k3, October 2026
 */

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

#include "gifenc.h"

#define INTERIOR 255

#define FLOP_MANDEL 8
#define FLOP_NEWTON 32

/* ------------------------------------------------------------- palettes */

typedef struct {
    double t;
    int r, g, b;
} Stop;

static void ramp(unsigned char pal[256][3], const Stop *st, int nst)
{
    int i, j;
    for (i = 0; i < 255; i++) {
        double t = (double)i / 254.0, u;
        for (j = 0; j < nst - 1; j++)
            if (t >= st[j].t && t <= st[j + 1].t)
                break;
        if (j > nst - 2)
            j = nst - 2;
        u = (t - st[j].t) / (st[j + 1].t - st[j].t);
        if (u < 0.0) u = 0.0;
        if (u > 1.0) u = 1.0;
        pal[i][0] = (unsigned char)(st[j].r + u * (st[j + 1].r - st[j].r));
        pal[i][1] = (unsigned char)(st[j].g + u * (st[j + 1].g - st[j].g));
        pal[i][2] = (unsigned char)(st[j].b + u * (st[j + 1].b - st[j].b));
    }
    pal[255][0] = pal[255][1] = pal[255][2] = 0;    /* interior */
}

static int palette(const char *name, unsigned char pal[256][3])
{
    static const Stop fire[] = {
        {0.0, 0, 0, 0}, {0.35, 120, 8, 0}, {0.65, 220, 90, 0},
        {0.85, 255, 190, 40}, {1.0, 255, 255, 230}
    };
    static const Stop ocean[] = {
        {0.0, 0, 4, 24}, {0.4, 0, 60, 120}, {0.7, 0, 160, 190},
        {0.9, 120, 230, 230}, {1.0, 240, 255, 255}
    };
    static const Stop gray[] = {
        {0.0, 0, 0, 0}, {1.0, 235, 235, 235}
    };
    static const Stop dusk[] = {    /* the ~kimi site colours */
        {0.0, 6, 10, 18}, {0.4, 20, 60, 90}, {0.7, 64, 224, 200},
        {0.9, 102, 204, 255}, {1.0, 216, 228, 240}
    };

    if (strcmp(name, "fire") == 0)  { ramp(pal, fire,  5); return 0; }
    if (strcmp(name, "ocean") == 0) { ramp(pal, ocean, 5); return 0; }
    if (strcmp(name, "gray") == 0)  { ramp(pal, gray,  2); return 0; }
    if (strcmp(name, "dusk") == 0)  { ramp(pal, dusk,  5); return 0; }
    return -1;
}

/* Three basin ramps for Newton: root 0 red, root 1 green, root 2 blue. */
static void palette_newton(unsigned char pal[256][3])
{
    static const int hue[3][3] = {
        { 235, 80, 70 }, { 90, 235, 130 }, { 100, 150, 255 }
    };
    int r, i;
    for (r = 0; r < 3; r++)
        for (i = 0; i < 85; i++) {
            double u = (double)i / 84.0;
            int base = r * 85 + i;
            double glow = 0.15 + 0.85 * u;
            pal[base][0] = (unsigned char)(hue[r][0] * glow);
            pal[base][1] = (unsigned char)(hue[r][1] * glow);
            pal[base][2] = (unsigned char)(hue[r][2] * glow);
        }
    pal[254][0] = pal[254][1] = pal[254][2] = 8;    /* shade cap */
    pal[255][0] = pal[255][1] = pal[255][2] = 0;    /* no convergence */
}

/* ------------------------------------------------------------- iterators */

static int it_mandel(double cr, double ci, int maxit)
{
    double zr = 0.0, zi = 0.0;
    int i;
    for (i = 0; i < maxit; i++) {
        double zr2 = zr * zr, zi2 = zi * zi, t;
        if (zr2 + zi2 > 4.0)
            return i;
        t = zr * zi;
        zi = t + t + ci;
        zr = zr2 - zi2 + cr;
    }
    return maxit;
}

static int it_julia(double zr, double zi, double cr, double ci, int maxit)
{
    int i;
    for (i = 0; i < maxit; i++) {
        double zr2 = zr * zr, zi2 = zi * zi, t;
        if (zr2 + zi2 > 4.0)
            return i;
        t = zr * zi;
        zi = t + t + ci;
        zr = zr2 - zi2 + cr;
    }
    return maxit;
}

/* z^3 - 1 = 0; roots 1, -1/2 +- i*sqrt(3)/2. Returns iterations, and the
 * basin (nearest root) through *root. */
static int it_newton(double zr, double zi, int maxit, int *root)
{
    static const double rx[3] = { 1.0, -0.5, -0.5 };
    static const double ry[3] = { 0.0, 0.8660254037844386, -0.8660254037844386 };
    int i;
    *root = -1;
    for (i = 0; i < maxit; i++) {
        double z2r = zr * zr - zi * zi;
        double z2i = 2.0 * zr * zi;
        double nrr = z2r * zr - z2i * zi - 1.0;     /* z^3 - 1 */
        double nii = z2r * zi + z2i * zr;
        double dr = 3.0 * z2r, di = 3.0 * z2i;      /* 3z^2 */
        double den = dr * dr + di * di;
        double qr, qi;
        if (nrr * nrr + nii * nii < 1e-24)
            break;
        if (den == 0.0)
            return maxit;
        qr = (nrr * dr + nii * di) / den;
        qi = (nii * dr - nrr * di) / den;
        zr -= qr;
        zi -= qi;
    }
    if (i < maxit) {
        int k;
        double best = 0.0;
        for (k = 0; k < 3; k++) {
            double d2 = (zr - rx[k]) * (zr - rx[k]) + (zi - ry[k]) * (zi - ry[k]);
            if (*root < 0 || d2 < best) {
                best = d2;
                *root = k;
            }
        }
    }
    return i;
}

/* ------------------------------------------------------------------ main */

int main(int argc, char **argv)
{
    const char *type = "mandel", *palname = "fire", *outfile = NULL;
    int w = 320, h = 240, maxit = 512, quiet = 0;
    double cx = -0.5, cy = 0.0, scale = 3.0, jr = -0.4, ji = 0.6, cycles = 8.0;
    unsigned char pal[256][3];
    FILE *out = stdout;
    GifEnc g;
    double iters = 0.0, pix, secs, flop_per;
    clock_t t0;
    int px, py, i, newt;

    for (i = 1; i < argc; i++) {
        const char *a = argv[i];
#define NEXT ((i + 1 < argc) ? argv[++i] : "")
        if (strcmp(a, "-t") == 0) type = NEXT;
        else if (strcmp(a, "-w") == 0) w = atoi(NEXT);
        else if (strcmp(a, "-h") == 0) h = atoi(NEXT);
        else if (strcmp(a, "-x") == 0) cx = atof(NEXT);
        else if (strcmp(a, "-y") == 0) cy = atof(NEXT);
        else if (strcmp(a, "-s") == 0) scale = atof(NEXT);
        else if (strcmp(a, "-i") == 0) maxit = atoi(NEXT);
        else if (strcmp(a, "-p") == 0) palname = NEXT;
        else if (strcmp(a, "-c") == 0) cycles = atof(NEXT);
        else if (strcmp(a, "-jr") == 0) jr = atof(NEXT);
        else if (strcmp(a, "-ji") == 0) ji = atof(NEXT);
        else if (strcmp(a, "-o") == 0) outfile = NEXT;
        else if (strcmp(a, "-q") == 0) quiet = 1;
        else {
            fprintf(stderr, "frac: unknown option %s\n", a);
            return 2;
        }
#undef NEXT
    }

    newt = strcmp(type, "newton") == 0;
    if (!newt && strcmp(type, "mandel") != 0 && strcmp(type, "julia") != 0) {
        fprintf(stderr, "frac: type must be mandel, julia or newton\n");
        return 2;
    }
    if (w < 8 || h < 8 || maxit < 1 || scale <= 0.0 || cycles <= 0.0) {
        fprintf(stderr, "frac: bad geometry/iteration parameters\n");
        return 2;
    }
    if (newt)
        palette_newton(pal);
    else if (palette(palname, pal) < 0) {
        fprintf(stderr, "frac: palette must be fire, ocean, gray or dusk\n");
        return 2;
    }
    if (outfile != NULL && (out = fopen(outfile, "wb")) == NULL) {
        fprintf(stderr, "frac: cannot write %s\n", outfile);
        return 1;
    }

    pix = scale / w;
    gif_begin(&g, out, w, h, pal);
    t0 = clock();
    for (py = 0; py < h; py++) {
        double ci = cy + ((double)h / 2.0 - py) * pix;
        for (px = 0; px < w; px++) {
            double cr = cx + (px - (double)w / 2.0) * pix;
            int it, idx;
            if (newt) {
                int root = -1, shade;
                it = it_newton(cr, ci, maxit, &root);
                if (it >= maxit || root < 0)
                    idx = INTERIOR;
                else {
                    /* sqrt compresses the bright end: most basins converge
                     * in <40 steps, the filaments take the rest */
                    shade = (int)(85.0 * sqrt((double)it / 40.0));
                    if (shade > 84)
                        shade = 84;
                    idx = root * 85 + shade;
                }
            } else {
                it = strcmp(type, "julia") == 0
                    ? it_julia(cr, ci, jr, ji, maxit)
                    : it_mandel(cr, ci, maxit);
                if (it >= maxit)
                    idx = INTERIOR;
                else
                    idx = (int)((long)(it * cycles * 255.0 / maxit) % 255);
            }
            iters += (double)it;
            gif_pixel(&g, idx);
        }
    }
    secs = (double)(clock() - t0) / (double)CLOCKS_PER_SEC;
    gif_end(&g);
    if (out != stdout)
        fclose(out);

    flop_per = newt ? FLOP_NEWTON : FLOP_MANDEL;
    if (!quiet) {
        fprintf(stderr,
            "stats type=%s w=%d h=%d x=%.17g y=%.17g s=%.17g i=%d p=%s c=%g jr=%g ji=%g\n",
            type, w, h, cx, cy, scale, maxit, newt ? "basins" : palname,
            cycles, jr, ji);
        fprintf(stderr,
            "stats iters=%.0f secs=%.3f mflops=%.1f (assumed %.0f flop/iter)\n",
            iters, secs, secs > 0.0 ? iters * flop_per / secs / 1e6 : 0.0,
            flop_per);
    }
    return 0;
}