2

Suppose I have an irrational number like \sqrt{3}. As it is irrational, it has no decimal representation. So when you try to express it with a IEEE 754 double, you will introduce an error.

A decimal representation with a lot of digits is:

1.7320508075688772935274463415058723669428052538103806280558069794519330169088
  00037081146186757248575675...

Now, when I calculate \sqrt{3}, I get 1.732051:

#include <stdio.h> // printf
#include <math.h>   // needed for sqrt

int main() {
    double myVar = sqrt (3);
    printf("as double:\t%f\n", myVar);
}

According to Wolfram|Alpha, I have an error of 1.11100... × 10^-7.

Is there any way I can calculate the error myself?

(I don't mind switching to C++, Python or Java. I could probably also use Mathematica, if there is no simple alternative)

Just to clarify: I don't want a solution that works only for sqrt{3}. I would like to get a function that gives me the error for any number. If that is not possible, I would at least like to know how Wolfram|Alpha gets more values.

My try

While writing this question, I found this:

#include <stdio.h> // printf
#include <math.h>  // needed for sqrt
#include <float.h> // needed for higher precision

int main() {
    long double r = sqrtl(3.0L);
    printf("Precision: %d digits; %.*Lg\n",LDBL_DIG,LDBL_DIG,r);
}

With this one, I can get the error down to 2.0 * 10^-18 according to Wolfram|Alpha. So I thought this might be close enough to get a good estimation of the error. I wrote this:

#include <stdio.h> // printf
#include <math.h>  // needed for sqrt
#include <float.h>

int main() {
    double myVar = sqrt (3);
    long double r = sqrtl(3.0L);
    long double error = abs(r-myVar) / r;
    printf("Double:\t\t%f\n", myVar);
    printf("Precision:\t%d digits; %.*Lg\n",LDBL_DIG,LDBL_DIG,r);
    printf("Error:\t\t%.*Lg\n", LDBL_DIG, error);
}

But it outputs:

Double:     1.732051
Precision:  18 digits; 1.73205080756887729
Error:      0

How can I fix that to get the error?

2
  • So are you asking two questions in one here? How you calculate the error in your answer, and also what's wrong with the second block of code?
    – Mike
    Mar 7, 2013 at 16:22
  • @Mike: Yes, I ask two questions. I have made the second one more precise. So if the second one is answered, the first one is automatically answered. If the first one gets an answer, I don't need an answer for the second one. Mar 7, 2013 at 16:36

6 Answers 6

3

What every Programmer should know about Floating Point Arithmetic by Goldberg is the definite guide you are looking for.

https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/02Numerics/Double/paper.pdf

0
1

printf rounds doubles to 6 places when you use %f without a precision.

e.g.

double x = 1.3;
long double y = 1.3L;
long double err = y - (double) x;
printf("Error %.20Lf\n", err);

My output: -0.00000000000000004445

If the result is 0, your long double and double are the same.

2
  • In line three: Why do you cast x to double? It is a double, so (double) doesn't change anything, does it? Mar 9, 2013 at 8:01
  • @moose: Just to be explicit about what's happening.
    – teppic
    Mar 9, 2013 at 11:26
1

One way to obtain an interval that is guaranteed to contain the real value of the computation is to use interval arithmetic. Then, comparing the double result to the interval tells you how far the double computation is, at worst, from the real computation.

Frama-C's value analysis can do this for you with option -all-rounding-modes.

double Frama_C_sqrt(double x);

double sqrt(double x)
{
  return Frama_C_sqrt(x);
}

double y;

int main(){
  y = sqrt(3.0);
}

Analyzing the program with:

frama-c -val t.c -float-normal -all-rounding-modes
[value] Values at end of function main:
      y ∈ [1.7320508075688772 .. 1.7320508075688774]

This means that the real value of sqrt(3), and thus the value that would be in variable y if the program computed with real numbers, is within the double bounds [1.7320508075688772 .. 1.7320508075688774].

Frama-C's value analysis does not support the long double type, but if I understand correctly, you were only using long double as reference to estimate the error made with double. The drawback of that method is that long double is itself imprecise. With interval arithmetic as implemented in Frama-C's value analysis, the real value of the computation is guaranteed to be within the displayed bounds.

7
  • I remember a math package for C about 25+ years ago that used interval arithmetic. I wonder why that concept seems to have gone by the wayside?
    – supercat
    Sep 3, 2013 at 16:53
  • @supercat In the context of static analysis, the popular and modern thing to do is work on relational domains, where the values of variables are assigned not only ranges but relations. As a general technique for computing approximation errors, it may have fallen out of fashion for the same kind of reason: it give safe but over-approximated bounds for expressions where sub-expressions are related, one useless but typical example being x - x. Sep 3, 2013 at 17:46
  • Certainly I can appreciate that there are many cases where a pessimistic range-based evaluation would yield useless results, and I can see that x-x is a particularly nice example [without a means of ensuring that the two x values are equivalent, it expands to the range (min-max)..(max-min)]. Still, in many cases the goal is to know whether a result can be guaranteed to be above or below some value. Knowing that the entire range is above, or the entire range is below, may eliminate the need for more refined analysis.
    – supercat
    Sep 3, 2013 at 18:21
  • BTW, am I the only guy who really misses 80-bit floating-point math? Even if such numbers generally got padded out to 16 bytes, they'd still for many purposes be better than the .NET Decimal type.
    – supercat
    Sep 3, 2013 at 18:24
  • @supercat You might like this series of posts to verify that a IOCCC winning program does not have undefined behavior (posts in order from bottom to top): blog.frama-c.com/index.php?tag/donut And on the subject of 80-bit FP, I too find it extremely convenient and I am glad that history was such that we now have it in addition to standard single- and double-precision, instead of having only these. Sep 3, 2013 at 18:28
0

You have a mistake in printing Double: 1.732051 here printf("Double:\t\t%f\n", myVar);

The actual value of double myVar is

1.732050807568877281 //18 digits

so 1.732050807568877281-1.732050807568877281 is zero

3
0

According to the C standard printf("%f", d) will default to 6 digits after the decimal point. This is not the full precision of your double.

It might be that double and long double happen to be the same on your architecture. I have different sizes for them on my architecture and get a non-zero error in your example code.

0

You want fabsl instead of abs when calculating the error, at least when using C. (In C, abs is integer.) With this substitution, I get:

Double:     1.732051
Precision:  18 digits; 1.73205080756887729
Error:      5.79643049346087304e-17

(Calculated on Mac OS X 10.8.3 with Apple clang 4.0.)

Using long double to estimate the errors in double is a reasonable approach for a few simple calculations, except:

  • If you are calculating the more accurate long double results, why bother with double?
  • Error behavior in sequences of calculations is hard to describe and can grow to the point where long double is not providing an accurate estimate of the exact result.
  • There exist perverse situations where long double gets less accurate results than double. (Mostly encountered when somebody constructs an example to teach students a lesson, but they exist nonetheless.)

In general, there is no simple and efficient way to calculate the error in a floating-point result in a sequence of calculations. If there were, it would be effectively a means of calculating a more accurate result, and we would use that instead of the floating-point calculations alone.

In special cases, such as when developing math library routines, the errors resulting from a particular sequence of code are studied carefully (and the code is redesigned as necessary to have acceptable error behavior). More often, error is estimated either by performing various “experiments” to see how much results fluctuate with varying inputs or by studying general mathematical behavior of systems.

You also asked “I would like to get a function that gives me the error for any number.” Well, that is easy, given any number x and the calculated result x', the error is exactly x'x. The actual problem is you probably do not have a description of x that can be used to evaluate that expression easily. In your example, x is sqrt(3). Obviously, then, the error is sqrt(3) – x, and x is exactly 1.732050807568877193176604123436845839023590087890625. Now all you need to do is evaluate sqrt(3). In other words, numerically evaluating the error is about as hard as numerically evaluating the original number.

Is there some class of numbers you want to perform this analysis for?

Also, do you actually want to calculate the error or just a good bound on the error? The latter is somewhat easier, although it remains hard for sequences of calculations. For all elementary operations, IEEE 754 requires the produced result to be the result that is nearest the mathematically exact result (in the appropriate direction for the rounding mode being used). In round-to-nearest mode, this implies that each result is at most 1/2 ULP (unit of least precision) away from the exact result. For operations such as those found in the standard math library (sine, logarithm, et cetera), most libraries will produce results within a few ULP of the exact result.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.