How can I calculate the inaccurate range when working with floating-point

[nguyenndk]

I'm newcomer in this forum. And I'm not a native speaker so please easy with me if I'm not good in English.

I'm writing a program to check a shape is Hexagon or not (Hexagon is a shape that have 6 equal edge .. you can check it at http://en.wikipedia.org/wiki/Hexagon)
//------------------------------------------------------------------
- Input data is
First Line: n (number of hexagon)
2->n+1 Line: 6 point(x,y) (floating-point number)
- Output data is
A string with n character (Y or N - describe if it's hexagon)

Example:
<HEXAGONS.INP>
100
-3 1 6 6.19615 0 6.19615 9 1 0 -4.19615 6 -4.19615
0 6 0 -4 6 6 6 -4 -1 1 9 1
<HEXAGONS.OUT>
YN

//------------------------------------------------------------------
- I've made a short program to solve this .. but I counter a problem.
When I working with floating point and use some variable also some math expression. In the end .. the result of me have some inaccurate around 0.000002
with this input("x1 y1 x2 y2 .. x6 y1" of 6 point): -3 1 6 6.19615 0 6.19615 9 1 0 -4.19615 6 -4.19615

- To solve this, I add a AcceptedErrorRange value in my code and compare with it before print out. I manage to solve my problem with it but I can't tell the the value of AcceptedErrorRange for many different "Input Data".

//------------------------------------------------------------------
Here is my code : It print out the result to file <HEXAGONS.OUT> if you have file <HEXAGONS.INP> .. And all the things print in console is for debug only (you can set "#define InDebugMode 0" if you want to turn it off)

Note: I compile it in Linux environment with g++ (Maybe it still have some error in other environment because of it's library)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fstream>
using namespace std;

//---------- Define -------------
/* 
#define fInputHexagons(a,i) { \
	fin >> a##.Node[i].x; \
	fin >> a##.Node[i].y; \
	}

#define fOutputHexagons(a,i) { \
	printf(" -> x: %f - y: %f \n", a.Node[i].x,a.Node[i].y);   \
	}
*/

#define InDebugMode 0 /*1->Enable ; 0->Disable*/

#define Max_Hexagons 100

#define mSwap(a,b,type) { type temp; temp=a; a=b; b=temp; }

#define mComparePoint(p1,p2)   (p1.x==p2.x && p1.y==p2.y)?0:1  /*if (p1==p2) -> false; else -> true*/

#define AcceptedErrorRange 0.001 /*The Value of Error in Result that can accept (Because can't be exactly accurate when using real number(floating point))*/

//---------- Struct & Var -------

ifstream fin("HEXAGONS.INP");
ofstream fout("HEXAGONS.OUT");
int N_Hexagons = 1;
const int N_Node = 6;
char result[100]={0};

struct Point
{
	double x; double y;
};

struct Hexagon
{
	Point Node[6];
};

//----------- Function ----------
void fInputHexagons(Hexagon *hex)
{
	for(int i=0;i<N_Hexagons;i++)
		for(int j=0;j<N_Node;j++)
		{	
			fin >> hex[i].Node[j].x;	
			fin >> hex[i].Node[j].y;	
		}
}

void fOutputHexagons(Hexagon *hex)
{
	printf(">> Print Out \n");
	for(int i=0;i<N_Hexagons;i++)
	{
		for(int j=0;j<N_Node;j++)
			printf("Hexagons %d -> Node %d -> x: %3.3f - y: %3.3f \n",i,j,hex[i].Node[j].x,hex[i].Node[j].y);
		printf("\n");
	}
}

double fDistance(Point a, Point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}

void fSort(double arr[],int s) // Selection Sort
{
	int m;
	for(int i=0;i<s-1;i++)
	{
		m=i;
		for(int j=i+1;j<s;j++)
			if(arr[j]>arr[m])
				m=j;
		if(m!=i)
			mSwap(arr[m],arr[i],double);
	}
}

int fIsHexagonCheck_1(double arr[]) /*Check First Times: The hexagon must have one longest distance; 2 equal long distance; 2 equal short distance; 1 zero distance (to itself)*/
{
	if(arr[0]>arr[1] && arr[1]==arr[2] && arr[2]>arr[3] && arr[3]==arr[4] && arr[5]==0) 
		return 1;
	else
		return 0;
}

Point fFindCenterPoint(double max_distance,double arr_distance[], Hexagon hex) /*Return the Center Point of the Hexagon*/
{
	Point returnPoint;
	for(int i=0;i<N_Node;i++)
		if(arr_distance[i]==max_distance)
		{
			if(InDebugMode)	printf("\n i=%d \n  arr_distance[i]-> %f \n  max_distance -> %f \n",i,arr_distance[i],max_distance);			
			
			returnPoint.x = (hex.Node[0].x + hex.Node[i].x)/2;
			returnPoint.y = (hex.Node[0].y + hex.Node[i].y)/2;

			if(InDebugMode)	printf("\n CenterPoint => (%f,%f) \n",returnPoint.x,returnPoint.y);			
		}	
	return returnPoint;
}

int fIsHexagonCheck_2(Point x,Hexagon h) /*Check Second Times: The distance from any Node of Hexagon to Center Point must be equal*/
{	
	double d1,d2,ErrorRange;
	for(int i=0;i<N_Node;i++)
	{
		d1 = fDistance(x,h.Node[i]);
		if(InDebugMode)		printf("\n Distance from Center Point to Node[%d] -> %f",i,d1);
		for(int j=0;j<N_Node;j++)
		{
			d2 = fDistance(x,h.Node[j]);
			if(InDebugMode)		printf("\n\t Distance from Center Point to Node[%d] -> %f",j,d2);						
			ErrorRange = sqrt((d1-d2)*(d1-d2));
			if(InDebugMode)		printf("\n\t => Error Range Between Node[%d] <=> Node[%d] : %f",i,j,ErrorRange);
			if(ErrorRange>AcceptedErrorRange)
				return 0;
		}
	}
	return 1;
}

//---------- Main Program --------
int main(int argc,char **argv)
{	
	//--- Declare Local Variables ---
	int i,j;
	double aDistance[Max_Hexagons][N_Node]={0};
	
	// --- Input Material ---
	fin >> N_Hexagons ;
	Hexagon *a = (Hexagon*)malloc((N_Hexagons+1)*sizeof(Hexagon));	

	fInputHexagons(a);
	if(InDebugMode) fOutputHexagons(a);
	
	//--- Process ---
	for(i=0;i<N_Hexagons;i++)	
	{
		if(InDebugMode) printf("\n");
		for(j=0;j<N_Node;j++)
		{
			aDistance[i][j]=fDistance(a[i].Node[0],a[i].Node[j]);
			if(InDebugMode) printf("fDistance from Node %d -> %d = %3.3f\n",i,j,fDistance(a[0].Node[i],a[0].Node[j]));
		}
	}
	
	if(InDebugMode)	
	{
		printf("\n- Before using qsort\n");
		for(i=0;i<N_Hexagons;i++)
		{			
			printf("Hexagon %d\n",i);
			for(j=0;j<N_Node;j++)
				printf("%f ",aDistance[i][j]);
			printf("\n");
		}
	}
	
	//Backup the array of distance that've not sorted yet!
	double aDistance_Backup[Max_Hexagons][N_Node];
	for(i=0;i<N_Hexagons;i++)
		for(j=0;j<N_Node;j++)
			aDistance_Backup[i][j] = aDistance[i][j];

	for(i=0;i<N_Hexagons;i++)
		fSort(aDistance[i],N_Node);
		//qsort(aDistance[i],N_Node,sizeof(double),fCompare); /*Can't use qsort in this case because it's double type*/
	
	if(InDebugMode)
	{
		printf("\n- After using qsort\n");
		for(i=0;i<N_Hexagons;i++)
		{
			printf("Hexagon %d\n",i);
			for(j=0;j<N_Node;j++)
				printf("%f ",aDistance[i][j]);
			printf("\n");
		}		
	}

	for(i=0;i<N_Hexagons;i++)
	{
		if(fIsHexagonCheck_1(aDistance[i]))
		{
			if(InDebugMode)	printf("Hexagon %d maybe is the real Hexagon\n",i);
			Point centerPoint = fFindCenterPoint(aDistance[i][0],aDistance_Backup[i],a[i]);
			if(mComparePoint(centerPoint,a[i].Node[0]))
			{
				if(InDebugMode)	printf("\nCenter Point (%f,%f) ",centerPoint.x,centerPoint.y);
				if(fIsHexagonCheck_2(centerPoint,a[i]))
				{
					if(InDebugMode)	printf("\n\n*** Hexagon %d \"IS\" the real Hexagon ***\n",i);	
					result[i]='Y';
				}
				else
				{
					if(InDebugMode)	printf("\n\n*** Hexagon %d \"IS NOT\" be the real Hexagon ***\n",i);				
					result[i]='N';
				}
			}
			else if(InDebugMode)
				printf("\nThere are some error when find Center Point\n");
		}
		else
			result[i]='N';

	}
		
	for(int i=0;i<N_Hexagons;i++)
		fout << (char)result[i];
	fout << endl;

	//--- Close Object ---
	fin.close();
	fout.close();
	return 0;
}

//------------------------------------------------------------------

Do someone have an idea to help me, please?

Thanks a lot.

[jsmith]

Read this:

http://docs.sun.com/source/806-3568/ncg_goldberg.html

Search for machine epsilon

And then try to avoid using floating point arithmetic.

(In your above code, your mComparePoint macro is useless since you can't really
compare floating point numbers for strict equality like you can integers.)

[turbozedd]

I assume you meant "Regular Hexagons".

For "equality" of floats you need to do something like fabs(a-b) < eps ? true : false
For a point you could test with the float method twice or define a distance function and use that e.g. dist(p1,p2) < eps ? true : false

eps is a tolerance supplied by you (.000001)

You may also want to look at: http://www.cplusplus.com/reference/std/limits/numeric_limits/
It shows how to find your machine epsilon.



jsmith: It is impossible to avoid floating point arithmetic in numerical methods (well, not quite, but inadvisable).

[jsmith]

Eh, well more useful than IEEE floating point would be fixed precision, but at that point
you could just as easily use a bigint library and multiply all the inputs by 10^k, k being
the number of digits of precision you need.

Speed could be a concern, although from the little I've used GMP it seems pretty fast.

[nguyenndk]

Thanks you two very much. I'll use it to solve this problem.
However, in some other problem, it's hard to solve with this solution .. such as when we use floating-point continuously (like a variable store some value when we recursion, back-tracking or while loop).I think the longer our program run the less precision it have because we use the old value to calculate the new one .. So, at that time, how can we find the right epsilon-value?

-----------------------------------------------------------------
turbozedd: Yes, This is "Regular Hexagons"

jsmith: Like turbozedd said, it's very hard to avoid floating in many cases .. such as when I need to compute a square root or somethings like that .. using big number very difficult.

If I've time .. maybe I try to build some function that use 2 or more double (or something 32bit else) and connect them to nx32bit type, then I'll build some basic function on it to increase the precision. I wonder if it have some library to do that (I mean a big floating-point)

[jsmith]

The link I posted above talks about error introduced by various calculations.

Specifically, look at the section on cancellation and catastrophic cancellation. This section should
show you that while there are ways to rewrite expressions to increase their resulting accuracy, these
techniques do not apply for all input values, thus the formula you have to use depends on the inputs.

If you really, really, really need a lot of precision you will either
1) go insane trying to figure out the right formulas to use for the inputs, or
2) take my advice and use integer math (with a bigint library).

