/*
 * spectre_ca.c — Spectre Face Cellular Automata (WASM core)
 *
 * Compiled with: clang --target=wasm32 -O2 -nostdlib -Wl,--no-entry
 *                -Wl,--export-all -Wl,--initial-memory=134217728
 *                -o spectre_ca.wasm spectre_ca.c
 *
 * Exports a flat memory layout that JS reads directly for rendering.
 * All tiling, adjacency, and CA stepping happen in WASM.
 */

#define NULL ((void*)0)

/* ── math imports from JS (hardware-precision sin/cos) ───────── */

extern double js_sin(double x);
extern double js_cos(double x);

static double floor_(double x) {
    int i = (int)x;
    return (double)(x < (double)i ? i - 1 : i);
}

/* ── bump allocator in WASM linear memory ────────────────────── */

extern unsigned char __heap_base;
static unsigned char *heap_ptr;

static void heap_reset(void) {
    heap_ptr = &__heap_base;
}

static void *alloc(int bytes) {
    unsigned long addr = (unsigned long)heap_ptr;
    addr = (addr + 7) & ~7UL;
    heap_ptr = (unsigned char *)(addr + bytes);
    return (void *)addr;
}

/* zero-fill allocation */
static void *zalloc(int bytes) {
    void *ptr = alloc(bytes);
    unsigned char *p = (unsigned char *)ptr;
    for (int i = 0; i < bytes; i++) p[i] = 0;
    return ptr;
}

/* ── constants ───────────────────────────────────────────────── */

#define MAX_TILES     50000
#define POINTS_PER    14
#define MAX_CHILDREN  8
#define MAX_DEPTH     5
#define PI            3.14159265358979323846

/* ── types ───────────────────────────────────────────────────── */

typedef struct { double x, y; } Point;
typedef struct { double a, b, c, d, e, f; } Affine;

typedef struct {
    float px[POINTS_PER];         /* 56 bytes  offset 0   */
    float py[POINTS_PER];         /* 56 bytes  offset 56  */
    int   neighbor[POINTS_PER];   /* 56 bytes  offset 112 */
    int   alive;                  /* 4 bytes   offset 168 */
    int   nextAlive;              /* 4 bytes   offset 172 */
    int   age;                    /* 4 bytes   offset 176 */
    int   liveFaces;              /* 4 bytes   offset 180 */
} Tile;                           /* total: 184 bytes      */

/* ── globals ─────────────────────────────────────────────────── */

static Tile   *tiles;
static int     tileCount;
static int     generation;

static int    cfgDepth      = 3;
static int    cfgBirthFaces = 3;
static int    cfgSurviveMin = 2;
static int    cfgSurviveMax = 3;
static int    cfgFaceWeight = 1;
static double cfgDensity    = 0.42;
static int    cfgRoot       = 0;
static double cfgSeed       = 1.0;

static double fitScale, fitOx, fitOy;
static double screenW = 1200.0, screenH = 800.0;

/* ── affine helpers ──────────────────────────────────────────── */

static Affine affIdent(void) { return (Affine){1,0,0, 0,1,0}; }

static Affine affMul(Affine A, Affine B) {
    return (Affine){
        A.a*B.a + A.b*B.d, A.a*B.b + A.b*B.e, A.a*B.c + A.b*B.f + A.c,
        A.d*B.a + A.e*B.d, A.d*B.b + A.e*B.e, A.d*B.c + A.e*B.f + A.f
    };
}

static Point affApply(Affine M, Point p) {
    return (Point){ M.a*p.x + M.b*p.y + M.c, M.d*p.x + M.e*p.y + M.f };
}

static Affine affRot(double angle) {
    double c = js_cos(angle), s = js_sin(angle);
    return (Affine){c, -s, 0, s, c, 0};
}

static Affine affTrans(double tx, double ty) { return (Affine){1,0,tx, 0,1,ty}; }

static Affine affTransTo(Point from, Point to) {
    return affTrans(to.x - from.x, to.y - from.y);
}

static double radians(double deg) { return deg * PI / 180.0; }

/* ── spectre base shape ──────────────────────────────────────── */

