用C语言在直方图中绘制随机数答案

【问题标题】：Plotting random numbers in a historgram in C language用C语言在直方图中绘制随机数
【发布时间】：2020-09-05 20:57:29
【问题描述】：

我使用“C 书中的数字食谱”中的代码来生成一个统一的随机数，该随机数被命名为函数float ran1(long *idum)。它在0 and 1 之间产生一个统一的随机数。我正在尝试对这些统一的随机数进行分类并制作直方图。我编写了以下代码，但是当我绘制 bin 与频率的关系图时，我没有得到均匀分布。

谁能帮我知道问题出在哪里？

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

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)


float ran1(long *idum);

int main(){
    
    double delta, x;
    int i , ibin ;
    long idum ;
    int histog[500] ; 
    
    
    delta=0.02 ;
    
    idum=-1 ;// necessary to generate the uniform random number, must be negative integer
    
    FILE *file1;
    
    file1=fopen("guassians.dat","w") ;
    
    
    for (i=0 ; i<=10000 ; ++i) {
        
        x=ran1(&idum) ; // here ran1(&idum) returns a uniform random number
      
        
        ibin=round(x/delta);
        
        if ( abs(ibin) < 250 ) {
            
            histog[ibin]=histog[ibin] + 1 ;
        }
    }
    
    for (ibin=-250 ; ibin <= 250 ; ++ibin){
        
        fprintf(file1, "%0.10lf \t %d \n", ibin*delta , histog[ibin]);
        
    }
    

        
    fclose(file1);
        
}



float ran1(long *idum){
    
    int j ;
    long k;
    static long iy=0;
    static long iv[NTAB];
    float temp ;
    
    if (*idum <= 0 || !iy) {
        
        if (-(*idum) < 1) *idum=1 ;
        else *idum = -(*idum) ;
        for (j=NTAB+7 ; j>=0 ; j--) {
            k=(*idum)/IQ ;
            *idum=IA*(*idum-k*IQ)-IR*k ;
            if (*idum < 0) *idum += IM;
            if (j < NTAB) iv[j] = *idum;
        }
        iy=iv[0] ;
    }
    k=(*idum)/IQ ;
    *idum=IA*(*idum-k*IQ)-IR*k ;
    if (*idum < 0) *idum += IM ;
    j=iy/NDIV ;
    iy=iv[j] ;
    iv[j] = *idum ;
    if ((temp=AM*iy) > RNMX) return RNMX ;
    else return temp;
        
}

【问题讨论】：

ran1 在哪里声明和定义？为什么它需要 long 参数作为指针传递？
Mohammed Alhissi，答案到达后更改问题的重要部分是糟糕的礼仪。帖子回滚。而是What should I do when someone answers my question?
@Bob__我已经修改了代码。

标签： c random histogram

【解决方案1】：

至少这些问题

int histog[500] 未初始化。
用负值索引histog[ibin] 是未定义的行为。

【讨论】：

感谢您的建议。我希望绘图围绕 0 对称。因此，我认为负索引是有意义的。
@MohammedAlhissi 使用int histog[500]={0}，代码可以使用索引[0...499]。
那么如何让直方图 bin 覆盖负值呢？

【解决方案2】：

我发现我应该使用指针而不是数组。现在，下面的代码生成[0:1]（不包括）之间的随机数，这些随机数是均匀分布的，

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

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)


float ran1(long *idum){

int j ;
long k;
static long iy=0;
static long iv[NTAB];
float temp ;

if (*idum <= 0 || !iy) {
    
    if (-(*idum) < 1) *idum=1 ;
    else *idum = -(*idum) ;
    for (j=NTAB+7 ; j>=0 ; j--) {
        k=(*idum)/IQ ;
        *idum=IA*(*idum-k*IQ)-IR*k ;
        if (*idum < 0) *idum += IM;
        if (j < NTAB) iv[j] = *idum;
    }
    iy=iv[0] ;
}
k=(*idum)/IQ ;
*idum=IA*(*idum-k*IQ)-IR*k ;
if (*idum < 0) *idum += IM ;
j=iy/NDIV ;
iy=iv[j] ;
iv[j] = *idum ;
if ((temp=AM*iy) > RNMX) return RNMX ;
else return temp;
    
}

int main(){

double delta ,x ;
int i , ibin, maxbin;
long idum=-1 ;
int *p ;


maxbin=100 ;

delta=1.0/maxbin ;

p= calloc(maxbin,sizeof(int)) ;

FILE *myfile1 ;

myfile1 = fopen("gs.dat","w") ;

for (i=0; i<500000 ; ++i) {

    x=ran1(&idum) ;

    ibin=round(x/delta);

    if ( ibin >= 0 && ibin< 100){

        *(p+ibin)=*(p+ibin) + 1 ;
     }

  }

for (ibin=0 ; ibin < 100 ; ++ibin){

       fprintf(myfile1, "%lf \t %lf \n", ibin*delta , *(p+ibin)/(500000*delta));
 }


free(p);

fclose(myfile1) ;
}

【讨论】：