# Re: rational to floating point

## Peter-Lawrence.Montgomery@cwi.nl (Peter L. Montgomery)14 Mar 2003 11:52:26 -0500

From comp.compilers

Related articles
[2 earlier articles]
Re: rational to floating point nmm1@cus.cam.ac.uk (2003-03-14)
Re: rational to floating point thant@acm.org (Thant Tessman) (2003-03-14)
Re: rational to floating point joachim_d@gmx.de (Joachim Durchholz) (2003-03-14)
Re: rational to floating point ajo@andrew.cmu.edu (Arthur J. O'Dwyer) (2003-03-14)
Re: rational to floating point gah@ugcs.caltech.edu (Glen Herrmannsfeldt) (2003-03-14)
Re: rational to floating point tmk@netvision.net.il (2003-03-14)
Re: rational to floating point Peter-Lawrence.Montgomery@cwi.nl (2003-03-14)
Re: rational to floating point francis@thibault.org (John Stracke) (2003-03-14)
| List of all articles for this month |

 From: Peter-Lawrence.Montgomery@cwi.nl (Peter L. Montgomery) Newsgroups: comp.compilers Date: 14 Mar 2003 11:52:26 -0500 Organization: CWI, Amsterdam References: 03-03-035 Keywords: arithmetic Posted-Date: 14 Mar 2003 11:52:26 EST

Thant Tessman <thant@acm.org> writes:
>The question is: Under what conditions will a rational number produce
>an infinite stream of digits for a given base? ...
>
>For the curious, a description and source for my interpreter can be
>found at:
>
> http://www.standarddeviance.com

When b = 10 (decimal), the rational number n/d = 1/6 lacks a finite
floating point expansion, even though gcd(d, b) = 2 > 1.
The problem is that d has a prime factor 3 which does not divide b.

Start by reducing n/d to lowest terms: if gcd(n, d) = g and
g > 1, replace n by n/g and d by d/g. Now try something like

bexp = 0; // Answer is b^bexp * (n/d)
do {
g = gcd(b, d);
d = d/g;
n = n*(b/g);
bexp--;
} while (g > 1);

This, given n/d = 1/6 and b = 10, will detect gcd(b, d) = g = 2.
Multiply numerator and denominator by b/g = 5, getting 5/30.
Remove the factor 10 from the denominator, getting (5/3)*10^(-1).
Now, since gcd(3, 10) = 1, terminate. The denominator 3 is not
equal to 1, so 1/6 lacks a finite expansion.

If you start with 3/8, you get 15/4 * 10^(-1), (75/2) * 10^(-2),
and (375/1) * 10^(-3). The denominator is now 1, so 3/8 has
a finite decimal expansion .375. Observe that the final numerator
375 is much larger than the original numerator or denominator.

If you merely want to detect finiteness, not produce the expansion,
you can try something like

d = d/gcd(d, n); // Lowest terms denominator
g = gcd(d, b);
while (g != 1) {
d = d/g;
g = gcd(d, g*g);
}
test whether g == 1

Using g*g in the last gcd reduces the number of iterations when a high
power of a prime divides the original d.
--
During the Korean train firebombing, there was a MAD DASH from the SUB.
Spell MAD DASH SUB backwards -- what names emerge?
Peter-Lawrence.Montgomery@cwi.nl Home: San Rafael, California
Microsoft Research and CWI

Post a followup to this message