Decay-based count

/* Qratror Labs https://qrator.net/ */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <unistd.h>
#include <math.h>


static inline uint64_t rdtsc(void) {
  union {
    uint64_t tsc_64;
    struct {
      uint32_t lo_32;
      uint32_t hi_32;
    };
  } tsc;
  asm volatile("rdtsc" :
               "=a" (tsc.lo_32),
               "=d" (tsc.hi_32));
  return tsc.tsc_64;
}

/* global config and other global stuff */
struct {
  double tau;
  uint64_t tVanish;
  uint16_t *rho_table; /* NB: tau must be no more than (1<<16) */

#if !defined(DECAY_MODEL)
  double powers[64];
  double beta_1;
#endif
} g;

int ratedet_init(double tau);
void ratedet_finalize();

int main(int argc, char **argv) {
  int nIters = 100000000;
  int nRounds = 1;
  switch(argc) {
  case 3:
    nRounds = atoi(argv[2]);
  case 2:
    nIters = atoi(argv[1]);
  case 1:
    break;
  default:
    printf("usage: %s [nIters] [nRounds]\n", argv[0]);
    exit(1);
  }

  if(ratedet_init(1<<16))
    exit(1);

  printf("tau: %lf", g.tau);
#ifdef DECAY_MODEL
  printf(", tVanish: %lu; %u...%u", g.tVanish, g.rho_table[0], g.rho_table[g.tVanish-1]);
#else
#ifdef MAXDELTA
  printf(", maxDelta: %u, accuracy: %lg", MAXDELTA, exp(-MAXDELTA/g.tau)/g.beta_1);
#else
  printf(", maxDelta: <unused>");
#endif
#endif
  printf("\n");

  uint64_t t = rdtsc();
  sleep(1);
  t = rdtsc() - t;
  const double freq = (double)t;

  uint64_t minT = ~0UL;
  double sT = 0.;

  int iter, round;
  for(round=0; round<nRounds; ++round) {
#ifdef DECAY_MODEL
    uint64_t currTime = g.tVanish;
    uint64_t value = 0;
#else
    uint64_t currTime = 0;
    uint64_t accTime = 0;
    double accValue = 0.;
#endif
    t = rdtsc();
    uint64_t tStep = 1000;
    for(iter=0; iter<nIters; ++iter) {
      /* FIXME: rdtsc() is too slow for this test */
#if 0
      currTime = rdtsc();
#else
      currTime = currTime + tStep;
      tStep = (tStep + 123) & 0xfff;
#endif
#ifdef DECAY_MODEL
      /* update */
      {
        const int64_t dt = labs(value - currTime);
        const int64_t offset = dt > g.tVanish ? g.tVanish : dt;
        const uint64_t maxt = value < currTime ? currTime : value;
        value = maxt + g.rho_table[offset];
      }
#else
      /* update */
      {
        uint64_t delta = currTime - accTime;
        accTime = currTime;

#ifdef MAXDELTA
        if(delta >= MAXDELTA)
          {
            accValue = g.beta_1;
          }
        else
#endif
          {
            const double *ppow = g.powers;
            for (; delta != 0; ++ppow, delta >>= 1) {
              if (delta & 1)
                accValue *= *ppow;
            }
            accValue += g.beta_1;
          }
      }
#endif
    }
    t = rdtsc() - t;

    if(0 == round) {
      /* This should be approximately the same for both versions */
      const double total =
#ifdef DECAY_MODEL
        exp((value-currTime)/(double)g.tau)
#else
        accValue/g.beta_1
#endif
        ;
      printf("total value: %lf\n", total);
    }

    sT += t;
    if(t < minT) minT = t;
  }

  const double meanT = sT/nRounds;
  printf("freq: %g MHz, iter time (best/mean): %lf/%lf cycles, rate (best/mean): %g/%g Miter/s\n",
         freq*1e-6,
         (double)minT/nIters, (double)meanT/nIters,
         freq/(double)minT*nIters*1e-6, freq/(double)meanT*nIters*1e-6);

  ratedet_finalize();
  return 0;
}

static double rho(double x) {
  return log(1. + exp(-x));
}

int ratedet_init(double tau) {
  g.tau = tau;
  g.tVanish = - tau * log(exp(.5/tau) - 1.);
  g.rho_table = (uint16_t*)malloc(g.tVanish * sizeof(uint16_t));
  if(!g.rho_table) return -1;
  uint64_t i;
  for(i=0; i<g.tVanish; ++i)
    g.rho_table[i] = .5 + tau * rho((double)i/tau);

#if !defined(DECAY_MODEL)
  const double beta = exp(-1./(double)tau);
  g.beta_1 = 1 - beta;
  double tmp = beta;
  for(i=0; i<64;++i) {
    g.powers[i] = tmp;
    tmp *= tmp;
  }
#endif

  return 0;
}

void ratedet_finalize() {
  if(g.rho_table)
    free(g.rho_table);
}