Search code examples
c++floating-accuracygmparbitrary-precision

Keeping accuracy when taking decimal to power of integer


My code is as follows (I have simplified it for ease of reading, sorry for the lack of functions):

#include <stdio.h> 
#include <string.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <time.h>
#include <stdlib.h>
#include <sstream>
#include <gmpxx.h>
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;

/* Set parameter t, the power in the energy function */

while(t<1001){

n=2;

/*set parameter n, the number of points going onto the sphere */

while(n<51){

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

  b=0;
  Time=0.0;

  /* Set parameter b, just a loop to distribute points many times (100) */

  while(b<100){

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

    /* Points have now been distributed randomly and normalised so they sit on 
       unit sphere */

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

    /*Energy has now been calculated for the system of points as a summation 
      function this is where accuracy is lost */

    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;

    /* Distributing points m-1 times and choosing the best random distribution */

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

    /* End of random distribution loop, the loop for a<m */

    En[b]=Energy;

    b=b+1;

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

    Time = Time + seconds;

  } 

  /* End of looping the entire body of the program, used to get an average reading */

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

  n=n+1;
  }

  /* End of n loop, here n increases so I get outputs for n from 2 to 50 for each t */

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

/* End of t loop, t changes to previously decided values to estimate Tammes's problem
   would like t to be as large as possible but t>200 makes energy calculations lose 
   accuracy */

  return 0;

} /* End of main function and therefore program. In original as seen by following link 
     below the code will use gradient flow algorithm before end of b, n and t loops to 
     minimise the energy function and therefore get accurate solutions. */

Every time I run the code for t>200 the energy output loses accuracy (as it is raised to a high power), I have been told I need to use arbitrary precision integers and to get the GMP library. I have done this and have managed to get the code run with the GMP library in my scope, but I don't really get what I am supposed to alter.

Do I alter t or energy (and Energy) or Distance or all three (/four)?? I don't really understand what I am supposed to change, but I am reading up now how to do it from the manual.

Note:My original question was here, but I thought that had really been answered and this warranted a new one. I shall accept the answer there when this actually works: Losing accuracy for large integers (pow?)

I have altered my code (shown below), but I am just coming up with Segmentation fault 11 as soon as I initialise the En[b]. I would really appreciate it if the comments were a little more in-depth as to what I am to do. Thanks for all the help so far, A.

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

int main()
{

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

  mpf_t Energy,energy,Power,D,En[501];

  mpf_set_default_prec(1024);

  mpf_init(Power);
  mpf_init(D);

  clock_t t1,t2;
  t1=clock();

  t=1000;

/* Set parameter t, the power in the energy function */

while(t<1001){

n=2;

/*set parameter n, the number of points going onto the sphere */

while(n<51){

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

  b=0;
  Time=0.0;

  /* Set parameter b, just a loop to distribute points many times (100) */

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

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

    /* Points distributed randomly and normalised so they sit on unit sphere */

    mpf_init (Energy);

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

        mpf_set_d(D,Distance);

        mpf_pow_ui(Power,D,t);      
        mpf_ui_div(Power,1.0,Power);
        mpf_add(Energy,Energy,Power);

      }
    }

    cout << "Energy=" << Energy << "\n";

    /*Energy calculated as a summation function this is where accuracy is lost */

    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;

    /* Distributing points m-1 times and choosing the best random distribution */

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

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

      mpf_init(energy);

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

        mpf_set_d(D,Distance);

        mpf_pow_ui(Power,D,t);      
        mpf_ui_div(Power,1.0,Power);
        mpf_add(energy,energy,Power);
        }
      }

      cout << "energy=" << energy << "\n";

      if(energy<Energy)
        for(i=1;i<=n;i++){
          for(j=1;j<=3;j++){
            mpf_set(Energy,energy);
            y[i][j]=x[i][j];
          }
        }
      else
        for(i=1;i<=n;i++){
          for(j=1;j<=3;j++){
            mpf_set(energy,Energy);
            x[i][j]=y[i][j];
          }
        }

      a=a+1;
    }

    /* End of random distribution loop, the loop for a<m */

    cout << "Energy=" << Energy << "\n";

    mpf_init(En[b]);
    mpf_set(En[b],Energy);

    for(i=0;i<=b;i++){
      cout << "En[" << i << "]=" << En[i] << "\n";
    }

    b=b+1;

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

    Time = Time + seconds;

  } 

  /* End of looping the entire body of the program, used to get an average reading */

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

  n=n+1;
  }

  /* End of n loop, here n increases so I get outputs for n from 2 to 50 for each t */

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

/* End of t loop, t changes to previously decided values to estimate Tammes's problem
   would like t to be as large as possible but t>200 makes energy calculations lose 
   accuracy */

  return 0;

} /* End of main function and therefore program. In original as seen by following link 
     below the code will use gradient flow algorithm before end of b, n and t loops to 
     minimise the energy function and therefore get accurate solutions. */

