Рационалы перечислимы. Например, этот код находит k-е рациональное в открытом интервале 0..1 с упорядочением, что {n1, d1} до {n2, d2}, если (d1<d2 || (d1==d2 && n1<n2)), предполагая, что {n,d} является взаимным.
RankedRational[i_Integer?Positive] :=
Module[{sum = 0, eph = 1, den = 1},
While[sum < i, sum += (eph = EulerPhi[++den])];
Select[Range[den - 1], CoprimeQ[#, den] &][[i - (sum - eph)]]/den
]
In[118]:= Table[RankedRational[i], {i, 1, 11}]
Out[118]= {1/2, 1/3, 2/3, 1/4, 3/4, 1/5, 2/5, 3/5, 4/5, 1/6, 5/6}
Теперь я хотел бы генерировать случайные рациональности, учитывая верхнюю оценку сортировки знаменателя равномерно, так что при достаточно больших знаменателях рациональности будут равномерно распределены по единичному интервалу.
Интуитивно можно было бы выбрать среди всех рациональных значений с помощью малых знаменателей с равными весами:
RandomRational1[maxden_, len_] :=
RandomChoice[(Table[
i/j, {j, 2, maxden}, {i,
Select[Range[j - 1], CoprimeQ[#, j] &]}] // Flatten), len]
Можно ли генерировать случайные рациональности с этим распределением более эффективно, не конструируя их все? Это не займет много места, чтобы эта таблица стала огромной.
In[197]:= Table[RankedRational[10^k] // Denominator, {k, 2, 10}]
Out[197]= {18, 58, 181, 573, 1814, 5736, 18138, 57357, 181380}
Или, может быть, можно эффективно генерировать рациональные решения с ограниченным знаменателем, имеющим другое "чувственно подобранное" распределение?
РЕДАКТИРОВАТЬ. Это код Mathematica, который запускает прием-отклонение, предложенное btilly.
Clear[RandomFarey];
RandomFarey[n_, len_] := Module[{pairs, dim = 0, res, gcds},
Join @@ Reap[While[dim < len,
gcds = cfGCD[pairs = cfPairs[n, len - dim]];
pairs = Pick[pairs, gcds, 1];
If[pairs =!= {},
dim += [email protected][res = pairs[[All, 1]]/pairs[[All, 2]]]];
]][[2, -1]]
]
Следующие скомпилированные пары функций целых чисел {i,j} такие, что 1<=i < j<=n:
cfPairs =
Compile[{{n, _Integer}, {len, _Integer}},
Table[{i, RandomInteger[{i + 1, n}]}, {i,
RandomChoice[2 (n - Range[n - 1])/(n (n - 1.0)) -> Range[n - 1],
len]}]];
и следующая скомпилированная функция вычисляет gcd. Он предполагает, что вход представляет собой пару целых положительных чисел.
cfGCD = Compile[{{prs, _Integer, 1}}, Module[{a, b, p, q, mod},
a = prs[[1]]; b = prs[[2]]; p = Max[a, b]; q = Min[a, b];
While[q > 0, mod = Mod[p, q]; p = q; q = mod]; p],
RuntimeAttributes -> Listable];
Тогда
In[151]:= data = RandomFarey[12, 10^6]; // AbsoluteTiming
Out[151]= {1.5423084, Null}
In[152]:= cdf = CDF[EmpiricalDistribution[data], x];
In[153]:= Plot[{cdf, x}, {x, 0, 1}, ImageSize -> 300]
