1

I'm running into a numerical issue in a large C programming project. (This is statistical research, not homework for a class). One step involves calculating sqrt(x^2 + y) - x, which I need to be positive, but sometimes I get sqrt(x^2 + y) - x < 0 even when x > 0 and y > 0. Example:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(){
  float x = 1210.9088134,
        y = 0.0062529947;

  printf("x =\t\t\t%70.50f\n", x);
  printf("y =\t\t\t%70.50f\n", y);

  printf("x*x =\t\t\t%70.50f\n", x*x);
  printf("x*x + y =\t\t%70.50f\n", x*x + y);
  printf("sqrt(x*x + y) =\t\t%70.50f\n", sqrt(x*x + y));

  printf("sqrt(x*x + y) - x =\t%70.50f\n", sqrt(x*x + y) - x);
}

My output:

x =                        1210.90881347656250000000000000000000000000000000000000
y =                           0.00625299476087093353271484375000000000000000000000
x*x =                   1466300.12500000000000000000000000000000000000000000000000
x*x + y =               1466300.12500000000000000000000000000000000000000000000000
sqrt(x*x + y) =            1210.90880127282912326336372643709182739257812500000000
sqrt(x*x + y) - x =          -0.00001220373337673663627356290817260742187500000000

This output is littered with strange behavior. Highlights:

  1. I assigned y 0.0062529947 but it prints out as 0.00625299476087093353271484375.
  2. x*x + y prints out as the same value as x*x.
  3. sqrt(x*x + y) - x < 0.

Why are 1-3 happening?

I should mention: I ran this example on both a 64-bit Mac OX 10.9.4 machine with gcc version:

$ gcc --version
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 5.1 (clang-503.0.40) (based on LLVM 3.4svn)
Target: x86_64-apple-darwin13.3.0
Thread model: posix

and a 64-bit CentOS server with gcc version:

$ gcc --version
gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-4)
Copyright (C) 2010 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Also, compilation returned no errors or warnings on either machine:

$ gcc -lm -Wall -pedantic -ansi test.c -o test
landau
  • 5,636
  • 1
  • 22
  • 50
  • 2
    1. float point accuracy 2. y is too much smaller than x*x 3. same as question 1 – leezii Sep 23 '14 at 05:00
  • http://stackoverflow.com/questions/588004/is-floating-point-math-broken?rq=1 – Mat Sep 23 '14 at 05:00
  • 1
    `y` is less than the difference between two adjacent floats, in the vicinity of `x * x`. You could push this situation back somewhat by using `double` instead of `float`, but this problem is unavoidable. `float` and `double` have limited precision, they cannot represent all real numbers. – M.M Sep 23 '14 at 05:01
  • 1
    You may want to read http://en.wikipedia.org/wiki/IEEE_floating_point – Mine Sep 23 '14 at 05:01
  • Think about how a computer represents a real number internally and you will realize why there is no such thing as storing exactly `0.0062529947` for example. – M.M Sep 23 '14 at 05:02
  • 1
    Look at [**What Every Computer Scientist Should Know About Floating-Point ...**](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) That will help answer your question. – David C. Rankin Sep 23 '14 at 05:03
  • 1
    The previous wikipedia link should eventually lead you to this [en.wikipedia.org/wiki/Single-precision_floating-point](http://en.wikipedia.org/wiki/Single-precision_floating-point_format) – user3386109 Sep 23 '14 at 05:07

0 Answers0