Keeping accuracy when taking decimal to power of integer

Question

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. */

when you have very long loops, it's easier to break down the loop body in functions to make the overall algorithm more apparent. You should at least make the indentation correct, or mark the ends of loop with a comment (like `// end of while(a>m)`). — didierc, Apr 01 '13 at 09:57
@n.m. The full program is given in the link, however I realised it was really long and might annoy people, so I shortened it and took all of the maths out so you could still see the pow(Distance,t). The one above only distributes points on a sphere randomly, many times, however the one in the link distributes them and then minimises the energy function to find optimal configurations. I want t to tend to infinity and therefore give solutions to Tammes's problem, however when t gets large all accuracy is lost. — adrem7, Apr 01 '13 at 10:04
@didierc Sorry for the long loops, I know you should use functions now, but I don't know how to use them and it would take time to learn and change a code that works (I agree it is messy though). I shall edit with comments added in to help out. — adrem7, Apr 01 '13 at 10:06
It is not enough to bring a library into scope, you need to actually use it. As in, insert types and operations defined by the library into your program. — n. m. could be an AI, Apr 01 '13 at 10:24
Never write pow(x, 2). That is a slow (and inaccurate in many cases) way of squaring. — Marc Glisse, Apr 01 '13 at 10:24
@n.m. Yes I realise that I am currently not using it, the issue is I don't know what to use it on... — adrem7, Apr 01 '13 at 10:33
@MarcGlisse Are they not the same? I don't really see how this makes a difference... — adrem7, Apr 01 '13 at 10:33
What I essentially want to do is raise a double to an integer power where the integer power can get arbitrarily large and the output will not lose accuracy by it doing so. I therefore need to know where and how to alter my current code to do this, and I am trying to figure both out at the moment, but struggling. — adrem7, Apr 01 '13 at 10:34
If you want to do infinite-precision calculations with Enegry, make Energy an infinite-precision type. — n. m. could be an AI, Apr 01 '13 at 11:20
But the energy energy depends on the distance and t so don't I need to change them all?? Bearing in mind I have never used functions before I am finding this very difficult, please excuse my ignorance. — adrem7, Apr 01 '13 at 11:27
If. like you are saying, the main precision loss is connected to Energy taken into Nth power, you probably need to change just that. — n. m. could be an AI, Apr 01 '13 at 12:02
@n.m. Energy is a function (mathematically speaking, not programming wise) of Distance, and so I think the loss of accuracy occurs by taking Energy=pow(Distance,t), therefore I guess Distance is losing the accuracy and Energy might just be continuing that loss. — adrem7, Apr 01 '13 at 12:08
You cn't exponentiate mpf_class or mpf_t how am I supposed to raise a decimal to an integer if the decimal isn't allowed to be in the first place?? — adrem7, Apr 01 '13 at 13:48
Did you miss `mpf_pow_ui` in the doc? Note that even if it didn't exist, you can always emulate it with a loop and multiplications. — Marc Glisse, Apr 01 '13 at 16:11
@MarcGlisse I have managed to do it now, it was an issue with how I was defining my variables, you are right, I did need to use mpg_pow_ui and did that in the second piece of code posted above. Do you have any input on why I am getting seg faults when initialising variables by any chance? — adrem7, Apr 01 '13 at 16:19
Most of your variables didn't get an mpf_init. You should use mpf_class, and the get_mpf_t() member function will let you use mpf_pow_ui. — Marc Glisse, Apr 01 '13 at 16:22
I have added the appropriate mpf_init commands, but could you just run by me why I should use mpf_class? Doesn't mpf_pow_ui work fine for mpf_t?? — adrem7, Apr 01 '13 at 16:39

score 0 · Accepted Answer · answered Apr 01 '13 at 20:49

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. */

Keeping accuracy when taking decimal to power of integer

1 Answers1

Linked