“Problem happens when z is very close to zj but their distance is still greater than EPS, since at this situation the formula in hj contains subtraction of two very close (double) float numbers, which causes heavy loss of accuracy."
It seemed that the double-precision floating-point format was being challenged.
Double precision only gives 15-17 significant decimal digits precision.
(See wikipedia)
Sign bit: 1 bit
Exponent: 11 bits
Significant precision: 53 bits (52 explicitly stored) -> 15-17 significant decimal digits
However, we could have a more professional way to circumvent/reduce "subtractive cancellation".
See the example given below and it can be a way to do some modifications.
The codes:
#include <iostream>
#include <limits.h>
#include <float.h>
#include <iomanip>
#include <cmath>
using namespace std;
typedef std::numeric_limits< double > dbl;
double scale=1e16;
int main(){
double a = 1.111111112355343;
double b = 1.111111112344554;
cout<< setprecision(30)<< a << endl << setprecision(30) << b <<endl;
double expn = floor(log10(abs(a-b))); // Assuming that the two numbers are very close to each other and both are greater than 1. The difference is approaching to the limitation of significant digits.
double a0 = a*pow(10,-expn-1)-floor(a*pow(10,-expn-1));
double b0 = b*pow(10,-expn-1)-floor(b*pow(10,-expn-1));
double a1 = round(a0*pow(10,15+expn+1)); // get rid of rubbish
double b1 = round(b0*pow(10,15+expn+1)); // get rid of rubbish
cout<< setprecision(30)<< a1 << endl << setprecision(30) << b1 <<endl;
cout<<"The difference between two number : "<<endl<< setprecision(20)<< a-b <<endl;
cout<<"Reference Value :"<<endl
<<" 92686996014.459171378255630735" <<endl;
cout<<"Calculations : "<<endl;
cout<<"Direct Calc : "<<setprecision(30)<<1/(a-b)<<endl;
cout<<"Optimised Calc : "<<setprecision(30)<<1e15/(a1-b1)<<endl;
return 0;
}
The output:
1.11111111235534298913307793555
1.11111111234455406382437558932
55343
44554
The difference between two number :
1.0788925308702346229e-11
Reference Value: 92686996014.459171378255630735
Calculations :
Direct Calc : 92687637682.8190765380859375
Optimised Calc : 92686996014.45916748046875
Hope it could be useful.
Cheers,
Hui