/*
 * bench.c - microbenchmarks for the PowerPC 604e under AIX 4.3.3.
 *
 * Portable C89, libm only. Prints CSV to stdout.
 *
 *   ./bench alu    integer / floating-point / branch microbenchmarks
 *   ./bench mem    memory hierarchy sweep
 *
 * ALU rows:   op,work,seconds,work_per_sec
 * MEM rows:   op,size_bytes,seconds,rate
 *               seqread/seqwrite -> rate is MB/s
 *               chase            -> rate is ns per access
 *
 * Integer and basic FP benchmarks are dependent chains (they measure
 * latency); the transcendental benchmarks are independent (throughput).
 *
 * Author: DeepSeek (deepseek-v4.1-flash), October 2026.
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>

static double   gsink;            /* written once, prevents dead code */
static unsigned utbl[256];
static unsigned btbl[256];
static double   dtbl[256];
static int      g_reps = 5;
static double   g_target = 0.20;

static double now(void)
{
    return (double)clock() / (double)CLOCKS_PER_SEC;
}

static void init_tables(void)
{
    unsigned s = 12345u;
    int i;
    srand(1);
    for (i = 0; i < 256; i++) {
        s = s * 1103515245u + 12345u;
        utbl[i] = (s >> 7) | 1u;                 /* odd -> no /0, no *0 */
        btbl[i] = (s >> 3);
        s = s * 1103515245u + 12345u;
        dtbl[i] = 1.0 + (double)((s >> 8) & 0xffffu) / 65536.0; /* [1,2) */
    }
}

/* ---- integer chains (operand from a table, dependent on previous) ---- */

#define INT_CHAIN(fname, OP) \
    static void fname(long n) { unsigned x = 1u; long i; \
        for (i = 0; i < n; i++) x = x OP utbl[i & 255]; \
        gsink += (double)x; }

INT_CHAIN(b_iadd, +)
INT_CHAIN(b_isub, -)
INT_CHAIN(b_imul, *)
INT_CHAIN(b_idiv, /)
INT_CHAIN(b_imod, %)
INT_CHAIN(b_iand, &)
INT_CHAIN(b_ior,  |)
INT_CHAIN(b_ixor, ^)

static void b_ishift(long n)
{
    unsigned x = 1u; long i;
    for (i = 0; i < n; i++) x = (x << 1) ^ (x >> 31) ^ utbl[i & 255];
    gsink += (double)x;
}

/* 64-bit (long long is 64-bit on AIX 4.3 PowerPC) */
static void b_lladd(long n)
{
    unsigned long long x = (unsigned long long)1; long i;
    for (i = 0; i < n; i++) x = x + (unsigned long long)utbl[i & 255];
    gsink += (double)(x & 0xffff);
}

static void b_llmul(long n)
{
    unsigned long long x = (unsigned long long)1; long i;
    for (i = 0; i < n; i++)
        x = x * (unsigned long long)utbl[i & 255] + (unsigned long long)1;
    gsink += (double)(x & 0xffff);
}

/* ---- double chains ---- */

#define DBL_CHAIN(fname, OP) \
    static void fname(long n) { double x = 1.5; long i; \
        for (i = 0; i < n; i++) x = x OP dtbl[i & 255]; \
        gsink += x; }

DBL_CHAIN(b_dadd, +)
DBL_CHAIN(b_dsub, -)
DBL_CHAIN(b_dmul, *)
DBL_CHAIN(b_ddiv, /)

/* ---- float chains ---- */
#define FLT_CHAIN(fname, OP) \
    static void fname(long n) { float x = 1.5f; long i; \
        for (i = 0; i < n; i++) x = x OP (float)dtbl[i & 255]; \
        gsink += (double)x; }

FLT_CHAIN(b_fadd, +)
FLT_CHAIN(b_fmul, *)
FLT_CHAIN(b_fdiv, /)

/* ---- transcendental chains (convergent, so values stay in range) ---- */
static void b_dsqrt(long n){ double x=1.5; long i; for(i=0;i<n;i++) x=sqrt(x); gsink+=x; }
static void b_dlog (long n){ double x=1.0; long i; for(i=0;i<n;i++) x=log(x+2.0); gsink+=x; }
static void b_dexp (long n){ double x=0.5; long i; for(i=0;i<n;i++) x=exp(x-1.0); gsink+=x; }
static void b_dsin (long n){ double x=1.0; long i; for(i=0;i<n;i++) x=sin(x); gsink+=x; }
static void b_dcos (long n){ double x=1.0; long i; for(i=0;i<n;i++) x=cos(x); gsink+=x; }
static void b_dpowy(long n){ double x=1.5; long i; for(i=0;i<n;i++) x=pow(x,1.0000001); gsink+=x; }

/* ---- branches ---- */
static void b_brpred(long n){ unsigned x=0; long i; for(i=0;i<n;i++){ if(i<2000000000u) x+=3; else x-=1; } gsink+=x; }
static void b_brrand(long n){ unsigned x=0; long i; for(i=0;i<n;i++){ if(btbl[i&255]&1u) x+=3; else x-=1; } gsink+=x; }

typedef struct {
    const char *name;
    double wpi;              /* work units per iteration */
    void (*fn)(long);
} Bench;

