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

#define RAND_MAX      32767       
#define MAX_RANDOM    2147483647  
#define MAX_DOUBLE    2147483647.0  

plural long ma[56];
int ia, ib;

plural double prand()
{

  register plural long ra;
  plural double r;
 
  if ( ++ia == 56 ) ia = 1;
  if ( ++ib == 56 ) ib = 1;
  ra = ma[ia] - ma[ib];
  if ( ra < 0 ) ra += MAX_RANDOM;
  ma[ia] = ra;
  r = ra / MAX_DOUBLE;

  return r;

}

/**********************************************************************/

#define MAXSEED 10000

init_random( Seed )
int Seed;
{
  int               i, k, ii;
  plural long       mj, mk;

  /* Sequential */
  srand(Seed);

  /* Parallel */
  Seed = Seed % MAXSEED;
  for (i=0; i<Seed; i++) 
  {
     mj = p_random();
  }

  mj = p_random();
  ma[55] = mj;
  mk = 1;
  for ( i = 1; i <= 54; i++ )
  {
    ii = ( 21 * i ) % 55;
    ma[ii] = mk;
    mk = mj - mk;
    if ( mk < 0 ) mk += MAX_RANDOM;
    mj = ma[ii];
  }
  for ( k = 1; k <= 4; k++ )
  {
    for ( i = 1; i <= 55; i++ );
    {
      ma[i] -= ma[1+(i+30) % 55];
      if ( ma[i] < 0 ) ma[i] += MAX_RANDOM;
    }
  }
  ia = 0;
  ib = 31;
}



double d_rand()
{
 double temp;

 temp = (double)rand() / (double)RAND_MAX; 
 if (temp == 1.0) temp = 0.0;

 if ( (temp >1.0) || (temp < 0 ) ) {
     perror("ERROR: d_rand is bad");
     printf("d_rand = %lf \n",temp);
     fflush(stdout);
     exit(-2);
 }

  printf("rand = %f\n", temp);

 return(temp);

}

