Bisecting Floating Point Numbers

http://squishythinking.com/julia/2014/02/22/bisecting-floats/

by jwmerrill • 11 years ago

65 26 11 years ago

adrianN• 11 years ago

I like numerical methods like this because they sit neatly in the space between high level algorithmic ideas and low level implementation details.

The library CGAL has many clever algorithms with floating point numbers. Its strength lies in the mixture of numerical approximation and exact computation it does to solve problems.

jacobolus2• 11 years ago

Any specific algorithms in CGAL you found particularly clever?

castratikron• 11 years ago

Floats are evil.

Here's some more fun you can have:

int main() { float a = 0x7FFFFFFF; printf("%f %f %f\n", a, a - 64, a - 65); return 0; }

ska2• 11 years ago

Float may be evil, but if so they are often a necessary evil.

chowells2• 11 years ago

I dispute "often". I honestly can't think of cases other than certain scientific calculations where they are more appropriate than fixed-point numbers.

The fact that scientists had so much input on the early design of hardware and programming languages means that floats are often a performance hack that exchanges semantics for speed. But that's hardly necessary in most applications.

And if hardware and language support was better for fixed point numbers, floats would be seen as the bizarre special case that's only appropriate in rare cases.

david-given• 11 years ago

Floats are evil.

Here's a program:

  int main(int argc, char* argv[]) {
    double x = 0.75;
    double t1 = 0.50000025;
    int i, j;
    int n8 = 93 * 100;
    for (i = 0;; i++) {
      for (j = 0; j<n8; j++)
        x = sqrt(exp(log(x) / t1));
      printf("%d %f\n", i, x);
    }
  }

Watch what happens to the performance as i passes 4632.

The full saga is here; AFAICT it's a pathological edge case in glibc. https://plus.google.com/118274112180368285367/posts/ZPL5LTK9...

sparky_z2• 11 years ago

This doesn't work if the output of the function you're checking is a noisy black box. For example, if fn() creates a finite element model, analyses it with an external software package, and post-processes the results to get an output value, you're gonna need a pretty high tolerance for your stopping condition.

anon42• 11 years ago

Well, he did say if the function is smooth and you have a reasonably small interval to start from.

jwmerrill• 11 years ago

If your function isn't a mathematical function (i.e. its output is not completely determined by its arguments), then it's a little confusing to figure out what isolating a root would mean. But even in that case, bisect_root terminates, because it always halves the interval at each step, even if your input procedure is returning nonsense. Witness

    julia> bisect_root((x) -> rand() - 1/2, 1.0, 2.0)
    (1.1980511149266218,1.198051114926622)

    julia> bisect_root((x) -> rand() - 1/2, 1.0, 2.0)
    (1.595600613503659,1.5956006135036591)

sparky_z2• 11 years ago

"...it's a little confusing to figure out what isolating a root would mean. "

Well, to take an example from some of my past work as a structural engineer, say you wanted to find out how much wind load a structure could sustain until it buckled under its own weight. So you create a finite element model of the structure that applies the wind load and then ramps up the gravity until the structure fails. The input is the magnitude of the wind pressure. The output is the total gravity at failure (in units of g). subtract 1 from the output, and you have a function where the root is the level of wind load you're looking for.

You're point about it always terminating when you use the bisection method is well taken, however. I was using a hybrid that included some Newton's method iterations to reduce the total number of steps, since each "function evaluation" took a lot of computer time.

lindstorm2• 11 years ago

Good! Instead of (x1 + x2)/2, it is better to do x1 + (x2 - x1)/2 [source:kdbuo]

lindstorm2• 11 years ago

Great post!

It may interest you that William Kahan has written quite extensively on the pitfalls of iterative algorithms on floating point numbers; his discussion of Muller's pathological recurrence on pp. 14+ of this paper seems particularly relevant.

http://www.cs.berkeley.edu/~wkahan/Mindless.pdf

Bisecting Floating Point Numbers

Related Stories