Rational approximations

[CBFalconer]

Here is an amusing toy, which may be of especial interest to the
Forthians who use rational fractions to represent real numbers.
The handling of criterion may also be interesting to some.

/* 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);
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 */

 

[Logan Shaw]

Here is an amusing toy, which may be of especial interest to the
Forthians who use rational fractions to represent real numbers.
The handling of criterion may also be interesting to some.

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

Actually, doubles are already rational numbers. The denominator
is always power of 2.

Now the numbers they approximate, that's a different story...

- Logan

 

[jacob navia]

CBFalconer a écrit :

Here is an amusing toy, which may be of especial interest to the
Forthians who use rational fractions to represent real numbers.
The handling of criterion may also be interesting to some.

/* 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);
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 */

Excuse me Chuck but I think you invoke UB in the first iteration
when we have:

num = 1
value = PI
we have then:
1/PI = 0.31830 ....
(int) 0.31830 == ZERO

Then you do

num/approx and you have a division by zero...

jacob

 

[James Dow Allen]

/* Find best rational approximation to a double */

One could try to argue they solve a slightly different problem,
but continued fractions provide optimal approximations, are
easy to code, and yield fascinating results, eg.
e = 2.71828... is {2;1,2,1,1,4,1,1,6,1,1,8,1,1,10,...] in
continued fraction form, and
the Golden Ratio is [1;1,1,1,1,1,1,1,...].

(This c.f. form of Golden Ratio means it is the real number
hardest to approximate by a rational! The Golden Ratio
bus scheduling algorithm depends on that property IIRC,
though most comp-sci algorithms called "Golden Ratio
algorithm" don't.)

As is so often the case these days, Wikipedia seems to
offer a good discussion:
http://en.wikipedia.org/wiki/Continued_fraction

value = 4 * atan(1.0);

Can't remember how to spell M_PI, eh? Neither can I, though
I usually just spell it prosaically: 3.141592635897

James Dow Allen

 

[Richard Bos]

Can't remember how to spell M_PI, eh?

No, he can remember that M_PI is not in ISO C.

Richard

 

[CBFalconer]

CBFalconer a écrit :

Here is an amusing toy, which may be of especial interest to the
Forthians who use rational fractions to represent real numbers.
The handling of criterion may also be interesting to some.

/* 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);
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 */

Excuse me Chuck but I think you invoke UB in the first iteration
when we have:

num = 1
value = PI
we have then:
1/PI = 0.31830 ....
(int) 0.31830 == ZERO

Then you do

num/approx and you have a division by zero...

You are right. Sloppy of me. I added the line:

if (0 == (int)approx) continue;

after the "apprex = ..." line.

--
"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
More details at: <http://cfaj.freeshell.org/google/>
Also see <http://www.safalra.com/special/googlegroupsreply/>

 

[dcorbit]

Reminiscient of:
http://www.acm.org/pubs/tog/GraphicsGems/gemsv/ch1-4/rat.h
http://www.acm.org/pubs/tog/GraphicsGems/gemsv/ch1-4/rat.c

 

[Peter J. Holzer]

["Followup-To:" header set to comp.programming.]

Can't remember how to spell M_PI, eh? Neither can I, though
I usually just spell it prosaically: 3.141592635897

I just found a spare "5". Would you like to buy it?

hp

 

[Keith Thompson]

Reminiscient of:
http://www.acm.org/pubs/tog/GraphicsGems/gemsv/ch1-4/rat.h
http://www.acm.org/pubs/tog/GraphicsGems/gemsv/ch1-4/rat.c

What is?

Please provide context.

 

[dcorbit]

What is?

C. B. Falconer's algorithm is like the graphics gems algorithm
(continued fraction to find best numerator/denomiator to fit a float).
The gem is somewhat more generalized. CB's algorithm is more compact.

 

[CBFalconer]

C. B. Falconer's algorithm is like the graphics gems algorithm
(continued fraction to find best numerator/denomiator to fit a float).
The gem is somewhat more generalized. CB's algorithm is more compact.

Whew. I'll say it is (more compact). Gem?