Search code examples
cmath.hfloorceil

The output of ceil() and floor() in C language is odd

I am doing my homework to realize a C programe and during my work there is one step that I need to get the interger part of the double type numbers. So I choose ceil() or floor() in to realize this. But the output is unpredictable and far from expected.

The whole program is the following :

 /*
 ************************************************************************
  Includes
 ************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <time.h>
#include <math.h>

/* Solvent Particle Propersities*/
typedef struct 
  {
    double vx,vy,rx,ry ;  /*  velocity AND position               */
    double    cell;                /*  index of grid  */
   }solvent_p;

  /* Cell Propersities*/
 typedef struct 
  {
  double vcm_x,vcm_y ;  /*  center of mass velocity  */
  int    number;        /*  number of solvent in the cell  */
  double roation_angle; /*  roation_angle of the cell */
  }cell_p;

/* periodic boundary condition judging and adjusting fuction PBC */

 void PBC(solvent_p *sol)
 {
  double L = 20.0;  // box size 20

  if(sol->rx >20) sol->rx=sol->rx-L;
  if(sol->rx < 0) sol->rx=sol->rx+L;
  if(sol->ry >20) sol->ry=sol->ry-L;
  if(sol->ry < 0) sol->ry=sol->ry+L;

 }

int main(int argc, char **argv)
{

   // Randome setup generates random numbers from GSL functions 
   const gsl_rng_type * T;
   gsl_rng * r;
   T = gsl_rng_default;
   gsl_rng_default_seed = ((unsigned long)(time(NULL))); //设seed值为当前时间 
   r = gsl_rng_alloc (T);

   solvent_p solvent[4000];
   int i,j,k,index=0;
   cell_p grid[400];
   double alpha=90.0; //roation angle

  /*  Iniinitializing solvent  
   *  All vx=0
   *  half vy = sqrt(3) and half vy=-sqrt(3) to make momentum zero and another requirement is the overall energy is 6000 equals energy of temperature=1 with 4000 solvent 3NkT/2 ,assuming mass of solvent = 1 ,which is all a test quantity 
  *  rx and ry are random number generated by GSL library
  *  cell=20*(ry+rx) is an index of which cell the solvent is in
  */
  for(i=0;i<=3999;i++)
  {
         solvent[i].vx=0.0;

         if(i<=1999)
              solvent[i].vy=sqrt(3);
         else
              solvent[i].vy=-sqrt(3);


      solvent[i].rx =20.0 * gsl_rng_uniform_pos(r);

      solvent[i].ry =20.0 * gsl_rng_uniform_pos(r);

     //printf("%f \t %f \n",solvent[i].rx,solvent[i].ry);

    solvent[i].cell=20*(floor(solvent[i].ry))+floor(solvent[i].rx)+1;
    }

     // grid setting up
     for (i=0;i<=399;i++)
   {
        grid[i].vcm_x=0;
        grid[i].vcm_y=0;
        grid[i].number=0;
        grid[i].roation_angle=0.0;
   }


        /*  Begining Simulation Work
         * 
         *  Fist process is preparing the system to equilibrium for 10000 processes
          * 
          *  the whole process involving two steps steaming and collision and the two steps are conducted one by one in our simulation 
          *  time duration for steaming is 0.1 which is assigned for Molecular Dynamics and time duration for collision is ignored
          * 
          *  
          */


       for(i=0;i<=9999;i++)
 {
       double temp;
       double delta_t_MD=0.1; //time step

        temp=gsl_rng_uniform_pos(r);
        double rand_rx = (temp < 0.5) ? temp : ((-1) * temp ); // randomshift rx;

        temp=gsl_rng_uniform_pos(r);
        double rand_ry = (temp < 0.5) ? temp : ((-1) * temp ); // randomshift ry


       //steaming
      for(j=0;j<=3999;j++)
      {
            //printf("%d \n",j);
            printf("1 :%d \t  %f \t %f \n",j,solvent[j].rx,solvent[j].ry);
            solvent[j].rx=solvent[j].rx+solvent[j].vx * delta_t_MD;
            solvent[j].ry=solvent[j].ry+solvent[j].vy * delta_t_MD;

            printf("2: %d \t  %f \t %f \n",j,solvent[j].rx,solvent[j].ry);

    //randomshift
    solvent[j].rx=solvent[j].rx+rand_rx;
    solvent[j].ry=solvent[j].ry+rand_ry;
    printf("3:  %d \t  %f \t %f \n",j,solvent[j].rx,solvent[j].ry);
    // periodic boundary condition
    PBC(&solvent[j]);
    printf("4: %d \t  %f \t %f \n",j,solvent[j].rx,solvent[j].ry);
    solvent[j].cell=20*(ceil(solvent[j].ry)-1)+ceil(solvent[j].rx);
    printf("5: %d \t  %f \t %f \t %f \t %f \n",j,solvent[j].rx,solvent[j].ry,ceil(solvent[j].rx),ceil(solvent[j].ry));
    index = (int)(solvent[j].cell);
    //printf("%d \t %d \t %f \t %f \t %f \n",j,index,solvent[j].cell,ceil(solvent[j].rx),ceil(solvent[j].ry));

    if ((index>=0) &&(index<=400))
    {
    grid[index].vcm_x=grid[index].vcm_x+solvent[i].vx;
    grid[index].vcm_y=grid[index].vcm_y+solvent[i].vy;
    grid[index].number=grid[index].number+1;
   }
}


// caculating vcm

for (k=0;k<=399;k++)
{
    if (grid[k].number >= 1)
    {
        grid[k].vcm_x=grid[k].vcm_x/grid[k].number;
        grid[k].vcm_y=grid[k].vcm_y/grid[k].number; 
    }

    double temp;
    temp=gsl_rng_uniform_pos(r);
    grid[k].roation_angle = (temp < 0.5) ? alpha : ((-1) * alpha ); 
}



//collsion


 }

gsl_rng_free (r); // free RNG

return 0;

}

