/*
 * cc randread.c -o ~/bin/randread -O2 -lm
 *
 * randread device [bufsize:512 [range%:90 [nprocs:32]]]
 *
 * requires TSC
 */
#include <sys/types.h>
#include <sys/sysctl.h>
#include <sys/stat.h>
#include <sys/file.h>
#include <sys/mman.h>
#include <sys/errno.h>
#include <sys/wait.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <assert.h>
#include <machine/atomic.h>
#include <machine/cpufunc.h>

typedef struct pdata {
        int64_t counter;
        int64_t lotime;
        int64_t hitime;
        int64_t tsc_total1;
        int64_t tsc_total2;
        int64_t unused00;
        int64_t unused01;
        int     unused02;
        int     reset;
} pdata_t;

int
main(int ac, char **av)
{
    char *buf;
    size_t bytes = 512;
    off_t limit;
    int fd;
    int i;
    int loops;
    int nprocs = 32;
    double range = 90.0;
    volatile pdata_t *pdata;
    int64_t tsc1;
    int64_t tsc2;
    int64_t delta;
    int64_t tscfreq = 0;
    int64_t lotime;
    int64_t hitime;
    size_t tscfreq_size = sizeof(tscfreq);

    sysctlbyname("hw.tsc_frequency", &tscfreq, &tscfreq_size, NULL, 0);
    assert(tscfreq != 0);

    if (ac < 2 || ac > 5) {
        fprintf(stderr, "%s <device> [bufsize:512 [range%:90 [nprocs:32]]]\n",
                av[0]);
        exit (1);
    }

    if (ac >= 3) {
        bytes = (size_t)strtoul(av[2], NULL, 0);
        if (bytes < 512 || (bytes ^ (bytes - 1)) != ((bytes << 1) - 1)) {
            fprintf(stderr, "bytes must be a power of 2 >= 512\n");
            exit (1);
        }
    }
    buf = malloc(bytes);

    if (ac >= 4) {
        range = strtod(av[3], NULL);
    }

    if (ac >= 5) {
        nprocs = strtol(av[4], NULL, 0);
        if (nprocs < 0 || nprocs > 512) {
            fprintf(stderr, "absurd nprocs (%d)\n", nprocs);
            exit(1);
        }
    }

    fd = open(av[1], O_RDONLY);
    if (fd < 0) {
        fprintf(stderr, "open %s: %s\n", av[1], strerror(errno));
        exit (1);
    }

    lseek(fd, 0L, 2);
    limit = lseek(fd, 0L, 1);
    limit = (off_t)((double)limit * range / 100.0);
    limit &= ~(off_t)(bytes - 1);
    printf("device %s bufsize %zd limit %4.3fGB nprocs %d\n",
        av[1], bytes, (double)limit / (1024.0*1024.0*1024.0), nprocs);

    pdata = mmap(NULL, nprocs * sizeof(*pdata), PROT_READ|PROT_WRITE,
                    MAP_SHARED|MAP_ANON, -1, 0);

    for (i = 0; i < nprocs; ++i) {
        if (fork() == 0) {
            close(fd);
            fd = open(av[1], O_RDONLY);
            srandomdev();
            pdata += i;

            tsc2 = rdtsc();
            pdata->lotime = 0x7FFFFFFFFFFFFFFFLL;

            for (;;) {
                long pos;

                if (pdata->reset) {
                        pdata->counter = 0;
                        pdata->tsc_total1 = 0;
                        pdata->tsc_total2 = 0;
                        pdata->lotime = 0x7FFFFFFFFFFFFFFFLL;
                        pdata->hitime = 0;
                        pdata->reset = 0;
                }

                pos = random() ^ ((long)random() << 31);
                pos &= 0x7FFFFFFFFFFFFFFFLLU;
                pos = (pos % limit) & ~(off_t)(bytes - 1);
                lseek(fd, pos, 0);
                read(fd, buf, bytes);
                tsc1 = tsc2;
                tsc2 = rdtsc();
                delta = tsc2 - tsc1;
                ++pdata->counter;
                pdata->tsc_total1 += delta;
                pdata->tsc_total2 += delta * delta;
                if (pdata->lotime > delta)
                        pdata->lotime = delta;
                if (pdata->hitime < delta)
                        pdata->hitime = delta;
            }
        }
    }

    tsc2 = rdtsc();
    loops = 0;

    for (;;) {
        int64_t count;
        int64_t total1;
        int64_t total2;
        double v;
        double lo;
        double hi;
        double s1;
        double s2;
        double stddev;

        sleep(1);
        lotime = pdata[0].lotime;
        hitime = pdata[0].hitime;
        total1 = 0;
        total2 = 0;
        count = 0;

        for (i = 0; i < nprocs; ++i) {
                count += pdata[i].counter;
                total1 += pdata[i].tsc_total1;
                total2 += pdata[i].tsc_total2;
                if (lotime > pdata[i].lotime)
                        lotime = pdata[i].lotime;
                if (hitime < pdata[i].hitime)
                        hitime = pdata[i].hitime;
                pdata[i].reset = 1;
        }
        tsc1 = tsc2;
        tsc2 = rdtsc();
        delta = tsc2 - tsc1;
        v = count * ((double)delta / (double)tscfreq);
        lo = (double)lotime / (double)tscfreq;
        hi = (double)hitime / (double)tscfreq;

        s1 = ((double)total2 - (double)total1 * (double)total1 / (double)count) / ((double)count - 1);
        if (s1 < 0.0)
                stddev = -sqrt(-s1);
        else
                stddev = sqrt(s1);
        stddev = stddev / (double)tscfreq;      /* normalize to 1 second units */

        if (loops) {
                printf("%6.0f/s avg=%6.2fuS bw=%-6.2fMB/s "
                       "lo=%-3.2fuS, hi=%-3.2fuS stddev=%3.2fuS\n",
                       v,
                       1e6 * nprocs / v,
                       (double)count * bytes / 1e6 / ((double)delta / (double)tscfreq),
                       lo * 1e6,
                       hi * 1e6,
                       stddev * 1e6);
        }
        ++loops;
    }
    return 0;
}