TOJ 3653 第k互素数

题目链接：http://www.tzcoder.cn/acmhome/problemdetail.do?method=showdetail&id=3653

描述

请你写一个程序，计算与m和n都互为素数的数中第k大的数，假如X与Y互为素数，那么X和Y的最大公约数必为1。

输入

第一行包含一个整数T代表测试数据的组数
每组数据有一行，分别包含m，n，k，其中（0<m，n，k<=10^9）

输出

对应每组测试数据，每行输出该组组号和第k个与m，n互为素数的整数

思路：
对于1e9这样的范围，一个一个算出与m，n互质的数显然不现实，思考O(log(n))的方法。设p(i)为小于等于i的正整数中与m,n互质数的个数，可以发现p(i)<=p(i+1)，这个函数是递增的，可以用二分来找出第k个与n，m互质的数。二分的下界为1，关于上界可以由素数定理（Π(x)~x/ln(x)）大致求出小于3e10的质数大约有1,243,550,985个,大于1e9，因此设上界为3e10。
接下来考虑计算p(i)的方法，参考n/2可以求出在1~n中可以被2整除的数量，因此可以求出n，m所有的质因子并利用容斥求出p(i)（加上所有i除以奇数个质数乘积的商然后减去所有i除以偶数个质数乘积的商，最后得到小于等于i中不与n，m互质的个数）。
考虑方案的时间复杂度。设N为质数个数，计算p(i)的复杂度为O(2^(N))，考虑N的最大值，nmax=2×3×5×7×11×13×17×19×23=223,092,870，mmax=29×31×37×41×43=58,642,669，Nmax=14。
总复杂度为O(log2(3e10)*2^N)(去掉了计算质数的复杂度)在1<=N<=14的条件下最大值约为395,255，可行。
一开始用递归一直超时，后来发现才在求质数组合的时候多了很多重复计算……

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll Ms=3e10+5;
int q[15],cnt;
void cook(int x)
{
    int n=sqrt(x);
    for(int i=2;i<=n;i++)
    {
        if(x%i==0)
        {
            q[cnt++]=i;
            x/=i;
            while(x%i==0)x/=i;
        }
        if(x<i)break;
    }
    if(x!=1)q[cnt++]=x;
}
ll M;
ll a[2][3500];
ll cul()
{
    int c[2]={1},it=0;
    a[0][0]=1;
    int st[2][3500]={};
    ll ans=0,tmp;
    for(int i=0;i<cnt;i++)
    {
        for(int j=0;j<c[it];j++)
        {
            for(int d=st[it][j];d<cnt;d++)
            {
                //ff++;
                if(a[it][j]%q[d])
                {
                    tmp=a[it][j]*q[d];
                    if(tmp>M)continue;
                    a[!it][c[!it]]=tmp;
                    st[!it][c[!it]]=d+1;
                    c[!it]++;
                    ans+=it?-M/tmp:M/tmp;
                }
            }
        }
        c[it]=0;
        it=!it;
        //cout<<c[it]<<endl;
        if(c[it]==0)break;
    }
    return ans;
}
int main()
{
    int t;
    scanf("%d",&t);
    for(int no=1;no<=t;no++)
    {
        int n,m,k;
        scanf("%d%d%d",&n,&m,&k);
        cnt=0;
        cook(n);cook(m);
        sort(q,q+cnt);
        cnt=unique(q,q+cnt)-q;
        //for(int i=0;i<cnt;i++)cout<<q[i]<<" ";cout<<endl;
        ll L=1,R=Ms;
        while(L<R)
        {
            M=(L+R)>>1;
            ll tmp=M-cul();
            if(tmp<k)L=M+1;
            else if(tmp==k)R=M;
            else R=M-1;
        }
        //cout<<ff<<endl;
        printf("Case %d: %lld\n",no,L);
    }
    return 0;
}

//1 223092870 58642669 1000000000

若有什么错误，欢迎指正^ _ ^ 。