440,649 Members | 2,142 Online
Need help? Post your question and get tips & solutions from a community of 440,649 IT Pros & Developers. It's quick & easy.

# Bit fiddling calculating fraction

 P: n/a Hello, I'm writing a function that should do the following: /** * Calculate and return fraction of valueA where max fractions is 31. * param valueA A five bit value, 0-31. * param valueB The fraction, a five bit value, 0-31. * */ ushort f ( ushort valueA, ushort valueB) { return valueA * (valueB / 31); //the paranthesis are only here for the looks of it. } I'm looking for a way to get approximatly the same result by doing some bit fiddling with the bit patterns of the values optimizing performance. Anyone got any ideas? Table lookup? Jan 5 '07 #1
22 Replies

 P: n/a

 P: n/a Sorry for being unclear. This is maybe better? ushort f ( ushort valueA, ushort valueB) { float fraction = valueB / 31.0; return (ushort)(valueA * fraction); } Only that I want to avoid using floating point operations and division/multiplication. Jan 5 '07 #3

 P: n/a do*********@gmail.com writes: I'm writing a function that should do the following: /** * Calculate and return fraction of valueA where max fractions is 31. * param valueA A five bit value, 0-31. * param valueB The fraction, a five bit value, 0-31. * */ ushort f ( ushort valueA, ushort valueB) { return valueA * (valueB / 31); //the paranthesis are only here for the looks of it. } C.l.c Guidelines suggest posting a complete program. We have to guess what ushort is. This is not as easy as you probably think since what you have written above makes no sense if we assume typedef unsigned short ushort; I'm looking for a way to get approximatly the same result by doing some bit fiddling with the bit patterns of the values optimizing performance. If we believe the block comment about the possible values in valueB, then return valueB == 31; will do the same as the return statement you wrote. If you put the parentheses in by mistake (thinking that it made no difference) then a divide is probably the way. Of course, it is more likely that you want to be dividing by 32. In that case a shift will be faster but the compiler will probably do that for you anyway. Write the simples clearest code and optimise it for speed only if you find it to be a problem. -- Ben. Jan 5 '07 #4

 P: n/a do*********@gmail.com a écrit : Hello, I'm writing a function that should do the following: /** * Calculate and return fraction of valueA where max fractions is 31. * param valueA A five bit value, 0-31. * param valueB The fraction, a five bit value, 0-31. * */ ushort f ( ushort valueA, ushort valueB) { return valueA * (valueB / 31); //the paranthesis are only here for the looks of it. } I'm looking for a way to get approximatly the same result by doing some bit fiddling with the bit patterns of the values optimizing performance. Anyone got any ideas? Table lookup? Any division by a signed constant can be replaced by a much cheaper multiplication and some shifts. For instance, for the code above, the lcc-win32 compiler will not generate a division by 31 but the following code: movzwl 8(%ebp),%edi ; valueA-->edi movzwl 12(%ebp),%eax ; valueB->eax movl \$-2078209981,%edx movl %eax,%ecx ; save valueB for later imull %edx ; multiply valueB * -2078209981 addl %ecx,%edx ; add valueA to high word sarl \$4,%edx ; shift right 4 sarl \$31,%ecx ; shift left 31 subl %ecx,%edx ; subtract the result to valueB movl %edx,%eax ; high word is result of division by 31 imull %eax,%edi ; multiply by valueA movzwl %di,%eax ; cast to unsigned short and set return value You see? No division. I would say you leave the optimization to the compiler. jacob Jan 5 '07 #5

 P: n/a 2; Those are the three basic techniques. Any look attractive to you? Dave. Jan 6 '07 #6

 P: n/a do*********@gmail.com wrote: > I'm writing a function that should do the following: /** * Calculate and return fraction of valueA where max fractions is 31. * param valueA A five bit value, 0-31. * param valueB The fraction, a five bit value, 0-31. * */ ushort f ( ushort valueA, ushort valueB) { return valueA * (valueB / 31); //the paranthesis are only here for the looks of it. } I'm looking for a way to get approximatly the same result by doing some bit fiddling with the bit patterns of the values optimizing performance. Anyone got any ideas? Table lookup? Try the following: /* Find best rational approximation to a double */ /* by C.B. Falconer, 2006-09-07. Released to public domain */ #include #include #include #include #include #include int main(int argc, char **argv) { int num, approx, bestnum, bestdenom; int lastnum = 500; double error, leasterr, value, criterion, tmp; char *eptr; value = 4 * atan(1.0); if (argc 2) lastnum = strtol(argv[2], NULL, 10); if (lastnum <= 0) lastnum = 500; if (argc 1) { tmp = strtod(argv[1], &eptr); if ((0.0 >= tmp) || (tmp INT_MAX) || (ERANGE == errno)) { puts("Invalid number, using PI"); } else value = tmp; } criterion = 2 * value * DBL_EPSILON; puts("Usage: ratvalue [number [maxnumerator]]\n" "number defaults to PI, maxnumerator to 500"); printf("Rational approximation to %.*f\n", DBL_DIG, value); for (leasterr = value, num = 1; num < lastnum; num++) { approx = (int)(num / value + 0.5); if (0 == (int)approx) continue; error = fabs((double)num / approx - value); if (error < leasterr) { bestnum = num; bestdenom = approx; leasterr = error; printf("%8d / %-8d = %.*f with error %.*f\n", bestnum, bestdenom, DBL_DIG, (double)bestnum / bestdenom, DBL_DIG, leasterr); if (leasterr <= criterion) break; } } return 0; } /* main, ratapprx */ -- Chuck F (cbfalconer at maineline dot net) Available for consulting/temporary embedded and systems. Jan 6 '07 #7

 P: n/a CBFalconer >I'm writing a function that should do the following:/** * Calculate and return fraction of valueA where max fractions is 31. * param valueA A five bit value, 0-31. * param valueB The fraction, a five bit value, 0-31. * */ushort f ( ushort valueA, ushort valueB) { return valueA * (valueB / 31); //the paranthesis are only here forthe looks of it.}I'm looking for a way to get approximatly the same result by doing somebit fiddling with the bit patterns of the values optimizingperformance.Anyone got any ideas? Table lookup? Try the following: /* Find best rational approximation to a double */ If this comment is an accurate description, then I can't see how it can help. The OP is starting with an exact, known, rational so getting close is not the issue. To the OP: Are you sure that simple integer arithmetic is not fast enough? I.e. return (valueA * valueB + 31/2) / 31; (and I m still having trouble with 31 being the correct divisor!) -- Ben. Jan 6 '07 #8

 P: n/a "CBFalconer" >I'm writing a function that should do the following:/** * Calculate and return fraction of valueA where max fractions is 31. * param valueA A five bit value, 0-31. * param valueB The fraction, a five bit value, 0-31. * */ushort f ( ushort valueA, ushort valueB) { return valueA * (valueB / 31); //the paranthesis are only here forthe looks of it.}I'm looking for a way to get approximatly the same result by doing somebit fiddling with the bit patterns of the values optimizingperformance.Anyone got any ideas? Table lookup? Try the following: /* Find best rational approximation to a double */ /* by C.B. Falconer, 2006-09-07. Released to public domain */ #include #include #include #include #include #include int main(int argc, char **argv) { int num, approx, bestnum, bestdenom; int lastnum = 500; double error, leasterr, value, criterion, tmp; char *eptr; value = 4 * atan(1.0); if (argc 2) lastnum = strtol(argv[2], NULL, 10); if (lastnum <= 0) lastnum = 500; if (argc 1) { tmp = strtod(argv[1], &eptr); if ((0.0 >= tmp) || (tmp INT_MAX) || (ERANGE == errno)) { puts("Invalid number, using PI"); } else value = tmp; } criterion = 2 * value * DBL_EPSILON; puts("Usage: ratvalue [number [maxnumerator]]\n" "number defaults to PI, maxnumerator to 500"); printf("Rational approximation to %.*f\n", DBL_DIG, value); for (leasterr = value, num = 1; num < lastnum; num++) { approx = (int)(num / value + 0.5); if (0 == (int)approx) continue; error = fabs((double)num / approx - value); if (error < leasterr) { bestnum = num; bestdenom = approx; leasterr = error; printf("%8d / %-8d = %.*f with error %.*f\n", bestnum, bestdenom, DBL_DIG, (double)bestnum / bestdenom, DBL_DIG, leasterr); if (leasterr <= criterion) break; } } return 0; } /* main, ratapprx */ I believe there is an O(log N) continued fraction algorithm to do this. Jan 6 '07 #9

 P: n/a "Ben Bacarisse" >>I'm looking for a way to get approximatly the same result by doing somebit fiddling with the bit patterns of the values optimizingperformance.Anyone got any ideas? Table lookup? If this comment is an accurate description, then I can't see how it can help. The OP is starting with an exact, known, rational so getting close is not the issue. What part of the OP's words "approximately the same result" didn't make it by your filter? Jan 6 '07 #10

 P: n/a "David T. Ashley" CBFalconer >>>I'm looking for a way to get approximatly the same result by doing somebit fiddling with the bit patterns of the values optimizingperformance.Anyone got any ideas? Table lookup? If this comment is an accurate description, then I can't see how itcan help. The OP is starting with an exact, known, rational so gettingclose is not the issue. What part of the OP's words "approximately the same result" didn't make it by your filter? Eh? What did I miss? The OP starts with a known rational valueA * valueB over 31 and CBF posts code to find a rational close to a given double. How does it help? The original rational even has a prime denominator so it can't even be reduced -- it's as good as he'll get (if a rational is what he wants). -- Ben. Jan 6 '07 #11

 P: n/a "Ben Bacarisse" Eh? What did I miss? The OP starts with a known rational valueA * valueB over 31 and CBF posts code to find a rational close to a given double. How does it help? The original rational even has a prime denominator so it can't even be reduced -- it's as good as he'll get (if a rational is what he wants). Oh, OK, I may have missed some of the nuances of the original post. In the larger scheme of things, just let me make the observation that these two problems are functionally equivalent. a)Finding a best rational approximation, with a certain maximum numerator and denominator, to a constant that is irrational. b)Finding a best rational approximation, with a certain maximum numerator and denominator, to a constant that is rational but that can't be expressed exactly without exceeding that certain maximum numerator and denominator. For example, the conversion constant from MPH to KPH is 1.609344 (I think this comes from the NIST or someone who should know). This constant is rational, but can't be expressed with numerator and denominator not exceeding 255 (for example). Using the continued fraction technique to find the best approximation with numerator and denominator not exceeding 255 (output pasted in at end) yields 103/64 as the best approximation. 103/64 is 1.609375 (pretty close). (By the way, it is sheer coincidence that the denominator came out to be a power of 2). Chuck's code has merit and relevance precisely in this situation. With Chuck's program, one would input the rational constant (1.609344) and get an approximation that could be expressed with smaller numerator and denominator. The "double" in Chuck's program may be in fact rational but just not expressible under the constraints on the numerator and denominator. My observation was that the math driving Chuck's program seemed to be ad hoc and probably O(N). That won't work well for approximations in, say, 256 bits. But the idea is sound. "Rational" and "rational under computational constraints" are two different things. Dave. ---------- C:\Documents and Settings\dashley>cfbrapab 1.609344 255 255 ------------------------------------------------------------------------------ MAJOR MODE: Finding closest rational number(s) under the constraints. ------------------------------------------------------------------------------ RI_IN Numerator: 25,146 ( 5 digits) ------------------------------------------------------------------------------ RI_IN Denominator: 15,625 ( 5 digits) ------------------------------------------------------------------------------ K_MAX: 255 ( 3 digits) ------------------------------------------------------------------------------ H_MAX: 255 ( 3 digits) ------------------------------------------------------------------------------ approx_num(1): 103 ( 3 digits) ------------------------------------------------------------------------------ approx_den(1): 64 ( 2 digits) ------------------------------------------------------------------------------ dap_num(1): 1, ( 109 digits) 609,375,000,000,000,000,000,000,000, 000,000,000,000,000,000,000,000,000, 000,000,000,000,000,000,000,000,000, 000,000,000,000,000,000,000,000,000 ------------------------------------------------------------------------------ dap_den(1): 1, ( 109 digits) 000,000,000,000,000,000,000,000,000, 000,000,000,000,000,000,000,000,000, 000,000,000,000,000,000,000,000,000, 000,000,000,000,000,000,000,000,000 ------------------------------------------------------------------------------ error_num(1): 31 ( 2 digits) ------------------------------------------------------------------------------ error_den(1): 1,000,000 ( 7 digits) ------------------------------------------------------------------------------ Jan 6 '07 #12

 P: n/a "David T. Ashley" >Eh? What did I miss? The OP starts with a known rationalvalueA * valueB over 31 and CBF posts code to find a rational closeto a given double. How does it help? The original rational even hasa prime denominator so it can't even be reduced -- it's as good ashe'll get (if a rational is what he wants). Oh, OK, I may have missed some of the nuances of the original post. The OP was rather short on nuance, I thought. I was simply defending my honour -- you seemed to imply I could not read properly! -- Ben. Jan 7 '07 #13

 P: n/a do*********@gmail.com wrote: I'm writing a function that should do the following: /** * Calculate and return fraction of valueA where max fractions is 31. * param valueA A five bit value, 0-31. * param valueB The fraction, a five bit value, 0-31. * */ ushort f ( ushort valueA, ushort valueB) { return valueA * (valueB / 31); //the paranthesis are only here for the looks of it. } I'm looking for a way to get approximatly the same result by doing some bit fiddling with the bit patterns of the values optimizing performance. Anyone got any ideas? Table lookup? Table look up is correct. In the comments you are suggesting that valueA and valueB are basically 5 bit numbers. If this is the case, then the table look up is at most 1024 entries. Furthermore, the table results can be stored in a single char. So clearly, that's the fastest possible answer. (If you think a 1K table is not overly burdensome of your resources.) Ok, if table look up is not for you, but you still have those range restrictions, and you just want an approximation then you can do the following: return (34 * valueA * valueB) >10; Hopefully your compiler can work out that 34 = 2*(16+1) which can be done all with some shifts and two shifts and an add. In such a small range, the above calculations use no more than 15 bits for any intermediate, which means that it will work portably on any machine with a correct C compiler implementation. If you are interested in the whole range on, say, 32 bits, I can offer you the following x86 assembly code: ;/* Input uint32_t valueA: eax, uint32_t valueB: ebx */ mul ebx ;/* eax = A * B */ cmp eax, 08d3dcb1bh ;/* "off by 1" threshold */ adc eax, 0ffffffffh ;/* compensate for off by 1 */ mov edx, 084210843h ;/* K = ceil(2^36 / 31) */ mul edx ;/* edx = A * B * K / 2^32 */ shr edx, 4 ;/* edx = A * B * (1/31) */ ;/* Output value: edx (eax destroyed) */ There is no way to re-express the above in pure C code, except to just stick with your original code, and hope the compiler is good enough to generate this code (not outside the realm of possibility, I am pretty sure gcc can do this, or something similar.) -- Paul Hsieh http://www.pobox.com/~qed/ http://bstring.sf.net/ Jan 7 '07 #14

 P: n/a "David T. Ashley" wrote: "CBFalconer" >/* Find best rational approximation to a double *//* by C.B. Falconer, 2006-09-07. Released to public domain */#include #include #include #include #include #include int main(int argc, char **argv){ int num, approx, bestnum, bestdenom; int lastnum = 500; double error, leasterr, value, criterion, tmp; char *eptr; value = 4 * atan(1.0); if (argc 2) lastnum = strtol(argv[2], NULL, 10); if (lastnum <= 0) lastnum = 500; if (argc 1) { tmp = strtod(argv[1], &eptr); if ((0.0 >= tmp) || (tmp INT_MAX) || (ERANGE == errno)) { puts("Invalid number, using PI"); } else value = tmp; } criterion = 2 * value * DBL_EPSILON; puts("Usage: ratvalue [number [maxnumerator]]\n" "number defaults to PI, maxnumerator to 500"); printf("Rational approximation to %.*f\n", DBL_DIG, value); for (leasterr = value, num = 1; num < lastnum; num++) { approx = (int)(num / value + 0.5); if (0 == (int)approx) continue; error = fabs((double)num / approx - value); if (error < leasterr) { bestnum = num; bestdenom = approx; leasterr = error; printf("%8d / %-8d = %.*f with error %.*f\n", bestnum, bestdenom, DBL_DIG, (double)bestnum / bestdenom, DBL_DIG, leasterr); if (leasterr <= criterion) break; } } return 0;} /* main, ratapprx */ I believe there is an O(log N) continued fraction algorithm to do this. I can believe that. I don't know it. I wrote the above to find the PI approximations, and generalized it. I didn't try to patent it and had no reason to try for speed-ups. -- Chuck F (cbfalconer at maineline dot net) Available for consulting/temporary embedded and systems. Jan 7 '07 #15

 P: n/a "CBFalconer" >I believe there is an O(log N) continued fraction algorithm to do this. I can believe that. I don't know it. I wrote the above to find the PI approximations, and generalized it. I didn't try to patent it and had no reason to try for speed-ups. Would this help you believe it more? http://www.dtashley.com/howtos/2007/...approximation/ : ) Dave Ashley Jan 7 '07 #16

 P: n/a Here is one I used to write a C++ version a long time ago: http://www.ics.uci.edu/~eppstein/numth/frap.c Jan 9 '07 #18

 P: n/a "user923005"

 P: n/a In article

 P: n/a In article -- dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131 home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/ Jan 9 '07 #21

 P: n/a "Dik T. Winter"

 P: n/a "David T. Ashley" In article

### This discussion thread is closed

Replies have been disabled for this discussion.