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?

 

[Barry]

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?

It isn't clear what you are doing. Based on the comments this
returns A or zero. Also, the parantheses do make a difference.

 

[dos.fishing]

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.

 

[Ben Bacarisse]

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.

 

[jacob navia]

(e-mail address removed) 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

 

[David T. Ashley]

Hello,

I'm writing a function that should do the following:

/**
* Calculate and return fraction of valueA where max fractions is 31.

Etc.

If I'm understanding your problem correctly, it is the standard rational
approximation problem with the denominator fixed.

Generally, you seek to approximate f(x) = r * x by g(x) = (h * x) div k;
where r and x are real, and h and k are integers.

You generally seek to choose h/k as close as possible to r.

There are three basic strategies.

#1: You can choose k as a power of 2 so that the mapping is

g(x) = (h * x)/2^q

where of course the division is handled by shifting.

#2: You can use an arbitrary h and k, which normally doesn't work out badly
since most processors have efficient integer multiplication and division
instructions.

Choosing h/k as close as possible to an arbitrary r isn't as easy as it
sounds. There is an O(log N) technique using number theory and continued
fractions. I'll be glad to send you a technical paper if you e-mail me
directly at (e-mail address removed). For example, the best approximation to PI with 16
bit h,k is 65298 / 20785. Finding such an approximation in 32 or 64 bits is
very hard if you don't know how to do it O(log N).

#3: You can sometimes use addition to implement the multiplication. For
example, if h is 11 and k is 4:

two_times = x + x;
four_times = two_times + two_times;
eight_times = four_times + four_times;
eleven_times = x + two_times + eight_times;
result = eleven_times >> 2;

Those are the three basic techniques.

Any look attractive to you?

Dave.

 

[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 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 <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("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 */

 

[Ben Bacarisse]

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

 

[David T. Ashley]

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

 

[David T. Ashley]

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?

 

[Ben Bacarisse]

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

 

[David T. Ashley]

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

 

[Ben Bacarisse]

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!

 

[websnarf]

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

 

[CBFalconer]

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

 

[David T. Ashley]

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/01/best_rational_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]

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 "intermediate 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/01/best_rational_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]

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