Asked by Sordin
on 31 May 2019

I am trying to generate two random numbers and such that their ratio is an irrational number. I understand that all numbers stored on a computer are rational, so one cannot have a truly irrational number in a simulation/experiment. So, is there any way of simulating the situation where the ratio of the two numbers has a very long repeat period?

Answer by John D'Errico
on 1 Jun 2019

Edited by John D'Errico
on 1 Jun 2019

Accepted Answer

Perhaps now I am starting to understand what your question comes down to.

The ratio of x1 and x2, if they are doubles, will never have a period longer than 16 digits, because that fraction, as a double precision number, only has roughly 16 significant digits. There simply are not more than 16 digits to which you can attribute any meaning in a double. Worse, as I point out in my initial comment to your question, ALL double precisionn numbers are rational, and non-repeating decimals.

If they are fully general integers, then the ratio of two integers can in theory have as long a period to repeat as you wish. You will just need to pick two numbers (mainly you need to pick the divisor) such that it has a long period. But that is not something that is easily controlled. In general, the numerator of a rational fraction does not contribute to the period of the decimal form, although it may reduce the period. The numerator (n) of a fraction n/d cannot increase the period beyond that of 1/d, however, it may yield a shorter period.

For example, the number 1/7 has an infinitely repeating period of length 6. However, the number 7/7 is an integer, so it has only repeating zeros after the decimal point. Or, the number 1/21 is a repeating decimal with period 6 (the "047619" repeats) but 7/21 repeats with period 1, since it reduces to 1/3.

So, what controls the period of a fraction? A quick link that suggests the solution is:

The problem is that you don't get a truly simple formula to use. But in general, consider the fraction 1/d. If you want the longest potential period, then you need to consider prime numbers for d. Why? Because the period is a function of the totient function of d. Euler's Totient function (I offer a totient calculator in my vpi toolbox. But the totient is easy to compute as long as you know the prime factorization of your number.)

It is usually written using the name phi. phi(d) is the number of integers m less than d, such that gcd(d,m)==1. Clearly the totient of a prime number is 1 less than the number itself. We can see how it works, by some examples:

The length of the repeating portion of a rational fraction n/d will be a factor of phi(d). So we can see that the MAXIMUM length of the period of a repeating fraction is d-1, since the largest possible value for the totient function of a number is 1 less than that number.

Consider 1/7. It has a period of 6. And phi(7)=6. And 6 is a factor of 6. So it works.

digits 50

vpa(sym(1)/7)

ans =

0.14285714285714285714285714285714285714285714285714

How about 1/13? It also has a period of 6. And phi(13)=12. 6 is clearly a factor of 12.

vpa(sym(1)/13)

ans =

0.076923076923076923076923076923076923076923076923077

So not all prime numbers have maximum length periods, because the requirement is merely that the period length be a FACTOR of totient(d). And it is not that easy to predict the period length either. Here is a good one, for a moderately small denominator.

digits 200

vpa(sym(1)/61)

ans =

0.016393442622950819672131147540983606557377049180327868852459016393442622950819672131147540983606557377049180327868852459016393442622950819672131147540983606557377049180327868852459016393442622950819672

Of ourse, phi(61) is 60, since 61 is a prime. And 1/61 has a repeating period of 60, with this repeating fragment:

016393442622950819672131147540983606557377049180327868852459

Is there an easy way to identify that fact? The trick can be to use a powermod utility. Lets see, the factors of 60 are given by my aliquotparts function:

aliquotparts(60)

ans =

1 2 3 4 5 6 10 12 15 20 30

The period (p) of 1/61 is then the SMALLEST power of 10, such that mod(10^p,61) == 1. So now, we test each power of 10.

p = aliquotparts(60)

p =

1 2 3 4 5 6 10 12 15 20 30

for p_i = p

powermod(10,p_i,61)

end

ans =

10

ans =

39

ans =

24

ans =

57

ans =

21

ans =

27

ans =

14

ans =

58

ans =

50

ans =

13

ans =

60

And of course, we have:

powermod(10,60,61)

ans =

1

So the smallest power of 10 such that (10^p_i-1) is divisible by 61 is in fact 60. Therefore the period of 1/61 has repeating part of length 60. In fact, we don't need to actually test 60 there, because it is provably true. That is, Fermat's theorem (sometimes called Fermat's Little Theorem)

tells us that for ANY prime p, if p does not dive a, then this relation always holds fpr prime p, and integer a > 1:

mod(a^(p-1) , p) == 1

So as long as p is not 2 or 5, then we know that

powermod(10,p-1,p) == 1

must hold true. The conclusion is that any prime that is not 2 or 5 will have a perion of at most p-1, when the fraction is represented in decimal digits.

So, do you want a fraction with a HUGE repeating period? Pick a large prime denominator. Then run a quick test.

d = nextprime(1e8)

d =

100000007

totient(d)

ans =

100000006

p = aliquotparts(totient(d))

p =

1 2 491 982 101833 203666 50000003

for p_i = p

powermod(10,p_i,d)

end

ans =

10

ans =

100

ans =

5539135

ans =

14400485

ans =

24154401

ans =

46828351

ans =

100000006

powermod(10,d-1,d)

ans =

1

If any of those intermediate results had been 1, then the period would be less than the maximum possible. For example, I claimed before that the period of 1/13 is 6? We can predict that would be true by this simple computation:

powermod(10,6,13)

ans =

1

And that is the smallest power of 10 that leaves a residue of 1 modulo 13, among the divisors of 12.

Anyway, the period of fractions of the form n/100000007 have period length 100000006. That is, just over 100 million decimal digits before it starts to repeat. Here are the first 200 digits of that number:

vpa(sym(1)/100000007)

ans =

0.0000000099999993000000489999965700002400999831930011764899176457057648005964639582475229226733954128623210996375230253733882238628243296022969278392150512549464121537511492374195533806312633558115650931904435

Why does the above test work? Not hard to prove, actually. It is actually more interesting to prove that for prime divisor d, the maximum value of the period is d-1. For this, a quick look at a theorem that goes back as far as Fermat will suffice.

One of the very pretty things (at least to me) is that this question, about something as prosaic as decimal fractions, actually uses some interesting results from number theory to solve.

What use this is to you? Your problem, not mine. ;-)

Sign in to comment.

Answer by James Tursa
on 31 May 2019

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 3 Comments

## Walter Roberson (view profile)

Direct link to this comment:https://se.mathworks.com/matlabcentral/answers/464911-simulation-of-irrational-numbers#comment_710490

## John D'Errico (view profile)

Direct link to this comment:https://se.mathworks.com/matlabcentral/answers/464911-simulation-of-irrational-numbers#comment_710503

## Sordin (view profile)

Direct link to this comment:https://se.mathworks.com/matlabcentral/answers/464911-simulation-of-irrational-numbers#comment_710508

Sign in to comment.