The Daily WTF: Curious Perversions in Information Technology

2007-02-01 Reply Admin

Stephen Touset:

int bits;
for (bits = 0; number != 0; number >>= 1; bits++) ;

Dammit. I suck.

int bits;
for (bits = 0; number != 0; number >>= 1, bits++) ;

2007-02-01 Reply Admin

Here's another fun way:

unsigned int WTFunction(unsigned int n) { unsigned int x = n; x |= x>>1; x |= x>>2; x |= x>>4; x |= x>>8; x |= x>>16; x ^= x>>1; return (x < n ? x << 1 : x); }

There are other easy ways to do this that involve bit twiddling and optimize better, but one can't give away every secret.

2007-02-01 Reply Admin

Okay, okay. Don't everyone beg all at once:

unsigned int reverse(unsigned int v) { v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1); // swap consecutive pairs v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2); // swap nibbles ... v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4); // swap bytes v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8); // swap 2-byte long pairs v = ( v >> 16 ) | ( v << 16);

return v; }

unsigned int WTFunction(unsigned int v) { unsigned int x = reverse(v); x = reverse(x & -x); return x < v ? x << 1 : x; }

CAPTCHA: craaazy

2007-02-01 Reply Admin

whicker:
I had my chuckle for the day. You can't possibly be serious. The algorithm is wrong, period.

Yes, the algorithm is broken, but not the + 1 that Alex mentioned, period. You may only visit this site for a laugh, but that won't stop me from correcting wrong claims so that those who also want to learn something don't learn stuff that's the result of incorrect reasoning.

2007-02-01 Reply Admin

whicker:
Lighten up. That I was anticipating people giving partial credit for the poor implementation was what I found so humorous. I was doubting that the 'fencepost' thinking even came into play with the snippet of this article. Sheesh.

Actually, the last line of the snippet is close enough to mathematical correctness to suggest that the fencepost thinking came into play, although it doesn't help much because of all the other problems.

(Mathematically, this should be ceil(alog(range + 1) / alog(2)).)

The guy who wrote this had his logarithms right, too (at least mathematically - because of rounding errors this is useless here, of course).

In my opinion, the author of this code was either a code monkey who remembered the correct formula but failed to reproduce it properly (let alone think about the intricacies of floating-point calculations) or a mathematics/science type who has less of a clue about computer programming than (s)he'd like to admit.

2007-02-01 Reply Admin

100+ comments for this? GRG's problem statement and understanding of why the code fails are both wrong, and the list of 12 is stupid. Why for example should alog check for negative numbers if it's being passed abs()?

GRG's statement of the problem is to find n such that 2^(n-1) <= r < 2^n for given r, so n-1 <= log_2(r) < n, so n = log_2(r) + 1 where flooring is automatic; not log_2(r+1), which is what the code is doing. The code is bound to be wrong in all cases EXCEPT when r is one less than a power of 2. The 0.41% and the fact that GRG thinks 8 bits oughta be sufficient for 256 suggests GRG isnt stating the problem right, and that 1 less than the no. of bits is required, ie. a simple n = r ? log_2(r) : <default>, where <default> could be 1 or 0 depending on the application. This makes the code wrong when r is one less than power of 2.

Limits etc oughta be checked before, not on this line. And rounding off errors can be adjusted with a simple addition of a small enough number, n = r ? log(r)/log(2) + EPS : 0, where EPS probably depends upon int limits (perhaps -log(INT_MAX-1) + log(INT_MAX)).

2007-02-02 Reply Admin

Stephen Touset:
Dammit. I suck.
int bits;
for (bits = 0; number != 0; number >>= 1, bits++) ;

Not too much. But considering that we are talking about C here, it should be either:

/* Old style */
int bits;
for (bits = 0; number; number >>= 1, bits++) ;

or

/* Safe style */
int bits;
for (bits = 0; 0 != number; number >>= 1, bits++)
   ;

2007-02-02 Reply Admin

/* Binary search to find highest 1-bit */ unsigned wtf_function(unsigned u) { unsigned sz = sizeof(unsigned)*8/2, nbits = 0;

do { if(u >> (sz + nbits)) nbits += sz; } while(sz >>= 1); return nbits + 1; }

2007-02-02 Reply Admin

Hmm. This actually made me go back and check my own code. I'm using a similar function for a quick-and-dirty convert netmask to CIDR script.

my $bits = 32-log(2**32-q2i($netmask))/log(2);

I'm sure it could be improved. Either way, I checked this and under the architecture this will run on, it's accurate to the 44th bit, which surely should be sufficient for 32-bit netmasks. Of course, I don't have to deal with rounding (unless some netmask is really screwed up, and if it is, I should catch it in the output anyway).

