# Simulation of irrational numbers

17 views (last 30 days)

Show older comments

##### 3 Comments

John D'Errico
on 31 May 2019

Edited: John D'Errico
on 31 May 2019

Anyway, I'm a bit confused.

"...generate two random numbers x1 and x2 such that their ratio x1/x2 is an irrational number."

By definition, the ratio of two numbers is never irrational, unless one of the numbers themselves was already irrational.

As the answers have said, you can generate any arbitrary sequence of digits of any length. That can be viewed as a decimal number, but what does that get you?

x = 0.45623626246273674527567786785685795686969679757080978965456345354854564

x =

0.456236262462737

So, MATLAB stores the number as a double, which can then be shown to be represented as:

sum(2 .^ [-2 -3 -4 -6 -9 -10 -13 -15 -16 -17 -18 -19 -22 -23 -26 -28 -31 -33 -38 -39 -41 -42 -44 -46 -47 -48 -49 -52])

ans =

0.456236262462737

As a sum of fractions (that is, negative powers of 2), it is clearly a rational number. Or, we can write that as the fraction:

2054705461620089/4503599627370496

ans =

0.456236262462737

2054705461620089/4503599627370496 == x

ans =

logical

1

That reproduces the double precision number, and does so exactly. But no matter what I do with a double precision number, it is not the original very long decimal number I typed in. At best, it would seem you could use a symbolic form.

But where are you going with this? You cannot easily recover the original numbers that formed the fractional approximation. And as a random string of digits, what use are they?

### Accepted Answer

John D'Errico
on 1 Jun 2019

Edited: John D'Errico
on 1 Jun 2019

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. ;-)

##### 0 Comments

### More Answers (1)

James Tursa
on 31 May 2019

##### 0 Comments

### See Also

### Categories

### Products

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!