Simulation of irrational numbers

조회 수: 28 (최근 30일)
Sordin
Sordin 2019년 5월 31일
편집: John D'Errico 2019년 6월 1일
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?
  댓글 수: 3
John D'Errico
John D'Errico 2019년 5월 31일
편집: John D'Errico 2019년 5월 31일
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
Sordin 2019년 5월 31일
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)?

댓글을 달려면 로그인하십시오.

채택된 답변

John D'Errico
John D'Errico 2019년 6월 1일
편집: John D'Errico 2019년 6월 1일
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. ;-)

추가 답변 (1개)

James Tursa
James Tursa 2019년 5월 31일
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.

카테고리

Help CenterFile Exchange에서 Number Theory에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by