蒙特卡罗(可能是模拟退火?)单位球面上N个相互排斥点的方法c++

我需要在c++中创建一个算法来使用蒙特卡罗方法模拟球体上相互排斥的点。到目前为止，我得到的是:

#include <stdio.h> 
#include <string.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <time.h>
#include <stdlib.h>
using namespace std;
int main()
{
  int a,f,g,n,m,i,j,k,r,s;
  double p,q,Energy,energy,y[101][4],x[101][4],Length,Distance;
 clock_t t1,t2;
  t1=clock();
  /*  set the number of points */
  n=12;
  /* check that there are no more than 100 points */
  if(n>100){
    cout << n << " is too many points for me :-( n";
    exit(0);
  }
  /* reset the random number generator */
  srand((unsigned)time(0));  
  for (i=1;i<=n;i++){
    x[i][1]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
    x[i][2]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
    x[i][3]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
    Length=sqrt(pow(x[i][1],2)+pow(x[i][2],2)+pow(x[i][3],2));
    for (k=1;k<=3;k++){
      x[i][k]=x[i][k]/Length;
    }
  }
  /* calculate the energy */
  Energy=0.0;
  for(i=1;i<=n;i++){
    for(j=i+1;j<=n;j++){
      Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                    +pow(x[i][3]-x[j][3],2));
      Energy=Energy+1.0/Distance;
    }
  }
  /* Save Original Points */
  for(i=1;i<=n;i++){
    y[i][1]=x[i][1];
    y[i][2]=x[i][2];
    y[i][3]=x[i][3];
  }
  /* Loop for random points m times*/
  m=10;
  if (m>100){
    cout << "The m="<< m << " loop is inefficient...lessen m n";
    exit(0);
  }
  a=1;
  while(a<m){
    /* assign random points */
    for (i=1;i<=n;i++){
      x[i][1]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[i][2]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[i][3]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      Length=sqrt(pow(x[i][1],2)+pow(x[i][2],2)+pow(x[i][3],2));
      for (k=1;k<=3;k++){
        x[i][k]=x[i][k]/Length;
      }
    }
    /* calculate the energy */
    energy=0.0;
    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                      +pow(x[i][3]-x[j][3],2));
        energy=energy+1.0/Distance;
      }
    }
    if(energy<Energy)
      for(i=1;i<=n;i++){
        for(j=1;j<=3;j++){
          Energy=energy;
          y[i][j]=x[i][j];
        }
      }
    else
      for(i=1;i<=n;i++){
        for(j=1;j<=3;j++){
          energy=Energy;
          x[i][j]=y[i][j];
        }
      }
    a=a+1;
  }
  /* Output the best random energy */
  cout << "Energy=" << Energy << "n";
  m=10;
  a=1;
  while(a<m){
    /* Choose random point to move */
    g=(rand() % n)+1;
    /* Choose a p small to give q in a range -p <= q <= p */
    p=0.1;
    /* q is how much I am moving the random point by */
    q=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0*p;
    /* Move the point by q */
    for(j=1;j<=3;j++){
      x[g][j]=((x[g][j])+q);
    }
    /* Bring it back onto sphere */
    Length=sqrt(pow(x[g][1],2)+pow(x[g][2],2)+pow(x[g][3],2));
    for (k=1;k<=3;k++){
      x[g][k]=x[g][k]/Length;
    }
    /* Calculate the new energy */
    energy=0.0;
    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                         +pow(x[i][3]-x[j][3],2));
        energy=energy+1.0/Distance;
      }
    }
    /* Choose best energy and therefore best point */
    if (energy<Energy)
      Energy=energy,x[g][1]=y[g][1],x[g][2]=y[g][2],x[g][3]=y[g][3];
    else
      energy=Energy,y[g][1]=x[g][1],y[g][2]=x[g][2],y[g][3]=x[g][3];
    a=a+1;  
  }
   /* Output the best single shift energy */
  cout << "Energy=" << Energy << "n";
  /* Set fail count to 0 */
  s=0;
  f=0;
  r=1;
  **p=0.1;**
  /* Maximum distance to move the random point */
  while (**p>0.00001**) {
    /* Number of loops to do */
    while (**r<3000**) {
      g=(rand() % n)+1;
      /* q is how much I am moving the random point by -p<=q<=p*/
      q=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0*p;
      /* Move the point by q */
      for(j=1;j<=3;j++){
        x[g][j]=((x[g][j])+q);
      }
      /* Bring it back onto sphere */
      Length=sqrt(pow(x[g][1],2)+pow(x[g][2],2)+pow(x[g][3],2));
      for (k=1;k<=3;k++){
        x[g][k]=x[g][k]/Length;
      }
      /* Calculate the new energy */
      energy=0.0;
      for(i=1;i<=n;i++){
        for(j=i+1;j<=n;j++){
          Distance=sqrt(pow(y[i][1]-y[j][1],2)+pow(y[i][2]-y[j][2],2)
                        +pow(y[i][3]-y[j][3],2));
          energy=energy+1.0/Distance;
        }
      }
      /* Choose best energy and therefore best point */
      if (energy<Energy)
        Energy=energy,x[g][1]=y[g][1],x[g][2]=y[g][2],x[g][3]=y[g][3],s=s+1;
      else
        energy=Energy,y[g][1]=x[g][1],y[g][2]=x[g][2],y[g][3]=x[g][3],f=f+1;
      r=r+1;
    }
    **/* Calculate percentage fails */
    if ((100.0*(f/r))>50.0)
      p=(p-0.00001);
    else
      p=p;**
    r=0;  
  }
  cout << "Overall Success Rate = " << ((s*1.0)/((s+f)*1.0))*100 << "%" << "n";
  cout << "Energy=" << fixed << setprecision(10) << Energy << "n";

  ofstream Bestpointssofar ("Bestpointssofar");
  for(i=1;i<=n;i++){
    Bestpointssofar << y[i][1] << " " <<   y[i][2] << " " << y[i][3] << "n";
  }
  Bestpointssofar.close(); 
  t2=clock();
    float diff ((float)t2-(float)t1);
    float seconds = diff / CLOCKS_PER_SEC;
    cout << fixed << setprecision(2) << "Run time: " << seconds << "(s)" << "n";
    return 0;
}

我认为这是可以的(注意我基本上是在尽量减少能量函数)，但我想使它更准确/使它运行得更快。要做到这一点，我认为我应该改变p的值，while循环条件或如何在代码末尾改变p。(所有这些都在*…*我试图鼓励他们让你明白我的意思。大约是代码的四分之三)。我已经坐了几个小时试图改变这些情况，但没有任何效果。对于n=12(球面上的12个点)，我的能量应该是49.16525306，但实际上我只能得到50.5到54.0之间的能量。我知道这是相对较好的，但我希望它更准确(即使它确实需要一些时间)。如果可能的话，我也希望成功率能提高(我的总体成功率绝对令人震惊)。

谢谢,a .

(注意:如果你想运行代码，你必须去掉两个*。