static const Point SPECTRE[14] = {
    {0, 0}, {1.0, 0.0}, {1.5, -0.8660254037844386},
    {2.366025403784439, -0.36602540378443865},
    {2.366025403784439, 0.6339745962155614},
    {3.366025403784439, 0.6339745962155614},
    {3.866025403784439, 1.5},
    {3.0, 2.0},
    {2.133974596215561, 1.5},
    {1.6339745962155614, 2.3660254037844393},
    {0.6339745962155614, 2.3660254037844393},
    {-0.3660254037844386, 2.3660254037844393},
    {-0.866025403784439, 1.5},
    {0.0, 1.0}
};

/* ── substitution rules ──────────────────────────────────────── */

enum {
    L_DELTA=0, L_GAMMA, L_THETA, L_LAMBDA, L_XI, L_PI, L_SIGMA, L_PHI, L_PSI,
    L_NULL
};

static const int SUPER_RULES[9][8] = {
    { L_XI,  L_DELTA, L_XI,   L_PHI, L_SIGMA, L_PI,  L_PHI,    L_GAMMA },
    { L_PI,  L_DELTA, L_NULL, L_THETA,L_SIGMA,L_XI,  L_PHI,    L_GAMMA },
    { L_PSI, L_DELTA, L_PI,   L_PHI, L_SIGMA, L_PI,  L_PHI,    L_GAMMA },
    { L_PSI, L_DELTA, L_XI,   L_PHI, L_SIGMA, L_PI,  L_PHI,    L_GAMMA },
    { L_PSI, L_DELTA, L_PI,   L_PHI, L_SIGMA, L_PSI, L_PHI,    L_GAMMA },
    { L_PSI, L_DELTA, L_XI,   L_PHI, L_SIGMA, L_PSI, L_PHI,    L_GAMMA },
    { L_XI,  L_DELTA, L_XI,   L_PHI, L_SIGMA, L_PI,  L_LAMBDA, L_GAMMA },
    { L_PSI, L_DELTA, L_PSI,  L_PHI, L_SIGMA, L_PI,  L_PHI,    L_GAMMA },
    { L_PSI, L_DELTA, L_PSI,  L_PHI, L_SIGMA, L_PSI, L_PHI,    L_GAMMA },
};

static const int T_RULES[7][3] = {
    {60,3,1}, {0,2,0}, {60,3,1}, {60,3,1}, {0,2,0}, {60,3,1}, {-120,3,3}
};

/* Per-level transforms: JS recomputes transforms at each supertile level
   using the evolving superQuad. We precompute all levels. */
static Affine LEVEL_TRANSFORMS[MAX_DEPTH][8];
static int    numLevels;

static void computeTransformsFromQuad(Point quad[4], Affine out[8]) {
    Affine reflect = {-1,0,0, 0,1,0};
    out[0] = affIdent();
    int totalAngle = 0;
    Affine rot = affIdent();
    Point tquad[4];
    for (int i = 0; i < 4; i++) tquad[i] = quad[i];

    for (int r = 0; r < 7; r++) {
        totalAngle += T_RULES[r][0];
        if (T_RULES[r][0] != 0) {
            rot = affRot(radians((double)totalAngle));
            for (int i = 0; i < 4; i++) tquad[i] = affApply(rot, quad[i]);
        }
        Point fromPt = tquad[T_RULES[r][2]];
        Point toPt   = affApply(out[r], quad[T_RULES[r][1]]);
        Affine slide  = affTransTo(fromPt, toPt);
        out[r + 1] = affMul(slide, rot);
    }
    for (int i = 0; i < 8; i++)
        out[i] = affMul(reflect, out[i]);
}

static void buildAllLevelTransforms(int depth) {
    /* Start with base quad */
    Point quad[4];
    quad[0] = SPECTRE[3]; quad[1] = SPECTRE[5];
    quad[2] = SPECTRE[7]; quad[3] = SPECTRE[11];

    numLevels = depth;

    for (int level = 0; level < depth; level++) {
        Affine Ts[8];
        computeTransformsFromQuad(quad, Ts);

        /* Store transforms for this level */
        for (int i = 0; i < 8; i++)
            LEVEL_TRANSFORMS[level][i] = Ts[i];

        /* Compute superQuad for next level (matches JS) */
        Point nextQuad[4];
        nextQuad[0] = affApply(Ts[6], quad[2]);
        nextQuad[1] = affApply(Ts[5], quad[1]);
        nextQuad[2] = affApply(Ts[3], quad[2]);
        nextQuad[3] = affApply(Ts[0], quad[1]);

        for (int i = 0; i < 4; i++) quad[i] = nextQuad[i];
    }
}

/* ── tiling generation ───────────────────────────────────────── */

