Math is simply 0

This is a discussion on Math is simply 0 within the C Programming forums, part of the General Programming Boards category; OK, I admit right now that I am wrong. That said, I can find absolutly no reason when this should ...

  1. #1
    Registered User crepincdotcom's Avatar
    Join Date
    Oct 2003
    Posts
    94

    Math is simply 0

    OK, I admit right now that I am wrong.

    That said, I can find absolutly no reason when this should be. The function excerpted below should compute the Force F force each set of 2 particels in a system. Yet, no matter how many ways I've tried it, the line
    Code:
    F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));
    always returns 0! I tried it with "G" instead of the 00....0667, with G as both a #declare and a local variable. I tried it with and without the type casts. I simple do not understand.

    Here is the function:

    Code:
    void ComputeForces( void )
    {
      int i,j,fxsign,fysign; //signs: 0 +, 1 -
      double r,F,theta,dx,dy;
    	float G=0.0000000000667;
    
      ForceSum.i=0;
      ForceSum.j=0;
    
      for (i=0;i<NUMPARTS;i++)
      {
        for (j=0;j<NUMPARTS;j++)
        {
          if (j!=i) 
          {
    				fxsign=0;
    				fysign=0;
    				if ((Particles[j].x - Particles[i].x) < 0) fxsign++;
    				if ((Particles[j].y - Particles[i].y) < 0) fysign++;
    				dx = abs(Particles[i].x - Particles[j].x);
            dy = abs(Particles[i].y - Particles[j].y);
            theta = atan(dy/dx); //theta in radians
            r = sqrt((dx * dx) + (dy * dy)); //distance
    				F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));
    				printf("F= %f, G= %f m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);
            //dx = cos(rad(theta)) * r;
            //dy = sin(rad(theta)) * r;
    				dx = cos(rad(theta)) * F;
            dy = sin(rad(theta)) * F;
    				//if (i==2) printf("compareing to %d, r - %f, force x %f, force y %f, theta %f\n",j,r,dx,dy,theta); //###DEBUG
    				if (fxsign==0)
    					ForceSum.i += dx;
    				else
    					ForceSum.i -= dx;
    
    				if (fysign==0)
            	ForceSum.j += dy;
    				else
    					ForceSum.j -= dy;
    				//if (i==2) printf("forcesum %f, %f\n",ForceSum.i,ForceSum.j); //###DEBUG
          }
        }
        MoveForce[i].i = ForceSum.i;
        MoveForce[i].j = ForceSum.j;
      }
    }
    And here is the whole source:
    Code:
    /* nbody.c by Jack Carrozzo */
    
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    
    #define NUMPARTS 3 //how many particles
    #define GRID_X   600
    #define GRID_Y   400
    #define MASS     1000
    #define TIMESTEP 0.04 //seconds per step
    #define rad(d)   (4.1415926 * d) / 180
    #define DEBUG    0
    
    struct PartType
    {
      int    x; //x coord for body
      int    y; //y coord
      int    m; //mass
      double vi;//i-hat of velocity vector
      double vj;//j-hat "                "
    };
    
    struct ForceType
    {
      double i;
      double j;
    };
    
    struct PartType Particles[NUMPARTS];
    struct ForceType ForceSum;
    struct ForceType MoveForce[NUMPARTS]; //where forces are stored until move()
    int ImageGrid[GRID_X][GRID_Y];
    int Steps,ERR;
    
    void InitParticles( void )
    {
      int i,j;
    
    	srand(time()*getpid());
      for (i=0;i<NUMPARTS;i++)
      {
        Particles[i].x = ((double)rand()/((double)pow(2,31))) * GRID_X;
        Particles[i].y = ((double)rand()/((double)pow(2,31))) * GRID_Y;
        Particles[i].m = MASS;
        Particles[i].vi = 0;
        Particles[i].vj = 0;
      }
    
    	for (i=0;i<GRID_X;i++)
    	{
    		for (j=0;j<GRID_Y;j++)
    		{
    			ImageGrid[i][j]=0;
    		}
    	}
    }
    
    void ComputeForces( void )
    {
      int i,j,fxsign,fysign; //signs: 0 +, 1 -
      double r,F,theta,dx,dy;
    	float G=0.0000000000667;
    
      ForceSum.i=0;
      ForceSum.j=0;
    
      for (i=0;i<NUMPARTS;i++)
      {
        for (j=0;j<NUMPARTS;j++)
        {
          if (j!=i) 
          {
    				fxsign=0;
    				fysign=0;
    				if ((Particles[j].x - Particles[i].x) < 0) fxsign++;
    				if ((Particles[j].y - Particles[i].y) < 0) fysign++;
    				dx = abs(Particles[i].x - Particles[j].x);
            dy = abs(Particles[i].y - Particles[j].y);
            theta = atan(dy/dx); //theta in radians
            r = sqrt((dx * dx) + (dy * dy)); //distance
    				F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));
    				printf("F= %f, G= %f m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);
            //dx = cos(rad(theta)) * r;
            //dy = sin(rad(theta)) * r;
    				dx = cos(rad(theta)) * F;
            dy = sin(rad(theta)) * F;
    				//if (i==2) printf("compareing to %d, r - %f, force x %f, force y %f, theta %f\n",j,r,dx,dy,theta); //###DEBUG
    				if (fxsign==0)
    					ForceSum.i += dx;
    				else
    					ForceSum.i -= dx;
    
    				if (fysign==0)
            	ForceSum.j += dy;
    				else
    					ForceSum.j -= dy;
    				//if (i==2) printf("forcesum %f, %f\n",ForceSum.i,ForceSum.j); //###DEBUG
          }
        }
        MoveForce[i].i = ForceSum.i;
        MoveForce[i].j = ForceSum.j;
      }
    }
    
    void MoveParticles( void )
    {
      int i;
      double A; //acceleration
    
      for (i=0;i<NUMPARTS;i++)
      {
        //F=MA, A=F/M
        A = (MoveForce[i].i / Particles[i].m);
        Particles[i].vi += (A * TIMESTEP);
        A = (MoveForce[i].j / Particles[i].m);
    		//###DEBUG!!!
    		A *= 25;
        Particles[i].vj += (A * TIMESTEP);
    
    		Particles[i].x += (int)(Particles[i].vi);
    		Particles[i].y += (int)(Particles[i].vj);
    
    		//if (i==2) printf("part 2, y %f, f %f, %f.\n",Particles[i].vj,MoveForce[i].i,MoveForce[i].j);
      }
      Steps++;
    }
    
    void DataOut( void )
    {
      int i,j;
    	char filename[100];
    
    	if (DEBUG)
    	{
    		printf("\n   [+] Time: %d [+]\n",Steps);
    		for (i=0;i<NUMPARTS;i++)
    		{
    			printf("----------------------\n");
    			printf("Particle Num     : %d\n",i);
    			printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
    			printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
    			printf("----------------------\n");
    		}
    	}
    
    	CreateBitmap(GRID_X,GRID_Y);
    
    	for (i=0;i<NUMPARTS;i++)
      {
        for (j=0;j<NUMPARTS;j++)
        {
    			SetPixel(i, j, 0, 0, 0); //blackout image
    		}
    	}
    
    	for (i=0;i<NUMPARTS;i++)
    	{
    		if (Particles[i].x >= 0 && Particles[i].y >= 0 && Particles[i].x <= GRID_X && Particles[i].y <= GRID_Y)
    		{
    			SetPixel((Particles[i].x), (Particles[i].y), 255, 255, 255);
    			SetPixel((Particles[i].x + 1), (Particles[i].y), 255, 255, 255);
    			SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 255, 255, 255);
    			SetPixel((Particles[i].x), (Particles[i].y -1), 255, 255, 255);
    
    			/*if (i==0)
    			{
    				SetPixel((Particles[i].x), (Particles[i].y), 255, 100, 100);
    				SetPixel((Particles[i].x + 1), (Particles[i].y), 255, 100, 100);
    				SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 255, 100, 100);
    				SetPixel((Particles[i].x), (Particles[i].y -1), 255, 100, 100);
    			}
    			if (i==1)
    			{
    				SetPixel((Particles[i].x), (Particles[i].y), 100, 255, 100);
    				SetPixel((Particles[i].x + 1), (Particles[i].y), 100, 255, 100);
    				SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 100, 255, 100);
    				SetPixel((Particles[i].x), (Particles[i].y -1), 100, 255, 100);
    			}
    			if (i==2)
    			{
    				SetPixel((Particles[i].x), (Particles[i].y), 100, 100, 255);
    				SetPixel((Particles[i].x + 1), (Particles[i].y), 100, 100, 255);
    				SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 100, 100, 255);
    				SetPixel((Particles[i].x), (Particles[i].y -1), 100, 100, 255);
    			}*/
    		}
    	}
    
    	sprintf(filename,"img/step%d.bmp",Steps);
    
    	WriteBitmapToFile(&filename);
    
    	DeleteBitmap();
    }
    
    //JACK: then take ImageGrid to setpoint(), REMEMBER TO zero IMAGE GRID!!
    
    
    
    //go to bottom to see main
    
    // bitmap.c
    // --------
    // author: Levi Lansing
    // email: levi_is@hotmail.com
    // last updated: 10/15/03
    //
    // DESCRIPTION:
    // - a simple set of bitmap opperations
    // - outputs a standard windows bitmap file
    
    //#include <stdio.h>
    
    typedef int BOOL;
    #define TRUE	1
    #define FALSE	0
    
    typedef unsigned short WORD;
    typedef unsigned long DWORD;
    typedef unsigned char BYTE;
    
    
    //bitmap headers
    typedef struct tag_BITMAPFILEHEADER {
    	WORD     bfType;
    	DWORD      bfSize;
    	WORD     bfReserved1;
    	WORD     bfReserved2;
    	DWORD      bfOffBits;
    } BITMAPFILEHEADER;
    
    typedef struct tag_BITMAPINFOHEADER{
    	DWORD      biSize;
    	DWORD      biWidth;
    	DWORD      biHeight;
    	WORD     biPlanes;
    	WORD     biBitCount;
    	DWORD      biCompression;
    	DWORD      biSizeImage;
    	DWORD      biXPelsPerMeter;
    	DWORD      biYPelsPerMeter;
    	DWORD      biClrUsed;
    	DWORD      biClrImportant;
    	DWORD      bmiColors;
    } BITMAPINFOHEADER;
    
    
    
    //some global variables to keep track of the bitmap info
    //don't mess with this stuff
    DWORD b_height;
    DWORD b_width;
    DWORD b_widthbytes;
    BYTE* b_pBits = NULL;
    
    // CALL THIS FIRST!
    // creates a bitmap and allocates the memory for it
    void CreateBitmap(int width,int height)
    {
    	if(b_pBits)
    		free(b_pBits);
    
    	b_width = width;
    	b_height = height;
    	b_widthbytes = b_width*3;
    	if(b_width%4)
    		b_widthbytes += 4-(b_width%4);
    	b_pBits = (BYTE*)malloc(b_widthbytes*height);
    }
    
    //frees the allocated memory for the bitmap
    void DeleteBitmap()
    {
    	if(b_pBits)
    		free(b_pBits);
    	b_pBits = NULL;
    }
    
    //sets the color of a single pixel in the bitmap
    // colors red,green,and blue accept values between 0 (dark) and 255 (light) only!
    void SetPixel(int x, int y, BYTE red, BYTE green, BYTE blue)
    {
    	if(!b_pBits)
    	{
    		printf("Error: you forgot to create the bitmap first!\n");
    		return;
    	}
    	if(x>=b_width || y>=b_height || x<0 || y<0)
    	{
    		printf("pixel (%d,%d) is out of bounds- bitmap is only %dx%d!\n",x,y,b_width,b_height);
    		return;
    	}
    
    	BYTE* pBits = b_pBits + y*b_widthbytes + x*3;
    	pBits[0] = (BYTE) blue;
    	pBits[1] = (BYTE) green;
    	pBits[2] = (BYTE) red;
    }
    
    BOOL WriteBitmapToFile(char* filename)
    {
    	if(!b_pBits)
    	{
    		printf("Error: you forgot to create the bitmap first!\n");
    		return;
    	}
    
    	BITMAPFILEHEADER bmfh;
    	BITMAPINFOHEADER bmi;
    	int i;
    
    	FILE* fp = fopen(filename,"wb");
    	if(!fp)
    		return FALSE;
    
    	bmfh.bfType = 0x4D42;
    	bmfh.bfSize = 40+14+b_widthbytes*b_height;
    	bmfh.bfReserved1 = 0;
    	bmfh.bfReserved2 = 0;
    	bmfh.bfOffBits = 14+40;
    
    	//write the file header
    	fwrite(&bmfh.bfType,2,1,fp);
    	fwrite(&bmfh.bfSize,4,1,fp);
    	fwrite(&bmfh.bfReserved1,2,1,fp);
    	fwrite(&bmfh.bfReserved2,2,1,fp);
    	fwrite(&bmfh.bfOffBits,4,1,fp);
    
    	bmi.biSize = 40;
    	bmi.biBitCount = 24;
    	bmi.biClrImportant = 0;
    	bmi.biClrUsed = 0;
    	bmi.biCompression = 0;
    	bmi.biHeight = b_height;
    	bmi.biPlanes = 1;
    	bmi.biSizeImage = b_width*b_height*3;
    	bmi.biWidth = b_width;
    	bmi.biXPelsPerMeter = 3780;
    	bmi.biYPelsPerMeter = 3780;
    	bmi.bmiColors = 0;
    
    	//write the bitmap header
    	fwrite(&bmi.biSize,40,1,fp);
    
    	//bitmap files are upside down! flip it over!
    	BYTE* pBits = b_pBits + (b_height-1)*b_widthbytes;
    	for(i=0; i<b_height; i++)
    	{
    		fwrite((void*)pBits,b_widthbytes,1,fp);
    		pBits = pBits - b_widthbytes;
    	}
    	fclose(fp);
    
    	return TRUE;
    }
    
    int main( int argc, char *argv[] )
    {
    	int i, iters;
    
    	printf("Simple NBODY solver by Jack Carrozzo.\n\n");
    
    	if (argc!=2)
    	{
    		printf("Usage: ./nbody <iterations\n");
    		exit(1);
    	}
    
    	iters = atoi(argv[1]);
    
    	InitParticles();
    	printf("\n   [+] Time: %d [+]\n",Steps);
    		for (i=0;i<NUMPARTS;i++)
    		{
    			printf("----------------------\n");
    			printf("Particle Num     : %d\n",i);
    			printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
    			printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
    			printf("----------------------\n");
    		}
    
    	DataOut();
    
    	for (i=0;i<iters;i++)
    	{
    		ComputeForces();
    		MoveParticles();
    		DataOut();
    	}
    
    	//Final
    	printf("\n   [+] Time: %d [+]\n",Steps);
    		for (i=0;i<NUMPARTS;i++)
    		{
    			printf("----------------------\n");
    			printf("Particle Num     : %d\n",i);
    			printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
    			printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
    			printf("----------------------\n");
    		}
    
    	printf("Done.\n");
    
      return 0;
    }
    Thank you very much, sorry I'm stupid,

    -Jack Carrozzo
    jack [{@}] crepinc.com

  2. #2
    Registered User
    Join Date
    Mar 2004
    Posts
    536
    Quote Originally Posted by crepincdotcom
    OK, I admit right now that I am wrong.

    That said, I can find absolutly no reason when this should be. The function excerpted below should compute the Force F force each set of 2 particels in a system. Yet, no matter how many ways I've tried it, the line
    Code:
    F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));
    always returns 0!
    Maybe it's not zero. Try using %e or %g for the print statement:

    Code:
      printf("F= %e, G= %e m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);
    Regards,

    Dave

  3. #3
    & the hat of GPL slaying Thantos's Avatar
    Join Date
    Sep 2001
    Posts
    5,681
    As a test try setting G to a larger number, say 1.

    Off the cuff without really looking I wouldn't be surpised if the number was too small which grav forces are between particles (unless the particals are say the earth and moon)

  4. #4
    .
    Join Date
    Nov 2003
    Posts
    307
    float doesn't "see" precision beyond 6 digits. limits.h will show you FLT_DIG which tells you how many digits of precision float variables have on your system.

    use double.

  5. #5
    & the hat of GPL slaying Thantos's Avatar
    Join Date
    Sep 2001
    Posts
    5,681
    He is using a double Jim

  6. #6
    Registered User
    Join Date
    Mar 2004
    Posts
    536
    Quote Originally Posted by Thantos
    He is using a double Jim
    But the he's using %f to print numbers whose magnitudes are are less than 1.0e-6. So it shows 0.000000.


    with %e everywhere:

    Code:
    F= 3.268262e-10, G= 6.670000e-11 m1= 1000, m2= 1000, r= 451.756572
    
       [+] Time: 2 [+]
    ----------------------
    Particle Num     : 0
    Particle Position: (328,185)
    Particle Velocity: (4.218655e-13,-8.258777e-14)
    ----------------------
    ----------------------
    Particle Num     : 1
    Particle Position: (38,375)
    Particle Velocity: (4.923984e-13,-1.053544e-13)
    ----------------------
    ----------------------
    Particle Num     : 2
    Particle Position: (428,147)
    Particle Velocity: (-2.067952e-29,4.846761e-30)
    Regards,

    Dave

  7. #7
    & the hat of GPL slaying Thantos's Avatar
    Join Date
    Sep 2001
    Posts
    5,681
    Well that has other problems also since the exponent part of a double has more bits then a float. So basically printf will be decoding it improperly

  8. #8
    Registered User
    Join Date
    Mar 2004
    Posts
    536
    Quote Originally Posted by Thantos
    Well that has other problems also since the exponent part of a double has more bits then a float. So basically printf will be decoding it improperly
    printf() works with %e, %f, %g whether it's given a double or float (since float arguments are always promoted to double whenever printf() is called).

    And I don't propose to determine whether the answers are correct; I just wanted to point out how he can get some output from his program so that he can check the results.

    Regards,

    Dave
    Last edited by Dave Evans; 05-16-2005 at 09:33 AM.

  9. #9
    Registered User crepincdotcom's Avatar
    Join Date
    Oct 2003
    Posts
    94
    Thanks Guys, so that fixes the displayed problem. But somewhere in the I still have type issues, because even after several iterations, all the ForceSums are still zero or my personal favorite, negative zero. Any ideas?

    By the way Thantos, that value for G is the gravitational constant, and when multiplyed by the masses of the 2 bodies it becomes sufficient enough to putt the bodies together.

    Thanks again,

    -Jack Carrozzo
    jack [{@}] crepinc.com
    -Jack C
    jack {at} crepinc.com
    http://www.crepinc.com

  10. #10
    ATH0 quzah's Avatar
    Join Date
    Oct 2001
    Posts
    14,826
    You mean to say that when you keep multiplying a number by a decimal less than one, it eventually goes to zero?! No! Tell me it's not so!

    Quzah.
    Hope is the first step on the road to disappointment.

  11. #11
    Registered User
    Join Date
    Mar 2004
    Posts
    536
    Quote Originally Posted by quzah
    You mean to say that when you keep multiplying a number by a decimal less than one, it eventually goes to zero?! No! Tell me it's not so!

    Quzah.
    That wasn't the case for the original problem (see my previous post), but it can happen. It's called floating point underflow:

    Code:
    #include <stdio.h>
    
    int main()
    {
      double x = 1.0;
    
      int i;
    
      for (i = 0; i < 10000; i++) {
        x *= 0.1;
        if (x == 0.0) {
          printf("for i = %d, x = %e\n", i, x);
          break;
        }
      }
      printf("After the loop, i = %d, x = %e\n", i, x);
      return 0;
    }

    Regards,

    Dave
    Last edited by Dave Evans; 05-17-2005 at 08:05 AM.

  12. #12
    & the hat of GPL slaying Thantos's Avatar
    Join Date
    Sep 2001
    Posts
    5,681
    I think quzah was being sarcastic

  13. #13
    Registered User
    Join Date
    Mar 2004
    Posts
    536
    Quote Originally Posted by Thantos
    I think quzah was being sarcastic
    Ok by me, but I can see how the original poster might have missed it. (By the way, if you use x = -1.0 in my example, you do get to see negative zero.)

    Regards,

    Dave

  14. #14
    Registered User crepincdotcom's Avatar
    Join Date
    Oct 2003
    Posts
    94
    So to recap, if I simply use a double rather than float it should work?

    Sorry to not have figured this out myself.

    Thanks guys, interesting discussion too ;-)
    -Jack C
    jack {at} crepinc.com
    http://www.crepinc.com

  15. #15
    and the hat of wrongness Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    32,530
    > 0.0000000000667
    As you're probably aware now, floating point numbers have a number of decimal digits of precision (in the case of doubles, it's about 15).

    In this particular number, you're using just 3 decimal digits (all good so far).

    > (double)(Particles[i].m) * (double)(Particles[j].m)))
    Imagine the partial result of this is huge, I mean something like
    123450000000000.000

    Now here comes the grief
    > 0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m)))
    In attempting to normalise your small number to be on the same level as your large number (so the calculation can be done), it gets moved more than 15 decimal digits away from its original value, and as a result, winds up being zero!.

    > F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));
    Since it's just a series of multiplies and divides, then rearranging the expression might keep more significant digits overlapping, thus giving a better answer.
    Say something like
    F = (((0.0000000000667 * (double)(Particles[i].m) / r * (double)(Particles[j].m) / r));


    http://cch.loria.fr/documentation/IE...M/goldberg.pdf
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.
    I support http://www.ukip.org/ as the first necessary step to a free Europe.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Math
    By knightjp in forum A Brief History of Cprogramming.com
    Replies: 16
    Last Post: 04-01-2009, 05:36 PM
  2. Help with C++ Math
    By aonic in forum C++ Programming
    Replies: 4
    Last Post: 01-29-2005, 03:40 AM
  3. Basic Math Problem. Undefined Math Functions
    By gsoft in forum C Programming
    Replies: 1
    Last Post: 12-28-2004, 02:14 AM
  4. Math Header?
    By Rune Hunter in forum C++ Programming
    Replies: 26
    Last Post: 09-17-2004, 06:39 AM
  5. toughest math course
    By axon in forum A Brief History of Cprogramming.com
    Replies: 12
    Last Post: 10-28-2003, 09:06 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21