結果

問題 No.2979 直角三角形の個数
ユーザー Diana773
提出日時 2024-12-11 00:22:48
言語 C++14
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 1,557 ms / 4,000 ms
コード長 4,920 bytes
コンパイル時間 1,945 ms
コンパイル使用メモリ 179,128 KB
実行使用メモリ 11,648 KB
最終ジャッジ日時 2024-12-11 00:22:56
合計ジャッジ時間 7,809 ms
ジャッジサーバーID
（参考情報） 		judge4 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 26
権限があれば一括ダウンロードができます

ソースコード

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

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
static long long N, B;
struct SmallData {
    long long total_contrib;  
    long long total_pairs;    
    vector<long long> uvals;  
} small_data;
static vector<int> mu; 
// u(m,s)=2(m+s)(2m+s)
inline __int128 calc_u128(long long m, long long s) {
    __int128 x = m+s;
    __int128 y = 2*m + s;
    __int128 val = 2 * x * y;
    return val;
}
inline long long calc_u(long long m, long long s) {
    __int128 val = calc_u128(m,s);
    if(val > (__int128)LLONG_MAX) return LLONG_MAX;
    return (long long)val;
}
long long enumerate_small_store() {
    long long ans = 0;
    long long t=(long long)B/2.0;
    long long limit = (long long)(floor(powl(t,0.5L)))+2;
    vector<long long> vals;
    for (long long m=1; m<=limit; m++) {
        long long low=1,high=1000000000;
        while(low<high) {
            long long mid=(low+high+1)>>1;
            __int128 val = calc_u128(m,mid);
            if(val<=B) low=mid; else high=mid-1;
        }
        long long Smax = low;
        if (Smax < 1) continue;
        for (long long s=1; s<=Smax; s++) {
            long long uv = calc_u(m,s);
            long long contrib = (long long)(N/uv);
            ans += contrib;
            vals.push_back(uv);
        }
    }
    sort(vals.begin(), vals.end());
    small_data.total_contrib = ans;
    small_data.total_pairs = (long long)vals.size();
    small_data.uvals = move(vals);
    return ans;
}
long long count_pairs_U_lessEqB(long long U) {
    if(U<=0) return 0;
    auto it = upper_bound(small_data.uvals.begin(), small_data.uvals.end(), U);
    return (long long)(it - small_data.uvals.begin());
}
long long count_pairs_anyU(long long U) {
    if(U<=B) {
        return count_pairs_U_lessEqB(U);
    }
    long long tt = (long long)U/2.0;
    long long limit = (long long)(floor(powl(tt,0.5L)))+2;
    long long cnt=0;
    for (long long m=1; m<=limit; m++) {
        long long low=1,high=1000000000;
        while(low<high) {
            long long mid=(low+high+1)>>1;
            __int128 val = calc_u128(m,mid);
            if(val<=U) low=mid; else high=mid-1;
        }
        long long Smax = low;
        if(Smax>=1) cnt+=Smax;
    }
    return cnt;
}
void linear_sieve_mu(int maxd) {
    mu.clear(); mu.resize(maxd+1,0);
    vector<int> minf(maxd+1,0), primes;
    mu[1]=1;
    for (int i=2;i<=maxd;i++){
        if(!minf[i]){
            minf[i]=i;
            primes.push_back(i);
            mu[i]=-1;
        }
        for (auto p : primes) {
            if((i64)p*i>maxd) break;
            minf[p*i]=p;
            if(i%p==0) {
                mu[p*i]=0;
                break;
            } else {
                mu[p*i]=-mu[i];
            }
        }
    }
}
inline int mu_of_d(i64 d) {
    if(d<(int)mu.size()) return mu[(int)d];
    // 对于本问题应该不会越界
    // 万一发生，也可以退化为因子分解
    // 这里直接因子分解
    int factors=0;
    i64 x=d;
    for(i64 p=2;p*p<=x;p++){
        if(x%p==0) {
            factors++;
            x/=p;
            if(x%p==0)return 0;
        }
    }
    if(x>1) factors++;
    return (factors%2==0)?1:-1;
}
// G函数
i64 G_func(i64 n){
    if(n<12)return 0LL;
    i64 ans=0;
    long long nd=(long long)n;
    long long V=cbrt(nd);
    
    i64 B2=n/V;
    long long tmp=4*B2+1;
    long long X=((long long)sqrt(tmp)-3)/4;
    //i64 X=(i64)floor(Xd);
    
    for(i64 x=1;x<=X;x++){
        long long Y=((long long)sqrt((x*x+2*B2))-3*x)/2;
        for(i64 y=1;y<=Y;y++){
            ans+=n/(2*(x+y)*(2*x+y));
        }
    }
    
    vector<i64> cnt(V+1);
    for(i64 i=1;i<=V;i++){
        i64 U=n/i;
        long long X2=((long long)sqrt(4*U+1)-3)/4;
        for(i64 x=1;x<=X2;x++){
            cnt[i]+=((long long)sqrt(x*x+2*U)-3*x)/2;
        }
    }
    for(i64 i=1;i<V;i++){
        ans+=i*(cnt[i]-cnt[i+1]);
    }
    return ans;
}
// f'(X)= sum_{odd d≥1} mu(d)*G(X/(d*d))
i64 f_prime(i64 X) {
    if(X<1)return 0;
    i64 limit=(i64)floor(sqrtl((long long)X));
    i64 ans=0;
    for(i64 d=1; d<=limit; d+=2) {
        int mu_val = mu_of_d(d);
        if(mu_val==0) continue;
        i64 div = X/(d*d);
        if(div<1) break;
        i64 g_val=G_func(div);
        ans += mu_val*g_val;
    }
    return ans;
}
// f(N)=f'(N)-f'(N/2)
i64 f_func(i64 N) {
    i64 fpN=f_prime(N);
    i64 fpN2=f_prime(N/2);
    return fpN - fpN2;
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cin>>N;
    {
        i64 maxd=(i64)(floorl(sqrtl((long long)N)))+10;
        if(maxd<1) maxd=1;
        if(maxd>10000000) maxd=10000000; //可根据N大小适当调整上限
        linear_sieve_mu((int)maxd);
    }
    // 枚举u≤B
    // 验证给定值 G(10000)=9696 f(10000)=4858
    // 本程序不额外输出验证，只输出f(N)结果
    i64 result = f_func(N);
    cout<<result<<"\n";
    return 0;
}
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
0