testing whether a double is a whole number

[David Marsh]

The program calculates the continued fraction representation of the input:

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

int main(int argc, char* argv[])
{
double diff, n, r, i;

if(argc != 2) exit(EXIT_FAILURE);

n = strtod(argv[1], NULL);

printf("\n [");

while(1)
{
i = (int)n;
printf(" %.f", i);
diff = n - i;
r = 1 / diff;
if(r is a whole number) /* ??? */
{
printf(" %.f", r);
break;
}
n = r;
}

printf("]\n");

return 0;
}

In the while() loop, I want to test whether r is a whole number. If it
is, it is the last denominator of the continued fraction and I'm done.
Incredibly, I was not able to do this in any straighforward way. For
example, if(r == (int)r) didn't work. I finally kludged it by writing a
function that counts the significant decimal digits in a real number,
but there has got to be a better way. How would you solve the problem?

(And, yes, I realize that if the input is irrational the program doesn't
terminate.)

David

 

[Walter Roberson]

The program calculates the continued fraction representation of the input:

r = 1 / diff;
if(r is a whole number) /* ??? */

In the while() loop, I want to test whether r is a whole number. If it
is, it is the last denominator of the continued fraction and I'm done.
Incredibly, I was not able to do this in any straighforward way. For
example, if(r == (int)r) didn't work. I finally kludged it by writing a
function that counts the significant decimal digits in a real number,
but there has got to be a better way. How would you solve the problem?

You don't. Your entire sequence is subject to increasing round-off
error. Even if r -appeared- to be a whole number, you wouldn't
be able to tell whether you had calculated the exact continued
fraction, or had instead merely ended up rounding off to that
value.

The algorithm you are using to calculate the continued fraction
only works when there is indefinite precision.

Consider, for example, something as simple as the input
value (double) 1.0 / (double) 3.0 . No matter how many bits
of precision you have, if you are using a binary representation
then you have the binary repeating fraction 01 (i.e.,
1/3 is binary 0.0101010101...) That is [using fraction notation]
1/4 + 1/16 + 1/64 + 1/256 + 1/1024 + 1/4096 + ...
and after 2N mantisa bits, your value is different from 1/3 by
the value 1/(3 * 2**(2N)).

If you expand this logic, then because of limited precision
on input, then even if you didn't lose precision during
the calculation, only the fractions which are sums of
negative powers of 2 are going to come out properly,
and *every* fraction that involves a prime greater than 2
would merely come out expressed in terms of negative
powers of 2, to the limit of the number of mantisa bits.


So... you cannot do any meaningful accurate continued-fraction
calculation with the algorithm you give.

The closest you can get is to lose a bit of precision,
understand that wrong answers will be given sometimes,
and change the "if" to check to see whether r is within
a tolerance of a whole number. For example,

if(r is a whole number) /* ??? */

would become something like

if ( abs(r - trunc(r)) < delta )

and you get to choose the delta according to how much slop you
are willing to put up with.

 

[Martin Ambuhl]

The program calculates the continued fraction representation of the input:

In the while() loop, I want to test whether r is a whole number. If it
is, it is the last denominator of the continued fraction and I'm done.
Incredibly, I was not able to do this in any straighforward way. For
example, if(r == (int)r) didn't work. I finally kludged it by writing a
function that counts the significant decimal digits in a real number,
but there has got to be a better way. How would you solve the problem?

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

#define ALMOST_ZERO (1.e-10)
#define MAXV 20


int main(int argc, char *argv[])
{
double thenumber, copy;
double denominator[MAXV], n, d, t;
int pass;

if (argc != 2) {
fprintf(stderr, "usage: %s number\n", argv[0]);
exit(EXIT_FAILURE);
}
thenumber = strtod(argv[1], 0);
if (thenumber < 0) {
printf("changing sign of %.f\n", thenumber);
thenumber *= -1;
}
printf("[ ");
for (copy = thenumber, pass = 0; pass < 20; pass++) {
copy = modf(copy, &denominator[pass]);
printf("%.f ", denominator[pass]);
if (copy < ALMOST_ZERO)
break;
copy = 1. / copy;
}
printf("]\n");
if (pass == 20)
pass--;
if (pass == 0) {
n = denominator[0];
d = 1;
}
else {
for (n = 1, d = denominator[pass], pass--; pass > 0; pass--) {
t = d * denominator[pass] + n;
n = d;
d = t;
}
n += d * denominator[0];
}
printf("%.f/%.f = %.*g\n"
" the original number was %.*g\n", n, d, DBL_DIG,
n / d, DBL_DIG, thenumber);

return 0;
}

 

[David Marsh]

You don't. Your entire sequence is subject to increasing round-off
error. Even if r -appeared- to be a whole number, you wouldn't
be able to tell whether you had calculated the exact continued
fraction, or had instead merely ended up rounding off to that
value.

The algorithm you are using to calculate the continued fraction
only works when there is indefinite precision.

So... you cannot do any meaningful accurate continued-fraction
calculation with the algorithm you give.

The closest you can get is to lose a bit of precision,
understand that wrong answers will be given sometimes,
and change the "if" to check to see whether r is within
a tolerance of a whole number. For example,

OK, I think I understand the precision issue. The source of the
algorithm (and everything I know about continued fractions) is:
http://en.wikipedia.org/wiki/Continued_fraction (Calculating continued
fraction representations).

Martin Ambuhl's program uses the fudge factor approach. I compared the
output of his program with mine. In every case my output was comparable
to his. Some examples below:

input David's program Martin's program
---------------------------------------------

..3 0;3,3 0;3,2,1
..5 0;2 0;2
..75 0;1,3 0;1,3
..829 0;1,4,1,5,1,1,2,1,3 0;1,4,1,5,1,1,2,1,3
..999 0;1,999 0;1,998,1

Anyway, it has become a discussion of algorithms and CS, not standard C,
so it's off-topic here. Thanks for your comments. BTW, are you Canadian?

David