maths inside a torus (donut)

[AeonFlux1212]

Hi! I'm trying to verify if I have the correct formula to determine whether a (x,y,z) coordinate is "inside" a torus with a minor radius of 3 and a major radius of 5 (see http://en.wikipedia.org/wiki/Torus). Part of my output seems ok but I get some weird results too... here's the code and the output

#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

const long double r = 3;
const long double R = 5;

int main(int argc, char** argv)
{
	fstream myFile("test.txt", ios::out | ios::trunc); //ignore the fstream object there I just didn't take it out...
	
	long double negXMin = 0;
	long double negXMax = -10;
	
	long double negYMin = 0;
	long double negYMax = -10;
	
	long double posXMin = 10;
	long double posXMax = 0;
	
	long double posYMin = 10;
	long double posYMax = 0;
	
	long double zMin = 0;
	long double zMax = 0;
		
	if(myFile.is_open())
	{
		for(long double z = -5; z <= 5; z += 0.1)
		{
			for(long double y = -10; y <= 10; y += 0.1)
			{
				for(long double x = -10; x <= 10; x += 0.1)
				{
					if(pow(R - sqrt(pow(x, 2) + pow(y, 2)), 2) + pow(z, 2) <= pow(r, 2)) //here's "the formula"
					{
						if(x < negXMin) negXMin = x;
						if(x < 0 && x > negXMax) negXMax = x;
						
						if(y < negYMin) negYMin = y;
						if(y < 0 && y > negYMax) negYMax = y;
						
						if(x > 0 && x < posXMin) posXMin = x;
						if(x > posXMax) posXMax = x;
						
						if(y > 0 && y < posYMin) posYMin = y;
						if(y > posYMax) posYMax = y;
						
						if(z < zMin) zMin = z;
						if(z > zMax) zMax = z;
					}
				}
			}
		}
	}
	
	cout << negXMin << endl
		 << negXMax << endl
		 << negYMin << endl
		 << negYMax << endl
		 << posXMin << endl
		 << posXMax << endl
		 << posYMin << endl
		 << posYMax << endl
		 << zMin << endl
		 << zMax << endl;
		 
	return 0;
}

-8
-0.1
-8
-0.1
5.55112e-016
7.9
5.55112e-016
7.9
-3
2.9

negXMin and negXMax respectively stands for "negative X min" and "negative X max", same for negYMin and negYMax, and posXMin means "positive X min", etc. These values should outline the overall shape of the torus, and should be respectively

-8
-2
-8
-2
2
8
2
8
-3
3

I suspect that the pow() and sqrt() functions are responsible for some of the gibberish I get as result, but I don't understand why I would get -0.1 as negXMax and negYMax... it's as if my formula to test if the coordinates are inside is incorrect... which is frustrating because it should be ok... lol. Does anyone has any idea what's going on there??

Thanks,
AeonFlux1212

[TheIdeasMan]

Hi,

First, why long double? Ordinary double has 15 or 16 significant places of precision, long double is only 2 more at 17 or 18.

Don't ever use floating point values in a for loop, they aren't exact and will cause you problems. It is easy to calc how iterations you need by using start, end and step values, then convert that to an int.

Also, the pow function is very inefficient for squaring numbers because it uses a series function to do it, so just do x * x instead .

Have you tried using a debugger?

[AeonFlux1212]

I tried long double because I would get -7.9 instead of the two -8 I get in the output, with ordinary doubles, so I tried to see if long double would fix it.

I'm not sure I understand how I can avoid floats in the for loops. How does start, end and step work ?

No, I haven't tried a debugger. The three for loops together generate around 4 millions values so it would be a bit complicated to check them out I think :-)

[AeonFlux1212]

OK! Thanks to TheIdeasMan I've fixed the floats in the for-loops problem, and it's working magnificently! The only problem left is the value I get for the X and Y values of what should be the inner layer of the donut. Here's the new code and output :

#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

const long double r = 3;
const long double R = 5;

int main(int argc, char** argv)
{
	fstream myFile("test.txt", ios::out | ios::trunc); //ignore the fstream object there I just didn't take it out...
	
	long double negXMin = 0;
	long double negXMax = -10;
	
	long double negYMin = 0;
	long double negYMax = -10;
	
	long double posXMin = 10;
	long double posXMax = 0;
	
	long double posYMin = 10;
	long double posYMax = 0;
	
	long double zMin = 0;
	long double zMax = 0;
	
	long double x, y, z;
		
	if(myFile.is_open())
	{
		for(short k = -50; k <= 50; k += 1)
		{
			for(short j = -100; j <= 100; j += 1)
			{
				for(short i = -100; i <= 100; i += 1)
				{
					x = (long double)i / 10;
					y = (long double)j / 10;
					z = (long double)k / 10;
					
					if((R - sqrt(x*x + y*y))*(R - sqrt(x*x + y*y)) + z*z <= r*r) //here's "the formula"
					{
						if(x < negXMin) negXMin = x;
						if(x < 0 && x > negXMax) negXMax = x;
						
						if(y < negYMin) negYMin = y;
						if(y < 0 && y > negYMax) negYMax = y;
						
						if(x > 0 && x < posXMin) posXMin = x;
						if(x > posXMax) posXMax = x;
						
						if(y > 0 && y < posYMin) posYMin = y;
						if(y > posYMax) posYMax = y;
						
						if(z < zMin) zMin = z;
						if(z > zMax) zMax = z;
					}
				}
			}
		}
	}
	
	cout << negXMin << endl
		 << negXMax << endl
		 << negYMin << endl
		 << negYMax << endl
		 << posXMin << endl
		 << posXMax << endl
		 << posYMin << endl
		 << posYMax << endl
		 << zMin << endl
		 << zMax << endl;
		 
	return 0;
}