Solution

  • The code now looks like this, for those in future apparently you must really learn how to use the GMP Library which can be found here http://gmplib.org, most issue I have had were solved by all those helpful people in the comments, so check them out if you are having issues. Thanks.

    #include <stdio.h> 
    #include <string.h>
    #include <math.h>
    #include <iostream>
    #include <iomanip>
    #include <fstream>
    #include <time.h>
    #include <stdlib.h>
    #include <sstream>
    #include <gmpxx.h>
    using namespace std;
    #define PI 3.14159265358979323846
    
    int main()
    {
    
      int a,b,c,d,f,i,j,k,m,n,s,Success,Fails;
      double p,theta,phi,Time,Averagetime,Distance,Length,DotProdForce,
             Forcemagnitude,ForceMagnitude[201],Force[201][4],E[1000001],Epsilon[4],Ep,
             x[201][4],new_x[201][4],y[201][4],A[201],alpha[201][201],degree,bestalpha[501];
      unsigned long int t;
    
      mpf_t Energy,energy,Power,D,En[501];
    
      mpf_set_default_prec(1024);
    
      mpf_init(Power);
      mpf_init(D);
    
      clock_t t1,t2;
      t1=clock();
    
      t=1000;
    
    /* Set parameter t, the power in the energy function */
    
    while(t<1001){
    
    n=2;
    
    /*set parameter n, the number of points going onto the sphere */
    
    while(n<3){
    
    cout << "N=" << n << "\n";
    
      b=0;
      Time=0.0;
    
      /* Set parameter b, just a loop to distribute points many times (100) */
    
      while(b<2){
    
        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;
          }
        }
    
        for(i=1;i<=n;i++){
          for(j=1;j<=3;j++){
            cout << "x[" << i << "][" << j << "]=" << x[i][j] << "\n";
          }
        }
    
        /* Points distributed randomly and normalised so they sit on unit sphere */
    
        mpf_init (Energy);
    
        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));
    
            mpf_set_d(D,Distance);
    
            mpf_pow_ui(Power,D,t);      
            mpf_ui_div(Power,1.0,Power);
            mpf_add(Energy,Energy,Power);
    
          }
        }
    
        cout << "Energy=" << Energy << "\n";
    
        /*Energy calculated as a summation function this is where accuracy is lost */
    
        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;
    
        /* Distributing points m-1 times and choosing the best random distribution */
    
        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;
            }
          }
    
          for(i=1;i<=n;i++){
            for(j=1;j<=3;j++){
              cout << "x[" << i << "][" << j << "]=" << x[i][j] << "\n";
            }
          }
    
          mpf_init(energy);
    
          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));
    
            mpf_set_d(D,Distance);
    
            mpf_pow_ui(Power,D,t);      
            mpf_ui_div(Power,1.0,Power);
            mpf_add(energy,energy,Power);
            }
          }
    
          cout << "energy=" << energy << "\n";
    
          if(energy<Energy)
            for(i=1;i<=n;i++){
              for(j=1;j<=3;j++){
                mpf_set(Energy,energy);
                y[i][j]=x[i][j];
              }
            }
          else
            for(i=1;i<=n;i++){
              for(j=1;j<=3;j++){
                mpf_set(energy,Energy);
                x[i][j]=y[i][j];
              }
            }
    
          a=a+1;
        }
    
        /* End of random distribution loop, the loop for a<m */
    
        cout << "Energy=" << Energy << "\n";
    
        mpf_init(En[b]);
    
        mpf_set(En[b],Energy);
    
        for(i=0;i<=b;i++){
          cout << "En[" << i << "]=" << En[i] << "\n";
        }
    
            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];
    
        for(i=1;i<=b;i++){
          cout << "Best Angle[" << i << "]: " << bestalpha[b] << "\n"; 
        }
    
        b=b+1;
    
        t4=clock();
        float diff ((float)t4-(float)t3);
        float seconds = diff / CLOCKS_PER_SEC;
    
        Time = Time + seconds;
    
      } 
    
      /* End of looping the entire body of the program, used to get an average reading */
    
      t2=clock();
        float diff ((float)t2-(float)t1);
        float seconds = diff / CLOCKS_PER_SEC;
    
      n=n+1;
      }
    
      /* End of n loop, here n increases so I get outputs for n from 2 to 50 for each t */
    
        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=1001;
    }
    
    /* End of t loop, t changes to previously decided values to estimate Tammes's problem
       would like t to be as large as possible but t>200 makes energy calculations lose 
       accuracy */
    
      return 0;
    
    } /* End of main function and therefore program. In original as seen by following link 
         below the code will use gradient flow algorithm before end of b, n and t loops to 
         minimise the energy function and therefore get accurate solutions. */