莫比乌斯反演解决euler71

2017-09-09

题目

给你n，列出所有以1~n作为分母的最简分数按照从小到大排序。

当n=1000000时，求在3/7左侧第一个数的分子

（这个问题的思考过程中顺便解决了求解在x/y左侧一共有多少个数字的问题）

题目链接：https://projecteuler.net/problem=71

题解

一开始题目看错了，看成了求3/7左侧一共有多少个数。

最简单的想法是，对于1~n找到每个数作为分母的时候，小于3/7的分子的取值范围是[1,m]

那么在这n个循环的过程中，每次都在求解[1,m]有多少个数与作为分母的数互质

根据莫比乌斯反演的思想，参考贾志鹏线性筛：
$$
\sum_{d|n} \mu (d) = [n=1] \\
\begin{aligned}
\sum_{i=1}^m{[gcd(i,n)=1]} & = \sum_{i=1}^m{\sum_{d|gcd(i,n)}{\mu (d)}} \\
& = \sum_{i=1}^m{\sum_{d|i \space and \space d|n}{\mu (d)}} \\
& = \sum_{d|n}{\sum_{i=1 \space and \space d|i}^{m}{\mu (d)}} \\
& = \sum_{d|n}{\left\lfloor m/d \right\rfloor \mu (d)}
\end{aligned}
$$
每次循环里的一次操作只要考虑约数即可，在$O(\sqrt n)$解决，这个代码就呼之欲出了。但并没有达到题目要求解的答案。对于题目要求的球左侧第一个数字，我们可以对于每个分母，求出小于当前数字3/7的和分母互质的最大分子，这个求的方法就是利用求互质数字的个数二分求解。

这样总的时间复杂度就是$O(n \sqrt n logn)$，所以运行时间达到了48s…代码如下：


namespace euler71{
    int pri[N] , mu[N+5] , tot;
    bool vis[N+5];
    void mubius_init(int n)
    {
        clr(vis,0); tot = 0;
        mu[0] = 0 , mu[1] = 1;
        for(int i=2 ; i<=n ; i++){
            if(!vis[i]){
                pri[tot++] = i;
                mu[i] = -1;
            }
            for(int j=0 ; j<tot ; j++){
                if((LL)pri[j]*i>n) break;
                vis[pri[j]*i] = true;
                if(i%pri[j] == 0){
                    mu[i*pri[j]] = 0;
                    break;
                }
                else mu[i*pri[j]] = -mu[i];
            }
        }
    }
    //找到1~m中有多少个数和n互质
    int num_of_relative_prime(int m , int n)
    {
        int p = (int) sqrt(n+0.5);
        int res = 0;
        for(int d=1 ; d<=p ; d++){
            if(n%d == 0){
                res += mu[d]*(m/d);
                if(n/d != d){
                    res += mu[n/d]*(m/(n/d));
                }
            }
        }
        return res;
    }
    //找到m往左数第一个跟n互质的数字
    int bin_search(int m , int n)
    {
        int cnt = num_of_relative_prime(m, n);
       // cout<<"bin_search: "<<m<<" "<<n<<" "<<cnt<<endl;
        int lf=1 , rg=m , mid , ans=0;
        while(lf<=rg){
            mid = (lf+rg)>>1;
            if(num_of_relative_prime(mid, n) == cnt) ans = mid , rg = mid-1;
            else lf = mid+1;
        }
        return ans;
    }
	//本题要求的答案解决的接口，分母为1~n组合出的最简分数有序序列，输出x/y左侧第一个数的分子
    int solution(int x , int y , int n)
    {
        mubius_init(n);
        double mx = 0; 
        int numerator = -1; //分子
        int fraction = -1;
        for(int i=2 ; i<=n ; i++){
            int m = i*x/y;
            int val = bin_search(i==y?m-1:m , i);
            
            if(val<0) continue;
            if(val*1.0 / i > mx) mx = val*1.0 / i , numerator = val , fraction = i;
        }
        cout<<numerator<<"/"<<fraction<<endl;
        return numerator;
    }
    //求所有以1~n作为分母的小于1的最简分数小于x/y的数目
    LL solve(int x , int y , int n)
    {
        mubius_init(n);
        LL res = 0;
        for(int i=1 ; i<=n ; i++){
            int m = i*x/y;
            res = res+num_of_relative_prime(m , i);
           // cout<<i<<" "<<m<<" "<<res<<endl;
        }
        return res - 1; //自己本身被算了一遍
    }
}
int main()
{
    cout<<euler71::solution(3,7,1000000)<<endl;
   // cout<<euler71::solve(3,7,1000000)<<endl;
	return 0;
}
//428570/999997 , 对于euler71需要的答案是分子428570，running 48.224 seconds

但是这个做法显然太复杂又麻烦了，换个简单的方法就是，外层的n的循环对于每个分母都能找到一个最大的分子，无所谓他们是否互质，不互质就大家除以公约数即可，反正大小一定保证是最大的解。

namespace euler71{
    LL gcd(LL x , LL y){return y==0?x:gcd(y,x%y);}
    LL solve(int x , int y , int n){
        LL x1 = 0 , x2 = 1; //x1/x2
        for(int i=1 ; i<=n ; i++){
            LL m = (i*x-1)/y;
            if(x1*i < x2*m)
                x1 = m , x2 = i; 
           // cout<<i<<" "<<m<<" "<<res<<endl;
        }
        return x1/gcd(x1,x2); //自己本身被算了一遍
    }
}
int main()
{
    cout<<euler71::solve(3,7,1000000)<<endl;
	return 0;
}
//428570/999997 , 对于euler71需要的答案是分子428570，running 48.224 seconds