c++ - 2个给定数字之间的分数密度

Question

我正在尝试对一个简单的Fraction类进行一些分析，并且我想要一些数据将该类型与doubles进行比较。

问题

知道我正在寻找一种使2个数字之间的分数密度更高的方法。分数基本上是2个整数(例如pair< long, long>)，并且s和t之间的密度是该范围内可表示数字的数量。它必须是在O(1)中完成的精确或非常好的近似，或者非常快。

为了简化一点，假设我想要s和t之间的所有数字(而不是分数)a/b，其中0 <= s <= a/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>的结果(以及 here如何获得 double 的密度):

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

Answer 1

密度介于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 ÷ π²) × n² + 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

我们将获得相同的值(value)。对于较大的时间间隔，所有密度的总和可能会更接近于实际值，因为舍入误差会相互抵消，但是线性插值的结果可能会变得不太准确。对于更大的M值，计算出的密度应与实际值收敛。

的近似逼近质量

使用以下简化公式:

Φ(n) = (3 ÷ π²) × n²

结果几乎总是小于实际值，但对于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