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("Rationa l approximation to %.*f\n", DBL_DIG, value);
for (leasterr = value, num = 1; num < lastnum; num++) {
approx = (int)(num / value + 0.5);
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 */

--
Chuck F (cbfalconer at maineline dot net)
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home .att.net>

[Logan Shaw]

CBFalconer wrote:

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("Rationa l approximation to %.*f\n", DBL_DIG, value);
for (leasterr = value, num = 1; num < lastnum; num++) {
approx = (int)(num / value + 0.5);
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 */

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]

CBFalconer wrote:

/* 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,1 0,...] 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]

"James Dow Allen" <jd*********@ya hoo.comwrote:

CBFalconer wrote:

value = 4 * atan(1.0);

Can't remember how to spell M_PI, eh?

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

Richard

[CBFalconer]

jacob navia wrote:

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("Rationa l approximation to %.*f\n", DBL_DIG, value);
for (leasterr = value, num = 1; num < lastnum; num++) {
approx = (int)(num / value + 0.5);
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 */

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.c om, 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/googlegroupsrep ly/>

[dcorbit]

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

[Keith Thompson]

dc*****@connx.c om writes:

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

What is?

Please provide context.

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[Peter J. Holzer]

["Followup-To:" header set to comp.programmin g.]
On 2006-09-08 07:43, James Dow Allen <jd*********@ya hoo.comwrote:

CBFalconer wrote:

> 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

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

hp
--
_ | Peter J. Holzer | Wieso sollte man etwas erfinden was nicht
|_|_) | Sysadmin WSR | ist?
| | | hj*@hjp.at | Was sonst wäre der Sinn des Erfindens?
__/ | http://www.hjp.at/ | -- P. Einstein u. V. Gringmuth in desd

[dcorbit]

Keith Thompson wrote:

dc*****@connx.c om writes:

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

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.

Please provide context.

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.