# include <math.h>
# define RAND_MAX 2147483647

/*================================================================*/
/*                                                                */
/*       Random Number Generator:  R 2 5 0                        */
/*                                                                */
/*       program version 1.0 for C                                */
/*       Dieter W. Heermann                                       */
/*================================================================*/

int init_r250(seed, m_f_ptr)
   int seed;
   int *m_f_ptr;
   {

    int  i,tmp, dummy, one ;
    int  *start;

    start = m_f_ptr;
    srand(seed);
    one    = 1;

    /* warm up the usual random number generator */
    for (i=0; i< 100; i++)
       {dummy = rand();
        }

    /* now draw the 250 (251)initial bit sequences */
    for (i=0 ; i<251; i++)
       {*m_f_ptr++ = rand();
       }

    /* now orthogonalize as best as we can  */
    m_f_ptr = start;
    for (i=0; i < 30; i++)
       { tmp      = *m_f_ptr;
	 *m_f_ptr = tmp | one;
	 one      = one << 1;
	 m_f_ptr++;
       }

    return 0;
  }

int r250 (n, x_ptr, m_f_ptr, save)
   int    n;
   float  *x_ptr;
   int    *m_f_ptr;
   int    save;
   {

    int ind ;
    int j, min,k,ll;
    float *ran_ptr;

    ind     = save;
    ll      = n + 250;
    ran_ptr = x_ptr;
    j       = 1;

    do 
       { ll = ll - (251 - ind);
         if ( ll < 147 )  min = ll;
            else         min = 147;
         for(k=ind; k < min+1; k++)
           {
              *(m_f_ptr+ k)   = *(m_f_ptr +k) ^ *(m_f_ptr +k+103);
	      *ran_ptr = *(m_f_ptr+k);
	      *ran_ptr = *ran_ptr / RAND_MAX;
	      ran_ptr++;
              j++;
            }
         ind = k;
         if ( ll < 250 )   min = ll;
           else           min = 250;

         for( k= ind; k < min +1; k++)
            {
              *(m_f_ptr+ k)   = *(m_f_ptr +k) ^ *(m_f_ptr +k-147);
	      *ran_ptr = *(m_f_ptr+k);
	      *ran_ptr = *ran_ptr / RAND_MAX;
	      ran_ptr++;
              j++;
            }
          if ( k == 251 )   {k=1;ind = 1;};
        }  while( (ind == 1) & ( j < n));

   return k;
  }