It's a hack, but then again, I'm not a coder. ;)

2007-02-02 Reply Admin

It's not the value 1 that can be represented in 0 bits. I'm saying that a range consisting of a single number can be represented in 0 bits. That is because you don't need any more information than the minimum value, consequently no data needs to be stored, hence 0 bits. It's presence is sufficient to encode the only possible value.

Ancient_Hacker · 2007-02-02 Reply Admin

Well, swell, but when I run this little test program:

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

void f() { float top, bot, ans, realb; int i, b, intans; b = 1; i = 0;

while( i < 30 ) { realb = (float) b; top = (float) log( realb ); bot = (float) log( 2.0 ); ans = top / bot; intans = (int) ans; if ( i == intans ) printf("loop: %3d, EQUALS %3d %25.20f\n", i, intans, ans ); else printf("loop: %3d, DOES NOT EQUAL %3d %25.20f\n", i, intans, ans ); b = b + b; i++; } }

int main( int argc, char * argv ) { f(); }

I get these results:

loop: 0, EQUALS 0 0.00000000000000000000 loop: 1, EQUALS 1 1.00000000000000000000 loop: 2, EQUALS 2 2.00000000000000000000 loop: 3, EQUALS 3 3.00000000000000000000 loop: 4, EQUALS 4 4.00000000000000000000 loop: 5, EQUALS 5 5.00000000000000000000 loop: 6, EQUALS 6 6.00000000000000000000 loop: 7, EQUALS 7 7.00000000000000000000 loop: 8, EQUALS 8 8.00000000000000000000 loop: 9, EQUALS 9 9.00000000000000000000 loop: 10, EQUALS 10 10.00000000000000000000 loop: 11, EQUALS 11 11.00000000000000000000 loop: 12, EQUALS 12 12.00000000000000000000 loop: 13, DOES NOT EQUAL 12 12.99999904632568400000 loop: 14, EQUALS 14 14.00000000000000000000 loop: 15, DOES NOT EQUAL 14 14.99999904632568400000 loop: 16, EQUALS 16 16.00000000000000000000 loop: 17, EQUALS 17 17.00000000000000000000 loop: 18, EQUALS 18 18.00000000000000000000 loop: 19, EQUALS 19 19.00000000000000000000 loop: 20, EQUALS 20 20.00000000000000000000 loop: 21, EQUALS 21 21.00000000000000000000 loop: 22, EQUALS 22 22.00000000000000000000 loop: 23, EQUALS 23 23.00000000000000000000 loop: 24, EQUALS 24 24.00000000000000000000 loop: 25, EQUALS 25 25.00000000000000000000 loop: 26, DOES NOT EQUAL 25 25.99999809265136700000 loop: 27, EQUALS 27 27.00000000000000000000 loop: 28, EQUALS 28 28.00000000000000000000 loop: 29, EQUALS 29 29.00000000000000000000

... and don't suggest taking ceil() !!

Ancient_Hacker · 2007-02-02 Reply Admin

Perhaps I should clarify things a bit.

As mentioned in the original posting, this general method was used about 54 times in different places throughout the code, with slight variations-- I only posted one of the instances. Sometimes the range was calculated correctly, as being hi - low + 1, othertimes not. Sometimes there was a "+1" or even "+ 2" inside the argument to log() to perhaps fix that error, sometimes it was a fudge factor to get the answer "more correct". Quite often there were commented out "WRITE(6,100) ANSWER" lines, showing the code was problematical until the right fudge factors were introduced. I apologize for all the red herrings spun off by my omitting this information.

There's been an alarming number of responses that gloss over and miss the whole point:

Although, in math theory, the formula is correct, IN A REAL ACTUALIZABLE COMPUTER, THE METHOD IS INCORRECT.

There is no way to guarantee the results of a log() function, calculated to a finite number of bits, then divided by log(2.0), is going to be exactly correct. In a distressing number of cases, the answer will be a bit or two off, which when either rounded UP or DOWN to an integer, is going to be off by an unacceptable amount (up to 100% off).

This is not my random speculation, I posted the results a few msgs above, showing unacceptable errors in three out of 20 cases.

IT'S BASICALLY INCORRECT

Also, there were some of the 54 instances that forgot to take abs() of the difference, which could result in zero or a negative number passed to log(), neither one of which will give acceptable results.

Alchymist · 2007-02-02 Reply Admin