static double *rawPx, *rawPy;
static int rawCount;

static void emitTile(Affine xform) {
    if (rawCount >= MAX_TILES) return;
    int base = rawCount * POINTS_PER;
    for (int i = 0; i < POINTS_PER; i++) {
        Point p = affApply(xform, SPECTRE[i]);
        rawPx[base + i] = p.x;
        rawPy[base + i] = p.y;
    }
    rawCount++;
}

/* depth counts down from totalDepth to 0.
   Level index = totalDepth - depth (outermost level uses highest-numbered transforms).
   JS builds: level0 from baseQuad, level1 from superQuad0, level2 from superQuad1...
   Flatten walks: outermost (level N-1) → ... → level 0 → emit base shapes.
   So at depth D, we use LEVEL_TRANSFORMS[numLevels - depth]. */
static void flattenRecurse(int label, int depth, Affine xform) {
    if (rawCount >= MAX_TILES) return;
    if (depth == 0) {
        if (label == L_GAMMA) {
            emitTile(xform);
            Affine gamma2 = affMul(xform, affMul(affTrans(SPECTRE[8].x, SPECTRE[8].y), affRot(PI / 6.0)));
            emitTile(gamma2);
            return;
        }
        emitTile(xform);
        return;
    }
    if (label < 0 || label > L_PSI) return;
    const int *children = SUPER_RULES[label];
    int levelIdx = depth - 1;
    for (int i = 0; i < MAX_CHILDREN; i++) {
        if (children[i] == L_NULL) continue;
        flattenRecurse(children[i], depth - 1, affMul(xform, LEVEL_TRANSFORMS[levelIdx][i]));
    }
}

/* ── edge adjacency (int32 keys, compact hash table) ─────────── */

#define EDGE_HASH_SIZE (1 << 20)  /* 1048576 slots */
#define EDGE_HASH_MASK (EDGE_HASH_SIZE - 1)

typedef struct {
    int ax, ay, bx, by;  /* int32 keys (screen coords * 1000) */
    int tileIndex;
    int faceIndex;
    int occupied;
    int _pad;             /* align to 32 bytes */
} EdgeSlot;               /* 32 bytes each */

static EdgeSlot *edgeTable;

static unsigned int edgeHash(int ax, int ay, int bx, int by) {
    unsigned int h = (unsigned int)ax * 73856093U;
    h ^= (unsigned int)ay * 19349663U;
    h ^= (unsigned int)bx * 83492791U;
    h ^= (unsigned int)by * 52428047U;
    h ^= h >> 16;
    h *= 0x45d9f3bU;
    return h & EDGE_HASH_MASK;
}

/* Build adjacency from raw double-precision world coords (rawPx/rawPy),
   NOT from the float32 screen coords. This preserves the precision needed
   for edge matching, just like the V1 JS version (which uses doubles). */
static void buildAdjacencyFromRaw(void) {
    for (int i = 0; i < EDGE_HASH_SIZE; i++)
        edgeTable[i].occupied = 0;

    for (int t = 0; t < tileCount; t++)
        for (int f = 0; f < POINTS_PER; f++)
            tiles[t].neighbor[f] = -1;

    for (int t = 0; t < tileCount; t++) {
        int base = t * POINTS_PER;
        for (int f = 0; f < POINTS_PER; f++) {
            int f2 = (f + 1) % POINTS_PER;
            /* Use raw double world coords, rounded to *1000 (matching V1) */
            double dax = rawPx[base + f]  * 1000.0;
            double day = rawPy[base + f]  * 1000.0;
            double dbx = rawPx[base + f2] * 1000.0;
            double dby = rawPy[base + f2] * 1000.0;
            int ax = (int)(dax + (dax >= 0.0 ? 0.5 : -0.5));
            int ay = (int)(day + (day >= 0.0 ? 0.5 : -0.5));
            int bx = (int)(dbx + (dbx >= 0.0 ? 0.5 : -0.5));
            int by = (int)(dby + (dby >= 0.0 ? 0.5 : -0.5));

            int cax, cay, cbx, cby;
            if (ax < bx || (ax == bx && ay <= by)) {
                cax = ax; cay = ay; cbx = bx; cby = by;
            } else {
                cax = bx; cay = by; cbx = ax; cby = ay;
            }

            unsigned int slot = edgeHash(cax, cay, cbx, cby);
            for (int probe = 0; probe < 128; probe++) {
                unsigned int idx = (slot + probe) & EDGE_HASH_MASK;
                if (!edgeTable[idx].occupied) {
                    edgeTable[idx].ax = cax;
                    edgeTable[idx].ay = cay;
                    edgeTable[idx].bx = cbx;
                    edgeTable[idx].by = cby;
                    edgeTable[idx].tileIndex = t;
                    edgeTable[idx].faceIndex = f;
                    edgeTable[idx].occupied = 1;
                    break;
                }
                if (edgeTable[idx].ax == cax && edgeTable[idx].ay == cay &&
                    edgeTable[idx].bx == cbx && edgeTable[idx].by == cby) {
                    tiles[t].neighbor[f] = edgeTable[idx].tileIndex;
                    tiles[edgeTable[idx].tileIndex].neighbor[edgeTable[idx].faceIndex] = t;
                    break;
                }
            }
        }
    }
}

