Simulation of irrational numbers

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?

[Walter Roberson]

Irrational numbers are irrational in all numeric bases.

Floating point numbers are stored as the ratio of two integers. The ratio of two floating point numbers is again the ratio of two integers. No matter what the denominator is, you can select a base in which the value was an integer together with at most one digit after the point.

Therefore you cannot force a "very long repeat period" unless you define a representation base.

[John D'Errico]

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?

[Sordin]

Hi John D'Errico,

I am simulating a physical situation where the ratio of the period of two waves is not rational (or rather, it has a very long repeat period), hence the combined wave does not appear to be periodic.

Is it possible to write a code such that the ratio of and has a repeat period longer than the length of the simulation (so that it appears as if it was irrational)?

Answer 1

[John D'Errico]

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:

https://math.stackexchange.com/questions/377683/length-of-period-of-decimal-expansion-of-a-fraction

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

https://en.wikipedia.org/wiki/Euler%27s_totient_function

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)

http://mathworld.wolfram.com/FermatsLittleTheorem.html

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

Answer 2

[James Tursa]

The random number generators that come with MATLAB have very long periods. You can simply use them to form your ratio and it will also have a very long period.