程序无缘无故溢出，请达人帮忙！

greentime

程序有时运行得好好的，可是有时无缘无故出现变量E，M溢出的情况，而且这个溢出的双精度变量E只用到加法运算，屏幕上显示此时其他的变量都很正常，不理解啊，程序是用devcpp编译的，程序用到rand()和sleep()函数
那个我昨晚把long和double改成int和float了，溢出问题依然存在

这是我的代码，写得比较搓，如果有耐心，可以看一下：
/*亚铁磁模型的伊辛模拟，哈密顿量为H=-sum_ij[J_AB*spin_A_i*spin_B_j]-sum_ij[J_A*spin_A_i*spin_A_j
-sum_ij[J_B*spin_B_i*spin_B_j],A和B相间分布，模拟系统的T_M,T_A,居里点*/
#include<stdio.h>
#include<math.h>
#include<time.h>
#include<stdlib.h>
main()
{
      /*定义系统的物理量*/
      int L, N, x, y, t, i, flip, flip_tmp, t_tmp;
      int up, down, left, right;
      float J_AB, J_A, J_B, kT, M, M_A, M_B, dM, A_A, A_B, A;
      float h_bar, uB, g_A, g_B, rnd, sum_ner, sum_far;
      float spin_A_up, spin_A_down, spin_B_up, spin_B_zero, spin_B_down, spin_tmp;
      float E, dE;
      FILE *fp;

      /*初始化物理量，把一些物理量单位化*/
      spin_A_up = 0.5;
      spin_A_down = -0.5;
      spin_B_up = 1.0;
      spin_B_down = -1.0;
      spin_B_zero = 0.0;
      fp=fopen("end.txt","w");
      
      L = 60;
      N = L * L;
      x = y = L;
      J_AB = -0.6;
      J_A = 1.0;
      J_B = 0.1;
      kT = 0.4;
      g_B = 2.0;
      g_A = 2.18;
      h_bar = uB = 1.0;
      flip=0;
      
      double spin[x][y];
      
      M = M_A = M_B = 0.0;
      A = A_A = A_B = 0.0;
      E = 0.0;
      
      printf("L=%d, N=%d\nJ_AB=%f, J_A=%f, J_B=%f, g_A=%f, g_B=%f\n",L,N,J_AB,J_A,J_B,g_A,g_B);
      
      /*初始化自旋的分布，采用随机分布*/
      srand((time(0)-1212300000));
      for(x=0;x<L;x++) {
          for(y=0;y<L;y++) {
              rnd = rand()/(RAND_MAX+0.0);
              if(((x%2)||(y%2))&&((x+y)%2)) {
                  if(rnd <= (1.0/3.0)) spin[x][y] = spin_B_down;
                  else if(rnd >= (2.0/3.0)) spin[x][y] = spin_B_up;
                  else spin[x][y] = spin_B_zero;
                  M_B = M_B - g_B * spin[x][y];
                  A_B = A_B + spin[x][y];
              }
              else {
                   if(rnd < 0.5)   spin[x][y] = spin_A_down;
                   else spin[x][y] = spin_A_up;
                   M_A = M_A - g_A * spin[x][y];
                   A_A = A_A + spin[x][y];
              }
          }
      }
      
      /*计算初始的能量*/
      for(x=0;x<L;x++) {
          if(x==L-1) right = 0;
          else right = x+1;
          if(x==0) left = L-1;
          else left = x-1;
          for(y=0;y<L;y++) {
              if(y==L-1) up = 0;
              else up = y+1;
              if(y==0) down = L-1;
              else down = y-1;
              sum_ner = spin[x][up]+spin[x][down]+spin[left][y]+spin[right][y];
              sum_far = spin[right][up]+spin[right][down]+spin[left][down]+spin[left][up];
              if(abs(spin[x][y])==0.5) E = E - 0.5*J_AB*spin[x][y]*sum_ner - 0.5*J_A*spin[x][y]*sum_far;
              else E = E - 0.5*J_AB*spin[x][y]*sum_ner - 0.5*J_B*spin[x][y]*sum_far;
          }
      }
      printf("Initial Energy E=%f, M_A=%f, M_B=%f, M=%f\n",E,M_A,M_B,M=M_A+M_B);
      
      /*用Metropolis算法计算系统的变化*/
      for(t=0;t<1000;t++) {
          sleep(100);
          srand(time(0)-1212300000);
          fprintf(fp,"%-f     %-f    %-d    %-d\n",E/N,M/N,flip,t);
          flip=0;
          for(i=0;i<N;i++) {
              x = (unsigned int)(1.0*L*rand()/(RAND_MAX+0.0));
              y = (unsigned int)(1.0*L*rand()/(RAND_MAX+0.0));
              if(x==L-1) right=0;
              else right=x+1;
              if(x==0) left=L-1;
              else left=x-1;
              if(y==L-1) up=0;
              else up=y+1;
              if(y==0) down=L-1;
              else down=y-1;
              sum_ner=spin[x][up]+spin[x][down]+spin[left][y]+spin[right][y];
              sum_far=spin[right][up]+spin[right][down]+spin[left][down]+spin[left][up];
              if(abs(spin[x][y])==0.5) {
                  spin_tmp=-spin[x][y];
                  dE=-0.5*J_AB*(spin_tmp-spin[x][y])*sum_ner-0.5*J_A*(spin_tmp-spin[x][y])*sum_far;
                  dM=-g_A*(spin_tmp-spin[x][y]);
              }
              else {
                   rnd=(1.0+rand())/(RAND_MAX+1.0);
                   if(abs(spin[x][y])==0) {
                       if(rnd<0.5) spin_tmp=spin_B_down;
                       else spin_tmp=spin_B_up;
                   }
                   else {
                        if(rnd<0.5) spin_tmp=-spin[x][y];
                        else spin_tmp=spin_B_zero;
                   }
                   dE=-0.5*J_AB*(spin_tmp-spin[x][y])*sum_ner-0.5*J_B*(spin_tmp-spin[x][y])*sum_far;
                   dM=-g_B*(spin_tmp-spin[x][y]);
              }
              if(dE<0) {
                  spin[x][y]=spin_tmp;
                  E = E + dE;
                  M = M + dM;
                  flip++;
              }
              else {
                   rnd=(1.0+rand())/(RAND_MAX+1.0);
                   if(rnd<exp(-dE/kT)) {
                       spin[x][y]=spin_tmp;
                       E = E + dE;
                       M = M + dM;
                       flip++;
                   }
              }
          }
          printf("E=%f, M=%f, flip=%d, t=%d, dE=%f, dM=%f\n",E,M=abs(M),flip,t,dE,dM);
      }
      getch();
}

