/*
 * peak.c - double-precision floating-point throughput probes (PowerPC 604e).
 *
 * DeepSeek's bench.c measured FP latency with dependent chains; these
 * kernels measure throughput: NACC independent accumulator chains, wide
 * enough to hide the 604e's 3-cycle FMA latency. Theoretical peak is one
 * FMA retired per cycle = 2 FLOP/cycle = 500 MFLOPS at 250 MHz.
 *
 *   peak [target_seconds] [tag]   calibrate + run all kernels, CSV out
 *   peak bits                     FMA-contraction bit test and exit
 *
 * FMA contraction note: xlc -qarch=604 contracts a*b+c into fmadd (single
 * rounding); that is exactly what we want the FMA kernel to show. The
 * "bits" mode proves whether the contraction happened: with
 *   a = 1 + 2^-27, b = 1 - 2^-27, c = -2^-53
 * the unfused evaluation (volatile temp) gives 1 - 2^-53 while a fused
 * fmadd gives 1 - 2^-52.
 *
 * Portable C89, no libraries beyond libm.
 * Build on obsidian:  /usr/vac/bin/xlc -O2 -qarch=604 -qtune=604 -o peak peak.c -lm
 *
 * kimi-k3, October 2026
 */

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

#define NACC 8                  /* independent chains; 3 would do, 8 is sure */

typedef double (*kernel_t)(long);

/* ------------------------------------------------------------- kernels */
/* Each returns the folded accumulators so no loop is dead code. Multiplier
 * 1.0000001 grows by e^10 over 1e8 iterations: no overflow, no denormals. */

static double k_fma(long reps)
{
    double a0 = 1.0,        a1 = 1.0 + 1e-12, a2 = 1.0 + 2e-12, a3 = 1.0 + 3e-12;
    double a4 = 1.0 + 4e-12, a5 = 1.0 + 5e-12, a6 = 1.0 + 6e-12, a7 = 1.0 + 7e-12;
    const double m = 1.0000001, c = 1e-9;
    long i;
    for (i = 0; i < reps; i++) {
        a0 = a0 * m + c; a1 = a1 * m + c; a2 = a2 * m + c; a3 = a3 * m + c;
        a4 = a4 * m + c; a5 = a5 * m + c; a6 = a6 * m + c; a7 = a7 * m + c;
    }
    return a0 + a1 + a2 + a3 + a4 + a5 + a6 + a7;
}

static double k_fma_dep(long reps)      /* one chain: the latency reference */
{
    double a = 1.0;
    const double m = 1.0000001, c = 1e-9;
    long i, n = reps * NACC;            /* same statement count as the rest */
    for (i = 0; i < n; i++)
        a = a * m + c;
    return a;
}

static double k_mul(long reps)
{
    double a0 = 1.0,        a1 = 1.0 + 1e-12, a2 = 1.0 + 2e-12, a3 = 1.0 + 3e-12;
    double a4 = 1.0 + 4e-12, a5 = 1.0 + 5e-12, a6 = 1.0 + 6e-12, a7 = 1.0 + 7e-12;
    const double m = 1.0000001;
    long i;
    for (i = 0; i < reps; i++) {
        a0 *= m; a1 *= m; a2 *= m; a3 *= m;
        a4 *= m; a5 *= m; a6 *= m; a7 *= m;
    }
    return a0 + a1 + a2 + a3 + a4 + a5 + a6 + a7;
}

static double k_add(long reps)
{
    double a0 = 1.0,        a1 = 1.0 + 1e-12, a2 = 1.0 + 2e-12, a3 = 1.0 + 3e-12;
    double a4 = 1.0 + 4e-12, a5 = 1.0 + 5e-12, a6 = 1.0 + 6e-12, a7 = 1.0 + 7e-12;
    const double c = 1e-9;
    long i;
    for (i = 0; i < reps; i++) {
        a0 += c; a1 += c; a2 += c; a3 += c;
        a4 += c; a5 += c; a6 += c; a7 += c;
    }
    return a0 + a1 + a2 + a3 + a4 + a5 + a6 + a7;
}

static double k_div(long reps)
{
    double a0 = 1.0,        a1 = 1.0 + 1e-12, a2 = 1.0 + 2e-12, a3 = 1.0 + 3e-12;
    double a4 = 1.0 + 4e-12, a5 = 1.0 + 5e-12, a6 = 1.0 + 6e-12, a7 = 1.0 + 7e-12;
    const double d = 1.0000001;
    long i;
    for (i = 0; i < reps; i++) {
        a0 /= d; a1 /= d; a2 /= d; a3 /= d;
        a4 /= d; a5 /= d; a6 /= d; a7 /= d;
    }
    return a0 + a1 + a2 + a3 + a4 + a5 + a6 + a7;
}

