Subj : Re: ratio approximation algorithm
To   : comp.programming
From : CBFalconer
Date : Thu Aug 18 2005 11:32 pm

Roger Willcocks wrote:
> "Mark Maglana" <mmaglana@gmail.com> wrote in message
>> Roger Willcocks wrote:
>
>>> The trick is to notice that the numerator (a.c) and the
>>> denominator (b.d) can only take on a small number of distinct
>>> values (for instance c >= a), so one can precompute a 6105-entry
>>> list of products, and find the minimum err = bdx - ac in 6105^2
>>> comparisons
>>
>> That one suggestion alone improved my search times by leaps and
>> bounds. Why the heck didn't I think of that earlier??? :-)
>>
.... snip code ...
>>
>> What I'm actually trying to do here is find as many combinations
>> within a given deadline (usually 15 seconds). With the above
>> method, what used to be a search result of 2 or 3 items has jumped
>> to 200+!
>>
>> Anyway, I guess the above method still makes 11,325^2 comparisons.
>> I'm thinking about using binary search for the innermost loop to
>> further improve the time. But I'm still trying to figure out how
>> to modify it such that it looks for a range rather than a
>> particular number. Am I making any sense or do I sound like a
<< raving idiot at the moment?
> 
> You only need to check the pairs where a >=c and b >= d. That
> should speed things up a bit.
.... snip ...

The numerator and denominator are each of the form a*b, where one
may be unity.  That reduces to finding a rational approximation to
whatever you are looking for, and then possibly factoring the
numerator and denominator.  I published code here to do just that a
while ago, the only difference being that it was looking for
approximations to PI, and did not attempt to factor the results. 
It did guarantee that there would be no common factors.

Here is that simple code again.  Modify it to input the value you
wish to approximate.

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

int main(int argc, char **argv)
{
   int           num, approx, bestnum, bestdenom;
   int           lastnum = 500;
   double        error, leasterr, pi, criterion;

   pi = 4 * atan(1.0);
   criterion = 2 * pi * DBL_EPSILON;
   if (argc > 1) lastnum = strtol(argv[1], NULL, 10);
   if (lastnum <= 0) lastnum = 500;
   printf("Usage: ratpi [maxnumerator]\n"
          "Rational approximation to PI = %.*f\n", DBL_DIG, pi);
   for (leasterr = pi, num = 3; num < lastnum; num++) {
      approx = (int)(num / pi + 0.5);
      error  = fabs((double)num / approx - pi);
      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 */

-- 
"If you want to post a followup via groups.google.com, don't use
 the broken "Reply" link at the bottom of the article.  Click on 
 "show options" at the top of the article, then click on the 
 "Reply" at the bottom of the article headers." - Keith Thompson

.