Sorry it is some extent long so I did not put in first.And something did not finished but the program framework is set up.

my programe is like this:

    printf("4: %d \t  %f \t %f \n",j,solvent[j].rx,solvent[i].ry);
    solvent[j].cell=20*(floor(solvent[j].ry))+floor(solvent[j].rx)+1;
    printf("5: %d \t  %f \t %f \t %f \t %f \n",j,solvent[j].rx,solvent[i].ry,floor(solvent[j].rx),floor(solvent[j].ry));

And although something as I wanted something more is wrong and below is I choose some parts of the output:

4: 3993       3.851240   17.047031 
5: 3993       3.851240   17.047031   3.000000    9.000000

although floor(solvent[j].rx) is right but floor(solvent[j].ry) is totally wrong.

And the final result of my program is 

 Segmentation fault (core dumped)


  ------------------
  (program exited with code: 139)

How to fix this one? Is there anything I use was wrong?

And to further test the problem I tried ceil() function in and the program and result is like this

program:

    printf("4: %d \t  %f \t %f \n",j,solvent[j].rx,solvent[i].ry);
    solvent[j].cell=20*(ceil(solvent[j].ry)-1)+ceil(solvent[j].rx);
    printf("5: %d \t  %f \t %f \t %f \t %f \n",j,solvent[j].rx,solvent[i].ry,ceil(solvent[j].rx),ceil(solvent[j].ry));

result:

2: 3999       14.571205      2.837654 
4: 3999       14.571205      2.837654 
5: 3999       14.571205      2.837654    15.000000   15.000000

So use the nearest output as an example to illustrate my desire result is:

a= 14.571205, ceil(a)=15.00 b= 2.837654 , ceil(b)=3.00 not the 15.000 in the output

And the problem worsening is that when I just use a and b to test ceil() everything seems perfect:

program:

 #include <stdio.h>
 #include <math.h>

 int main()
{
        double a= 14.571205;  
        double b= 2.837654 ;  

        printf("%f \t %f \n",ceil(a),ceil(b));
        return 0;
 }

output: 15.000000 3.000000

I am using GCC in Linux Ubuntu.

==============================================================================

The problem has been solved.

The real fatal problem is here

    if ((index>=0) &&(index<=400))
    {
    grid[index].vcm_x=grid[index].vcm_x+solvent[i].vx;
    grid[index].vcm_y=grid[index].vcm_y+solvent[i].vy;
    grid[index].number=grid[index].number+1;
   }
}

Both solvent[i].vy and solvent[i].vx should be changed i for j.

Just as @Jon and @Blckknght @Eric Lippert have pointed out.

And I obvious get in the wrong trap and mislead @ouah and @Blckknght and In fact, the output sentences do have problem but they do not caused the prorame to crash only the out of boundary will do that.

Thank you ALL!

And I do like to share Eric Lippert 's comment which is powerful and insightful:

More generally, the problem you have is called "select isn't broken" by experienced programmers. That is, you have done something completely wrong, you cannot figure out what you did wrong, and so you conclude that a piece of basic, easy infrastructure written by experts and tested extensively is wrong. I assure you, floor and ceil do exactly what they are supposed to do. They are not wrong. The problem lies somewhere in your code, not in the standard library. – Eric Lippert

The final result

Solution

  • Your array is declared as:

    solvent_p solvent[4000];

    but you have this loop:

     for(i=0;i<=9999;i++)

    with this function call inside:

     printf("1 :%d \t  %f \t %f \n",j,solvent[j].rx,solvent[i].ry);

    which means you are accessing out-of-bounds array elements.

    EDIT:

    OP test case has been edited to fix the out-of-bound accesses.

    My second suggestion is to check index value (which is used to index grid array) never exceeds the number of elements of grid array after this assignment: index = (int)(solvent[j].cell).