/* ── pseudo-random ───────────────────────────────────────────── */

/* Match V1's hash: hash(n) = fract(sin(n * 127.1 + 311.7) * 43758.5453123) */
static double hashRand(double n) {
    double x = js_sin(n * 127.1 + 311.7) * 43758.5453123;
    x = x - floor_(x);
    if (x < 0.0) x += 1.0;
    return x;
}

/* ── public API ──────────────────────────────────────────────── */

void setConfig(int depth, int birthFaces, int surviveMin, int surviveMax,
               int faceWeight, double density, int root, double seed,
               double width, double height) {
    cfgDepth      = depth < 1 ? 1 : (depth > MAX_DEPTH ? MAX_DEPTH : depth);
    cfgBirthFaces = birthFaces < 0 ? 0 : (birthFaces > 14 ? 14 : birthFaces);
    cfgSurviveMin = surviveMin < 0 ? 0 : (surviveMin > 14 ? 14 : surviveMin);
    cfgSurviveMax = surviveMax < 0 ? 0 : (surviveMax > 14 ? 14 : surviveMax);
    if (cfgSurviveMin > cfgSurviveMax) {
        int tmp = cfgSurviveMin; cfgSurviveMin = cfgSurviveMax; cfgSurviveMax = tmp;
    }
    cfgFaceWeight = faceWeight < 1 ? 1 : (faceWeight > 4 ? 4 : faceWeight);
    cfgDensity    = density < 0.0 ? 0.0 : (density > 1.0 ? 1.0 : density);
    cfgRoot       = root < 0 ? 0 : (root > 3 ? 3 : root);
    cfgSeed       = seed;
    screenW       = width;
    screenH       = height;
}

int buildTiling(void) {
    heap_reset();

    /* allocate working memory */
    tiles     = (Tile *)zalloc(MAX_TILES * (int)sizeof(Tile));
    rawPx     = (double *)alloc(MAX_TILES * POINTS_PER * (int)sizeof(double));
    rawPy     = (double *)alloc(MAX_TILES * POINTS_PER * (int)sizeof(double));
    edgeTable = (EdgeSlot *)zalloc(EDGE_HASH_SIZE * (int)sizeof(EdgeSlot));

    buildAllLevelTransforms(cfgDepth);

    int rootLabels[] = { L_DELTA, L_GAMMA, L_PSI, L_SIGMA };
    rawCount = 0;
    flattenRecurse(rootLabels[cfgRoot], cfgDepth, affIdent());
    tileCount = rawCount;

    /* bounds */
    double minX = 1e18, minY = 1e18, maxX = -1e18, maxY = -1e18;
    for (int i = 0; i < tileCount * POINTS_PER; i++) {
        if (rawPx[i] < minX) minX = rawPx[i];
        if (rawPy[i] < minY) minY = rawPy[i];
        if (rawPx[i] > maxX) maxX = rawPx[i];
        if (rawPy[i] > maxY) maxY = rawPy[i];
    }
    double boundsW = maxX - minX; if (boundsW < 1e-6) boundsW = 1.0;
    double boundsH = maxY - minY; if (boundsH < 1e-6) boundsH = 1.0;

    double margin = 28.0;
    fitScale = (screenW - margin*2) / boundsW;
    { double sy = (screenH - margin*2) / boundsH; if (sy < fitScale) fitScale = sy; }
    fitOx = (screenW - boundsW * fitScale) * 0.5 - minX * fitScale;
    fitOy = (screenH - boundsH * fitScale) * 0.5 - minY * fitScale;

    for (int t = 0; t < tileCount; t++) {
        int base = t * POINTS_PER;
        for (int i = 0; i < POINTS_PER; i++) {
            tiles[t].px[i] = (float)(rawPx[base + i] * fitScale + fitOx);
            tiles[t].py[i] = (float)(rawPy[base + i] * fitScale + fitOy);
        }
        tiles[t].alive = hashRand((double)t * 0.173 + 10.7) < cfgDensity ? 1 : 0;
        tiles[t].age = tiles[t].alive ? 1 : 0;
        tiles[t].nextAlive = 0;
        tiles[t].liveFaces = 0;
    }

    buildAdjacencyFromRaw();
    generation = 0;
    return tileCount;
}