-8
-0.1
-8
-0.1
0.1
8
0.1
8
-3
3

I should be getting

-8
-2
-8
-2
2
8
2
8
-3
3

I'll be looking for a new donut formula... meanwhile if anyone has an idea of the problem please let me know...

Thanks,
AeonFlux1212

[TheIdeasMan]

Hi,

You can format the output std::cout to show as many decimal places as you want, also the default is not to show trailing zero's.

http://www.cplusplus.com/reference/iomanip/setprecision/

How does start, end and step work ?

It is pretty easy:

const double StartX = -10.0;
const double EndX = 10.0;
const Step = 0.1;

double x = 1000.0;

NumberOfXIterations = static_cast<int>( (EndX -StartX) / Step ); // is 200

for (int a = 0; a < NumberOfXIterations ; ++a) {
     x = StartX + (a * Step) ;
}

No, I haven't tried a debugger. The three for loops together generate around 4 millions values so it would be a bit complicated to check them out I think :-)

You don't have to test with 4 million values, and you can set conditional breakpoints.

With the FP values, the inaccuracy of them causes problems on line 45 and 48. When you increment by what you think is 0.1, it is actually something like 0.1 + 5e-16. So once you increment x until you think it is zero, it is actually 5e-16. Now with the logic on line 45:

if(x > 0 && x < posXMin) posXMin = x;

x > 0 is true, and x < 10 is true, so now posXMin is set to 5e-16. With the next iteration, x is now 0.1+ 5e-16, so the x < 5e-16 part is now false, and posXMin remains at 5e-16 for the rest of the program.

To fix this you need to incorporate a precision value into your comparisons:

const double precison = 0.01; // set this to whatever you think is reasonable

// greater than zero
bool IsGreaterThanZero(double x, const double precision) {
if( (x > precision)  )  { 
   return true;
}
else {
    return false;
}
}

This is more concisely written with the ternary operator:

const double precison = 0.01; // set this to whatever you think is reasonable

// greater than zero
bool IsGreaterThanZero(const double x, const double precision) {
     return (x > precision) ? true : false;
}

For equality with zero:

//  equality with zero
bool IsEqualZero(double x, const double precision) {
     return (std::abs(x) < precision) ? true : false;
}

For comparison with a number:

//  comparison with a number
bool IsEqualNumber(const double x, const double number, const double precision) {
     return (std::abs(x -number) < precision) ? true : false;
}

I haven't compiled any of this code, hopefully haven't pulled any clangers!

Hope this makes sense, and helps you out.

[helios]

(x ? true : false) == x.

The formula is correct, it's just that, for example, (cos(90°) * R, sin(90°) * R, 0) is inside the torus, so if you ignore the other coordinates, the minimum x and y that are above or below zero will always be just below or above zero.

[AeonFlux1212]

@IdeasMan

This is more concisely written with the ternary operator:

const double precison = 0.01; // set this to whatever you think is reasonable

// greater than zero
bool IsGreaterThanZero(const double x, const double precision) {
     return (x > precision) ? true : false;
}

For equality with zero:

//  equality with zero
bool IsEqualZero(double x, const double precision) {
     return (std::abs(x) < precision) ? true : false;
}

For comparison with a number:

//  comparison with a number
bool IsEqualNumber(const double x, const double number, const double precision) {
     return (std::abs(x -number) < precision) ? true : false;
}

Thanks a lot for this, I'll keep that in mind.

Now with the logic on line 45:

if(x > 0 && x < posXMin) posXMin = x;

x > 0 is true, and x < 10 is true, so now posXMin is set to 5e-16. With the next iteration, x is now 0.1+ 5e-16, so the x < 5e-16 part is now false, and posXMin remains at 5e-16 for the rest of the program.

I understand. The only problem is that, according to the torus formula which is checked in my first if statement, I should be getting no values between -2 and 2 (for both x and y), because that's where the 'hole' in the donut should be... that's why I don't understand why I get such a small value. I think I fixed the FP problem (check my last post), but now I get 0.1 (the increment value) instead of 5e-16...

Anyways thanks a lot for your help! It's a math problem now I guess... I don't see anything wrong with the code...

[AeonFlux1212]

@helios

My God you're right... I just figured that out...



Thanks again for the input guys, @helios and @TheIdeasMan, this is solved... :-)

AeonFlux1212

[ne555]

> This is more concisely written with the ternary operator:
> return (x > precision) ? true : false;
Or simply return x > precision;