/* time_fadd.c  try to measure time of double  A = A + B;         */
/*              roughly time of one floating point add            */
/*              using double difference and minimum and stability */

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

#define dabs(a) (a)<0.0?-(a):(a)
void do_count(int * count_check, int rep, double * B);

int main(void)
{
  double t1, t2, t3, t4, tmeas, t_min, t_prev, ts, tavg;
  double A, B, Q;
  int stable;
  int i, j;
  int count_check, outer;
  int rep, min_rep;


  t_min = 10.0;    /* 10.0 seconds typical minimum measurement time */
  Q  = 5.0;        /* 5.0 typical approximate percentage stability */
  min_rep = 32;    /* minimum of 32 typical */
  outer = 100000;  /* some big number */

  printf("time_fadd.c \n");
  printf("min time %g seconds, min stability %g percent, outer loop=%d\n",
         t_min, Q, outer);


  stable = 5; /* max tries */
  t_prev = 0.0;
  for(rep=min_rep; rep<100000; rep=rep+rep) /* increase until good measure */
  {
    A = 0.0;
    B = 0.1;
    t1 = (double)clock()/(double)CLOCKS_PER_SEC;
    for(i=0; i<outer; i++) /* outer control loop */
    {
      count_check = 0;
      for(j=0; j<rep; j++)   /* inner control loop */
      {
         do_count(&count_check, rep, &B);
      }
    }
    t2 = (double)clock()/(double)CLOCKS_PER_SEC;
    if(count_check != rep) printf("bad count_check_1 %d \n", count_check);

    A = 0.0;
    t3 = (double)clock()/(double)CLOCKS_PER_SEC;
    for(i=0; i<outer; i++) /* outer measurement loop */
    {
      count_check = 0;
      for(j=0; j<rep; j++)   /* inner measurement loop */
      {
         do_count(&count_check, rep, &B);
         A = A + B;   /* item being measured, approximately FADD time */
      }
    }
    t4 = (double)clock()/(double)CLOCKS_PER_SEC;
    if(count_check != rep) printf("bad count_check_2 %d \n", count_check);

    tmeas = (t4-t3) - (t2-t1); /* the double difference */
    printf("rep=%d, t measured=%g \n", rep, tmeas);

    if((t4-t3)<t_min) continue; /* need more rep */

    if(t_prev==0.0)
    {
      printf("tmeas=%g, t_prev=%g, rep=%d \n", tmeas, t_prev, rep);
      t_prev = tmeas;
    } 
    else /* compare to previous */
    {
      printf("tmeas=%g, t_prev=%g, rep=%d \n", tmeas, t_prev, rep);
      ts = 2.0*(dabs(tmeas-t_prev)/(tmeas+t_prev));
      tavg = 0.5*(tmeas+t_prev);
      if(100.0*ts < Q)  break; /* above minimum and stable */
      t_prev = tmeas;
    }
    stable--;
    if(stable==0) break;
    rep = rep/2; /* hold rep constant */
  } /* end loop increasing rep */

  /* stable? and over minimum */
  if(stable==0) printf("rep=%d  unstable \n", rep);
  if(tmeas<t_min) printf("time measured=%g, under minimum \n", tmeas);
  printf("raw time=%g, fadd time=%g, rep=%d, stable=%g\% \n\n", tmeas, 
         (tavg/(double)outer)/(double)rep, rep, ts);
  return 0;
} /* end time_fadd.c */

/* do_count to prevent dead code elimination  */
void do_count(int * count_check, int rep, double * B)
{
  (*count_check)++;
  /* could change B but probably don't have to. */
}