c++ - 失去大整数的准确性(pow？)

Question

我正在尝试最小化能量函数以找到解决 Tammes 问题的方法。我的代码使用梯度流方法工作，我可以更改变量 t 的值(能量函数的幂)来估计解决方案。

Tammes 问题是 t 趋于无穷大时的结果，因此显然我希望 t 尽可能大。

问题是我的解决方案(计算的能量)随着 t 的增加而失去准确性。

我认为这与 pow() 函数以及它使用如此大的整数进行计算这一事实有关。

我的代码如下:

#include <stdio.h> 
#include <string.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <time.h>
#include <stdlib.h>
#include <sstream>
using namespace std;
#define PI 3.14159265358979323846

int main()
{

  int a,b,c,d,f,i,j,k,m,n,s,t,Success,Fails;
  double p,theta,phi,Time,Averagetime,Energy,energy,Distance,Length,DotProdForce,
         Forcemagnitude,ForceMagnitude[201],Force[201][4],E[1000001],En[501],Epsilon[4],Ep,
         x[201][4],new_x[201][4],y[201][4],A[201],alpha[201][201],degree,bestalpha[501];

  clock_t t1,t2;
  t1=clock();

t=1;

while(t<1001){

n=2;

while(n<51){

cout << "N=" << n << "\n";

  b=1;
  Time=0.0;

  while(b<101){

    clock_t t3,t4;
    t3=clock();

    if(n>200){
      cout << n << " is too many points for me :-( \n";
      exit(0);
    }

    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;
      }
    }

    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/pow(Distance,t);
      }
    }

    for(i=1;i<=n;i++){
      y[i][1]=x[i][1];
      y[i][2]=x[i][2];
      y[i][3]=x[i][3];
    }

    m=100;

    if (m>100){
      cout << "The m="<< m << " loop is inefficient...lessen m \n";
      exit(0);
    }

    a=1;

    while(a<m){

      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;
        }
      }

      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/pow(Distance,t);
        }
      }

      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;
    }

    /* 
    ostringstream String1;
    String1 << "BestTamrandompoints_" << t << "_" << n;  //add number

    ofstream File1 (String1.str().c_str());
    File1 << "Energy=" << Energy << "\n";
    for(i=1;i<=n;i++){
      File1 << x[i][1] << " " <<   x[i][2] << " " << x[i][3] << "\n";
    }
    File1.close(); 
    */

    E[0]=Energy;

    a=1;
    s=0;
    f=0;
    Forcemagnitude=1.0;

    if(n>80 && t<11)
      p=0.005;
    else if(n>60 && t<11)
      p=0.01;
    else if(n>40 && t<11)
      p=0.02;
    else if(n>25 && t<11)
      p=0.05;
    else if(n>40 && t>11)
      p=0.002;
    else if(n>30 && t>11)
      p=0.005;
    else if(n>20 && t>11)
      p=0.01;
    else if(n>5 && t>11)
      p=0.025;
    else
      p=0.05;

    while(Forcemagnitude>0.000001 && a<1000001){

      for(i=1;i<=n;i++){ 
        Force[i][1]=0.0; 
        Force[i][2]=0.0; 
        Force[i][3]=0.0; 
      } 

      for(i=1;i<=n;i++){ 
        for(j=1;j<i;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));

          Force[i][1]=Force[i][1]+((x[i][1]-x[j][1])/(pow((Distance),3))); 
          Force[i][2]=Force[i][2]+((x[i][2]-x[j][2])/(pow((Distance),3))); 
          Force[i][3]=Force[i][3]+((x[i][3]-x[j][3])/(pow((Distance),3))); 
        } 

        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));

          Force[i][1]=Force[i][1]+((x[i][1]-x[j][1])/(pow((Distance),3))); 
          Force[i][2]=Force[i][2]+((x[i][2]-x[j][2])/(pow((Distance),3))); 
          Force[i][3]=Force[i][3]+((x[i][3]-x[j][3])/(pow((Distance),3)));
        }
      }

      for(i=1;i<=n;i++){

        DotProdForce=Force[i][1]*x[i][1]+Force[i][2]*x[i][2]+Force[i][3]*x[i][3];

        y[i][1]=x[i][1];
        y[i][2]=x[i][2];
        y[i][3]=x[i][3];

        Force[i][1]=Force[i][1]-DotProdForce*y[i][1];
        Force[i][2]=Force[i][2]-DotProdForce*y[i][2];
        Force[i][3]=Force[i][3]-DotProdForce*y[i][3];

        x[i][1] = y[i][1]+p*(Force[i][1]);
        x[i][2] = y[i][2]+p*(Force[i][2]);
        x[i][3] = y[i][3]+p*(Force[i][3]);

        Length=sqrt(pow(x[i][1],2)+pow(x[i][2],2)+pow(x[i][3],2));

        for (j=1;j<=3;j++){
          x[i][j]=x[i][j]/Length;
        }
      }

      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/pow(Distance,t);
        }
      }

      E[a]=energy;

      for(i=1;i<=n;i++){
        ForceMagnitude[i]=pow((pow(Force[i][1],2)+pow(Force[i][2],2)
                               +pow(Force[i][3],2)),0.5);
      }

      cout << "E[" << a << "]=" << E[a] << "\n";

      cout << fixed << setprecision(40) << "Energy=" << Energy << " energy=" << energy << "\n";

      if (energy<Energy)
        Energy=energy,s=s+1;
      else
        energy=Energy,f=f+1,p=(9.5*p)/10;

      for(i=1;i<=n-1;i++){
        if(ForceMagnitude[i]<ForceMagnitude[i+1])
          ForceMagnitude[i]=ForceMagnitude[i+1];
        else
          ForceMagnitude[i+1]=ForceMagnitude[i];
      }

      Ep=sqrt(pow(E[a]-E[a-1],2));

      if (Ep<0.00000000000000000000000000000000000000000000000000000000001){
        ForceMagnitude[n]=0.000001;
        cout << "BROKEN \n";
      }
      else
        ForceMagnitude[n]=ForceMagnitude[n];

      Forcemagnitude=ForceMagnitude[n];

      cout << fixed << setprecision(60) << "FM=" << Forcemagnitude << "\n";
      cout << "Energy=" << Energy << "\n";

      a=a+1;
    }

    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));

        degree=(180/PI);

        alpha[i][j]=degree*acos((2.0-pow(Distance,2))/2.0);
      }
    }

    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        cout << "alpha[" << i << "][" << j << "]=" << alpha[i][j] << "\n";
      }
    }

    for(i=1;i<=n-1;i++){
      for(j=i+1;j<=n-1;j++){
        if(alpha[i][j]>alpha[i][j+1])
          alpha[i][j]=alpha[i][j+1];
        else
          alpha[i][j+1]=alpha[i][j];
      }
    }

    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        cout << "alpha[" << i << "][" << j << "]=" << alpha[i][j] << "\n";
      }
    }

    for(i=1;i<=n-2;i++){
      if(alpha[i][n]>alpha[i+1][n])
        alpha[i][n]=alpha[i+1][n];
      else
        alpha[i+1][n]=alpha[i][n];
    }

    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        cout << "alpha[" << i << "][" << j << "]=" << alpha[i][j] << "\n";
      }
    }

    bestalpha[b]=alpha[n-1][n];

    En[b]=Energy;

    b=b+1;

    t4=clock();
    float diff ((float)t4-(float)t3);
    float seconds = diff / CLOCKS_PER_SEC;

    Time = Time + seconds;

  }

  Averagetime = Time/(b-1);

  cout << fixed << setprecision (4) << "Average Time: " << Averagetime << "(s) \n";
  cout << fixed << setprecision(10) << "Energy=" << Energy << "\n";
  cout << "s=" << s << " f=" << f << "\n";

  Success=Success+s;
  Fails=Fails+f;

  cout << "Successes=" << Success << " Failures=" << Fails << "\n";
  cout << setprecision(5) << "Success Rate=" << ((Success*1.0)/((Success*1.0)+(Fails*1.0)))*100.0 << "\n";

  t2=clock();
    float diff ((float)t2-(float)t1);
    float seconds = diff / CLOCKS_PER_SEC;

  for(i=1;i<=n;i++){
    for(j=1;j<=3;j++){
      y[i][j]=x[i][j];
    }
  }

  a=1;
  d=0;
  c=0;

  while(a<2){

  theta=x[a][2]/x[a][3];
  theta=d*PI+atan(theta);

  for(i=1;i<=n;i++){
    new_x[i][1]=x[i][1];
    new_x[i][2]=x[i][2]*cos(theta)-x[i][3]*sin(theta);
    new_x[i][3]=x[i][2]*sin(theta)+x[i][3]*cos(theta);
  }

  phi=new_x[a][1]/new_x[a][3];
  phi=c*PI+atan(phi);

  for(i=1;i<=n;i++){
    x[i][1]=new_x[i][1]*cos(phi)-new_x[i][3]*sin(phi);
    x[i][2]=new_x[i][2];
    x[i][3]=new_x[i][1]*sin(phi)+new_x[i][3]*cos(phi);
  }

  Epsilon[1]=sqrt(pow((x[a][1]),2));  
  Epsilon[2]=sqrt(pow((x[a][2]),2));
  Epsilon[3]=sqrt(pow((x[a][3]-1.0),2));  

  cout << "x[" << a << "][1]=" << x[a][1] << " x[" << a << "][2]=" << x[a][2] << " x[" << a << "][3]=" << x[a][3] << "\n";  

   if(Epsilon[1]<0.00001 && Epsilon[2]<0.00001 && Epsilon[3]<0.00001)
    a=a+1;
  else if(d==0 && c==0)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        c=1;
      }
    }
  else if(d==0 && c==1)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        d=1;
        c=0;
      }
    }
  else if(d==1 && c==0)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        c=1;
      }
    }
  else if(d==1 && c==1)
    break;

  }

  cout << "a=" << a << " d=" << d << " c=" << c << "\n";

  ostringstream String2;
  String2 << "TamPoints_" << t << "_" << n;  //add number

  ofstream File2 (String2.str().c_str());
  for(i=1;i<=n;i++){
    File2 << x[i][1] << " " <<   x[i][2] << " " << x[i][3] << "\n";
  }
  File2.close(); 

  ostringstream String3;
  String3 << "TamEnergies_" << t << "_" << n;  //add number

  ofstream File3 (String3.str().c_str());
  for(i=1;i<a;i++){
    File3 << fixed << setprecision (10) << E[i] << "\n";
  }
  File3.close(); 

  ofstream File4 ("mypoints");
  for(i=1;i<=n;i++){
    File4 << x[i][1] << " " <<   x[i][2] << " " << x[i][3] << "\n";
  }
  File4.close(); 

  ostringstream String4;
  String4 << "TamInfo_" << t << "_" << n;  //add number

  ofstream File5 (String4.str().c_str());
  File5 << "Iterations=" << a-1 << "\n";
  File5 << "Successes=" << s << " Failures=" << f << "\n";
  File5 << fixed << setprecision(20) << "Energy=" << Energy << "\n";
  File5 << fixed << setprecision(5) << "Total run time: " << seconds << "(s) \n";
  File5 << fixed << setprecision(5) << "Average run time: " << Averagetime << "(s) \n";
  File5.close(); 

  ostringstream String5;
  String5 << "Tam%Energies_" << t << "_" << n;  //add number

  ofstream File6 (String5.str().c_str());
  for(i=1;i<b;i++){
    File6 << fixed << setprecision(20) << En[i] << "\n";
  }
  File6.close();

  ostringstream String6;
  String6 << "TamAngles_" << t << "_" << n << ".2";  //add number

  ofstream File7 (String6.str().c_str());
  for(i=1;i<b;i++){
      File7 << fixed << setprecision(15) << bestalpha[i] << "\n";
  }
  File7.close();

  cout << fixed << setprecision(5) << "Run time: " << seconds << "(s)" << "\n";

  n=n+1;
  }

    if(t==1)
    t=2;
    else if(t==2)
    t=5;
    else if(t==5)
    t=10;
    else if(t==10)
    t=25;
    else if(t==25)
    t=50;
    else if(t==50)
    t=100;  
    else if(t==100)
    t=250;
    else if(t==250)
    t=500;
    else if(t==500)
    t=1000;
    else
    t=t+1;
}

  return 0;

}

注意:我很欣赏我的代码非常粗糙，但是我的论文很快就会完成，因此我不能花时间让它变得漂亮，我只需要它可以工作。对于缺少功能，我提前表示歉意，我对任何类型的编码都是新手，直到我觉得花所有时间整理会适得其反，才意识到要使用这些功能。我确定我的问题是能量计算:1/pow(Distance,t)，因此随着 t 的增加，它们会失去准确性。该代码适用于较小的 t(至少高达 100)，但对于 t=250，代码的准确性开始下降。

有什么办法可以快速解决这个问题吗？我听说过 bignum 图书馆，但我对此一无所知，今晚将继续阅读。提前感谢大家，A。

Answer 1

您走在正确的轨道上——您需要任意精度的整数。我用过 the GMP library在过去，但那是很久以前的事了……祝你好运!