Hi all,

I have encountered a peculiar problem when performing some simple arithmetic in a program I am developing. As the function of the program is pretty specific to my field I have included a smaller program which illustrates the same problem without the technicalities.

Given a function (in the mathematical sense of the word) e.g. f(x), and an interval on which the function is defined e.g. x1 to x2, the program divides this interval into n segments of equal width and searches for a zero crossing of the function, i.e. a point on the interval where there is a solution to the equation f(x) = 0. The program loops over each segment incrementing x from x1, evaluating f(x) as it goes. It stops when it detects a zero crossing and outputs the refined interval on which the solution lies. Here is the program...

//...problem code

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

int main(){

  int n = 1000; // number of segments
  double x, x1, x2, dx, fp, fc, xb1, xb2; 

  x1=-0.1; x2=0.1; dx=(x2-x1)/n; // inteval defined by x1 and x2, dx is width of segments
  fp=erf(x=x1); // error function defined in cmath
  
  for(int i=0; i<n; i++){ // loop over intervals
    fc=erf(x+=dx); // determine incremented value of erf()
    if(fc*fp<=0.0){ // sign change indicates root of function
      xb1=x-dx; // record bound and break from loop
      xb2=x;
      break;
    }
  }
  cout << xb1 << " " << xb2 << endl;
}

In this case f(x) is the erf() function defined in cmath, this has the property erf(0)=0. The problem arises when n is defined such that the program encounters x=0 (as above), when this happens x-dx is not 0.0 as I would expect, but is a very small value of order 1e-16 and whose exact value is dependent on the value of n. Any help on this matter would be very much appreciated.

Best wishes

It would help if you used code tags.
I mean, you almost made it, but you still fell short.
Consider using the "review post" feature to make sure what you think you're getting is what you're actually going to get.

It sounds like you are experiencing the joys of floating point math, where close enough is what counts. For example, how many threes behind the decimal point does it take to equal 1/3, etc? Similarly what value is 1 x 10^16? For all practical purposes it is zero, unless you have some fantastically small error tolerances. In other words 0.0 equals 0.0000003 if the delta is +/- 0.0001. When comparing two variables of type double and float using floating point math choose the delta you're comfortable with rather than looking for identity.

Comments
I concur

Hi, please accept my apologies for the posting mistake.

Thanks for your prompt reply, I have since done a little reading on floating point maths. Please correct me if I have misunderstood... From what I now understand, my problem arises where the program (see original post) increments x -- a comparatively large -ve number -- by a small +ve number dx. Eventually x becomes a small -ve number and the evaluation of x+= dx, amounts to the subtraction of two very nearly equal numbers. The result of this is a value whose only significant bits are the few lower order ones by which the two operands differed. The effect of these "rounding errors" is noticed only when the program encounters x+dx=0, because prior to this the rounding error << x+dx and its low order bits are right shifted away as the computer performs the arithmetic.

Is this the correct interpretation?

Sorry, I'm not familiar with right shifting and low order bits, etc. What you describe sounds plausable though. I just think of it as "rounding errors" associated with floating point math when working at extremely large and extremely small numbers or when trying to express rational numbers that don't have exact decimal equivalences, and I haven't really tried to get to the bottom of it by looking into gory details of binary math with bit shifting to the right or left, etc. Maybe someone else can affirm your research.

This question has already been answered. Start a new discussion instead.