5 ms·
Bisecting Floating Point Numbers
- adrianN 11y agoI 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.
- jacobolus 11y agoAny specific algorithms in CGAL you found particularly clever?
- castratikron 11y agoFloats 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; }
- ska 11y agoFloat may be evil, but if so they are often a necessary evil.
- chowells 11y agoI 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.
- zurn 11y agoHow big improvements do you feel improved hardware fixed point support could yield, and what form might the improvements take? Imagining for example that Intel decided to make a instruction set extension in their next CPU.
- jjuhl 11y agoI don't know the answer to your question. But if you can get access to a machine with a POWER CPU you could test it, since those CPUs contain hardware fixed-point arithmetic units.
- chowells 11y agoModern CPUs are already so complex and so fast that I can't imagine the potential speedup would be significant as compared to modern optimized fixed point libraries. Language support is a much bigger issue. Programmers are fundamentally lazy, and will always do the thing that takes the least work (for them, right now). So for fixed point numbers to really be used, that means that they need to be fully first class in the language. Literals need to exist, and with syntax no more complex than it is for other numerical literals. And of course, all operators must use the same syntax for both. This is hard to retrofit into a language. Java and C can't retrofit either in. C++, Python, and Ruby can do operators, but not literals. You have to get into pretty uncommon languages before support for additional numerical types becomes good. Haskell and Scala, for instance, can meet both of the requirements. Given that, I'd say any kind of widespread use of fixed point numbers in the cases where they're more appropriate than floats is at least 30 years away. Maybe 50.
- david-given 11y agoFloats 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/ZPL5LTK9gpm https://plus.google.com/118274112180368285367/posts/ZPL5LTK9...
- sparky_z 11y agoThis 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.
- anon4 11y agoWell, he did say if the function is smooth and you have a reasonably small interval to start from.
- jwmerrill 11y agoIf 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_z 11y 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.
- lindstorm 11y agoGood! Instead of (x1 + x2)/2, it is better to do x1 + (x2 - x1)/2 [source:kdbuo]
- lindstorm 11y agoGreat 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 http://www.cs.berkeley.edu/~wkahan/Mindless.pdf