test.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

    int status = 0;
    for (i = 0; i < N; ++i)
      {
        status |= (test_b[i] != test_a[i]);
      }
    gsl_test (status, "%s, random number state consistency",
              gsl_rng_name (r));
  }

  gsl_rng_free (r);
}


void
rng_parallel_state_test (const gsl_rng_type * T)
{
  unsigned long int test_a[N], test_b[N];
  unsigned long int test_c[N], test_d[N];
  double test_e[N], test_f[N];

  int i;

  gsl_rng *r1 = gsl_rng_alloc (T);
  gsl_rng *r2 = gsl_rng_alloc (T);

  for (i = 0; i < N; ++i)
    {
      gsl_rng_get (r1);         /* throw away N iterations */
    }

  gsl_rng_memcpy (r2, r1);              /* save the intermediate state */

  for (i = 0; i < N; ++i)
    {
      /* check that there is no hidden state intermixed between r1 and r2 */
      test_a[i] = gsl_rng_get (r1);     
      test_b[i] = gsl_rng_get (r2);
      test_c[i] = gsl_rng_uniform_int (r1, 1234);       
      test_d[i] = gsl_rng_uniform_int (r2, 1234);
      test_e[i] = gsl_rng_uniform (r1); 
      test_f[i] = gsl_rng_uniform (r2);
    }

  {
    int status = 0;
    for (i = 0; i < N; ++i)
      {
        status |= (test_b[i] != test_a[i]);
        status |= (test_c[i] != test_d[i]);
        status |= (test_e[i] != test_f[i]);
      }
    gsl_test (status, "%s, parallel random number state consistency",
              gsl_rng_name (r1));
  }

  gsl_rng_free (r1);
  gsl_rng_free (r2);

}

void
rng_read_write_test (const gsl_rng_type * T)
{
  unsigned long int test_a[N], test_b[N];

  int i;

  gsl_rng *r = gsl_rng_alloc (T);

  for (i = 0; i < N; ++i)
    {
      gsl_rng_get (r);  /* throw away N iterations */
    }

  { /* save the state to a binary file */
    FILE *f = fopen("test.dat", "wb"); 
    gsl_rng_fwrite(f, r);
    fclose(f);
  }

  for (i = 0; i < N; ++i)
    {
      test_a[i] = gsl_rng_get (r);
    }

  { /* read the state from a binary file */
    FILE *f = fopen("test.dat", "rb"); 
    gsl_rng_fread(f, r);
    fclose(f);
  }

  for (i = 0; i < N; ++i)
    {
      test_b[i] = gsl_rng_get (r);
    }

  {
    int status = 0;
    for (i = 0; i < N; ++i)
      {
        status |= (test_b[i] != test_a[i]);
      }
    gsl_test (status, "%s, random number generator read and write",
              gsl_rng_name (r));
  }

  gsl_rng_free (r);
}

void
generic_rng_test (const gsl_rng_type * T)
{
  gsl_rng *r = gsl_rng_alloc (T);
  const char *name = gsl_rng_name (r);
  unsigned long int kmax = 0, kmin = 1000;
  double sigma = 0;
  const unsigned long int ran_max = gsl_rng_max (r);
  const unsigned long int ran_min = gsl_rng_min (r);

  int status = rng_max_test (r, &kmax, ran_max);

  gsl_test (status,
            "%s, observed vs theoretical maximum (%lu vs %lu)",
            name, kmax, ran_max);

  status = rng_min_test (r, &kmin, ran_min, ran_max);

  gsl_test (status,
            "%s, observed vs theoretical minimum (%lu vs %lu)",
            name, kmin, ran_min);

  status = rng_sum_test (r, &sigma);

  gsl_test (status,
            "%s, sum test within acceptable sigma (observed %.2g sigma)",
            name, sigma);

  status = rng_bin_test (r, &sigma);

  gsl_test (status,
            "%s, bin test within acceptable chisq (observed %.2g sigma)",
            name, sigma);

  gsl_rng_set (r, 1);   /* set seed to 1 */
  status = rng_max_test (r, &kmax, ran_max);

  gsl_rng_set (r, 1);   /* set seed to 1 */
  status |= rng_min_test (r, &kmin, ran_min, ran_max);

  gsl_rng_set (r, 1);   /* set seed to 1 */
  status |= rng_sum_test (r, &sigma);

  gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
  status |= rng_max_test (r, &kmax, ran_max);

  gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
  status |= rng_min_test (r, &kmin, ran_min, ran_max);

  gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
  status |= rng_sum_test (r, &sigma);

  gsl_test (status, "%s, maximum and sum tests for non-default seeds", name);

  gsl_rng_free (r);
}

int
rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max)
{
  unsigned long int actual_uncovered;
  double expect_uncovered;
  int status;
  unsigned long int max = 0;
  int i;

  for (i = 0; i < N2; ++i)
    {
      unsigned long int k = gsl_rng_get (r);
      if (k > max)
        max = k;
    }

  *kmax = max;

  actual_uncovered = ran_max - max;
  expect_uncovered = (double) ran_max / (double) N2;

  status = (max > ran_max) || (actual_uncovered > 7 * expect_uncovered) ;

  return status;
}

int
rng_min_test (gsl_rng * r, unsigned long int *kmin, 
              unsigned long int ran_min, unsigned long int ran_max)
{
  unsigned long int actual_uncovered;
  double expect_uncovered;
  int status;
  unsigned long int min = 1000000000UL;
  int i;

  for (i = 0; i < N2; ++i)
    {
      unsigned long int k = gsl_rng_get (r);
      if (k < min)
        min = k;
    }

  *kmin = min;

  actual_uncovered = min - ran_min;
  expect_uncovered = (double) ran_max / (double) N2;

  status = (min < ran_min) || (actual_uncovered > 7 * expect_uncovered);

  return status;
}

int
rng_sum_test (gsl_rng * r, double *sigma)
{
  double sum = 0;
  int i, status;

  for (i = 0; i < N2; ++i)
    {
      double x = gsl_rng_uniform (r) - 0.5;
      sum += x;
    }

  sum /= N2;

  /* expect the average to have a variance of 1/(12 n) */

  *sigma = sum * sqrt (12.0 * N2);

  /* more than 3 sigma is an error */
  
  status = (fabs (*sigma) > 3 || fabs(*sigma) < 0.003);

  if (status) {
      fprintf(stderr,"sum=%g, sigma=%g\n",sum,*sigma);
  }

  return status;
}

#define BINS 17
#define EXTRA 10
int
rng_bin_test (gsl_rng * r, double *sigma)
{
  int count[BINS+EXTRA];
  double chisq = 0;
  int i, status;

  for (i = 0; i < BINS+EXTRA; i++)
      count[i] = 0 ;


  for (i = 0; i < N2; i++)
    {
      int j = gsl_rng_uniform_int (r, BINS);
      count[j]++ ;
    }

  chisq = 0 ;
  for (i = 0; i < BINS; i++)
    {
      double x = (double)N2/(double)BINS ;
      double d = (count[i] - x) ;
      chisq += (d*d) / x;
    }

  *sigma = sqrt(chisq/BINS) ;

  /* more than 3 sigma is an error */
  
  status = (fabs (*sigma) > 3 || fabs(*sigma) < 0.003);

  for (i = BINS; i < BINS+EXTRA; i++)
    {
      if (count[i] != 0)
        {
          status = 1 ;
          gsl_test (status, 
                    "%s, wrote outside range in bin test "
                    "(%d observed vs %d expected)",
                    gsl_rng_name(r), i, BINS - 1);
        }
    }

  return status;
}