static double k_sqrt(long reps)         /* libm call: measures what exists */
{
    double a0 = 1.0,        a1 = 1.0 + 1e-12, a2 = 1.0 + 2e-12, a3 = 1.0 + 3e-12;
    double a4 = 1.0 + 4e-12, a5 = 1.0 + 5e-12, a6 = 1.0 + 6e-12, a7 = 1.0 + 7e-12;
    const double c = 1e-9;
    long i;
    for (i = 0; i < reps; i++) {
        a0 = sqrt(a0 + c); a1 = sqrt(a1 + c); a2 = sqrt(a2 + c); a3 = sqrt(a3 + c);
        a4 = sqrt(a4 + c); a5 = sqrt(a5 + c); a6 = sqrt(a6 + c); a7 = sqrt(a7 + c);
    }
    return a0 + a1 + a2 + a3 + a4 + a5 + a6 + a7;
}

/* ------------------------------------------------------------- harness */

typedef struct {
    const char *name;
    kernel_t    fn;
    int         flop;           /* FLOP weight per statement: FMA = 2 */
} Kernel;

static const Kernel kernels[] = {
    { "fma",     k_fma,     2 },
    { "fma_dep", k_fma_dep, 2 },
    { "mul",     k_mul,     1 },
    { "add",     k_add,     1 },
    { "div",     k_div,     1 },
    { "sqrt",    k_sqrt,    1 }
};
#define NK (sizeof(kernels) / sizeof(kernels[0]))

static double guard;            /* sinks every kernel result; printed at end */

static double time_once(kernel_t fn, long reps)
{
    clock_t t0 = clock();
    guard += fn(reps);
    return (double)(clock() - t0) / (double)CLOCKS_PER_SEC;
}

/* AIX's clock() reports in CLOCKS_PER_SEC units but only advances at the
 * scheduler tick (CLK_TCK = 100 Hz on obsidian). A probe shorter than a few
 * ticks can read zero, so grow reps until the probe spans at least 50 ms
 * before trusting the scale factor. (First version skipped this and scaled
 * 20000 reps by 0.4/1e-6: 8e9 reps, found the hard way.) */
static long calibrate(kernel_t fn, double target)
{
    long reps = 20000;
    double t = 0.0;
    int tries;

    for (tries = 0; tries < 10; tries++) {
        t = time_once(fn, reps);
        if (t >= 0.05)
            break;
        reps *= 4;
    }
    if (t <= 0.0)
        t = 1e-6;
    reps = (long)((double)reps * (target / t));
    if (reps < 1000)
        reps = 1000;
    return reps;
}

/* ------------------------------------------------------------- bits test */

static void fma_bits(void)
{
    double a = 1.0 + 1.0 / 134217728.0;         /* 1 + 2^-27 */
    double b = 1.0 - 1.0 / 134217728.0;         /* 1 - 2^-27 */
    double c = -1.0 / 9007199254740992.0;       /* -2^-53  */
    volatile double t;                          /* forces a rounding point */
    double unfused, fused;
    union { double d; unsigned long u[2]; } x = { 0.0 }, y = { 0.0 };

    t = a * b;
    unfused = t + c;                            /* expect 1 - 2^-53 */
    fused = a * b + c;                          /* fmadd: expect 1 - 2^-52 */
    x.d = unfused;
    y.d = fused;
    printf("a=1+2^-27 b=1-2^-27 c=-2^-53\n");
    printf("unfused a*b+c: %08lx%08lx  %.17g\n", x.u[0], x.u[1], unfused);
    printf("fused   a*b+c: %08lx%08lx  %.17g\n", y.u[0], y.u[1], fused);
    if (x.u[0] == y.u[0] && x.u[1] == y.u[1])
        printf("identical: no FMA contraction by this compiler\n");
    else
        printf("differ: FMA contraction is active\n");
}

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

int main(int argc, char **argv)
{
    double target = 2.0;
    const char *tag = "run";
    size_t k;

    if (argc > 1 && strcmp(argv[1], "bits") == 0) {
        fma_bits();
        return 0;
    }
    if (argc > 1)
        target = atof(argv[1]);
    if (argc > 2)
        tag = argv[2];
    if (target < 0.2)
        target = 0.2;

    printf("# tag,kernel,reps,seconds,mstmt_s,mflop_s\n");
    guard = 0.0;
    for (k = 0; k < NK; k++) {
        long reps = calibrate(kernels[k].fn, target * 0.2);
        double best = 0.0, stmts;
        int r;
        for (r = 0; r < 3; r++) {
            double t = time_once(kernels[k].fn, reps);
            if (best == 0.0 || t < best)
                best = t;
        }
        stmts = (double)reps * NACC;
        printf("%s,%s,%ld,%.3f,%.1f,%.1f\n",
               tag, kernels[k].name, reps, best,
               stmts / best / 1e6,
               stmts * kernels[k].flop / best / 1e6);
    }
    printf("# guard %.17g\n", guard);
    return 0;
}