Bit fiddling calculating fraction

[dos.fishing]

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?

[Ben Bacarisse]

"David T. Ashley" <dt*@e3ft.comwr ites:

"Ben Bacarisse" <be********@bsb .me.ukwrote in message
news:87******** ****@bsb.me.uk. ..

>CBFalconer <cb********@yah oo.comwrites:

>>>>
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
performanc e.
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 "approximat ely 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.

[David T. Ashley]

"Ben Bacarisse" <be********@bsb .me.ukwrote in message
news:87******** ****@bsb.me.uk. ..

>
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\dashle y>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,00 0,000,
000,000,000,000 ,000,000,000,00 0,000,
000,000,000,000 ,000,000,000,00 0,000,
000,000,000,000 ,000,000,000,00 0,000
------------------------------------------------------------------------------
dap_den(1): 1, ( 109
digits)
000,000,000,000 ,000,000,000,00 0,000,
000,000,000,000 ,000,000,000,00 0,000,
000,000,000,000 ,000,000,000,00 0,000,
000,000,000,000 ,000,000,000,00 0,000
------------------------------------------------------------------------------
error_num(1): 31 ( 2
digits)
------------------------------------------------------------------------------
error_den(1): 1,000,000 ( 7
digits)
------------------------------------------------------------------------------

[Ben Bacarisse]

"David T. Ashley" <dt*@e3ft.comwr ites:

"Ben Bacarisse" <be********@bsb .me.ukwrote in message
news:87******** ****@bsb.me.uk. ..

>>
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.

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.

[websnarf]

do*********@gma il.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/

[CBFalconer]

"David T. Ashley" wrote:

"CBFalconer " <cb********@yah oo.comwrote in message

.... snip ...

>>
/* Find best rational approximation to a double */
/* by C.B. Falconer, 2006-09-07. Released to public domain */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <errno.h>

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("Rationa l 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)nu m / 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.
<http://cbfalconer.home .att.net>

[David T. Ashley]

"CBFalconer " <cb********@yah oo.comwrote in message
news:45******** *******@yahoo.c om...

>>
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

[dos.fishing]

Thanks everybody for your constructive information! =)

Jacob,
Really amazing how much compilers can optimize these days. Only wish I
knew more x86 assembler to fully understand your code snippet. Why is
it doing "multiply valueB * -2078209981"?
Your name rang I bell, I didn't connect at first. Really cool that one
of the creators of lcc-win32 is answering my question! =D

Dave,
Your 3 basic strategies looks usefull to me, except that I don't really
understand #2.
I can use # 3 to multiply my valueA with valueB and then apply #1 by
using 32 as denominator instead of 31 to do the division.

Chuck F
Interesting code snippet, but I'm more interested in performance than
precision.

Ben,

>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!)

You're right about the divisor. It's 32. I tried 31 only because i
thought it might simplyfy the operation since the values are 5bit
values.

return (valueA * valueB + 16) / 32
The row above looks like a good candidate together with Dave's methods.
Just curious do you guys think the compiler will optimize the division
away for me (and maybe do something about the multiplication also)?
Jacob?

Paul,
Yes, I'll probably use table lookup in one implementation to see if
it's worth the memory.

return (34 * valueA * valueB) >10;

I'm not so good at thinking in binary, could you please explain what
this expression does? How does the multiply by 34 and right shift by 10
correspond to division by 31 (or 32)?
About your assembler code snippet, still trying to figure out how it
works. ;o)

By the way, anyone got a good x86 assembler instruction howto site to
recommend?

Thanks all for your help! =)
//Kid

[user923005]

Here is one I used to write a C++ version a long time ago:
http://www.ics.uci.edu/~eppstein/numth/frap.c

[David T. Ashley]

"user923005 " <dc*****@connx. comwrote in message
news:11******** *************@s 34g2000cwa.goog legroups.com...

Here is one I used to write a C++ version a long time ago:
http://www.ics.uci.edu/~eppstein/numth/frap.c

This is correct. Dave Eppstein discovered or rediscovered the continued
fraction results about 7 years before I did. (I didn't look at the program
closely, but I assume this is how it works.)

My path was that I read Khinchin's book, and I believe Khinchin states that
the best approximation is either the convergent or some other type of
fraction (forget what it is called, maybe "intermedia te fraction"). Then I
looked at it closely ... and I realized that you could make a stronger
statement. You can say specifically that the number of interest is
bracketed on the left and right by the convergent and a specific one of the
intermediate fractions (but either can be on the left). This appears in the
documents I cited:

http://www.dtashley.com/howtos/2007/...approximation/

Khinchin was of course brilliant. I can only assume that he didn't make the
strongest statement possible because he wasn't that interested in the
problem.

Dave Eppstein got there 7 years before me, and I discovered what he had done
after I did what I did. I do suspect, of course, that neither of us was the
first.

This is relatively minor stuff for a number theorist, but major stuff for an
an embedded systems developer (me). It took me about half a week to figure
out what would have probably taken a mathematician 30 minutes.

[Dik T. Winter]

In article <SZ************ *************** ***@giganews.co m"David T. Ashley" <dt*@e3ft.comwr ites:
....

My path was that I read Khinchin's book, and I believe Khinchin states that
the best approximation is either the convergent or some other type of
fraction (forget what it is called, maybe "intermedia te fraction").

That is one of the possible terms. The classic book on continued fractions
is by O. Perron, but I do not know whether it has been translated in English.

I
looked at it closely ... and I realized that you could make a stronger
statement. You can say specifically that the number of interest is
bracketed on the left and right by the convergent and a specific one of the
intermediate fractions (but either can be on the left). This appears in the
documents I cited:

Given successive convergents h[n-2]/k[n-2], h[n-1]/k[n-1] and h[n]/k[n],
We have h[n] = a[n] * h[n-1] + h[n-2] and similar for k[n], where a[n]
is the n-th term in the expression of the continued fraction. When
we use integers in the range a[n]/2 .. a[n] rather than a[n] itself,
we get the fractions that are also "best". (Note: when a[n] is even
there are some rules that state whether a[n]/2 also can be used.)
What you can state is that h[n-1]/k[n-1] is on one side while the
next intermediate fractions and the next convergent are on the other side.

Khinchin was of course brilliant. I can only assume that he didn't make the
strongest statement possible because he wasn't that interested in the
problem.

I would think that also Khinchin will have stated this in his book.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, [phone]
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/