2个给定数字之间的分数密度
Density of fractions between 2 given numbers
我正在尝试对一个简单的 Fraction
class 进行一些分析,我想要一些数据将该类型与 doubles
.
进行比较
问题
我知道我正在寻找一些好的方法来获得 2 个数字之间分数的密度。分数基本上是 2 个整数(例如 pair< long, long>
),s
和 t
之间的密度是该范围内可表示数字的数量。它需要在 O(1) 或非常快的时间内完成精确或非常好的近似。
为了让它更简单一点,假设我想要 s 和 t 之间的所有数字(不是分数)a/b,其中 0 <= s <= a/b < t <= M , 且 0 <= a,b <= M (b > 0, a, b 为整数)
例子
如果我的分数是一种只能数到 6 (M = 6) 的数据类型,并且我希望密度介于 0 和 1 之间,答案将是 12。这些数字是:
0, 1/6, 1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5, 5/6.
我已经想到了
一种非常幼稚的方法是遍历所有可能的分数,然后计算那些不能简化的分数。类似于:
long fractionsIn(double s, double t){
long density = 0;
long M = LONG_MAX;
for(int d = 1; d < floor(M/t); d++){
for(int n = ceil(d*s); n < M; n++){
if( gcd(n,d) == 1 )
density++;
}
}
return density;
}
但是gcd()
非常慢,所以它不起作用。我也尝试做一些数学,但我什么也学不好。
解决方案
感谢@m69 的回答,我为 Fraction = pair<Long,Long>
:
编写了这段代码
//this should give the density of fractions between first and last, or less.
double fractionsIn(unsigned long long first, unsigned long long last){
double pi = 3.141592653589793238462643383279502884;
double max = LONG_MAX; //i can't use LONG_MAX directly
double zeroToOne = max/pi * max/pi * 3; // = approx. amount of numbers in Farey's secuence of order LONG_MAX.
double res = 0;
if(first == 0){
res = zeroToOne;
first++;
}
for(double i = first; i < last; i++){
res += zeroToOne/(i * i+1);
if(i == i+1)
i = nextafter(i+1, last); //if this happens, i might not count some fractions, but i have no other choice
}
return floor(res);
}
主要变化是 nextafter,这对大数字 (1e17) 很重要
结果
正如我在开头解释的那样,我试图将 Fractions
与 double
进行比较。这是 Fraction = pair<Long,Long>
的结果(以及 我如何获得双打密度):
Density between 0,1: | 1,2 | 1e6,1e6+1 | 1e14,1e14+1 | 1e15-1,1e15 | 1e17-10,1e17 | 1e19-10000,1e19 | 1e19-1000,1e19
Doubles: 4607182418800017408 | 4503599627370496 | 8589934592 | 64 | 8 | 1 | 5 | 0
Fraction: 2.58584e+37 | 1.29292e+37 | 2.58584e+25 | 2.58584e+09 | 2.58584e+07 | 2585 | 1 | 0
0到1之间的密度
如果表示分数的整数在0~M范围内,则0(含)和1(不含)之间的分数密度为:
M: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
0~(1): 1 2 4 6 10 12 18 22 28 32 42 46 58 64 72 80 96 102 120 128 140 150 172 180 200 212 230 242 270 278 308 ...
这是 OEIS 上的序列 A002088。如果向下滚动到公式部分,您会找到有关如何对其进行近似的信息,例如:
Φ(n) = (3 ÷ π2) × n2 + O[n × (ln n)2/3 × (ln ln n)4/3]
(遗憾的是,没有给出有关 O[x] 部分中涉及的常数的更多详细信息。请参阅下面关于近似值质量的讨论。)
跨范围分布
从0到1的区间包含了M以内可以表示的唯一分数总数的一半;例如这是 M = 15(即 4 位整数)时的分布:
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
72 36 12 6 4 2 2 2 1 1 1 1 1 1 1 1
总共 144 个独特的分数。如果查看不同 M 值的序列,您会发现此序列中的步骤会收敛:
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1: 1 1
2: 2 1 1
3: 4 2 1 1
4: 6 3 1 1 1
5: 10 5 2 1 1 1
6: 12 6 2 1 1 1 1
7: 18 9 3 2 1 1 1 1
8: 22 11 4 2 1 1 1 1 1
9: 28 14 5 2 2 1 1 1 1 1
10: 32 16 5 3 2 1 1 1 1 1 1
11: 42 21 7 4 2 2 1 1 1 1 1 1
12: 46 23 8 4 2 2 1 1 1 1 1 1 1
13: 58 29 10 5 3 2 2 1 1 1 1 1 1 1
14: 64 32 11 5 4 2 2 1 1 1 1 1 1 1 1
15: 72 36 12 6 4 2 2 2 1 1 1 1 1 1 1 1
不仅0和1之间的密度是分数总数的一半,1和2之间的密度是四分之一,2和3之间的密度接近十二分之一,依此类推。
随着 M 值的增加,分数在 0-1、1-2、2-3 ... 范围内的分布收敛为:
1/2, 1/4, 1/12, 1/24, 1/40, 1/60, 1/84, 1/112, 1/144, 1/180, 1/220, 1/264 ...
这个数列可以从1/2开始计算,然后:
0-1: 1/2 x 1/1 = 1/2
1-2: 1/2 x 1/2 = 1/4
2-3: 1/4 x 1/3 = 1/12
3-4: 1/12 x 2/4 = 1/24
4-5: 1/24 x 3/5 = 1/40
5-6: 1/40 x 4/6 = 1/60
6-7: 1/60 x 5/7 = 1/84
7-8: 1/84 x 6/8 = 1/112
8-9: 1/112 x 7/9 = 1/144 ...
您当然可以直接计算这些值中的任何一个,不需要中间的步骤:
0-1: 1/2
6-7: 1/2 x 1/6 x 1/7 = 1/84
(另请注意,分布序列的后半部分由 1 组成;这些都是除以 1 的整数。)
近似给定区间内的密度
使用 OEIS 页面上提供的公式,您可以计算或近似计算 0-1 区间内的密度,然后乘以 2,这是可以表示为分数的唯一值的总数。
给定两个值 s 和 t,然后可以计算并求和间隔 s ~ s+1、s+1 ~ s+2、... t-1 ~ t 中的密度,或使用插值以获得更快但不太精确的近似值。
示例
假设我们使用 10 位整数,能够表示 0 到 1023 之间的值。使用从 OEIS 页面链接的 this table,我们发现 0~1 之间的密度为 318452,分数总数为636904。
如果我们想求区间s~t = 100~105内的密度:
100~101: 1/2 x 1/100 x 1/101 = 1/20200 ; 636904/20200 = 31.53
101~102: 1/2 x 1/101 x 1/102 = 1/20604 ; 636904/20604 = 30.91
102~103: 1/2 x 1/102 x 1/103 = 1/21012 ; 636904/21012 = 30.31
103~104: 1/2 x 1/103 x 1/104 = 1/21424 ; 636904/21424 = 29.73
104~105: 1/2 x 1/104 x 1/105 = 1/21840 ; 636904/21840 = 29.16
四舍五入得出总和:
32 + 31 + 30 + 30 + 29 = 152
暴力算法给出了这个结果:
32 + 32 + 30 + 28 + 28 = 150
因此,对于这个低 M 值和只有 5 个值的小区间,我们偏离了 1.33%。如果我们在第一个值和最后一个值之间使用线性插值:
100~101: 31.53
104~105: 29.16
average: 30.345
total: 151.725 -> 152
我们会得出相同的值。对于较大的间隔,所有密度的总和可能会更接近真实值,因为舍入误差会相互抵消,但线性插值的结果可能会变得不太准确。对于更大的 M 值,计算的密度应该与实际值收敛。
Φ(n) 的近似质量
使用这个简化的公式:
Φ(n) = (3 ÷ π2) × n2
结果几乎总是小于实际值,但对于 n ≥ 182,它们在 1% 以内,对于 n ≥ 1880,它们在 0.1% 以内,对于 n ≥ 19494,它们在 0.01% 以内。我建议硬编码较低的范围(可以找到前 50,000 个值 here),然后从近似值足够好的地方开始使用简化的公式。
这是一个简单的代码示例,其中 Φ(n) 的前 182 个值被硬编码。分布序列的逼近似乎增加了与Φ(n)的逼近相似幅度的误差,因此应该可以得到一个像样的逼近。代码简单地迭代间隔 s~t 中的每个整数并对分数求和。为了加快代码速度并仍然获得良好的结果,您可能应该计算间隔中几个点的分数,然后使用某种非线性插值。
function fractions01(M) {
var phi = [0,1,2,4,6,10,12,18,22,28,32,42,46,58,64,72,80,96,102,120,128,140,150,172,180,200,212,230,242,270,278,308,
324,344,360,384,396,432,450,474,490,530,542,584,604,628,650,696,712,754,774,806,830,882,900,940,964,1000,
1028,1086,1102,1162,1192,1228,1260,1308,1328,1394,1426,1470,1494,1564,1588,1660,1696,1736,1772,1832,1856,
1934,1966,2020,2060,2142,2166,2230,2272,2328,2368,2456,2480,2552,2596,2656,2702,2774,2806,2902,2944,3004,
3044,3144,3176,3278,3326,3374,3426,3532,3568,3676,3716,3788,3836,3948,3984,4072,4128,4200,4258,4354,4386,
4496,4556,4636,4696,4796,4832,4958,5022,5106,5154,5284,5324,5432,5498,5570,5634,5770,5814,5952,6000,6092,
6162,6282,6330,6442,6514,6598,6670,6818,6858,7008,7080,7176,7236,7356,7404,7560,7638,7742,7806,7938,7992,
8154,8234,8314,8396,8562,8610,8766,8830,8938,9022,9194,9250,9370,9450,9566,9654,9832,9880,10060];
if (M < 182) return phi[M];
return Math.round(M * M * 0.30396355092701331433 + M / 4); // experimental; see below
}
function fractions(M, s, t) {
var half = fractions01(M);
var frac = (s == 0) ? half : 0;
for (var i = (s == 0) ? 1 : s; i < t && i <= M; i++) {
if (2 * i < M) {
var f = Math.round(half / (i * (i + 1)));
frac += (f < 2) ? 2 : f;
}
else ++frac;
}
return frac;
}
var M = 1023, s = 100, t = 105;
document.write(fractions(M, s, t));
将 Φ(n) 的近似值与前 50,000 个值的列表进行比较表明,添加 M÷4 是公式第二部分的可行替代;我还没有对更大的 n 值进行测试,因此请谨慎使用。
蓝色:简化的公式。红色:改进的简化公式。
分布近似的质量
将M=1023的结果与暴力算法的结果进行比较,实际误差很小,不会超过-7或+6,并且在区间205~206以上它们仅限于- 1~+1。然而,很大一部分范围(57~1024)每个整数的分数少于 100 个,而在 171~1024 区间每个整数只有 10 个或更少的分数。这意味着 -1 或 +1 的小误差和舍入误差会对结果产生很大影响,例如:
interval: 241 ~ 250
fractions/integer: 6
approximation: 5
total: 50 (instead of 60)
为了改进每个整数分数很少的区间的结果,我建议将上述方法与范围最后一部分的单独方法结合起来:
范围最后一部分的替代方法
如前所述,并在代码示例中实现,范围的后半部分,M÷2 ~ M,每个整数有 1 个小数。还有,区间 M÷3 ~ M÷2 有 2;区间 M÷4 ~ M÷3 有 4。这当然又是 Φ(n) 数列:
M/2 ~ M : 1
M/3 ~ M/2: 2
M/4 ~ M/3: 4
M/5 ~ M/4: 6
M/6 ~ M/5: 10
M/7 ~ M/6: 12
M/8 ~ M/7: 18
M/9 ~ M/8: 22
M/10 ~ M/9: 28
M/11 ~ M/10: 32
M/12 ~ M/11: 42
M/13 ~ M/12: 46
M/14 ~ M/13: 58
M/15 ~ M/14: 64
M/16 ~ M/15: 72
M/17 ~ M/16: 80
M/18 ~ M/17: 96
M/19 ~ M/18: 102 ...
在这些间隔之间,一个整数可以有不同数量的分数,具体取决于 M 的确切值,例如:
interval fractions
202 ~ 203 10
203 ~ 204 10
204 ~ 205 9
205 ~ 206 6
206 ~ 207 6
区间204~205位于区间的边缘,因为M÷5=204.6;它有 6 + 3 = 9 个分数,因为 M 模 5 是 3。如果 M 是 1022 或 1024 而不是 1023,它将有 8 个或 10 个分数。 (这个例子很简单,因为 5 是素数;见下文。)
同样,我建议使用 Φ(n) 的硬编码值来计算范围最后一部分的分数数。如果您使用上面列出的前 17 个值,这将涵盖每个整数少于 100 个小数的范围部分,这样会将舍入误差的影响降低到 1% 以下。前 56 个值将为您提供 0.1%,前 182 个值将为您提供 0.01%。
连同 Φ(n) 的值,您可以硬编码每个模值的边缘间隔的分数数,例如:
modulo: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
M/ 2 1 2
M/ 3 2 3 4
M/ 4 4 5 5 6
M/ 5 6 7 8 9 10
M/ 6 10 11 11 11 11 12
M/ 7 12 13 14 15 16 17 18
M/ 8 18 19 19 20 20 21 21 22
M/ 9 22 23 24 24 25 26 26 27 28
M/10 28 29 29 30 30 30 30 31 31 32
M/11 32 33 34 35 36 37 38 39 40 41 42
M/12 42 43 43 43 43 44 44 45 45 45 45 46
M/13 46 47 48 49 50 51 52 53 54 55 56 57 58
M/14 58 59 59 60 60 61 61 61 61 62 62 63 63 64
M/15 64 65 66 66 67 67 67 68 69 69 69 70 70 71 72
M/16 72 73 73 74 74 75 75 76 76 77 77 78 78 79 79 80
M/17 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
M/18 96 97 97 97 97 98 98 99 99 99 99 100 100 101 101 101 101 102
这与以下内容完全相同:(Sum of phi(k))
其中 m <= k <= M
其中 phi(k)
是 Euler Totient Function and with phi(0) = 1
(as defined by the problem). There is no known closed form for this sum. However there are many optimizations known as mentioned in the wiki link. This is known as the Totient Summatory Function in Wolfram. The same website also links to the series: A002088 并提供了一些渐近近似值。
推理是这样的:考虑 {1/M, 2/M, ...., (M-1)/M, M/M}
形式的值的数量。所有那些可以减少到更小值的分数都不会被计算在 phi(M)
中,因为它们不是相对质数。它们会出现在另一个totient的总和中。
例如,phi(6) = 12
并且您有 1 + phi(6)
,因为您还计算了 0
。
我正在尝试对一个简单的 Fraction
class 进行一些分析,我想要一些数据将该类型与 doubles
.
问题
我知道我正在寻找一些好的方法来获得 2 个数字之间分数的密度。分数基本上是 2 个整数(例如 pair< long, long>
),s
和 t
之间的密度是该范围内可表示数字的数量。它需要在 O(1) 或非常快的时间内完成精确或非常好的近似。
为了让它更简单一点,假设我想要 s 和 t 之间的所有数字(不是分数)a/b,其中 0 <= s <= a/b < t <= M , 且 0 <= a,b <= M (b > 0, a, b 为整数)
例子
如果我的分数是一种只能数到 6 (M = 6) 的数据类型,并且我希望密度介于 0 和 1 之间,答案将是 12。这些数字是:
0, 1/6, 1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5, 5/6.
我已经想到了
一种非常幼稚的方法是遍历所有可能的分数,然后计算那些不能简化的分数。类似于:
long fractionsIn(double s, double t){
long density = 0;
long M = LONG_MAX;
for(int d = 1; d < floor(M/t); d++){
for(int n = ceil(d*s); n < M; n++){
if( gcd(n,d) == 1 )
density++;
}
}
return density;
}
但是gcd()
非常慢,所以它不起作用。我也尝试做一些数学,但我什么也学不好。
解决方案
感谢@m69 的回答,我为 Fraction = pair<Long,Long>
:
//this should give the density of fractions between first and last, or less.
double fractionsIn(unsigned long long first, unsigned long long last){
double pi = 3.141592653589793238462643383279502884;
double max = LONG_MAX; //i can't use LONG_MAX directly
double zeroToOne = max/pi * max/pi * 3; // = approx. amount of numbers in Farey's secuence of order LONG_MAX.
double res = 0;
if(first == 0){
res = zeroToOne;
first++;
}
for(double i = first; i < last; i++){
res += zeroToOne/(i * i+1);
if(i == i+1)
i = nextafter(i+1, last); //if this happens, i might not count some fractions, but i have no other choice
}
return floor(res);
}
主要变化是 nextafter,这对大数字 (1e17) 很重要
结果
正如我在开头解释的那样,我试图将 Fractions
与 double
进行比较。这是 Fraction = pair<Long,Long>
的结果(以及
Density between 0,1: | 1,2 | 1e6,1e6+1 | 1e14,1e14+1 | 1e15-1,1e15 | 1e17-10,1e17 | 1e19-10000,1e19 | 1e19-1000,1e19
Doubles: 4607182418800017408 | 4503599627370496 | 8589934592 | 64 | 8 | 1 | 5 | 0
Fraction: 2.58584e+37 | 1.29292e+37 | 2.58584e+25 | 2.58584e+09 | 2.58584e+07 | 2585 | 1 | 0
0到1之间的密度
如果表示分数的整数在0~M范围内,则0(含)和1(不含)之间的分数密度为:
M: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
0~(1): 1 2 4 6 10 12 18 22 28 32 42 46 58 64 72 80 96 102 120 128 140 150 172 180 200 212 230 242 270 278 308 ...
这是 OEIS 上的序列 A002088。如果向下滚动到公式部分,您会找到有关如何对其进行近似的信息,例如:
Φ(n) = (3 ÷ π2) × n2 + O[n × (ln n)2/3 × (ln ln n)4/3]
(遗憾的是,没有给出有关 O[x] 部分中涉及的常数的更多详细信息。请参阅下面关于近似值质量的讨论。)
跨范围分布
从0到1的区间包含了M以内可以表示的唯一分数总数的一半;例如这是 M = 15(即 4 位整数)时的分布:
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
72 36 12 6 4 2 2 2 1 1 1 1 1 1 1 1
总共 144 个独特的分数。如果查看不同 M 值的序列,您会发现此序列中的步骤会收敛:
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1: 1 1
2: 2 1 1
3: 4 2 1 1
4: 6 3 1 1 1
5: 10 5 2 1 1 1
6: 12 6 2 1 1 1 1
7: 18 9 3 2 1 1 1 1
8: 22 11 4 2 1 1 1 1 1
9: 28 14 5 2 2 1 1 1 1 1
10: 32 16 5 3 2 1 1 1 1 1 1
11: 42 21 7 4 2 2 1 1 1 1 1 1
12: 46 23 8 4 2 2 1 1 1 1 1 1 1
13: 58 29 10 5 3 2 2 1 1 1 1 1 1 1
14: 64 32 11 5 4 2 2 1 1 1 1 1 1 1 1
15: 72 36 12 6 4 2 2 2 1 1 1 1 1 1 1 1
不仅0和1之间的密度是分数总数的一半,1和2之间的密度是四分之一,2和3之间的密度接近十二分之一,依此类推。
随着 M 值的增加,分数在 0-1、1-2、2-3 ... 范围内的分布收敛为:
1/2, 1/4, 1/12, 1/24, 1/40, 1/60, 1/84, 1/112, 1/144, 1/180, 1/220, 1/264 ...
这个数列可以从1/2开始计算,然后:
0-1: 1/2 x 1/1 = 1/2
1-2: 1/2 x 1/2 = 1/4
2-3: 1/4 x 1/3 = 1/12
3-4: 1/12 x 2/4 = 1/24
4-5: 1/24 x 3/5 = 1/40
5-6: 1/40 x 4/6 = 1/60
6-7: 1/60 x 5/7 = 1/84
7-8: 1/84 x 6/8 = 1/112
8-9: 1/112 x 7/9 = 1/144 ...
您当然可以直接计算这些值中的任何一个,不需要中间的步骤:
0-1: 1/2
6-7: 1/2 x 1/6 x 1/7 = 1/84
(另请注意,分布序列的后半部分由 1 组成;这些都是除以 1 的整数。)
近似给定区间内的密度
使用 OEIS 页面上提供的公式,您可以计算或近似计算 0-1 区间内的密度,然后乘以 2,这是可以表示为分数的唯一值的总数。
给定两个值 s 和 t,然后可以计算并求和间隔 s ~ s+1、s+1 ~ s+2、... t-1 ~ t 中的密度,或使用插值以获得更快但不太精确的近似值。
示例
假设我们使用 10 位整数,能够表示 0 到 1023 之间的值。使用从 OEIS 页面链接的 this table,我们发现 0~1 之间的密度为 318452,分数总数为636904。
如果我们想求区间s~t = 100~105内的密度:
100~101: 1/2 x 1/100 x 1/101 = 1/20200 ; 636904/20200 = 31.53
101~102: 1/2 x 1/101 x 1/102 = 1/20604 ; 636904/20604 = 30.91
102~103: 1/2 x 1/102 x 1/103 = 1/21012 ; 636904/21012 = 30.31
103~104: 1/2 x 1/103 x 1/104 = 1/21424 ; 636904/21424 = 29.73
104~105: 1/2 x 1/104 x 1/105 = 1/21840 ; 636904/21840 = 29.16
四舍五入得出总和:
32 + 31 + 30 + 30 + 29 = 152
暴力算法给出了这个结果:
32 + 32 + 30 + 28 + 28 = 150
因此,对于这个低 M 值和只有 5 个值的小区间,我们偏离了 1.33%。如果我们在第一个值和最后一个值之间使用线性插值:
100~101: 31.53
104~105: 29.16
average: 30.345
total: 151.725 -> 152
我们会得出相同的值。对于较大的间隔,所有密度的总和可能会更接近真实值,因为舍入误差会相互抵消,但线性插值的结果可能会变得不太准确。对于更大的 M 值,计算的密度应该与实际值收敛。
Φ(n) 的近似质量
使用这个简化的公式:
Φ(n) = (3 ÷ π2) × n2
结果几乎总是小于实际值,但对于 n ≥ 182,它们在 1% 以内,对于 n ≥ 1880,它们在 0.1% 以内,对于 n ≥ 19494,它们在 0.01% 以内。我建议硬编码较低的范围(可以找到前 50,000 个值 here),然后从近似值足够好的地方开始使用简化的公式。
这是一个简单的代码示例,其中 Φ(n) 的前 182 个值被硬编码。分布序列的逼近似乎增加了与Φ(n)的逼近相似幅度的误差,因此应该可以得到一个像样的逼近。代码简单地迭代间隔 s~t 中的每个整数并对分数求和。为了加快代码速度并仍然获得良好的结果,您可能应该计算间隔中几个点的分数,然后使用某种非线性插值。
function fractions01(M) {
var phi = [0,1,2,4,6,10,12,18,22,28,32,42,46,58,64,72,80,96,102,120,128,140,150,172,180,200,212,230,242,270,278,308,
324,344,360,384,396,432,450,474,490,530,542,584,604,628,650,696,712,754,774,806,830,882,900,940,964,1000,
1028,1086,1102,1162,1192,1228,1260,1308,1328,1394,1426,1470,1494,1564,1588,1660,1696,1736,1772,1832,1856,
1934,1966,2020,2060,2142,2166,2230,2272,2328,2368,2456,2480,2552,2596,2656,2702,2774,2806,2902,2944,3004,
3044,3144,3176,3278,3326,3374,3426,3532,3568,3676,3716,3788,3836,3948,3984,4072,4128,4200,4258,4354,4386,
4496,4556,4636,4696,4796,4832,4958,5022,5106,5154,5284,5324,5432,5498,5570,5634,5770,5814,5952,6000,6092,
6162,6282,6330,6442,6514,6598,6670,6818,6858,7008,7080,7176,7236,7356,7404,7560,7638,7742,7806,7938,7992,
8154,8234,8314,8396,8562,8610,8766,8830,8938,9022,9194,9250,9370,9450,9566,9654,9832,9880,10060];
if (M < 182) return phi[M];
return Math.round(M * M * 0.30396355092701331433 + M / 4); // experimental; see below
}
function fractions(M, s, t) {
var half = fractions01(M);
var frac = (s == 0) ? half : 0;
for (var i = (s == 0) ? 1 : s; i < t && i <= M; i++) {
if (2 * i < M) {
var f = Math.round(half / (i * (i + 1)));
frac += (f < 2) ? 2 : f;
}
else ++frac;
}
return frac;
}
var M = 1023, s = 100, t = 105;
document.write(fractions(M, s, t));
将 Φ(n) 的近似值与前 50,000 个值的列表进行比较表明,添加 M÷4 是公式第二部分的可行替代;我还没有对更大的 n 值进行测试,因此请谨慎使用。
蓝色:简化的公式。红色:改进的简化公式。
分布近似的质量
将M=1023的结果与暴力算法的结果进行比较,实际误差很小,不会超过-7或+6,并且在区间205~206以上它们仅限于- 1~+1。然而,很大一部分范围(57~1024)每个整数的分数少于 100 个,而在 171~1024 区间每个整数只有 10 个或更少的分数。这意味着 -1 或 +1 的小误差和舍入误差会对结果产生很大影响,例如:
interval: 241 ~ 250
fractions/integer: 6
approximation: 5
total: 50 (instead of 60)
为了改进每个整数分数很少的区间的结果,我建议将上述方法与范围最后一部分的单独方法结合起来:
范围最后一部分的替代方法
如前所述,并在代码示例中实现,范围的后半部分,M÷2 ~ M,每个整数有 1 个小数。还有,区间 M÷3 ~ M÷2 有 2;区间 M÷4 ~ M÷3 有 4。这当然又是 Φ(n) 数列:
M/2 ~ M : 1
M/3 ~ M/2: 2
M/4 ~ M/3: 4
M/5 ~ M/4: 6
M/6 ~ M/5: 10
M/7 ~ M/6: 12
M/8 ~ M/7: 18
M/9 ~ M/8: 22
M/10 ~ M/9: 28
M/11 ~ M/10: 32
M/12 ~ M/11: 42
M/13 ~ M/12: 46
M/14 ~ M/13: 58
M/15 ~ M/14: 64
M/16 ~ M/15: 72
M/17 ~ M/16: 80
M/18 ~ M/17: 96
M/19 ~ M/18: 102 ...
在这些间隔之间,一个整数可以有不同数量的分数,具体取决于 M 的确切值,例如:
interval fractions
202 ~ 203 10
203 ~ 204 10
204 ~ 205 9
205 ~ 206 6
206 ~ 207 6
区间204~205位于区间的边缘,因为M÷5=204.6;它有 6 + 3 = 9 个分数,因为 M 模 5 是 3。如果 M 是 1022 或 1024 而不是 1023,它将有 8 个或 10 个分数。 (这个例子很简单,因为 5 是素数;见下文。)
同样,我建议使用 Φ(n) 的硬编码值来计算范围最后一部分的分数数。如果您使用上面列出的前 17 个值,这将涵盖每个整数少于 100 个小数的范围部分,这样会将舍入误差的影响降低到 1% 以下。前 56 个值将为您提供 0.1%,前 182 个值将为您提供 0.01%。
连同 Φ(n) 的值,您可以硬编码每个模值的边缘间隔的分数数,例如:
modulo: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
M/ 2 1 2
M/ 3 2 3 4
M/ 4 4 5 5 6
M/ 5 6 7 8 9 10
M/ 6 10 11 11 11 11 12
M/ 7 12 13 14 15 16 17 18
M/ 8 18 19 19 20 20 21 21 22
M/ 9 22 23 24 24 25 26 26 27 28
M/10 28 29 29 30 30 30 30 31 31 32
M/11 32 33 34 35 36 37 38 39 40 41 42
M/12 42 43 43 43 43 44 44 45 45 45 45 46
M/13 46 47 48 49 50 51 52 53 54 55 56 57 58
M/14 58 59 59 60 60 61 61 61 61 62 62 63 63 64
M/15 64 65 66 66 67 67 67 68 69 69 69 70 70 71 72
M/16 72 73 73 74 74 75 75 76 76 77 77 78 78 79 79 80
M/17 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
M/18 96 97 97 97 97 98 98 99 99 99 99 100 100 101 101 101 101 102
这与以下内容完全相同:(Sum of phi(k))
其中 m <= k <= M
其中 phi(k)
是 Euler Totient Function and with phi(0) = 1
(as defined by the problem). There is no known closed form for this sum. However there are many optimizations known as mentioned in the wiki link. This is known as the Totient Summatory Function in Wolfram. The same website also links to the series: A002088 并提供了一些渐近近似值。
推理是这样的:考虑 {1/M, 2/M, ...., (M-1)/M, M/M}
形式的值的数量。所有那些可以减少到更小值的分数都不会被计算在 phi(M)
中,因为它们不是相对质数。它们会出现在另一个totient的总和中。
例如,phi(6) = 12
并且您有 1 + phi(6)
,因为您还计算了 0
。