[ 本帖最后由 greentime 于 2008-6-2 21:13 编辑 ]

r2r4

我那时的伊辛，找了别人的copy，修改了一下，交了作业，可惜全忘了

想得到答案，先把你的代码整理了，注释加全了，
有些你认为不是很必要的注释也加上，没准你的注释和代码就不一样，错误就在这里
不企图大家都懂物理意义这个，但程序运算还是能看明白的

一陀加减乘除，谁看的过来阿~~

r2r4

小建议，
能用double就不要用float，如果你认为这是精度问题

if(x==L-1) right = 0;

我相信这里你不是想和 L 比较,类似的比较不少（影响应该不大)

a = a +b;
可以写成 a += b;

if(((x%2)||(y%2))&&((x+y)%2)) {

可以修改一下

greentime

谢谢各位，我今晚再搞搞，调试得头都大了，要是不行再把详细注释的程序贴出来
我主要是嫌注释影响我调试阅读，呵呵

greentime

问题暂时解决了，我换到linux下用gcc编译，可以稳定运行，当然有小幅度的修改，发现gcc下面rand()的范围竟然是0~十亿多，OMG，windows下面也只有32767,差别太大了

77h2_eleven

原帖由 greentime 于 2008-6-3 14:29 发表
问题暂时解决了，我换到linux下用gcc编译，可以稳定运行，当然有小幅度的修改，发现gcc下面rand()的范围竟然是0~十亿多，OMG，windows下面也只有32767,差别太大了

没对随机数求余么？

sdu_lizhipeng

我也头大了。。。。。

lenovo

只想说，程序不会无缘无故的溢出。
就如这个世界没有无缘无故的爱，
也没有无缘无故的恨。

Sorehead

lenovo版主受到什么刺激了？