Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Decay-based count by Qrator
/* Copyright (c) 2017, Qrator Labs */
/* All rights reserved. */
/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions are met: */
/* 1. Redistributions of source code must retain the above copyright notice, this */
/* list of conditions and the following disclaimer. */
/* 2. Redistributions in binary form must reproduce the above copyright notice, */
/* this list of conditions and the following disclaimer in the documentation */
/* and/or other materials provided with the distribution. */
/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND */
/* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED */
/* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
/* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR */
/* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES */
/* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND */
/* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT */
/* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS */
/* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
/* The views and conclusions contained in the software and documentation are those */
/* of the authors and should not be interpreted as representing official policies, */
/* either expressed or implied, of the FreeBSD Project. */
#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);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.