static void run(const Bench *b)
{
    long n = 1024;
    double t0, t1, secs, best = 1e30;
    int r;

    for (;;) {
        t0 = now(); b->fn(n); t1 = now(); secs = t1 - t0;
        if (secs >= 0.02 || n > (1L << 28)) break;
        n *= 4;
    }
    if (secs > 0.0) {
        double f = g_target / secs;
        if (f > 1.0) n = (long)((double)n * f);
    }
    if (n < 1024) n = 1024;

    for (r = 0; r < g_reps; r++) {
        t0 = now(); b->fn(n); t1 = now(); secs = t1 - t0;
        if (secs < best) best = secs;
    }
    printf("%s,%.0f,%.6f,%.0f\n", b->name, b->wpi * (double)n, best,
           (b->wpi * (double)n) / best);
}

static void run_alu(void)
{
    static const Bench table[] = {
        { "iadd",      1.0, b_iadd   },
        { "isub",      1.0, b_isub   },
        { "imul",      1.0, b_imul   },
        { "idiv",      1.0, b_idiv   },
        { "imod",      1.0, b_imod   },
        { "iand",      1.0, b_iand   },
        { "ior",       1.0, b_ior    },
        { "ixor",      1.0, b_ixor   },
        { "ishift",    1.0, b_ishift },
        { "lladd",     1.0, b_lladd  },
        { "llmul",     1.0, b_llmul  },
        { "dadd",      1.0, b_dadd   },
        { "dsub",      1.0, b_dsub   },
        { "dmul",      1.0, b_dmul   },
        { "ddiv",      1.0, b_ddiv   },
        { "fadd",      1.0, b_fadd   },
        { "fmul",      1.0, b_fmul   },
        { "fdiv",      1.0, b_fdiv   },
        { "dsqrt",     1.0, b_dsqrt  },
        { "dlog",      1.0, b_dlog   },
        { "dexp",      1.0, b_dexp   },
        { "dsin",      1.0, b_dsin   },
        { "dcos",      1.0, b_dcos   },
        { "dpow",      1.0, b_dpowy  },
        { "brpred",    1.0, b_brpred },
        { "brrand",    1.0, b_brrand }
    };
    int i;
    for (i = 0; i < (int)(sizeof(table) / sizeof(table[0])); i++)
        run(&table[i]);
}

static void run_mem(void)
{
    static const long sizes[] = {
        4096L, 16384L, 32768L, 65536L, 131072L, 524288L,
        1048576L, 2097152L, 8388608L, 16777216L
    };
    int si;
    int nsz = (int)(sizeof(sizes) / sizeof(sizes[0]));

    for (si = 0; si < nsz; si++) {
        long bytes = sizes[si];
        long ints = bytes / 4;
        long i, r, passes;
        int *a = (int *)malloc((size_t)bytes);
        double t0, t1, best;

        if (a == (int *)0) { continue; }
        for (i = 0; i < ints; i++) a[i] = (int)i;

        passes = (32L * 1024L * 1024L) / bytes;
        if (passes < 1) passes = 1;

        best = 1e30;
        for (r = 0; r < g_reps; r++) {
            long p;
            unsigned s = 0;
            t0 = now();
            for (p = 0; p < passes; p++)
                for (i = 0; i < ints; i++) s += (unsigned)a[i];
            t1 = now();
            gsink += (double)s;
            if (t1 - t0 < best) best = t1 - t0;
        }
        printf("seqread,%ld,%.6f,%.1f\n", bytes, best,
               (double)passes * (double)bytes / best / 1e6);

        best = 1e30;
        for (r = 0; r < g_reps; r++) {
            long p;
            t0 = now();
            for (p = 0; p < passes; p++)
                for (i = 0; i < ints; i++) a[i] = (int)i;
            t1 = now();
            if (t1 - t0 < best) best = t1 - t0;
        }
        printf("seqwrite,%ld,%.6f,%.1f\n", bytes, best,
               (double)passes * (double)bytes / best / 1e6);

        free(a);
    }

    {
        static const long ns[] = { 4096L, 16384L, 65536L, 262144L,
                                   1048576L, 4194304L };
        int k;
        int nn = (int)(sizeof(ns) / sizeof(ns[0]));
        for (k = 0; k < nn; k++) {
            long N = ns[k];
            long *idx = (long *)malloc((size_t)N * sizeof(long));
            long i, r, reps, total;
            double t0, t1, best = 1e30;

            if (idx == (long *)0) continue;
            for (i = 0; i < N; i++) idx[i] = i;
            for (i = N - 1; i > 0; i--) {           /* Sattolo: one cycle */
                long j = (long)(rand() % (int)i);
                long t = idx[i]; idx[i] = idx[j]; idx[j] = t;
            }
            reps = 2000000L / N;
            if (reps < 1) reps = 1;
            total = reps * N;
            for (r = 0; r < g_reps; r++) {
                long chain = 0, it;
                t0 = now();
                for (it = 0; it < total; it++) chain = idx[chain];
                t1 = now();
                gsink += (double)chain;
                if (t1 - t0 < best) best = t1 - t0;
            }
            printf("chase,%ld,%.6f,%.2f\n", N * 8, best,
                   best / (double)total * 1e9);
            free(idx);
        }
    }
}

int main(int argc, char **argv)
{
    const char *mode = (argc > 1) ? argv[1] : "alu";

    init_tables();

    if (strcmp(mode, "mem") == 0)
        run_mem();
    else
        run_alu();

    fflush(stdout);
    return 0;
}