[turbozedd]

You need to use a fixed epsilon.

There are already implementations of quad precision if you want to mess around with them, but increasing the precision in that manner does not give any improvement (there are theorems about this) (see Higham, chapter 1 of Accuracy and Stability of Numerical Algorithms, http://books.google.com/books?id=epilvM5MMxwC&dq=higham+Accuracy+and+Stability+of+Numerical+Algorithms&printsec=frontcover&source=bl&ots=CGWnf2x1xS&sig=0KGxsMfbkN6xzFi-dRzTnbDgRjw&hl=en&ei=FzK5Sov8KYPWtgOHpIEg&sa=X&oi=book_result&ct=result&resnum=1#v=onepage&q=&f=false )

9 digits of precision is usually more than enough, double supports 15 digits in most implementations.
In the old days of single floats, people would compute an answer as a single and then use a specialized algorithm to "polish" the answer as a double, gaining a couple digits of accuracy. (note: this is not the same as the previous paragraph which is just doubling the word size everywhere)

Could you post the algorithm for deciding whether a set of 6 points is a hexagon? I'm not sure what your code is trying to do.



@jsmith:
Floating point has been developed by numerical analysts for 70 years (now IEEE std) because it is the best solution for the type of computing we do.

http://en.wikipedia.org/wiki/Floating_point

The advantage of floating-point representation over fixed-point (and integer) representation is that it can support a much wider range of values...when stored in the same space, floating-point numbers achieve their greater range at the expense of slightly less precision.

Unless you are calculating pi or something unlimited precision is not needed, but
1) Numerical analysts have already done alot of this-- unless you are do something cutting edge its not a problem
2) Software implementations of number systems are too slow and use too much memory.