[转帖]蒙特卡罗Ising算法C的程序

Doublen

[size=13.63636302947998px]偶然见到的的一个小程序，作者是newkitty，有兴趣地可以来看看。以下是原文：

[size=13.63636302947998px]最近学习蒙特卡罗Ising算法，想其看看参考的程序，但是怎么也找不到c的源程序，于是只好硬着头皮自己编了 。。。。。。终于做出一个。。。希望拿出来让大家看看是否正确，这样也可以给其他在学习的人参考一下。。。有不对的地方，麻烦大家多多指点，谢谢。。。

[size=13.63636302947998px]#include<stdio.h>
#include<stdlib.h>
#include<math.h>

[size=13.63636302947998px]#define L 40
double funiseed (double iseed)  //同余法产生随机数
{
double mult=1277,modulo,rmodulo;
modulo=rmodulo=pow(2,17);
iseed=iseed*mult;
if(iseed<modulo) ;
else iseed=iseed-modulo;
iseed=iseed/rmodulo;
return iseed;
}
main()
{
    int i,j,la[L+1][L+1],ip[L+1],im[L+1],count=0;
int ici,ien,mcs,mcsmax,n0,mcstep;
float w[9],m,a;
double jkt,iseed;
    //初始化数组    for(i=0;i<=L;i++) {
  for(j=0;j<=L;j++) {
  la[j] = -1;
  a=1./(L*L);
  }
}
//************************************
for(i=1;i<=L;i++) {
  ip=i+1;
  im=i-1;
}
ip[L]=1;
im[1]=L;
//***********************************
printf("输入样本步数：");
scanf("%d",&mcstep);
printf("输入n0的值：");
scanf("%d",&n0);
printf("输入jkt的值：");
scanf("%lf",&jkt);

//此方法用在2维Ising模型中

[size=13.63636302947998px]    for(j=-4;j<=4;j=j+2) {
  w[j+4]=1;
  if(j>0) w[j+4]=exp(-2*jkt*j);
}

[size=13.63636302947998px] //蒙特卡罗部分

[size=13.63636302947998px]     m=-(L*L);

[size=13.63636302947998px]  for(mcs=1;mcs<=mcstep+n0;mcs++)  {

[size=13.63636302947998px]     for(i=1;i<=L;i++) {
      for(j=1;j<=L;j++) {
                             
             ici=la[j];
    ien=la[ ip ][j]+la[ im ][j]+la[ ip[j] ]+la[ im[j] ];
             ien=ici*ien;

[size=13.63636302947998px]//  printf("put in the seed:");
//         scanf("%lf",&iseed);
             iseed=1+rand()%10;

[size=13.63636302947998px]    if (funiseed ( iseed )*10 < w[ien] ) {   //测试是否反转
     la[j]=-ici;
     m=m-2*ici;

[size=13.63636302947998px]              }
    }
    }

[size=13.63636302947998px]    if (mcs >= n0) {
     count=count+1;
     if (count=n0) {
                      count=0;
                      printf("after %d step m=%f\n",mcs-n0,m*a);   
     }  
    }
  }  

[size=13.63636302947998px]本帖转自52MC