結果
問題 | No.1322 Totient Bound |
ユーザー |
![]() |
提出日時 | 2021-03-31 15:21:37 |
言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
結果 |
RE
|
実行時間 | - |
コード長 | 3,038 bytes |
コンパイル時間 | 2,196 ms |
コンパイル使用メモリ | 207,584 KB |
最終ジャッジ日時 | 2025-01-20 01:51:37 |
ジャッジサーバーID (参考情報) |
judge5 / judge3 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
sample | AC * 3 |
other | AC * 31 RE * 5 |
ソースコード
#include <bits/stdc++.h>using namespace std;//template#define rep(i,a,b) for(int i=(int)(a);i<(int)(b);i++)#define ALL(v) (v).begin(),(v).end()using ll=long long int;const int inf = 0x3fffffff; const ll INF = 0x1fffffffffffffff; const double eps=1e-12;template<typename T>inline bool chmax(T& a,T b){if(a<b){a=b;return 1;}return 0;}template<typename T>inline bool chmin(T& a,T b){if(a>b){a=b;return 1;}return 0;}//endll mpow(ll x,ll t,ll m){ll res=1;while(t){if(t&1)res=(__int128_t(res)*x)%m;x=(__int128_t(x)*x)%m; t>>=1;} return res;}bool isp(ll n){if(n<2 or (n&1)==0)return (n==2);ll d=n-1; while((d&1)==0)d>>=1;vector<ll> seeds;if(n<(1<<30))seeds={2, 7, 61};else seeds={2, 325, 9375, 28178, 450775, 9780504};for(auto& x:seeds){if(n<=x)break;ll t=d,y=mpow(x,t,n);while(t!=n-1 and y!=1 and y!=n-1){y=(__int128_t(y)*y)%n; t<<=1;}if(y!=n-1 and (t&1)==0)return 0;} return 1;}mt19937 izyi(398);vector<ll> pollard(ll n){if(n<=1)return {};if(isp(n))return {n};if((n&1)==0){vector<ll> v=pollard(n>>1); v.push_back(2);return v;}for(ll x=2,y=2,d;;){ll c=izyi()%(n-2)+2;do{x=(__int128_t(x)*x+c)%n;y=(__int128_t(y)*y+c)%n;y=(__int128_t(y)*y+c)%n;d=__gcd(x-y+n,n);}while(d==1);if(d<n){vector<ll> lb=pollard(d),rb=pollard(n/d);lb.insert(lb.end(),ALL(rb)); return lb;}}}void picalc(const ll& N,vector<int>& lo,vector<ll>& hi){ll M=round(sqrt(N));lo.resize(M+1); hi.resize(M+1);rep(i,1,M+1)lo[i]=i-1,hi[i]=N/i-1;int cnt=0;for(int p=2;p<=M;p++){if(lo[p-1]==lo[p])continue;const ll n=N/p;const ll q=ll(p)*p;const int w=M/p;const ll L=min(M,N/q);for(int i=1;i<=w;i++)hi[i]-=hi[i*p]-cnt;const int t=min<int>(sqrt(n),L);for(int i=w+1;i<=t;i++)hi[i]-=lo[n/i]-cnt;for(int i=L,j=n/L;i>t;j++){int c=lo[j];while(j+1<=M and lo[j+1]==c)j++;c-=cnt;for(int e=max<int>(t,n/(j+1));i>e;i--)hi[i]-=c;}for(int i=M,j=M/p;j>=p;j--){const int c=lo[j]-cnt;for(int e=j*p;i>=e;i--)lo[i]-=c;}cnt++;}}ll N;const int MX=201010;vector<int> lo;vector<ll> hi;vector<ll> ps;ll pi(ll x){if(x<=round(sqrt(N)))return lo[x];else return hi[N/x];}void dfs(int i,ll phi,ll& res){res+=max(0LL,pi(N/phi)+isp(N/phi+1)-(i+1))+1;if(phi*ps[i]*ps[i]<=N)dfs(i,phi*ps[i],res);rep(k,i+1,ps.size()){if(phi*ps[k]*(ps[k]-1)>N)break;dfs(k,phi*(ps[k]-1),res);}}int main(){cin>>N;picalc(N,lo,hi);bitset<MX> isP;rep(p,2,MX)isP[p]=1;for(int p=2;p*p<MX;p++)if(isP[p]){for(int q=p*p;q<MX;q+=p)isP[q]=0;}rep(p,2,MX)if(isP[p])ps.push_back(p);ll res=1+pi(N)+isp(N+1);rep(i,0,ps.size()){if(ps[i]*(ps[i]-1)>N)break;dfs(i,ps[i]-1,res);}cout<<res<<'\n';return 0;}