#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#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[i][j] = -1;
  a=1./(L*L);
  }
 }
 //************************************
 for(i=1;i<=L;i++) {
  ip[i]=i+1;
  im[i]=i-1;
 }
 ip[L]=1;
 im[1]=L;
 //***********************************
 printf("輸入樣本步數：");
 scanf("%d",&mcstep);
 printf("輸入n0的值：");
 scanf("%d",&n0);
 printf("輸入jkt的值：");
 scanf("%lf",&jkt);
 
 //此方法用在2維Ising模型中
    for(j=-4;j<=4;j=j+2) {
  w[j+4]=1;
  if(j>0) w[j+4]=exp(-2*jkt*j);
 }
 //蒙特卡羅部分
     m=-(L*L);
  for(mcs=1;mcs<=mcstep+n0;mcs++)  {
     for(i=1;i<=L;i++) {
      for(j=1;j<=L;j++) {
                             
             ici=la[i][j];
    ien=la[ ip[i] ][j]+la[ im[i] ][j]+la[i][ ip[j] ]+la[i][ im[j] ];
             ien=ici*ien;
//  printf("put in the seed:");
//         scanf("%lf",&iseed);
             iseed=1+rand()%10;
    if (funiseed ( iseed )*10 < w[ien] ) {   //測試是否反轉
     la[i][j]=-ici;
     m=m-2*ici;
              }
    }
    }
    if (mcs >= n0) {
     count=count+1;
     if (count=n0) {
                      count=0;
                      printf("after %d step m=%f\n",mcs-n0,m*a);    
     }  
    } 
  }  
}