結果

問題 No.2620 Sieve of Coins
ユーザー Rubikun
提出日時 2024-01-27 00:38:28
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
MLE  
実行時間 -
コード長 7,743 bytes
コンパイル時間 2,329 ms
コンパイル使用メモリ 207,264 KB
最終ジャッジ日時 2025-02-18 23:56:01
ジャッジサーバーID
（参考情報） 		judge3 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample MLE * 1 -- * 4
other -- * 53
権限があれば一括ダウンロードができます

ソースコード

diff #
プレゼンテーションモードにする

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
template<class T>bool chmax(T &a, const T &b) { if (a<b) { a=b; return true; } return false; }
template<class T>bool chmin(T &a, const T &b) { if (b<a) { a=b; return true; } return false; }
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define mp make_pair
#define si(x) int(x.size())
const int mod=998244353,MAX=100000005,INF=1<<30;
//Dirichlet
vector<int> prime;//i番目の素数
bool is_prime[MAX+1];
ll mawaru[MAX+1];
void sieve(int n){
    for(int i=0;i<=n;i++){
        is_prime[i]=true;
        mawaru[i]=1;
    }
    
    is_prime[0]=is_prime[1]=false;
    
    for(int i=2;i<=n;i++){
        if(is_prime[i]){
            prime.push_back(i);
            is_prime[i]=false;
            mawaru[i]=i;
            for(int j=2*i;j<=n;j+=i){
                if(is_prime[j]){
                    is_prime[j]=false;
                    ll now=j;
                    while(now%i==0){
                        mawaru[j]*=i;
                        now/=i;
                    }
                }
            }
        }
    }
}
//K+1ぐらいまで篩う
// https://maspypy.com/dirichlet-%e7%a9%8d%e3%81%a8%e3%80%81%e6%95%b0%e8%ab%96%e9%96%a2%e6%95%b0%e3%81%ae%e7%b4%af%e7%a9%8d%e5%92%8c
template<typename T>
pair<vector<T>,vector<T>> Dirichlet_seki(vector<T> a,vector<T> b,vector<T> A,vector<T> B,ll N,int type){
    ll K,L;
    if(type==0){
        K=pow((double)N/(log2(N)+1),2.0/3);
        L=pow((double)N,1.0/3)*pow((log2(N)+1),2.0/3);
    }
    if(type==1){
        K=pow((double)N/(log2((log2(N)+1))+1),2.0/3);
        L=pow((double)N,1.0/3)*pow(log2((log2(N)+1))+1,2.0/3);
    }
    if(type==2){
        K=pow((double)N,2.0/3);
        L=pow((double)N,1.0/3);
    }
    chmax(K,(ll)(pow(N,1.0/2)));
    K+=2;
    L+=2;
    
    assert(si(a)==K+1);
    assert(si(A)==L+1);
    
    //cout<<N<<" "<<K<<" "<<L<<endl;
    
    vector<T> Asmall(K+1),Bsmall(K+1);
    for(ll i=1;i<=K;i++){
        Asmall[i]=Asmall[i-1]+a[i];
        Bsmall[i]=Bsmall[i-1]+b[i];
    }
    
    vector<T> c(K+1),C(L+1);
    if(type==0){
        for(ll i=1;i<=K;i++){
            for(ll j=1;i*j<=K;j++){
                c[i*j]+=a[i]*b[j];
            }
        }
    }
    if(type==1){
        c=b;
        for(ll p:prime){
            for(ll i=K/p;i>=1;i--){
                ll n=p*i;
                ll q=p,m=i;
                while(1){
                    c[n]+=a[q]*c[m];
                    if(m%p) break;
                    q*=p;
                    m/=p;
                }
            }
        }
    }
    if(type==2){
        c[1]=a[1]*b[1];
        for(ll p:prime){
            ll now=1;
            for(ll k=1;;k++){
                now*=p;
                if(now>K) break;
                ll xx=1,yy=now;
                for(ll s=0;s<=k;s++){
                    c[now]+=a[xx]*b[yy];
                    xx*=p;
                    yy/=p;
                }
            }
        }
        for(ll n=1;n<=K;n++){
            if(mawaru[n]==n) continue;
            c[n]=c[mawaru[n]]*c[n/mawaru[n]];
        }
    }
    
    for(ll i=1;i<=L;i++){
        ll X=N/i;
        ll M=sqrt(X);
        while(M*M>X) M--;
        while((M+1)*(M+1)<=X) M++;
        for(ll ii=1;ii<=M;ii++){
            ll z=X/ii;
            T kake;
            if(z<=K) kake=Bsmall[z];
            else kake=B[N/z];
            kake*=a[ii];
            
            C[i]+=kake;
        }
        for(ll jj=1;jj<=M;jj++){
            ll z=X/jj;
            if(M>=z) continue;
            T kake;
            if(z<=K) kake=Asmall[z];
            else kake=A[N/z];
            if(M<=K) kake-=Asmall[M];
            else kake-=A[N/M];
            kake*=b[jj];
            C[i]+=kake;
        }
    }
    
    return mp(c,C);
}
//type=0 N^(2/3)logN^1/3 , type=1 N^(2/3)loglogN^1/3, type=2 N^(2/3)
//type1のときは a が乗法的とする
template<typename T>
pair<vector<T>,vector<T>> Dirichlet_waru(vector<T> a,vector<T> c,vector<T> A,vector<T> C,ll N,int type){
    ll K,L;
    if(type==0){
        K=pow((double)N/(log2(N)+1),2.0/3);
        L=pow((double)N,1.0/3)*pow((log2(N)+1),2.0/3);
    }
    if(type==1){
        K=pow((double)N/(log2((log2(N)+1))+1),2.0/3);
        L=pow((double)N,1.0/3)*pow(log2((log2(N)+1))+1,2.0/3);
    }
    if(type==2){
        K=pow((double)N,2.0/3);
        L=pow((double)N,1.0/3);
    }
    chmax(K,(ll)(pow(N,1.0/2)));
    K+=2;
    L+=2;
    //cout<<K+1<<" "<<L+1<<endl;
    assert(si(a)==K+1);
    assert(si(A)==L+1);
    
    vector<T> Asmall(K+1),Bsmall(K+1);
    for(ll i=1;i<=K;i++){
        Asmall[i]=Asmall[i-1]+a[i];
        //Bsmall[i]=Bsmall[i-1]+b[i];
    }
    
    vector<T> b(K+1),B(L+1);
    if(type==0){
        for(ll i=1;i<=K;i++){
            b[i]=c[i]/a[1];
            for(ll j=1;i*j<=K;j++){
                c[i*j]-=a[j]*b[i];
            }
        }
    }
    if(type==1){
        b=c;
        for(ll p:prime){
            for(ll i=1;i<=K/p;i++){
                ll n=p*i;
                ll q=p,m=i;
                while(1){
                    b[n]-=a[q]*b[m];
                    if(m%p) break;
                    q*=p;
                    m/=p;
                }
            }
        }
    }
    if(type==2){
        b[1]=c[1]/a[1];
        for(ll p:prime){
            ll now=1;
            for(ll k=1;;k++){
                now*=p;
                if(now>K) break;
                ll xx=1*p,yy=now/p;
                ll sum=c[now];
                for(ll s=1;s<=k;s++){
                    sum-=a[xx]*b[yy];
                    xx*=p;
                    yy/=p;
                }
                b[now]=sum/a[1];
            }
        }
        for(ll n=1;n<=K;n++){
            if(mawaru[n]==n) continue;
            b[n]=b[mawaru[n]]*b[n/mawaru[n]];
        }
    }
    
    for(ll i=1;i<=K;i++){
        Bsmall[i]=Bsmall[i-1]+b[i];
    }
    
    for(ll i=1;i<=L;i++){
        ll X=N/i;
        ll M=sqrt(X);
        while(M*M>X) M--;
        while((M+1)*(M+1)<=X) M++;
        for(ll jj=1;jj<=M;jj++){
            ll z=X/jj;
            if(M>=z) continue;
            T kake;
            if(z<=K) kake=Asmall[z];
            else kake=A[N/z];
            if(M<=K) kake-=Asmall[M];
            else kake-=A[N/M];
            kake*=b[jj];
            C[i]-=kake;
        }
    }
    
    for(ll i=L;i>=1;i--){
        ll X=N/i;
        ll M=sqrt(X);
        while(M*M>X) M--;
        while((M+1)*(M+1)<=X) M++;
        for(ll ii=2;ii<=M;ii++){
            ll z=X/ii;
            T kake;
            if(z<=K) kake=Bsmall[z];
            else kake=B[N/z];
            kake*=a[ii];
            
            C[i]-=kake;
        }
        for(ll ii=1;ii<=1;ii++){
            ll z=X/ii;
            T kake;
            if(ii<=K) kake=a[ii];
            else kake=A[N/ii];
            
            if(z<=K){
                
            }else{
                B[N/z]=C[i]/kake;
            }
        }
    }
    
    return mp(b,B);
}
//type=0 N^(2/3)logN^1/3 , type=1 N^(2/3)loglogN^1/3, type=2 N^(2/3)
int main(){
    
    std::ifstream in("text.txt");
    std::cin.rdbuf(in.rdbuf());
    cin.tie(0);
    ios::sync_with_stdio(false);
    
    sieve(MAX-3);
    
    ll N;cin>>N;
    ll K,L;
    K=pow((double)N,2.0/3);
    L=pow((double)N,1.0/3);
    chmax(K,(ll)(pow(N,1.0/2)));
    K+=2;
    L+=2;
    vector<ll> c(K+1),C(L+1),a(K+1),A(L+1);
    for(ll i=1;i<=K;i++){
        a[i]=1;
    }
    for(ll i=1;i<=L;i++){
        A[i]=N/i;
    }
    for(ll i=1;i<=K;i++){
        ll x=round(floor(sqrtl(i)+1e-8));
        if(x*x==i) c[i]=1;
    }
    for(ll i=1;i<=L;i++){
        C[i]=round(floor(sqrtl(N/i)+1e-8));
    }
    
    auto res=Dirichlet_waru(c,a,C,A,N,2);
    
    cout<<res.se[1]<<endl;
    
}
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
0