听闻存在根号 3log 做法被顶礼膜拜所以来写一个根号 2log 做法震撼人心

luogu_gza · 2025-02-11 15:24:43 · 题解

根号 3log 也太菜了吧！！！！！

代码里面本来应该是杜教筛的部分我其实写的是线性筛，只要知道这是能做到 O(n^{\frac 23}+\sqrt n \log^2 n) 的即可。

stern-brocot tree，莫比乌斯反演，类欧几里得算法。

先来问一个问题，怎么求小于 \frac pq 的最简真分数的个数？

不难注意到这就是 \sum_{i=1}^{n}\sum_{j=1}^{n}[\frac ji \leq \frac pq][\gcd(i,j)=1]，做简单的莫比乌斯反演就有答案为 \sum_{d=1}^{n}\mu(d)\sum_{i=1}^{\lfloor \frac nd \rfloor} \lfloor \frac{pi}{q} \rfloor。

后半部分是一个类欧的形式，最终求也确实用类欧求。

然后我们在 stern-brocot tree 上做二分。

先来看怎么通过一个 01 序列（即每一位表示向左还是向右走）获取到一个分数。

有结论 1：如果当前这个分数是由 \frac ab 和 \frac cd 合并而来的，向左走的分数是由 \frac ab 和 \frac{a+c}{b+d} 合并而来，向右走则是由 \frac{a-c}{b-d} 和 \frac bd 合并而来。

结论 2：任意一个 \frac ab 在 stern-brocot tree 上走到它只需要转向 O(\log(\max(a,b)) 次。

我们就有一个很清晰的思路了，每次先判断这一次走的一条长链是走左边还是右边，然后再倍增求这一次走的长链走多远即可。

复杂度不知道是 3log 根号还是 2log 根号，反正考场上随便手搓了几组 n=10^7 的数据都跑了 200ms 以下。

孩子是 3log 的，被干爆了，怎么优化呢？倍增的时候不要每次从 \log(n) 开始倍增，而是先检查最大的 d 满足能跳 2^d 步，然后从 d 往下倍增。

这样就是可爱的 2log 做法了。

#include<bits/stdc++.h>
using namespace std;
namespace gza{
    #define int long long
//  #define R read()
    #define pc putchar
    #define pb push_back
    #define fo(i,a,b) for(int i=a;i<=b;i++)
    #define rep(i,a,b) for(int i=a;i>=b;i--) 
    #define m1(a,b) memset(a,b,sizeof a)
    #define MT int _=read();while(_--)
    namespace IO{
        inline int read()
        {
            int x=0,flag=0;
            char ch=getchar();
            while(!isdigit(ch)){if(ch=='-') flag=1;ch=getchar();}
            while(isdigit(ch)) x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
            if(flag) return -x;
            return x;
        }
        template<typename T> inline void write(T x)
        {
            if(x<0) pc('-'),x=-x;
            if(x>9) write(x/10);
            pc(x%10+'0'); 
        }
    }
    using namespace IO;

    struct mat{
        int a,b,c,d;
        mat(){a=b=c=d=0;}
        mat(int A,int B,int C,int D){a=A,b=B,c=C,d=D;}
        inline mat operator* (const int& A)const
        {
            if(A>0) return mat(a,b,c+a*A,d+b*A);
            return mat(a-A*c,b-A*d,c,d); 
        }
        inline pair<int,int> get()
        {
            return {a+c,b+d};
        }
    }now;
    int l[24],r[24];
    const int N=1e7+10,M=7e5+10;
    int n,k;
    int primes[M],cnt;
    bool vis[N];
    signed mu[N];
    inline void get_primes(int n)
    {
        mu[1]=1;
        fo(i,2,n)
        {
            if(!vis[i]) primes[++cnt]=i,mu[i]=-1;
            for(int j=1;j<=cnt&&i*primes[j]<=n;j++)
            {
                vis[i*primes[j]]=1;
                if(!(i%primes[j])) break;
                mu[i*primes[j]]=-mu[i];
            }
        }
        fo(i,1,n) mu[i]+=mu[i-1];
    }
    __int128 f(__int128 a,__int128 b,__int128 c,__int128 n)
    {
//      write(a),pc(' ');
//      write(b),pc(' ');
//      write(c),pc(' ');
//      write(n),pc(' ');
//      puts("");
        __int128 n1=(n+1),S1=n1*n/2;
        __int128 m=(a*n+b)/c,ac=a/c,bc=b/c,nm=n*m;
        if(m==0) return 0;
        if(!a) return n1*bc;
        if(a>=c&&b>=c) return (f(a%c,b%c,c,n)+S1*ac+n1*bc);
        __int128 lst=f(c,c-b-1,a,m-1);
        return nm-lst;
    }
    int calc(int p,int q)
    {
        int res=0;
        for(int l=1,r;l<=n;l=r+1)
        {
            r=n/(n/l);
            res+=(mu[r]-mu[l-1])*f(p,0,q,n/l); 
        }
//      cout<<p<<' '<<q<<' '<<res<<endl;
        return res;
    }
    void main(){
        get_primes(1e7);
        n=read(),k=read();
        now=mat(0,1,1,1);
        l[0]=1;
        r[0]=-1;
        fo(i,1,23) l[i]=l[i-1]+l[i-1];
        fo(i,1,23) r[i]=r[i-1]+r[i-1];
        int vv=calc(1,2);
        if(vv==k) return void(puts("1 2"));
        int flag=(vv<=k);//0:l 1:r
        while(1)
        {
            if(flag==0)
            {
                int sum=1;
                mat t=now;
                int d=23;
                fo(i,0,23)
                {
                    mat tmp=t*l[i];
                    pair<int,int> P=tmp.get();
                    int num=calc(P.first,P.second);
                    if(num<=k){d=i-1;break;}
                }
                rep(i,d,0)
                {
                    mat tmp=t*l[i];
                    pair<int,int> P=tmp.get();
                    int num=calc(P.first,P.second);
                    if(num>k) sum+=(1<<i),t=tmp;
                }
                rep(i,23,0) if(sum>>i&1) now=now*l[i];
            }
            else
            {
                int sum=1;
                mat t=now;
                int d=23;
                fo(i,0,23)
                {
                    mat tmp=t*r[i];
                    pair<int,int> P=tmp.get();
                    int num=calc(P.first,P.second);
                    if(num>=k){d=i-1;break;}
                }
                rep(i,d,0)
                {
                    mat tmp=t*r[i];
                    pair<int,int> P=tmp.get();
                    int num=calc(P.first,P.second);
                    if(num<k) sum+=(1<<i),t=tmp;
                }
                rep(i,23,0) if(sum>>i&1) now=now*r[i];
            }
            pair<int,int> P=now.get();
            int num=calc(P.first,P.second);
//          cout<<P.first<<' '<<P.second<<' '<<num<<endl;
            if(k==num) break;
            flag^=1;
        }
        pair<int,int> tmp=now.get();
        int x=tmp.first,y=tmp.second;
        write(x),pc(' '),write(y),puts("");
    }
}
signed main(){
//  freopen("fraction.in","r",stdin);
//  freopen("fraction.out","w",stdout); 
    gza::main(); 
}