void stepCA(void) {
    for (int t = 0; t < tileCount; t++) {
        int seen[16], counts[16], seenCount = 0;
        for (int f = 0; f < POINTS_PER; f++) {
            int neighborIdx = tiles[t].neighbor[f];
            if (neighborIdx < 0 || !tiles[neighborIdx].alive) continue;
            int found = -1;
            for (int s = 0; s < seenCount; s++)
                if (seen[s] == neighborIdx) { found = s; break; }
            if (found >= 0) counts[found]++;
            else if (seenCount < 16) { seen[seenCount] = neighborIdx; counts[seenCount] = 1; seenCount++; }
        }
        int total = 0;
        for (int s = 0; s < seenCount; s++) {
            int c = counts[s] > cfgFaceWeight ? cfgFaceWeight : counts[s];
            total += c;
        }
        tiles[t].liveFaces = total;
        tiles[t].nextAlive = tiles[t].alive
            ? (total >= cfgSurviveMin && total <= cfgSurviveMax ? 1 : 0)
            : (total == cfgBirthFaces ? 1 : 0);
    }
    for (int t = 0; t < tileCount; t++) {
        tiles[t].alive = tiles[t].nextAlive;
        if (tiles[t].alive) {
            int a = tiles[t].age + 1; tiles[t].age = a > 24 ? 24 : a;
        } else {
            int a = tiles[t].age - 1; tiles[t].age = a < 0 ? 0 : a;
        }
    }
    generation++;
}

void reseed(double seed) {
    cfgSeed = seed;
    for (int t = 0; t < tileCount; t++) {
        tiles[t].alive = hashRand((double)t * 1.618 + seed * 0.73 + 5.0) < cfgDensity ? 1 : 0;
        tiles[t].age = tiles[t].alive ? 1 : 0;
    }
    generation = 0;
}

void toggleTile(int index) {
    if (index < 0 || index >= tileCount) return;
    tiles[index].alive = tiles[index].alive ? 0 : 1;
    tiles[index].age = tiles[index].alive ? 4 : 0;
}

int   getTileCount(void)        { return tileCount; }
int   getGeneration(void)       { return generation; }
float getTilePx(int t, int i)   { return (t >= 0 && t < tileCount && i >= 0 && i < POINTS_PER) ? tiles[t].px[i] : 0.0f; }
float getTilePy(int t, int i)   { return (t >= 0 && t < tileCount && i >= 0 && i < POINTS_PER) ? tiles[t].py[i] : 0.0f; }
int   getTileAlive(int t)       { return (t >= 0 && t < tileCount) ? tiles[t].alive : 0; }
int   getTileAge(int t)         { return (t >= 0 && t < tileCount) ? tiles[t].age : 0; }
int   getTileLiveFaces(int t)   { return (t >= 0 && t < tileCount) ? tiles[t].liveFaces : 0; }
int   getTilesPtr(void)         { return (int)(unsigned long)tiles; }
int   getTileStride(void)       { return (int)sizeof(Tile); }
int   getRawPxPtr(void)         { return (int)(unsigned long)rawPx; }
int   getRawPyPtr(void)         { return (int)(unsigned long)rawPy; }

int getAliveCount(void) {
    int count = 0;
    for (int t = 0; t < tileCount; t++) if (tiles[t].alive) count++;
    return count;
}

int getFaceMin(void) {
    int lo = 999;
    for (int t = 0; t < tileCount; t++) if (tiles[t].liveFaces < lo) lo = tiles[t].liveFaces;
    return lo < 999 ? lo : 0;
}

int getFaceMax(void) {
    int hi = -1;
    for (int t = 0; t < tileCount; t++) if (tiles[t].liveFaces > hi) hi = tiles[t].liveFaces;
    return hi > -1 ? hi : 0;
}