Author |
Topic: Rational search. (Read 1366 times) |
|
Grimbal
wu::riddles Moderator
Uberpuzzler
Gender: 
Posts: 7527
|
 |
Rational search.
« on: May 26th, 2004, 6:30am » |
 Quote  Modify
|
You are given interval [a,b], a<b real.
Design a fast algorithm to find the simplest rational p/q in this interval.
"Simplest" means that any other p'/q' in [a,b] has p'>=p and q'>=q.
update:
Assume a, p, q>0
or, for the general case,
consider as simplest rational p/q such that |p'|>=|p| and q'>=q>0
|
| « Last Edit: May 29th, 2004, 10:32am by Grimbal » |
 IP Logged |
|
|
|
towr
wu::riddles Moderator
Uberpuzzler
Some people are average, some are just mean.
Gender:
Posts: 13730
|
|
Re: Rational search.
« Reply #1 on: May 26th, 2004, 2:23pm » |
Quote Modify
|
This won't be fast, but it's a start.. maybe..
q=0;
do
q++;
p = ceil(a*q);
until p/q < b
Maybe a 'binary search'-like approach might work.. (Of course we could f.i. encounter an interval that has 1/3 in it, but not p/4 or p/2 for some integer p, so there are some inherent problems)
Anyway, at least you know this problem isn't being ignored 
|
|
IP Logged |
Wikipedia, Google, Mathworld, Integer sequence DB
|
|
|
Eigenray
wu::riddles Moderator
Uberpuzzler
Gender:
Posts: 1948
|
|
Re: Rational search.
« Reply #2 on: May 29th, 2004, 12:28am » |
Quote Modify
|
(You should probably either add some absolute values, or require a,p,q > 0, or something.)
My first thought was to try something like this:
h1=0; h2=1;
k1=1; k2=0;
s=(int)r=(a+b)/2;
while(a*k2 > h2 || b*k2 < h2) {
t=h2; h2=s*h2+h1; h1=t;
t=k2; k2=s*k2+k1; k1=t;
s=(int)r=1/(r-s);
}
But that is not guaranteed to work. It is true, though, that if p/q is the simplest rational in [a,b], then we must have k2 [ge] q > k1, so maybe use towr's code from there.
|
|
IP Logged |
|
|
|
Barukh
Uberpuzzler
Gender:
Posts: 2276
|
|
Re: Rational search.
« Reply #3 on: May 29th, 2004, 10:20am » |
Quote Modify
|
Excellent problem, grimbal!
I would like to understand what Eigenray's solution is...
Idea: why don't use continued fractions? 
|
|
IP Logged |
|
|
|
Grimbal
wu::riddles Moderator
Uberpuzzler
Gender:
Posts: 7527
|
|
Re: Rational search.
« Reply #4 on: May 29th, 2004, 10:25am » |
Quote Modify
|
on May 29th, 2004, 12:28am, Eigenray wrote:| (You should probably either add some absolute values, or require a,p,q > 0, or something.) |
|
Yes. assume a, p, q > 0.
If you can solve that, it easy to solve the general case with any a<b with p/q having the smallest |p| and smallest q>0.
|
| « Last Edit: May 29th, 2004, 10:35am by Grimbal » |
IP Logged |
|
|
|
Grimbal
wu::riddles Moderator
Uberpuzzler
Gender:
Posts: 7527
|
|
Re: Rational search.
« Reply #5 on: May 29th, 2004, 12:32pm » |
Quote Modify
|
on May 29th, 2004, 12:28am, Eigenray wrote:
h1=0; h2=1;
k1=1; k2=0;
s=(int)r=(a+b)/2;
while(a*k2 > h2 || b*k2 < h2) {
t=h2; h2=s*h2+h1; h1=t;
t=k2; k2=s*k2+k1; k1=t;
s=(int)r=1/(r-s);
}
|
|
Excellent approach! In fact, it uses the trick I use in my solution, but with an additional optimization.
But as you noted, it is not doing what is wanted. Your code is approximating (a+b)/2 with a tolerance of (b-a)/2. If you have [a,b] = [10,20], it will return 15/1. 
But you are on the right track.
Continuous fractions are probably a good idea also.
|
|
IP Logged |
|
|
|
towr
wu::riddles Moderator
Uberpuzzler
Some people are average, some are just mean.
Gender:
Posts: 13730
|
|
Re: Rational search.
« Reply #6 on: May 30th, 2004, 7:33am » |
Quote Modify
|
Once you have q it's easy.. It's also easy to find an upper bound for q, namely [lceil]1/(b-a)[rceil], as there must be some p/q in any interval longer than 1/q.
|
|
IP Logged |
Wikipedia, Google, Mathworld, Integer sequence DB
|
|
|
Grimbal
wu::riddles Moderator
Uberpuzzler
Gender:
Posts: 7527
|
|
Re: Rational search.
« Reply #7 on: May 31st, 2004, 5:02pm » |
Quote Modify
|
Hint: you can find a rational between 2 other rationals p/q and r/s by computing (p+q)/(r+s)
Now, if you start from 0 = 0/1 and oo = 1/0, you can find 0/1 < 1/1 < 1/0.
Then, 0/1 < 1/2 < 1/1 and 1/1 < 2/1 < 1/0.
you can repeat the process as long as you want.
0/1
1/4
1/3
2/5
1/2
3/5
2/3
3/4
1/1
4/3
3/2
5/3
2/1
5/2
3/1
4/1
1/0
Etc. And Guess what? This way, you generate all rationals >0.
|
| « Last Edit: May 31st, 2004, 5:57pm by Grimbal » |
IP Logged |
|
|
|
Eigenray
wu::riddles Moderator
Uberpuzzler
Gender:
Posts: 1948
|
|
Re: Rational search.
« Reply #8 on: May 31st, 2004, 9:19pm » |
Quote Modify
|
on May 29th, 2004, 10:20am, Barukh wrote:I would like to understand what Eigenray's solution is...
Idea: why don't use [hidden]? |
|
The thing is, you're just guaranteed that nothing simpler than a given convergent has a smaller error. You're not guaranteed that the simplest fraction within a given error is a convergent, though. So if a convergent h/k has too large an error, you know that anything with a denominator [le] k does also. Unfortunately, denominators tend to increase exponentially I think, so only having to check from k+1 probably wouldn't speed things up by more than a multiplicative constant.
I guess a binary search on the Stern-Brocot Tree is the way to go here.
|
|
IP Logged |
|
|
|
Barukh
Uberpuzzler
Gender:
Posts: 2276
|
|
Re: Rational search.
« Reply #9 on: Jun 1st, 2004, 12:00am » |
Quote Modify
|
It becomes more and more interesting! 
Here’s an alogorithm based on the idea of continued fractions (in fact, it’s a code that compiles).
Code:#define FZERO(r) (r < 1.0e-15)
RationalSearch(double a, double b)
{
unsigned s[100];
int n;
unsigned a0, b0;
unsigned p, q;
/* Matching continued fractions */
for (n = 0, a = 1.0/a, b = 1.0/b; ; n++)
{
a0 = (unsigned)floor(1.0/a);
b0 = (unsigned)floor(1.0/b);
if (a0 != b0) break;
s[n] = a0;
a = 1.0/a - a0;
b = 1.0/b - b0;
if (FZERO(a) || FZERO(b)) break;
}
if (a0 != b0) s[n] = (a0 < b0) ? a0+1 : b0+1;
/* Rational point */
for (q = 1, p = s[n--]; n >= 0; n--)
{
unsigned t = p;
p = p * s[n] + q;
q = t;
}
printf("%d/%d\n", p, q);
} |
|
In words: Simple continued fractions (CF) are built for a and b simultaneously – until they don’t match, or a fraction is completed for one of the number. In the former case, 1 is added to the smallest number.
The argument goes as follows:
1. It is known that the convergents of CF give the best rational approximations with such a small denominator.
2. Also, it is clear that while the convergents for a and b match, there cannot be a rational with the smaller denominator in the range.
3. If a CF is completed before first non-match, then the corresponding end-point is the “best” rational itself.
4. When the first non-match occurs, I use the fact that the sequence of (finite) CFs is monotonic if only the last number changes.
Comments will be greatly appreciated...
|
| « Last Edit: Jun 2nd, 2004, 9:54am by Barukh » |
IP Logged |
|
|
|
Grimbal
wu::riddles Moderator
Uberpuzzler
Gender:
Posts: 7527
|
|
Re: Rational search.
« Reply #10 on: Jun 1st, 2004, 9:00am » |
Quote Modify
|
on May 31st, 2004, 9:19pm, Eigenray wrote:
The Stern-Brocot Tree! I knew the tree, but I never knew the name. I knew it should have a name, though.
|
|
IP Logged |
|
|
|
Barukh
Uberpuzzler
Gender:
Posts: 2276
|
|
Re: Rational search.
« Reply #11 on: Jun 3rd, 2004, 11:06pm » |
Quote Modify
|
The following entry in the PlanetMath Encyclopedia shows the close relation between Stern-Brocot Trees (SBT) and CF, so the two approaches are to produce the same result!
Specifically, if a rational number p/q has a CF-representation in the form [a0, a1, …, an], then its path in the SBT is as follows (using Grimbal’s representation):
Start at 1/1
Go down a0 times
Go up a1 times
…
Go up/down (depending on the parity of n) an-1 times.
If I am not missing something, it is easy now to estimate the efficiency of different algorithms:
Towr’s : q iterations
SBT : [sum] ak - 1 iterations
CF : n iterations
Let me give an example to clarify: let a = 3/302, b = 5/502. The simplest approximation p/q = 2/201 = [0, 100, 2]. Then, towr’s algorithm will take 201 steps, SBT – 101 steps, and CF – 2 steps.
|
| « Last Edit: Jun 3rd, 2004, 11:08pm by Barukh » |
IP Logged |
|
|
|
Grimbal
wu::riddles Moderator
Uberpuzzler
Gender:
Posts: 7527
|
|
Re: Rational search.
« Reply #12 on: Jun 4th, 2004, 5:53am » |
Quote Modify
|
Yes, but then, you could optimize the search. Instead of just going up or down one step with each iteration, we can compute how many steps you can go in one direction. This comes close to what you do to compute continued fractions, but you don't need to track back afterwards to compute the actual fraction.
You have to go up as many steps as necessary to exceed the minimum a, at which time you are either within the interval or above b. If not, you go down as many times as necessary to get below b, and so on.
How many times to go up? Given p/q and r/s as upper and lower limit of the search, you want n s.t. (p+nr)/(q+ns) >= a. This translates to n >= (aq-p)/(r-as), or n = ceil((aq-p)/(r-as)). If you want to keep p/q as lower bound, you have to stop one step before.
Applying all this gives
Code:
#include <stdio.h>
#include <math.h>
main(){
double a = 3.0/302.0;
double b = 5.0/502.0;
int p, q, r, s;
int n, iter;
printf("a = %0.10lf\n", a);
printf("b = %0.10lf\n", b);
if( a<=0 || b<=a ) exit(1);
iter = 0;
p = ceil(a) - 1; q=1;
r=1; s=0;
printf("%d: %d/%d = %0.10lf\n", iter++, p+r, q+s, (double)(p+r)/(q+s));
while( 1 ){
n = ceil((r-b*s)/(b*q-p)) - 1;
if( n<=0 ) break;
r += n*p;
s += n*q;
printf("%d: %d/%d = %0.10lf\n", iter++, p+r, q+s, (double)(p+r)/(q+s));
n = ceil((a*q-p)/(r-a*s)) - 1;
if( n<=0 ) break;
p += n*r;
q += n*s;
printf("%d: %d/%d = %0.10lf\n", iter++, p+r, q+s, (double)(p+r)/(q+s));
}
printf("result: %d/%d = %0.10lf\n", p+r, q+s, (double)(p+r)/(q+s));
}
|
|
output:
a = 0.0099337748
b = 0.0099601594
0: 1/1 = 1.0000000000
1: 1/101 = 0.0099009901
2: 2/201 = 0.0099502488
result: 2/201 = 0.0099502488
This gives 3 iterations, which includes the initial iteration of 0 steps up when computing the starting point.
|
| « Last Edit: Jun 4th, 2004, 5:55am by Grimbal » |
IP Logged |
|
|
|
|
RIDDLES SITE
WRITE MATH!
Home 
Help 
Search 
Members 
Login 
Register
Reply
Notify of replies
Send Topic
Print