# Thread: Math is simply 0

1. ## 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);
//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);
//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;

WORD     bfType;
DWORD      bfSize;
WORD     bfReserved1;
WORD     bfReserved2;
DWORD      bfOffBits;

DWORD      biSize;
DWORD      biWidth;
DWORD      biHeight;
WORD     biPlanes;
WORD     biBitCount;
DWORD      biCompression;
DWORD      biSizeImage;
DWORD      biXPelsPerMeter;
DWORD      biYPelsPerMeter;
DWORD      biClrUsed;
DWORD      biClrImportant;
DWORD      bmiColors;

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

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;

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;

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. 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. 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. 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. He is using a double Jim

6. 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. 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. 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

9. 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

10. 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.

11. 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

12. I think quzah was being sarcastic

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

15. > 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