/**
 * HPC Challenge Class 2 -- Global Random Access
 *
 * IBM Research
 * (C) Copyright IBM Corp. 2005, All Rights Reserved
 *
 */
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include <sys/time.h>
#include <time.h>

typedef unsigned long long u64Int;
typedef long long s64Int;

#ifndef N
  #define N         (21)
#endif
#define TableSize (1ULL<<N)
#define NUPDATE   (4ULL * TableSize)
#define POLY      0x0000000000000007ULL
#define PERIOD    1317624576693539401LL

shared u64Int *Table;
shared unsigned verifyfailed;

double    mysecond();
u64Int    starts (s64Int);
void      RandomAccessUpdate();
void      RandomAccessVerify();

/* *********************************************************************** */
/*                     main program                                        */
/* *********************************************************************** */
int main() 
{
  u64Int i;
  Table = (shared u64Int *)upc_all_alloc(TableSize,sizeof(u64Int));
  assert(Table != 0);
  double time = 0;
  double GUPs = 0;
  verifyfailed = 0;

  if(MYTHREAD == 0) 
	printf("Table size = %llu MBytes/CPU, %llu MB/total on %d threads\n",
	       TableSize*8/1024/1024/THREADS, 
	       TableSize*8/1024/1024,
               THREADS);

  /* ------------------- */
  /* begin timed portion */
  /* ------------------- */
  upc_barrier;
  time = mysecond();
  upc_forall (i = 0; i < TableSize; i++; i) Table[i] = i;
  upc_barrier;
  RandomAccessUpdate();
  upc_barrier;
  time = mysecond() - time;
  GUPs = ((double) 1e-9 * NUPDATE) / time;

  /* ------------------- */
  /* end timed portion   */
  /* ------------------- */
  if (MYTHREAD == 0) printf ("End timed portion: verifying\n");
  RandomAccessUpdate();                 /* do it again */
  upc_barrier;
  RandomAccessVerify();
  upc_barrier;

  if(MYTHREAD == 0) 
    {
      if (verifyfailed == 0) printf ("Verification: SUCCESS\n");
      else                   printf ("Verification *FAILED*, %d\n", verifyfailed);
      printf("Default number of updates (RECOMMENDED) = %llu\n", NUPDATE);
      printf("Real time used = %.6f seconds\n", time );
      printf("%.9f Billion(10^9) Updates    per second [GUP/s]\n", GUPs);
    }
  upc_free (Table);
  return 0;
}

/* *********************************************************************** */
/*                principal loop (that does all the work)                  */
/* *********************************************************************** */

void RandomAccessUpdate()
{
  u64Int i, ran = starts(NUPDATE/THREADS * MYTHREAD);
  upc_forall (i = 0; i < NUPDATE; i++; i) 
    {
      ran = (ran << 1) ^ (((s64Int) ran < 0) ? POLY : 0);
      Table[ran & (TableSize-1)] ^= ran;
    }
}

/* *********************************************************************** */
/*                check loop                                               */
/* *********************************************************************** */

void RandomAccessVerify()
{
  u64Int i;
  upc_forall (i = 0; i < TableSize; i++; i) 
    if (Table[i] != i)
      { verifyfailed++; printf ("Failed at index=%lld: %lld\n", i, Table[i]); return; }
}

/* ************************************************************ */
/* Utility routine to start random number generator at Nth step */
/* ************************************************************ */

u64Int starts(s64Int n)
{
  /* s64Int i, j; */
  int i, j;
  u64Int m2[64];
  u64Int temp, ran;

  while (n < 0)         n += PERIOD;
  while (n > PERIOD)    n -= PERIOD;

  if (n == 0)           return 0x1;

  temp = 0x1;
  for (i=0; i<64; i++)
    {
      m2[i] = temp;
      temp = (temp << 1) ^ ((s64Int) temp < 0 ? POLY : 0);
      temp = (temp << 1) ^ ((s64Int) temp < 0 ? POLY : 0);
    }

  for (i=62; i>=0; i--)
    if ((n >> i) & 1)
      break;

  ran = 0x2;
  while (i > 0)
    {
      temp = 0;
      for (j=0; j<64; j++)
        if ((ran >> j) & 1)
          temp ^= m2[j];
      ran = temp;
      i -= 1;
      if ((n >> i) & 1)
        ran = (ran << 1) ^ ((s64Int) ran < 0 ? POLY : 0);
    }

  return ran;
}

double mysecond()
{
  struct timeval tv;
  gettimeofday(&tv, 0);
  return tv.tv_sec + ((double) tv.tv_usec / 1000000);
}