Cosine redux
By jag on Jan 20, 2006
It's been quite a while since my last blog entry: I had a bit of a technology meltdown that I'm not quite done with yet :-(
I was doing some recreational hacking over the holidays that involved evaluating cosines. I ended up doing (once again!) an implementation of cosine (don't ask why). Given the confused flamage that my previous posts on cosine generated, I figure that showing some code would be useful. I would never have thought that cosine would generate so much email traffic. Yes, I know about Taylor series. No I wasn't trying to insult centuries of standard mathematical practice. Performance is often the art of cheating carefully.
So, here's my implementation of cosine, with its argument in turns (1 turn == 360 degrees):
public static float cost(float f) {
int bits = Float.floatToRawIntBits(f);
int mantissa = (bits&0x7FFFFF)|(1<<23);
int shift = (bits<<1>>>24)-127+9; // exponent, unbiased, with shift
if(shift>=32 || shift<=-32) return 1;
int fractionBits = shift>=0 ? mantissa<<shift : mantissa>>-shift;
int tableIndex = (fractionBits>>>(30-resultPrecision))&tableSizeMask;
switch(fractionBits>>>30) { // Quadrant is top two bits
case 0: return cosTable[tableIndex];
case 1: return -cosTable[tableSizeMask-tableIndex];
case 2: return -cosTable[tableIndex];
default/\*case 3\*/: return cosTable[tableSizeMask-tableIndex];
}
}
Lets go through this slowly:
-
int bits = Float.floatToRawIntBits(f);
Get the IEEE 754 bits - int mantissa = (bits&0x7FFFFF)|(1<<23);
The mantissa is the bottom 23 bits - to which the hidden bit must be prepended. - int shift = (bits<<1>>>24)-127+9;
Extract the exponent, correct for the exponent bias, then add a bias to move the binary point to the top of the word. - if(shift>=32 || shift<=-32) return 1;
If the shift is too large, the fraction bits would be zero, therefore the result is 1. - int fractionBits = shift>=0 ? mantissa<<shirt : mantissa>>-shift;
Shift the mantissa so that it's a fixed point number with the binary point at the top of the 32 bit int. The magic is in what's not here: because the argument is in turns, I get to ignore all of the integer bits (range reduction made trivial); and because it's the cosine function, which is symmetric about the origin, I get to ignore the sign. - int tableIndex = (fractionBits>>>(30-resultPrecision))&tableSizeMask;
The top two bits are the quadrant... extract the bits below that to derive a table index. - switch(fractionBits>>>30) {
One case for each quadrant. This switch could be eliminated by making the table 4 times larger. - case 0: return cosTable[tableIndex];
Yes! It's just a table lookup! Truly trivial. - case 1: return -cosTable[tableSizeMask-tableIndex];
...
Since this is just a table
lookup, the resulting approximation can be pretty jagged if the table is small. But it's easy to tune the
table size depending on accuracy needs. The smaller the table is, the higher the cache hit rate will
be, and the more likely it is that the whole table will fit in cache.
Table lookups are a very common way to implement mathematical functions, particularly the periodic ones like cosine. There are all kinds of elaborations. One of the most common for improving accuracy is to do some sort of interpolation between table elements (linear or cubic, usually).
Having taken a chunk of time off in August to hang out with family
(including one very active 90-year-old aunt's birthday party) I'm back at
work. This week I've been visiting customers and participating in
developer events. I spent Monday in Zurich visiting developers. The rest
of the week I've spent in St Petersburg (not Florida). We did a big
developer event on Thursday and ended up with hugely more folks attending
than we had expected. We had people sitting in the aisles in the theatre.
Peter Moore and his mates from 
Lots and lots of cool stuff... Not nearly enough time to write about it
all. The thrills started even before the keynote: the jazz group that
provided the opening music including one of my all time favorite
musicians: 
I was just looking through the 
Since I've been involved in real time programming and the Real Time
Specification for Java (
I can publish my blog to two sites at the same time! Both blogs.sun.com
and today.java.net. Bonus!
I need to stop talking to reporters. It's so easy for what results to get
misunderstood. I was not trying to say that all open source
projects are chaotic: there is a spectrum. Apache is at the very high end
of the scale, on average exhibiting excellent behaviour. Their governence
rules are very effective. Apache and the Linux kernel set the Gold
Standard. But they aren't all of the open source universe, and there's
some decidedly oddball behaviour that goes on. The problem is that it's
often the crazy behavior that becomes publicly visible and it tarnishes
everything. When I made the comment that got so misunderstood I was
talking about the perception of the open source community by
outsiders.
