2013-03-31 48 views
1

我想尽量减少能量函数来找到泰姆斯问题的解决方案。我的代码通过使用梯度流方法工作,我可以改变我的变量t(能量函数的功率)的值来估计解。失去大整数的精确度(pow?)

泰姆斯的问题是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(距离,t),并且随着t的增加它们将失去准确性。该代码适用于小t(至少高达100),但是对于t = 250,代码在准确性上开始失败。

有没有什么办法可以快速解决这个问题?我听说过这个高档图书馆,但我对此一无所知,今晚将在此上面阅读。感谢所有的提前,A.

+0

我知道你不是在问这个问题。但是,将事物声明为主要的两倍是非常接近于吃掉所有堆栈空间。 – carlosdc

+0

@carlosdc:据我所知,这是一种FORTRAN风格(给定1索引)。 – user7116

回答

1

你在正确的轨道 - 你需要任意的精度整数。过去我使用the GMP library,但那是很久以前的...祝你好运!

+0

因此,昨晚我得到了GMP库,并设法让它工作,现在我需要知道如何使用它并重写我的代码。有没有人知道我应该改变什么?只是吨,或能源,能源和距离?还是全部? – adrem7

+0

注意:对于那些稍后我使用'$ g ++ -o program“filename.cc”-lgmpxx -lgmp -I/folder/where/header/is -L/folder/where/library/is“进行编译的人以及将#include 放入我的范围。他们通常都会在您最初解压缩的文件夹中。如果你找不到它,这个库实际上被称为libgmpxx.dylib和gmpxx.4.dylib,并且头文件被称为gmpxx.h – adrem7

+0

我真的很困惑如何使用GMP,我建立了另一个问题,因为我认为它应该是它自己的一个,但我想知道你对我在那里做了什么有什么建议或评论?谢谢,答:http://stackoverflow.com/questions/15741101/using-gmp-library-to-alter-my-code-what-to-change – adrem7

0

如果你想使任意大,你需要一个bignum库。

对于可能工作取决于你如何大需要吨至是,尝试宣告你的号码为unsigned long long(如果他们将永远是正)速战速决或long long代替int(这里了解更多:http://msdn.microsoft.com/en-us/library/s3f49ktz(v=vs.80).aspx

+0

你的意思只是吨,或距离以及所有可能导致精度损失? – adrem7

+0

首先,我会尝试将它当作你当前声明为int的所有变量。然后,如果有效,您可以根据需要花时间弄清究竟哪些需要更改,哪些可以保留为整数。如果它不起作用,那么你不会浪费太多时间尝试。 – Cam

+0

虽然这没有奏效,但还是提示了一些提示,当我有更多时间时可以看一看。 – adrem7