fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. static int ix = 1, iy = 1, iz = 1;  /* 1--30000の任意の整数 */
  6.  
  7. /*
  8.  * 改良版乱数(Wichmann--Hillの乱数発生法)
  9.  */
  10. void init_rnd(int x, int y, int z)  /* x, y, z = 1..30000 */
  11. {
  12.   ix = x;
  13.   iy = y;
  14.   iz = z;
  15. }
  16.  
  17. double Rand(void)  /* 0 <= rnd() < 1 */
  18. {
  19.   double r;
  20.  
  21.   ix = 171 * (ix % 177) -  2 * (ix / 177);
  22.   iy = 172 * (iy % 176) - 35 * (iy / 176);
  23.   iz = 170 * (iz % 178) - 63 * (iz / 178);
  24.   if (ix < 0) ix += 30269;
  25.   if (iy < 0) iy += 30307;
  26.   if (iz < 0) iz += 30323;
  27.   r = ix / 30269.0 + iy / 30307.0 + iz / 30323.0;
  28.   while (r >= 1.0)
  29.     r = r - 1.0;
  30.  
  31.   return r;
  32. }
  33.  
  34. /*
  35.  * ボックス＝ミューラー法
  36.  */
  37. double NormalRandom(double mu, double sigma)
  38. {
  39.   double t, u, r;
  40.   static int sw = 0;
  41.   static double rs;
  42.  
  43.   if (sw == 0) {
  44.     t = sqrt(-2.0 * log(1.0 - Rand()));
  45.     u = 2.0 * M_PI * Rand();
  46.     r = t * cos(u);
  47.     rs = t * sin(u);
  48.     sw = 1;
  49.   } else {
  50.     r = rs;
  51.     sw = 0;
  52.   }
  53.  
  54.   return sigma * r + mu;
  55. }
  56.  
  57. int main(void)
  58. {
  59.   int i, n, ix, *histo;
  60.   double mu, sigma, x, s1 = 0.0, s2 = 0.0;
  61.   double *array;
  62.  
  63.   init_rnd((time(NULL) % 30000) + 1, ((time(NULL) + 23456) % 30000) + 1, ((time(NULL) + 45678) % 30000) + 1);  /* 1..30000 の任意の整数で初期化. */
  64. //  printf("個数 n = ");
  65. //  scanf("%d", &n);
  66.   n = 10000;
  67. //  printf("平均 μ = ");
  68. //  scanf("%lf", &mu);
  69.   mu = 1.5;
  70. //  printf("分散 σ = ");
  71. //  scanf("%lf", &sigma);
  72.   sigma = 0.5;
  73.  
  74.   if ((histo = (int *)malloc(sizeof(int) * n)) == NULL)
  75.     exit(1);
  76.   if ((array = (double *)malloc(sizeof(double) * n)) == NULL)
  77.     exit(1);
  78.  
  79.   /*
  80.    * 乱数はarray[]に保存する
  81.    * 度数分布を同時に求める
  82.    */
  83.   for (i = 0; i < n; i++) {
  84.     x = array[i] = NormalRandom(mu, sigma);
  85.     ix = (int)floor(2.0 * x) + 10;
  86.     if (ix >= 0 && ix < 20)
  87.       histo[ix]++;
  88.     s1 += x;
  89.     s2 += x * x;
  90.   }
  91.  
  92.   /*
  93.    * 平均と標準偏差の計算
  94.    */
  95.   for (i = 0; i < 20; i++)
  96.     printf("%4.1f -- %4.1f: %5.1f%%\n", 0.5 * (i - 10), 0.5 * (i - 9), 100.0 * histo[i] / n);
  97.   s1 /= n;
  98.   s2 = sqrt((s2 - n * s1 * s1) / (n - 1));
  99.   printf("平均 %g  標準偏差 %g\n", s1, s2);
  100.  
  101.   free(array);
  102.   free(histo);
  103.  
  104.   return 0;
  105. }
  106.  
Success #stdin #stdout 0.02s 1852KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
-5.0 -- -4.5:   0.0%
-4.5 -- -4.0:   0.0%
-4.0 -- -3.5:   0.0%
-3.5 -- -3.0:   0.0%
-3.0 -- -2.5:   0.0%
-2.5 -- -2.0:   0.0%
-2.0 -- -1.5:   0.0%
-1.5 -- -1.0:   0.0%
-1.0 -- -0.5:   0.0%
-0.5 --  0.0:   0.2%
 0.0 --  0.5:   2.4%
 0.5 --  1.0:  13.4%
 1.0 --  1.5:  35.0%
 1.5 --  2.0:  33.2%
 2.0 --  2.5:  13.7%
 2.5 --  3.0:   2.0%
 3.0 --  3.5:   0.1%
 3.5 --  4.0:   0.0%
 4.0 --  4.5:   0.0%
 4.5 --  5.0:   0.0%
平均 1.49171  標準偏差 0.501545
https://ideone.com/Ry2nq
language:
C (gcc 8.3)
created:
12 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!