Oliver:
Duston:
Could have been worse...could have been:
if (x == 1) { nbits = 1; } else { if (x == 2 || x == 3) { nbits = 2; } else { if (x >= 4 && x < 8) { . . .

But if that's the case, the function would be clearly understandable and easily verified as correct, much much better than the original WTF. So what if it takes log(N) operations to get the result for input N? We are talking about very important data here.

We can cut down on the number of comparisons, as well

if (x < 0) { // Some warning message or possibly sizeof(int) * 8} else if (x < 2) { nbits = 1; } else if (x < 4 ) { nbits = 2; } else if (x < 8) { . .

(or the FORTRAN equivalent).

How does this compare in speed to calculating multiple floating point operations using the formula given? I don't have the background to determine this? Anyone out there able to analyse it?

2007-02-02 Reply Admin

Sure enough,

 if (x<0) throw new NegativeRangeExeption();
 nbits = x>0 ? 1+(int)floor(log(x)/log(2.0)) : 0;

let's assume the compiler evals log(2.0) to a constant for us...

A typical example where unit testing would have been a good idea.

(1..12).each do |n|
  ((1<<n)..(2<<n-1)).each do |x|
    assert( nbits(x) == n+1 )
  end
end
</pre>

2007-02-02 Reply Admin

And you can use fancy-schmancy cmov to avoid the branch.

2007-02-02 Reply Admin

It's interesting to read how some of you guys have such a great distrust of floating point numbers, like if the processor could generate random results in the least significant bits.

If any junior programmers are reading this, please don't listen too much to the Dogmatic voices heard on this page. They are wrong. Floating point math can be used to get exact answers to certain questions, if you know what you are doing. It might be true that floating point solutions to the problem in this WTF were harder to read and harder to maintain. Bit it isn't so simple as: "floating point can never be used if exact results are required".

Consider that a typical floating point processor will return exact results for various floating point operations, if the results can be represented. This is quite a neat engineering feat, in my oppinion. There's some hard math behind the circuits that do for instance division or cosine. When the result cannot be represented exactly, there are specifications for how much off the result will be.

Finally, I want to present my C-language solution to the problem. You know that (at least intel-) floating point hardware has operations to convert from integer to floating point? Why not use that hardware:

int bits_reqd(unsigned int num) { double f=(double)num; //convert to double int exponent; frexp(f,&exponent); //get exponent return exponent; }

(Tested on an intel-processor and microsoft toolchain - gives exact results for all 32-bit unsigned ints)

2007-02-02 Reply Admin

Ancient_Hacker:
Well, swell, but when I run this little test program
... and don't suggest taking ceil() !!

Try this test program and output:

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

void main() { float num, den; int pow, b;

for(b = 1; b < 50; b++) { num = log(b); den = log(2); pow = ceil(num/den); printf("A range of %d numbers needs %d bits.\n", b, pow); } }

A range of 1 numbers needs 0 bits. A range of 2 numbers needs 1 bits. A range of 3 numbers needs 2 bits. A range of 4 numbers needs 2 bits. A range of 5 numbers needs 3 bits. A range of 6 numbers needs 3 bits. A range of 7 numbers needs 3 bits. A range of 8 numbers needs 3 bits. A range of 9 numbers needs 4 bits. A range of 10 numbers needs 4 bits. A range of 11 numbers needs 4 bits. A range of 12 numbers needs 4 bits. A range of 13 numbers needs 4 bits. A range of 14 numbers needs 4 bits. A range of 15 numbers needs 4 bits. A range of 16 numbers needs 4 bits. A range of 17 numbers needs 5 bits. A range of 18 numbers needs 5 bits. A range of 19 numbers needs 5 bits. A range of 20 numbers needs 5 bits. A range of 21 numbers needs 5 bits. A range of 22 numbers needs 5 bits. A range of 23 numbers needs 5 bits. A range of 24 numbers needs 5 bits. A range of 25 numbers needs 5 bits. A range of 26 numbers needs 5 bits. A range of 27 numbers needs 5 bits. A range of 28 numbers needs 5 bits. A range of 29 numbers needs 5 bits. A range of 30 numbers needs 5 bits. A range of 31 numbers needs 5 bits. A range of 32 numbers needs 5 bits. A range of 33 numbers needs 6 bits. A range of 34 numbers needs 6 bits. A range of 35 numbers needs 6 bits. A range of 36 numbers needs 6 bits. A range of 37 numbers needs 6 bits. A range of 38 numbers needs 6 bits. A range of 39 numbers needs 6 bits. A range of 40 numbers needs 6 bits. A range of 41 numbers needs 6 bits. A range of 42 numbers needs 6 bits. A range of 43 numbers needs 6 bits. A range of 44 numbers needs 6 bits. A range of 45 numbers needs 6 bits. A range of 46 numbers needs 6 bits. A range of 47 numbers needs 6 bits. A range of 48 numbers needs 6 bits. A range of 49 numbers needs 6 bits.

2007-02-02 Reply Admin

Or perhaps this would be a little better:

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

void main() { float num, den; int pow, b, i; b=1; for(i = 1; i < 32; i++) { num = log10(b); den = log10(2); pow = ceil(num/den); printf("A range of %d numbers needs %d bits.\n", b, pow); b=b+b; } }

A range of 1 numbers needs 0 bits. A range of 2 numbers needs 1 bits. A range of 4 numbers needs 2 bits. A range of 8 numbers needs 3 bits. A range of 16 numbers needs 4 bits. A range of 32 numbers needs 5 bits. A range of 64 numbers needs 6 bits. A range of 128 numbers needs 7 bits. A range of 256 numbers needs 8 bits. A range of 512 numbers needs 9 bits. A range of 1024 numbers needs 10 bits. A range of 2048 numbers needs 11 bits. A range of 4096 numbers needs 12 bits. A range of 8192 numbers needs 13 bits. A range of 16384 numbers needs 14 bits. A range of 32768 numbers needs 15 bits. A range of 65536 numbers needs 16 bits. A range of 131072 numbers needs 17 bits. A range of 262144 numbers needs 18 bits. A range of 524288 numbers needs 19 bits. A range of 1048576 numbers needs 20 bits. A range of 2097152 numbers needs 21 bits. A range of 4194304 numbers needs 22 bits. A range of 8388608 numbers needs 23 bits. A range of 16777216 numbers needs 24 bits. A range of 33554432 numbers needs 25 bits. A range of 67108864 numbers needs 26 bits. A range of 134217728 numbers needs 27 bits. A range of 268435456 numbers needs 28 bits. A range of 536870912 numbers needs 29 bits. A range of 1073741824 numbers needs 30 bits.

2007-02-02 Reply Admin

[image]

Ancient_Hacker · 2007-02-02 Reply Admin

Floating point math can be used to get exact answers to certain questions, if you know what you are doing.

Yes, that's true. And if my grandmother had wheels, she'd be a car.

You yourself seem to have not a clue re fp numbers.
To reiterate, fp numbers in every computer I've seen or heard of are stored as fixed length fractions times an exponent.

Now you might recall from math class that the Greeks waay back before Britney Spears, about 200BC, blew their minds by proving there are numbers, like sqrt(2) which CANNOT be represented by any fraction. Same with most log and trig functions. The results are TRANSCENDENTAL, an unending decimal expansion. You can't exactly represent a non-repeating decimal number in any computer. No how. No way. Just not possible.

So the basic concept of using log functions to calculate exact results is just plain wrong.

Ancient_Hacker · 2007-02-02 Reply Admin

I'm sure you're well-intentioned, but your example proves nothing. You only tested the powers of two, and on one computer. That proves nothing. To prove the correctness of this basic approach you'd have to try it on every computer in existence, in every IEEE and/or native rounding modes.

Using "ceil" just papers over the basic gotcha. All it would take to burst your bubble would be to have the ratio turn out to be one LSB above, instead of one below. pop

To repeat: you can't get exact output with inexact inputs. log() of anything is a transcendental number, impossible to exactly represent in a computer or anywhere else.

2007-02-02 Reply Admin

Answer to Ancient Hacker:

No, ofcourse sqrt(2) cannot be represented as a binary floating point number. So if your question is: "What is the square-root of 2", you will not get the correct result if you use typical computer-implemented fixed-length floating point numbers.

However, if your question is:

"what are the first 3 digits of the decimal expansion of the mathematical expression sqrt(2)"

...you will get perfectly accurate results to the question using floating point numbers in your calculations.

You seem to feel very strongly about this, and I am sorry if I have offended you.

Now, do you think the code I provided will ever fail, on any computer with binary 64-bit doubles? (and that's going to be almost all computers, decimal or hexadecimal floating point isn't used on typical modern computers).

Ancient_Hacker · 2007-02-02 Reply Admin

No, ofcourse sqrt(2) cannot be represented as a binary floating point number. So if your question is: "What is the square-root of 2", you will not get the correct result if you use typical computer-implemented fixed-length floating point numbers.

You're still not getting it. The problem isn't limited to "typical computer-implemented fixed-length floating point numbers". There is just NO WAY represent sqrt(2). Even if you put every bit of memory on the planet into your computer, you still couldn't represent sqrt(2).

"what are the first 3 digits of the decimal expansion of the mathematical expression sqrt(2) ..you will get perfectly accurate results to the question using floating point numbers in your calculations. "

No, you still don't get it. There are an infinite number of numbers where the sqrt() of them are 10.9999999999999XXXXX. A computer using 32 or 64 or 80 bit math is going to get the wrong answer, even to three digits.

When you need exact and predictable results, you can't use fp math, and especially not transcendental functions.

I'm sorry if many folks are little bit fuzzy about the concept of exact results. Such is the state of our craft.

Now, do you think the code I provided will ever fail, on any computer with binary 64-bit doubles?

Doesnt matter what I think, or what happens on any or all computers. The basic idea is foobar, applying patches like CEIL() is a kludge, and things like SetIEEERoundingMode( ROUND_UP ) will likely mess up the results.

Regards,

A_H

2007-02-02 Reply Admin

There is just NO WAY represent sqrt(2). Even if you put every bit of memory on the planet into your computer, you still couldn't represent sqrt(2).

See http://www.wolfram.com/products/mathematica/newin5/numberobjects.html for algebraic number representation. I think that's one of the non-typical ways Anders was alluding to.

And please stop generalizing in vain. It's true that one can't reason about FP implementation of every computer in the Universe, but once you rely on some standard, which in practical cases you do, you can reason quite a bit about the results of FP calculations. Read http://citeseer.ist.psu.edu/cache/papers/cs/23861/http:zSzzSzwww.math.rpi.eduzSzFacultyzSzHolmeszSzNumCompzSzMisczSzgoldberg.pdf/goldberg91what.pdf , try finding the words "Just as integer programs can be proven to be correct, so can floating-point programs, although what is proven in that case is that the rounding error of the result satisfies certain bounds."

2007-02-02 Reply Admin

There is just NO WAY represent sqrt(2). Even if you put every bit of memory on the planet into your computer, you still couldn't represent sqrt(2).

See http://www.wolfram.com/products/mathematica/newin5/numberobjects.html for algebraic number representation. I think that's one of the non-typical ways Anders was alluding to.

And please stop generalizing in vain. It's true that one can't reason about FP implementation of every computer in the Universe, but once you rely on some standard, which in practical cases you do, you can reason quite a bit about the results of FP calculations. Read http://citeseer.ist.psu.edu/cache/papers/cs/23861/http:zSzzSzwww.math.rpi.eduzSzFacultyzSzHolmeszSzNumCompzSzMisczSzgoldberg.pdf/goldberg91what.pdf , try finding the words "Just as integer programs can be proven to be correct, so can floating-point programs, although what is proven in that case is that the rounding error of the result satisfies certain bounds."

2007-02-02 Reply Admin

GRG's statement of the problem is wrong, and s/he needs algebraically the floor of log_2(range), which is not the number of bits required to hold the number. The code is doing log_2(range + 1) in an ideal world, and thus is wrong when range is one less than a power of 2. To get rid of floating point errors: (a) double precision must be used, ie. alog(range*1.0d0)/alog(2.0d0) (b) range=0 must be accounted for by using a default (c) int limits etc must be checked for, but this has to happen before the "twelve-on-a" line.

For all values of nonnegative ints within int limits, double precision should be precise enough without needing a fudge factor. If a fudge factor is still necessary, the code needs to do alog(range + fudge)/alog(2.0d0) where fudge is a small double precision number, and emphatically not 1.0d0.

Any other errors than above are GRG's imagination, and for what it's worth I dont think I'd trust the software in its new hands either :)

Ancient_Hacker · 2007-02-02 Reply Admin

See http://www.wolfram.com/products/mathematica/newin5/numberobjects.html

We were discussing explicit calculations in Fortran, not symbolic math in Mathematica.

It's true that one can't reason about FP implementation of every computer in the Universe, but once you rely on some standard, which in practical cases you do, you can reason quite a bit about the results of FP calculations.

Who is this "you"? Maybe Donald Knuth can predict the error in doing two alog() and a divide using IEEE math, but 99.999% of the hackers out there can't.

It's irresponsible to spread such statements.

2007-02-02 Reply Admin

Ancient_Hacker:
We were discussing explicit calculations in Fortran, not symbolic math in Mathematica.

I quote You quoting Anders: "The problem isn't limited to "typical computer-implemented fixed-length floating point numbers". There is just NO WAY.." I perceived that as an argument beyond Fortran.

Ancient_Hacker:
Who is this "you"? Maybe Donald Knuth can predict the error in doing two alog() and a divide using IEEE math, but 99.999% of the hackers out there can't.
It's irresponsible to spread such statements.

Of course, I'm not advocating use of floating point where not appropriate. I really suggest everybody think twice and consult their local Knuth for ways of avoiding FP even when it seems appropriate.

But it's also not that "When you need exact and predictable results, you can't use fp math". There are quite simple (understandable) precision guarantees for non-transcendental operations. Gaining insight on them won't do harm to a programmer.

I agree about the transcendental functions though.

2007-02-02 Reply Admin

First let me say that I agree that using floating point to do exact calculations is usually dubious, if not for any other reason, then for the reason that it will waste time for any person reading the code, trying to determine if it is in fact correct.

Ancient hacker wrote: "No, you still don't get it. There are an infinite number of numbers where the sqrt() of them are 10.9999999999999XXXXX. A computer using 32 or 64 or 80 bit math is going to get the wrong answer, even to three digits."

Oops, it seems you are right :-). (if not for sqrt(2) though, but for sqrt of other numbers).

But I'm not giving up that easily. Let's say that you have an embedded application, a high end video camera doing video compression. Now, let's say, that for some reason, we need to calculate the number of bits required to store a certain integer. If we notice that the floating point processor has an instruction for converting a 32-bit integer to a 64bit floating point number, that the exponent of the resulting floating point number is exactly the number of bits we need, and that this instruction is backed by dedicated hardware, do you think we should not use it? Even though it would make our product faster? Even though we have full control over the hardware we run on?

All I'm saying is that floating point for exact calculations is bad in 99.999% of all cases, but not in ALL cases.

Anyway, always nice with a little argument.

2007-02-02 Reply Admin

I just had to add:

Ancient hacker, you state:

"Even if you put every bit of memory on the planet into your computer, you still couldn't represent sqrt(2)."

No not as a floating point. Not as a fixed-point number. Not in any way?

Let's look at you're statement again:

"Even if you put every bit of memory on the planet into your computer, you still couldn't represent sqrt(2)."

1: It states that a certain entity can not be represented on a any computer. 2: It is itself stored in a computer. 3: It does, itself, represent the unrepresentable entity. Exactly.

2007-02-02 Reply Admin

replace you're with your.

2007-02-02 Reply Admin

I'm sure you're well-intentioned, but your example proves *nothing*. You only tested the powers of two, and on one computer. That proves *nothing*. To prove the correctness of this basic approach you'd have to try it on every computer in existence, in every IEEE and/or native rounding modes.

The real WTF here is that you're so fixated on floating point esoterica that you are blind to the fact that the posted algorithm calculates the right answer. You also don't appear to understand that the powers of 2 are the boundary conditions and if you get the right answer for those everything else must be correct by definition.

To repeat: you can't get exact output with inexact inputs.

Why do you continue to make statements that contradict the bare facts that have just been placed before you?

2007-02-02 Reply Admin

Floating point should never be used for exact things. It's an approximation.

While the resolution of a standard float or double is high, there are numerous operations you can do that will not only reduce resolution, but increase the error! Basically, if you use an approximation in a bunch of calculations, your approximation can get worse and worse and you can end up dropping bits.

Here's a document that's floating around the 'net - "What every computer scientist should know about floating point arithmetic". http://docs.sun.com/source/806-3568/ncg_goldberg.html

It describes in detail what happens in the floating point unit and how you can lose precision doing certain operations in certain ways. In other words, if you need the precision, you must exercise care in how you calculate to avoid inadvertently losing precision.

Ancient_Hacker · 2007-02-02 Reply Admin

"Even if you put every bit of memory on the planet into your computer, you still couldn't represent sqrt(2)."

1: It states that a certain entity can not be represented on a any computer. 2: It is itself stored in a computer. 3: It does, itself, represent the unrepresentable entity. Exactly.

There's a name for this, it's called "being disengenious".

You know full well that we're talking about representing sqrt(2) as a floating-point number. Yet you waste net bandwidth and clarity by quibbling about the obvious and irrelevant. It matters not one whit that we can write "sqrt(2)". What matters is whether we can represent its value as a floating-point number, and the answer is unequivocably NO.

Ancient_Hacker · 2007-02-02 Reply Admin

you are blind to the fact that the posted algorithm calculates the right answer

You're confusing algorithm and implementation.

The algorithm, by mathematical definition, calculates the right answer.

Any implementation on a real computer, might work, or it might not. And one test run for a very limited set of data on one computer, with one compiler, with one set of settings, proves nothing.

Ancient_Hacker · 2007-02-02 Reply Admin

If we notice that the floating point processor has an instruction for converting a 32-bit integer to a 64bit floating point number, that the exponent of the resulting floating point number is exactly the number of bits we need, and that this instruction is backed by dedicated hardware, do you think we should not use it?

Well unfortunately exponent fields are usually biased halfway up, so zero exponent is coded as +128. And there are still all those IBM mainframes (370,4xxx,Z) with 4-bit aligned mantissas, so the exponent goes up by jerks. But if you're sure the exponent is always correct, use it.

Regards,

A_H

2007-02-02 Reply Admin

Ancient_Hacker:
>you are blind to the fact that the posted algorithm *calculates the right answer*
You're confusing algorithm and implementation.

The algorithm, by mathematical definition, calculates the right answer.

Any implementation on a real computer, might work, or it might not. And one test run for a very limited set of data on one computer, with one compiler, with one set of settings, proves nothing.

What I find interesting is that a number of people have written actual code which they compiled and executed and found correct results. Ancient Hacker, however, proposes the chance of a possible failure, and yet has no tangible evidence. Can you, Ancient Hacker, find a compiler or architecture that does not come to the proper result with the code listed?

Spouting paranoid theories proves much, much less than these solid, tangible examples, which are backed by theories. How many test runs with what data sets one how many computers with how many compilers with how many different settings are needed, Ancient Hacker?

For the given situation, we know the entire data set: 1 through MAX_INT+1. (Since we are dealing with ranges of numbers, the smallest range would be where all the numbers are the same, or 1, and the largest would be where one number is 0 and the other is the largest number possible, MAX_INT, and the range would be MAX_INT-0+1 or MAX_INT+1.)

We also know that the vast majority of modern processors implement the IEEE floating point numbers.

Sure, there are still a variety of external variables not accounted for--so what, exactly, would you have us do to guarantee proper results?

Ancient_Hacker · 2007-02-02 Reply Admin

Okay, let's step back a bit.

One the one hand, we have a simple and quick algorithm that we KNOW works, can be proven to be correct, and will work correctly on each and every computer that has ever had a C compiler:

for( tot = 0; input > 0; input >>= 1, tot++ }

on the other hand, we have the log() algorithm, which cannot be proved to be correct, and is known to break with 32-bit floats. And can be mighty slow if you don't have log() on the chip.

Now which one should we choose? Hmmmm....

Ancient_Hacker · 2007-02-02 Reply Admin

How many test runs with what data sets one how many computers with how many compilers with how many different settings are needed, Ancient Hacker?

You're playing a losing game. There is no number of tests that can "prove" a shaky bit of code.

Not my proverb-- You might recall Dijkstra famous statement: that software testing cannot prove that there are no errors in your code; it can only demonstrate that they exist.

Proof-by-testing can never work. The counterexamples are legion.

How about we just use a reliable algorithm that's guaranteed to work?

2007-02-02 Reply Admin

You also don't appear to understand that the powers of 2 are the boundary conditions and if you get the right answer for those everything else must be correct by definition.

Um, I think the point A_H was trying to make was that when dealing with floating point values, what you think are boundary values aren't. And simple testing bears this out:

cat bits.c

#include <stdio.h>
#include <math.h>
#include <stdint.h>
 
int main()
{
   float num, den;
   int i, pow, res;
   res = -1;
   den = logf (2.0);
   for(i = 1; i < INT32_MAX; i++)
   {
      num = logf ((float)i);
      pow = (int)floorf (num/den);
      if (res != pow)
         printf("%d has %d bits.\n", i, pow+1);
      res = pow;
   }
   return 0;
}

gcc -W -Wall -o bits bits.c -lm ./bits 1 has 1 bits. 2 has 2 bits. 4 has 3 bits. 8 has 4 bits. 16 has 5 bits. 32 has 6 bits. 64 has 7 bits. 128 has 8 bits. 256 has 9 bits. 512 has 10 bits. 1024 has 11 bits. 2048 has 12 bits. 4096 has 13 bits. 8193 has 14 bits. 16384 has 15 bits. 32769 has 16 bits. 65536 has 17 bits. 131072 has 18 bits. 262144 has 19 bits. 524288 has 20 bits. 1048576 has 21 bits. 2097151 has 22 bits. 4194303 has 23 bits. 8388601 has 24 bits. 16777201 has 25 bits. 33554418 has 26 bits. 67108869 has 27 bits. 134217541 has 28 bits. 268435209 has 29 bits. 536870673 has 30 bits. 1073741889 has 31 bits. 2147480641 has 32 bits.

The boundary values for 14, 16, and 22-32 bits, are WRONG.

Results are similar if you use the ceil formula. Considering that IEEE floats have 23 bits of mantissa, it's not surprising that this implementation starts to be consistently off at that point. If you change to doubles with 52 bits of mantissa, it should be fine for 32-bit integers, but bit-fiddling to get the exact answer is still the faster, more correct way to go.

Why introduce error terms at all into a result that can be perfectly calculated with integer math?

2007-02-02 Reply Admin

No way to precisely represent sqrt(2)? How about representing it as sqrt(2)? Why do you have this impetus to turn everything into decimals?

2007-02-02 Reply Admin

As fun as it was reading three pages of posts discussing the finer points of floating point integers and IEEE standards, I felt kind of bad that no one has pointed out that tyhik has already posted the optimal solution to this problem that should have been obvious to anyone: a binary search.

My attempt is below: (note: uint = typedef unsigned int uint; here)

uint bitlen(uint val) {
uint step = 16, range = step;
  while(step >>= 1) range += (val >= 1 << range) ? step : -step;
  return 1 + range - (val < 1 << range);
}

... which resulted in roughly identical performance, but was notably more complex than tyhik's (superior) solution:

uint bitlen(uint val) {
uint step = 16, bits = 0;
  do if(val >> step + bits) { bits += step; } while(step >>= 1);
  return 1 + bits;
}

Now, compared to:

uint bitlen(uint val) {
  return 1 + (uint)(log(val) / log(2));
}

Either of our methods are anywhere from five to ten times faster than using logarithms, are easier to understand than natural logarithms (personal opinion), and are written in purely integeral math (not needing even multiplication or division so they will work on processors such as the 65c02), and could easily be extended to any bitsize your processor supports, or even to arrays of bits extending to any arbitrary number of bytes. Whereas you are always limited to ~50-60 bits with the double approach.

So this function is: faster, cleaner, safer and more extensible than any other solution posted here. tyhik is the winner of this thread :)

Sorry to interrupt, please feel free to continue discussing why ceil() should or should not ever be used for this specific purpose ;)

2007-02-02 Reply Admin

I apologize, before anyone points it out ...

Please use:

uint step = sizeof(uint) << 2

... to find the middle of the integer, as the size of an `unsigned int' is clearly not a given.

2007-02-03 Reply Admin

Could have been even worse:

switch (x) {
case 1:
  rbits = 1;
  break;
case 2:
case 3:
  rbits = 2;
  break;
case 4:
case 5:
case 6:
case 7:
  rbits = 3;
  break;
// ....

2007-02-03 Reply Admin

range >>= range;

This is machine dependend coding, depending wether the highest bit is on the left or right, should be checked in advance. Divide / 2 would make it independend.

2007-02-03 Reply Admin

I thought the real WTF was that the 54 separate instances of the code were not converted to a single function.

Copy once, shame on you. Copy 54 times, you're fired!

If it was a C-Fortran hybrid, it would have been really fun to define a C macro for the Fortran code and run the C preprocessor on the Fortran files. Hmm, sounds like a good job security move.

Reminds me of a student I once had in a C class... She did not realize that her primary programming language (a VB dialect) supported functions/subroutines. "I just like to to put blank lines and spaces between the different pieces of logic."

Paramedics were called shortly after the "introduction to pointers section."

cklam · 2007-02-03 Reply Admin

Ben Follis:
You could get the same effect by successively multiplying 2 until it is > the number and then counting the number of mults.

Hmmm ... if the code is nested inside many code loops you could run into performance issues. On the other hand, your solution would definitely work.

IMHO, the real WTF here is that the code is replicated so many times in the code project instead of having been isolated in a function.

cklam · 2007-02-03 Reply Admin

:
Assuming N is a positive integer, the number of bits is either
1 + floor( log(N) / log(2) )

or

ceil( log(N+1) / log(2) )

The former is better suited to bit-twiddling implementations though, so should be considered superior.

Why use floating point in the first place ? You convert from int to floating point back int int ending up dealing with all these nice floating point rounding adjustments.

And it is not really a mathematical problem in the first place - it is a problem concerning storage representation.

You do not have to use math for that (OK - Fortran cannot shift bits in the first place - although it will be able to - see Fortran-2008); just loop over the number and count bits ... that's all. May not be beautiful/elegant but it works.

And isolate the code in a function.

Conceded, there are not bit-shift operators possible

cklam · 2007-02-04 Reply Admin

Ancient_Hacker:
< .. Snip .. >
As mentioned in the original posting, this general method was used about 54 times in different places throughout the code, with slight variations-- I only posted one of the instances. Sometimes the range was calculated correctly, as being hi - low + 1, othertimes not. Sometimes there was a "+1" or even "+ 2" inside the argument to log() to perhaps fix that error, sometimes it was a fudge factor to get the answer "more correct". Quite often there were commented out "WRITE(6,100) ANSWER" lines, showing the code was problematical until the right fudge factors were introduced.

< .. Snip .. >

Whow. That means if the code it is refactored to use a single correct function then you have to check in 54 places for fudging in the code that follows the specific variation being replaced by the function.

They did not even use cut-and-paste - they re-invented the wheel 54 times. There is no excuse for that - that is just plain old stupid.

2007-02-04 Reply Admin

webrunner:
Remco Gerlich:
I hate people whining about floating point being "inexact." Even though it is technically true, a IEEE 754 float or double is the closest possible approximation in the given number of bits. On modern hardware, the error for a double is never greater than 1 part in 2**53.

"closest possible" means "inexact" when it's not possible to be exact.

Remember, in binary, no .1 for you!

What the hell sort of soup tastes like .1?

Twelve on a Line

Leave a comment on “Twelve on a Line”