如何快速计算指定区间内分数类型的密度以对比double类型

我最近在分析一个简单的Fraction类（本质是pair<long, long>形式的分子分母对），核心目标是对比它和double类型在不同区间的数值表示密度——也就是统计区间内可表示的最简分数数量，要求计算是精确值或高精度近似，且时间复杂度得是O(1)或者极快。

问题定义 #

我们要找所有满足0<=s<=a/b<t<=M、0<=a,b<=M（b>0，a、b为整数）的最简分数a/b。举个直观例子：当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。

最初的尝试（朴素方法） #

一开始我写了个遍历所有可能分数的代码，通过统计不可约分的数量来计算密度，但gcd()操作的开销太大，完全没法在大数场景下使用；同时尝试数学推导也没得到有效的简化公式。

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;
}

优化后的解决方案 #

后来感谢@m69的思路，我基于Farey序列的近似公式，编写了适配Fraction=pair<Long,Long>的代码，还加入了nextafter来处理超大数值（比如1e17）的精度问题：

//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);
}

最终对比结果 #

用这个方法得到的Fraction与double类型的密度对比数据如下：

内容的提问来源于stack exchange，提问作者Pernoctador