【转】数论板子
twh大神：
bool issqr(__int128_t x)//开方 { __int128_t y=(__int128_t)ceil(sqrt((long double)x)); for(;y*y<=x;++y); for(--y;y*y>x;--y); return y*y==x; } bool iscub(__int128_t x)//三次方 { __int128_t y=(__int128_t)ceil(pow((long double)x,1.0/3)); for(;y*y*y<=x;++y); for(--y;y*y*y>x;--y); return y*y*y==x; }
int tot,pr[maxx],vis[maxx]; void init()//素数表 O(n) { for(int i=2;i<maxx;i++) { if(!vis[i])pr[tot++]=i; for(int j=0;j<tot&&i*pr[j]<maxx;j++) { vis[i*pr[j]]=1; if(i%pr[j]==0)break; } } }
int d[maxx],pr[maxx],tot; void init() { for(int i = 2; i < maxx; ++i) { if(!d[i]) pr[tot++] = d[i] = i; for(int j = 0, k; (k = i * pr[j]) < maxx; ++j) { d[k] = pr[j]; if(d[i] == pr[j]) break; } } }
const int s=10;//Miller-Rabin inline LL mul_mod(LL a, LL b, LL m) O(n^(1/4)) { LL c = a*b-(LL)((long double)a*b/m+0.5)*m; return c<0 ? c+m : c; } LL fast_exp(LL a,LL x,LL m) { LL b = 1; while (x) { if (x & 1) b = mul_mod(b, a, m); a = mul_mod(a, a, m); x >>= 1; } return b; } bool MR(LL n) { if (!(n&1)) return n == 2; LL t = 0, u; for (u = n-1; !(u&1); u >>= 1) ++t; for (int i = 0; i < s; ++i) { LL a = rand()%(n-2)+2, x = fast_exp(a, u, n); for (LL j = 0, y; x != 1 && j < t; ++j, x = y) { y = mul_mod(x, x, n); if (y == 1 && x != n-1) return false; } if (x != 1) return false; } return true; }
Miller-Rabin Pollard-rho
（输入一个数n，输出小于n且与n互素的整数个数）
（大数欧拉，启发式分解）
#include <cstdio> #include <cstdlib> #include <ctime> #include <algorithm> typedef long long ll; const int s = 10, MAX_F = 70; ll cnt, f[MAX_F]; inline ll mul_mod(ll a, ll b, ll m) { ll c = a*b-(ll)((long double)a*b/m+0.5)*m; return c<0 ? c+m : c; } ll fast_exp(ll a, ll x, ll m) { ll b = 1; while (x) { if (x & 1) b = mul_mod(b, a, m); a = mul_mod(a, a, m); x >>= 1; } return b; } bool MR(ll n) { if (!(n&1)) return n == 2; ll t = 0, u; for (u = n-1; !(u&1); u >>= 1) ++t; for (int i = 0; i < s; ++i) { ll a = rand()%(n-2)+2, x = fast_exp(a, u, n); for (ll j = 0, y; x != 1 && j < t; ++j, x = y) { y = mul_mod(x, x, n); if (y == 1 && x != n-1) return false; } if (x != 1) return false; } return true; } inline ll abs(ll x) { return x<0 ? -x : x; } ll gcd(ll a, ll b) { return b ? gcd(b, a%b) : a; } ll PR(ll n, ll a) { ll x = rand()%n, y = x, k = 1, i = 0, d = 1; while (d == 1) { if ((x = (mul_mod(x, x, n)+a)%n) == y) return n; d = gcd(n, abs(y-x)); if (++i == k) { k <<= 1; y = x; } } return d; } void decomp(ll n) { if (n == 1) return; if (MR(n)) { f[cnt++] = n; return; } ll d = n, c = n-1; while (d == n) d = PR(n, c--); do { n /= d; } while (!(n%d)); decomp(d); decomp(n); } int main() { srand(time(0)); ll n; while (scanf("%lld", &n), n) { cnt = 0; decomp(n); std::sort(f, f+cnt); cnt = std::unique(f, f+cnt)-f; ll ans = n; for (int i = 0; i < cnt; ++i) { // printf("%lld\n", f[i]); ans = ans/f[i]*(f[i]-1); } printf("%lld\n", ans); } return 0; }
int _pow(LL a,LL n) { LL ret=1; while(n) { if(n&1) ret=ret*a%mod; a=a*a%mod; n>>=1; } return ret; } int inv(int x) { return _pow(x,mod-2); }
拓展gcd
void exgcd(ll a,ll b,ll &g,ll &x,ll &y) { if(!b) {g=a;x=1;y=0;} else {exgcd(b,a%b,g,y,x);y-=x*(a/b);} } void work(ll a , ll b , ll d ,ll& g , ll& x , ll& y) { exgcd(a,b,g,x,y); //此处的x，y接收 ax+by=gcd(a,b)的值 x *= d/g; //求ax+by=c 的解x int t = b/g; x = (x%t + t) % t; y = abs( (a*x - d) / b); } ll a,b,g,x,y,d; void solve() { S_3(a,b,d); work(a,b,d,g,x,y); /* if(x==0) { for(int i=1;x<=0;i++) { x=x+b/g*i; y=y-a/g*i; } }*/ }
逆元
void ex_gcd(LL a, LL b, LL &x, LL &y, LL &d) { if (!b) {d = a, x = 1, y = 0;} else{ ex_gcd(b, a % b, y, x, d); y -= x * (a / b); } } LL inv2(LL t, LL p) {//如果不存在，返回-1 LL d, x, y; ex_gcd(t, p, x, y, d); return d == 1 ? (x % p + p) % p : -1; }
错排 排列组合
const int mod=1e9+7; const int maxn=1e4+7; LL D[maxn], C[maxn][105]; int p[1000050]; void init() { D[0]=1;D[1]=0; for(int i=1;i<maxn;i++)//错排 D[i]=((i-1)*(D[i-1]+D[i-2])+mod)%mod; int i; for (p[0]=i=1;i<=200000;i++) p[i]=1ll*p[i-1]*i%mod; for (int i = 0; i < maxn; i++)//排列组合 { for (int j = 0; j <= min(105, i); j++) { if (j == 0) C[i][j] = 1; else C[i][j] = (C[i-1][j] + C[i-1][j-1]) % mod; } } }
C(n,m)
先init 再comb(n,m) n在下
namespace COMB { int F[N<<1], Finv[N<<1], inv[N<<1]; void init() { inv[1] = 1; for (int i = 2; i < N<<1; i++) { inv[i] = (LL)(MOD - MOD/i) * inv[MOD%i] % MOD; } F[0] = Finv[0] = 1; for (int i = 1; i < N<<1; i++) { F[i] = (LL)F[i-1] * i % MOD; Finv[i] = (LL)Finv[i-1] * inv[i] % MOD; } } int comb(int n, int m) { if (m < 0 || m > n) return 0; return (LL)F[n] * Finv[n-m] % MOD * Finv[m] % MOD; } } using namespace COMB;
C（n,m） 不mod
double ln[3100000]; int n,p; void init() { int maxc=3e6; for(int i=1;i<=maxc;i++) { ln[i]=log(i); } } db comb(int x,int y) { db ans=0; for(int i=1;i<=y;i++) { ans+=ln[x-i+1]-ln[i]; } return ans; } db calc(int x) { return x*exp(comb(n,p-1)-comb(n+x,p)); }
C(n,m) mod p (p为素数)
//C(m,n) m在上 n在下,p为素数 LL n,m,p; LL quick_mod(LL a, LL b) { LL ans = 1; a %= p; while(b) { if(b & 1) { ans = ans * a % p; b--; } b >>= 1; a = a * a % p; } return ans; } LL C(LL n, LL m) { if(m > n) return 0; LL ans = 1; for(int i=1; i<=m; i++) { LL a = (n + i - m) % p; LL b = i % p; ans = ans * (a * quick_mod(b, p-2) % p) % p; } return ans; } LL Lucas(LL n, LL m) { if(m == 0) return 1; return C(n % p, m % p) * Lucas(n / p, m / p) % p; } int main() { int T; scan_d(T); while(T--) { scan_d(n),scan_d(m),scan_d(p); print(Lucas(n,m)); } return 0; }
C(n,m) mod p (p可能为合数)
const int N = 200005; bool prime[N]; int p[N]; int cnt; void isprime() { cnt = 0; memset(prime,true,sizeof(prime)); for(int i=2; i<N; i++) { if(prime[i]) { p[cnt++] = i; for(int j=i+i; j<N; j+=i) prime[j] = false; } } } LL quick_mod(LL a,LL b,LL m) { LL ans = 1; a %= m; while(b) { if(b & 1) { ans = ans * a % m; b--; } b >>= 1; a = a * a % m; } return ans; } LL Work(LL n,LL p) { LL ans = 0; while(n) { ans += n / p; n /= p; } return ans; } LL Solve(LL n,LL m,LL P) { LL ans = 1; for(int i=0; i<cnt && p[i]<=n; i++) { LL x = Work(n, p[i]); LL y = Work(n - m, p[i]); LL z = Work(m, p[i]); x -= (y + z); ans *= quick_mod(p[i],x,P); ans %= P; } return ans; } int main() { int T; isprime(); cin>>T; while(T--) { LL n,m,P; cin>>n>>m>>P; n += m - 2; m--; cout<<Solve(n,m,P)<<endl; } return 0; }
FFT：
//china no.1 #pragma comment(linker, "/STACK:1024000000,1024000000") #include <vector> #include <iostream> #include <string> #include <map> #include <stack> #include <cstring> #include <queue> #include <list> #include <stdio.h> #include <set> #include <algorithm> #include <cstdlib> #include <cmath> #include <iomanip> #include <cctype> #include <sstream> #include <functional> #include <stdlib.h> #include <time.h> #include <bitset> #include <complex> using namespace std; #define pi acos(-1) #define PI acos(-1) #define endl '\n' #define srand() srand(time(0)); #define me(x,y) memset(x,y,sizeof(x)); #define foreach(it,a) for(__typeof((a).begin()) it=(a).begin();it!=(a).end();it++) #define close() ios::sync_with_stdio(0); cin.tie(0); #define FOR(x,n,i) for(int i=x;i<=n;i++) #define FOr(x,n,i) for(int i=x;i<n;i++) #define W while #define sgn(x) ((x) < 0 ? -1 : (x) > 0) #define bug printf("***********\n"); #define db double #define ll long long typedef long long LL; const int INF=0x3f3f3f3f; const LL LINF=0x3f3f3f3f3f3f3f3fLL; const int dx[]={-1,0,1,0,1,-1,-1,1}; const int dy[]={0,1,0,-1,-1,1,-1,1}; const int maxn=1e3+10; const int maxx=1<<17; const double EPS=1e-8; const double eps=1e-8; const int mod=1e9+7; template<class T>inline T min(T a,T b,T c) { return min(min(a,b),c);} template<class T>inline T max(T a,T b,T c) { return max(max(a,b),c);} template<class T>inline T min(T a,T b,T c,T d) { return min(min(a,b),min(c,d));} template<class T>inline T max(T a,T b,T c,T d) { return max(max(a,b),max(c,d));} template <class T> inline bool scan_d(T &ret){char c;int sgn;if (c = getchar(), c == EOF){return 0;} while (c != '-' && (c < '0' || c > '9')){c = getchar();}sgn = (c == '-') ? -1 : 1;ret = (c == '-') ? 0 : (c - '0'); while (c = getchar(), c >= '0' && c <= '9'){ret = ret * 10 + (c - '0');}ret *= sgn;return 1;} inline bool scan_lf(double &num){char in;double Dec=0.1;bool IsN=false,IsD=false;in=getchar();if(in==EOF) return false; while(in!='-'&&in!='.'&&(in<'0'||in>'9'))in=getchar();if(in=='-'){IsN=true;num=0;}else if(in=='.'){IsD=true;num=0;} else num=in-'0';if(!IsD){while(in=getchar(),in>='0'&&in<='9'){num*=10;num+=in-'0';}} if(in!='.'){if(IsN) num=-num;return true;}else{while(in=getchar(),in>='0'&&in<='9'){num+=Dec*(in-'0');Dec*=0.1;}} if(IsN) num=-num;return true;} void Out(LL a){if(a < 0) { putchar('-'); a = -a; }if(a >= 10) Out(a / 10);putchar(a % 10 + '0');} void print(LL a){ Out(a),puts("");} //freopen( "in.txt" , "r" , stdin ); //freopen( "data.txt" , "w" , stdout ); //cerr << "run time is " << clock() << endl; typedef complex<long double>Complex; /* * 进行FFT和IFFT前的反转变换。 * 位置i和 （i二进制反转后位置）互换 * len必须去2的幂 */ void reverse(Complex y[],int len) { int i,j,k; for(i = 1, j = len/2;i < len-1; i++) { if(i < j)swap(y[i],y[j]); //交换互为小标反转的元素，i<j保证交换一次 //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的 k = len/2; while( j >= k) { j -= k; k /= 2; } if(j < k) j += k; } } /* * 做FFT * len必须为2^k形式， * on==1时是DFT，on==-1时是IDFT */ void fft(Complex a[],int n,int on) { reverse(a,n); for(int i=1;i<n;i<<=1) { Complex wn=Complex(cos(PI/i),on*sin(PI/i)); for(int j=0;j<n;j+=(i<<1)) { Complex w=Complex(1,0); for(int k=0;k<i;k++,w=w*wn) { Complex x=a[j+k],y=w*a[j+k+i]; a[j+k]=x+y; a[j+k+i]=x-y; } } } if(on==-1) for(int i=0;i<n;i++) a[i].real()/=n; } struct node { int c[maxx]; }A,B,C; Complex a[maxx],b[maxx],c[maxx]; int main() { int n,x; W(scan_d(n)) { FOr(0,n,i) { scan_d(x); x+=20000; A.c[x]++; B.c[x+x]++; C.c[x+x+x]++; } FOR(0,maxx,i) { a[i]=A.c[i],b[i]=B.c[i]; } fft(a,maxx,1); fft(b,maxx,1); FOR(0,maxx,i) { c[i]=a[i]*(a[i]*a[i]-3.0l*b[i]); } fft(c,maxx,-1); FOR(0,maxx,i) { LL ans=((LL)(c[i].real()+0.5)+2*C.c[i])/6; if(ans>0) printf("%d : %lld\n",i-60000,ans); } } }
一个数因子数和：
#include <cmath> #include <cstdio> #include <cstdlib> #include <cstring> #include <iostream> #include <map> #include <algorithm> using namespace std; #define pi acos(-1) #define endl '\n' #define srand() srand(time(0)); #define me(x,y) memset(x,y,sizeof(x)); #define foreach(it,a) for(__typeof((a).begin()) it=(a).begin();it!=(a).end();it++) #define close() ios::sync_with_stdio(0); cin.tie(0); #define FOR(x,n,i) for(int i=x;i<=n;i++) #define FOr(x,n,i) for(int i=x;i<n;i++) #define W while #define sgn(x) ((x) < 0 ? -1 : (x) > 0) #define bug printf("***********\n"); typedef long long LL; const int INF=0x3f3f3f3f; const LL LINF=1e18+7; const int dx[]= {-1,0,1,0,1,-1,-1,1}; const int dy[]= {0,1,0,-1,-1,1,-1,1}; const int maxn=2005; const int maxx=1e5+100; const double EPS=1e-7; //const int mod=998244353; template<class T>inline T min(T a,T b,T c) { return min(min(a,b),c); } template<class T>inline T max(T a,T b,T c) { return max(max(a,b),c); } template<class T>inline T min(T a,T b,T c,T d) { return min(min(a,b),min(c,d)); } template<class T>inline T max(T a,T b,T c,T d) { return max(max(a,b),max(c,d)); } inline LL Scan() { LL Res=0,ch,Flag=0; if((ch=getchar())=='-')Flag=1; else if(ch>='0' && ch<='9')Res=ch-'0'; while((ch=getchar())>='0'&&ch<='9')Res=Res*10+ch-'0'; return Flag ? -Res : Res; } LL a,b,mod=9901; LL ans=1; LL tot,pr[maxx],vis[maxx]; void init()//素数表 O(nsqrt(n)) { for(int i=2;i<maxx;i++) { if(!vis[i])pr[tot++]=i; for(int j=0;j<tot&&i*pr[j]<maxx;j++) { vis[i*pr[j]]=1; if(i%pr[j]==0)break; } } } LL _pow(LL a,LL n) { LL ret=1; while(n) { if(n&1) ret=ret*a%mod; a=a*a%mod; n>>=1; } return ret; } LL inv(LL x) { return _pow(x,mod-2); } void solve() { for(int i=0;i<tot;i++) { int c=0; if(a%pr[i]==0) { while(a%pr[i]==0) { c++; a/=pr[i]; // cout<<c<<" "<<pr[i]<<endl; } if(pr[i]%mod==0) continue; else if(pr[i]%mod==1) { ans=(ans*(c*b+1))%mod; } else { LL ret=(_pow(pr[i],(c*b+1))-1+mod)%mod; ans=(ans*ret*inv(pr[i]-1))%mod; } } } if(a>1) { if(a%mod==0) return ; else if(a%mod==1) { ans=(ans*(b+1))%mod; } else { LL ret=(_pow(a,(b+1))-1+mod)%mod; ans=(ans*ret*inv(a-1))%mod; } } } int main() { init(); while(~scanf("%lld%lld",&a,&b)) { ans=1; solve(); cout<<ans<<endl; } }
baby-step A^x=B(mod c)
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #include<cmath> #define LL __int64 #define N 1000000 using namespace std; struct Node { int idx; int val; } baby[N]; bool cmp(Node n1,Node n2) { return n1.val!=n2.val?n1.val<n2.val:n1.idx<n2.idx; } int gcd(int a,int b) { return b==0?a:gcd(b,a%b); } int extend_gcd(int a,int b,int &x,int &y) { if(b==0) { x=1; y=0; return a; } int gcd=extend_gcd(b,a%b,x,y); int t=x; x=y; y=t-a/b*y; return gcd; } int inval(int a,int b,int n) { int e,x,y; extend_gcd(a,n,x,y); e=((LL)x*b)%n; return e<0?e+n:e; } int PowMod(int a,int b,int MOD) { LL ret=1%MOD,t=a%MOD; while(b) { if(b&1) ret=((LL)ret*t)%MOD; t=((LL)t*t)%MOD; b>>=1; } return (int)ret; } int BinSearch(int num,int m) { int low=0,high=m-1,mid; while(low<=high) { mid=(low+high)>>1; if(baby[mid].val==num) return baby[mid].idx; if(baby[mid].val<num) low=mid+1; else high=mid-1; } return -1; } int BabyStep(int A,int B,int C) { LL tmp,D=1%C; int temp; for(int i=0,tmp=1%C; i<100; i++,tmp=((LL)tmp*A)%C) if(tmp==B) return i; int d=0; while((temp=gcd(A,C))!=1) { if(B%temp) return -1; d++; C/=temp; B/=temp; D=((A/temp)*D)%C; } int m=(int)ceil(sqrt((double)C)); for(int i=0,tmp=1%C; i<=m; i++,tmp=((LL)tmp*A)%C) { baby[i].idx=i; baby[i].val=tmp; } sort(baby,baby+m+1,cmp); int cnt=1; for(int i=1; i<=m; i++) if(baby[i].val!=baby[cnt-1].val) baby[cnt++]=baby[i]; int am=PowMod(A,m,C); for(int i=0; i<=m; i++,D=((LL)(D*am))%C) { int tmp=inval(D,B,C); if(tmp>=0) { int pos=BinSearch(tmp,cnt); if(pos!=-1) return i*m+pos+d; } } return -1; } int main() { int A,B,C; while(scanf("%d%d%d",&A,&C,&B)!=EOF) { if(B>=C) { puts("Orz,I can’t find D!"); continue; } int ans=BabyStep(A,B,C); if(ans==-1) puts("Orz,I can’t find D!"); else printf("%d\n",ans); } return 0; }
莫比乌斯：
int tot,p[maxx],miu[maxx],v[maxx]; void Mobius(int n) { int i,j; for(miu[1]=1,i=2;i<=n;i++) { if(!v[i]) p[tot++]=i,miu[i]=-1; for(j=0;j<tot&&i*p[j]<=n;j++) { v[i*p[j]]=1; if(i%p[j]) miu[i*p[j]]=-miu[i]; else break; } } }
递推式
#include <cstdio> #include <cstring> #include <cmath> #include <algorithm> #include <vector> #include <string> #include <map> #include <set> #include <cassert> using namespace std; #define rep(i,a,n) for (long long i=a;i<n;i++) #define per(i,a,n) for (long long i=n-1;i>=a;i--) #define pb push_back #define mp make_pair #define all(x) (x).begin(),(x).end() #define fi first #define se second #define SZ(x) ((long long)(x).size()) typedef vector<long long> VI; typedef long long ll; typedef pair<long long,long long> PII; const ll mod=1e9+7; ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;} // head long long _,n; namespace linear_seq { const long long N=10010; ll res[N],base[N],_c[N],_md[N]; vector<long long> Md; void mul(ll *a,ll *b,long long k) { rep(i,0,k+k) _c[i]=0; rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod; for (long long i=k+k-1;i>=k;i--) if (_c[i]) rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod; rep(i,0,k) a[i]=_c[i]; } long long solve(ll n,VI a,VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+... // printf("%d\n",SZ(b)); ll ans=0,pnt=0; long long k=SZ(a); assert(SZ(a)==SZ(b)); rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1; Md.clear(); rep(i,0,k) if (_md[i]!=0) Md.push_back(i); rep(i,0,k) res[i]=base[i]=0; res[0]=1; while ((1ll<<pnt)<=n) pnt++; for (long long p=pnt;p>=0;p--) { mul(res,res,k); if ((n>>p)&1) { for (long long i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0; rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod; } } rep(i,0,k) ans=(ans+res[i]*b[i])%mod; if (ans<0) ans+=mod; return ans; } VI BM(VI s) { VI C(1,1),B(1,1); long long L=0,m=1,b=1; rep(n,0,SZ(s)) { ll d=0; rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod; if (d==0) ++m; else if (2*L<=n) { VI T=C; ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; L=n+1-L; B=T; b=d; m=1; } else { ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; ++m; } } return C; } long long gao(VI a,ll n) { VI c=BM(a); c.erase(c.begin()); rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod; return solve(n,c,VI(a.begin(),a.begin()+SZ(c))); } }; int main() { for (scanf("%I64d",&_);_;_--) { scanf("%I64d",&n); printf("%I64d\n",linear_seq::gao(VI{31,197,1255,7997},n-2)); } }
1-n 素数个数
//china no.1 #pragma comment(linker, "/STACK:1024000000,1024000000") #include <vector> #include <iostream> #include <string> #include <map> #include <stack> #include <cstring> #include <queue> #include <list> #include <stdio.h> #include <set> #include <algorithm> #include <cstdlib> #include <cmath> #include <iomanip> #include <cctype> #include <sstream> #include <functional> #include <stdlib.h> #include <time.h> #include <bitset> #include <complex> using namespace std; //#define pi acos(-1) #define PI acos(-1) #define endl '\n' #define srand() srand(time(0)); #define me(x,y) memset(x,y,sizeof(x)); #define foreach(it,a) for(__typeof((a).begin()) it=(a).begin();it!=(a).end();it++) #define close() ios::sync_with_stdio(0); cin.tie(0); #define FOR(x,n,i) for(int i=x;i<=n;i++) #define FOr(x,n,i) for(int i=x;i<n;i++) #define W while #define sgn(x) ((x) < 0 ? -1 : (x) > 0) #define bug printf("***********\n"); #define db double #define ll long long typedef long long LL; const int INF=0x3f3f3f3f; const LL LINF=0x3f3f3f3f3f3f3f3fLL; const int dx[]={-1,0,1,0,1,-1,-1,1}; const int dy[]={0,1,0,-1,-1,1,-1,1}; const int maxn=1e3+10; const int maxx=1e5+10; const double EPS=1e-8; const double eps=1e-8; const LL mod=1e6+3; template<class T>inline T min(T a,T b,T c) { return min(min(a,b),c);} template<class T>inline T max(T a,T b,T c) { return max(max(a,b),c);} template<class T>inline T min(T a,T b,T c,T d) { return min(min(a,b),min(c,d));} template<class T>inline T max(T a,T b,T c,T d) { return max(max(a,b),max(c,d));} template <class T> inline bool scan_d(T &ret){char c;int sgn;if (c = getchar(), c == EOF){return 0;} while (c != '-' && (c < '0' || c > '9')){c = getchar();}sgn = (c == '-') ? -1 : 1;ret = (c == '-') ? 0 : (c - '0'); while (c = getchar(), c >= '0' && c <= '9'){ret = ret * 10 + (c - '0');}ret *= sgn;return 1;} inline bool scan_lf(double &num){char in;double Dec=0.1;bool IsN=false,IsD=false;in=getchar();if(in==EOF) return false; while(in!='-'&&in!='.'&&(in<'0'||in>'9'))in=getchar();if(in=='-'){IsN=true;num=0;}else if(in=='.'){IsD=true;num=0;} else num=in-'0';if(!IsD){while(in=getchar(),in>='0'&&in<='9'){num*=10;num+=in-'0';}} if(in!='.'){if(IsN) num=-num;return true;}else{while(in=getchar(),in>='0'&&in<='9'){num+=Dec*(in-'0');Dec*=0.1;}} if(IsN) num=-num;return true;} void Out(LL a){if(a < 0) { putchar('-'); a = -a; }if(a >= 10) Out(a / 10);putchar(a % 10 + '0');} void print(LL a){ Out(a),puts("");} //freopen( "in.txt" , "r" , stdin ); //freopen( "data.txt" , "w" , stdout ); //cerr << "run time is " << clock() << endl; const int N = 5e6 + 2; bool np[N]; int prime[N], pi[N]; int getprime() { int cnt = 0; np[0] = np[1] = true; pi[0] = pi[1] = 0; for(int i = 2; i < N; i++) { if(!np[i]) prime[++cnt] = i; pi[i] = cnt; for(int j = 1; j <= cnt && i * prime[j] < N; j++) { np[i * prime[j]] = true; if(i % prime[j] == 0) break; } } return cnt; } const int M = 7; const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17; int phi[PM + 1][M + 1], sz[M + 1]; void init() { getprime(); sz[0] = 1; for(int i = 0; i <= PM; i++) phi[i][0] = i; for(int i = 1; i <= M; i++) { sz[i] = prime[i] * sz[i - 1]; for(int j = 1; j <= PM; j++) phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1]; } } int sqrt2(LL x) { LL r = (LL)sqrt(x - 0.1); while(r * r <= x) ++r; return int(r - 1); } int sqrt3(LL x) { LL r = (LL)cbrt(x - 0.1); while(r * r * r <= x) ++r; return int(r - 1); } LL get_phi(LL x, int s) { if(s == 0) return x; if(s <= M) return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s]; if(x <= prime[s]*prime[s]) return pi[x] - s + 1; if(x <= prime[s]*prime[s]*prime[s] && x < N) { int s2x = pi[sqrt2(x)]; LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2; for(int i = s + 1; i <= s2x; i++) ans += pi[x / prime[i]]; return ans; } return get_phi(x, s - 1) - get_phi(x / prime[s], s - 1); } LL getpi(LL x) { if(x < N) return pi[x]; LL ans = get_phi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1; for(int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; i++) ans -= getpi(x / prime[i]) - i + 1; return ans; } LL lehmer_pi(LL x) { if(x < N) return pi[x]; int a = (int)lehmer_pi(sqrt2(sqrt2(x))); int b = (int)lehmer_pi(sqrt2(x)); int c = (int)lehmer_pi(sqrt3(x)); LL sum = get_phi(x, a) + (LL)(b + a - 2) * (b - a + 1) / 2; for (int i = a + 1; i <= b; i++) { LL w = x / prime[i]; sum -= lehmer_pi(w); if (i > c) continue; LL lim = lehmer_pi(sqrt2(w)); for (int j = i; j <= lim; j++) sum -= lehmer_pi(w / prime[j]) - (j - 1); } return sum; } int main() { init(); LL n; while(~scanf("%I64d",&n)) { print(lehmer_pi(n)); } return 0; }
大数：
#include <iostream> #include <cstring> #include <stdio.h> using namespace std; #define DIGIT 4 //四位隔开,即万进制 #define DEPTH 10000 //万进制 #define MAX 25100 //题目最大位数/4，要不大直接设为最大位数也行 typedef int bignum_t[MAX+1]; int read(bignum_t a,istream&is=cin) { char buf[MAX*DIGIT+1],ch ; int i,j ; memset((void*)a,0,sizeof(bignum_t)); if(!(is>>buf))return 0 ; for(a[0]=strlen(buf),i=a[0]/2-1;i>=0;i--) ch=buf[i],buf[i]=buf[a[0]-1-i],buf[a[0]-1-i]=ch ; for(a[0]=(a[0]+DIGIT-1)/DIGIT,j=strlen(buf);j<a[0]*DIGIT;buf[j++]='0'); for(i=1;i<=a[0];i++) for(a[i]=0,j=0;j<DIGIT;j++) a[i]=a[i]*10+buf[i*DIGIT-1-j]-'0' ; for(;!a[a[0]]&&a[0]>1;a[0]--); return 1 ; } void write(const bignum_t a,ostream&os=cout) { int i,j ; for(os<<a[i=a[0]],i--;i;i--) for(j=DEPTH/10;j;j/=10) os<<a[i]/j%10 ; } int comp(const bignum_t a,const bignum_t b) { int i ; if(a[0]!=b[0]) return a[0]-b[0]; for(i=a[0];i;i--) if(a[i]!=b[i]) return a[i]-b[i]; return 0 ; } int comp(const bignum_t a,const int b) { int c[12]= { 1 } ; for(c[1]=b;c[c[0]]>=DEPTH;c[c[0]+1]=c[c[0]]/DEPTH,c[c[0]]%=DEPTH,c[0]++); return comp(a,c); } int comp(const bignum_t a,const int c,const int d,const bignum_t b) { int i,t=0,O=-DEPTH*2 ; if(b[0]-a[0]<d&&c) return 1 ; for(i=b[0];i>d;i--) { t=t*DEPTH+a[i-d]*c-b[i]; if(t>0)return 1 ; if(t<O)return 0 ; } for(i=d;i;i--) { t=t*DEPTH-b[i]; if(t>0)return 1 ; if(t<O)return 0 ; } return t>0 ; } void add(bignum_t a,const bignum_t b) { int i ; for(i=1;i<=b[0];i++) if((a[i]+=b[i])>=DEPTH) a[i]-=DEPTH,a[i+1]++; if(b[0]>=a[0]) a[0]=b[0]; else for(;a[i]>=DEPTH&&i<a[0];a[i]-=DEPTH,i++,a[i]++); a[0]+=(a[a[0]+1]>0); } void add(bignum_t a,const int b) { int i=1 ; for(a[1]+=b;a[i]>=DEPTH&&i<a[0];a[i+1]+=a[i]/DEPTH,a[i]%=DEPTH,i++); for(;a[a[0]]>=DEPTH;a[a[0]+1]=a[a[0]]/DEPTH,a[a[0]]%=DEPTH,a[0]++); } void sub(bignum_t a,const bignum_t b) { int i ; for(i=1;i<=b[0];i++) if((a[i]-=b[i])<0) a[i+1]--,a[i]+=DEPTH ; for(;a[i]<0;a[i]+=DEPTH,i++,a[i]--); for(;!a[a[0]]&&a[0]>1;a[0]--); } void sub(bignum_t a,const int b) { int i=1 ; for(a[1]-=b;a[i]<0;a[i+1]+=(a[i]-DEPTH+1)/DEPTH,a[i]-=(a[i]-DEPTH+1)/DEPTH*DEPTH,i++); for(;!a[a[0]]&&a[0]>1;a[0]--); } void sub(bignum_t a,const bignum_t b,const int c,const int d) { int i,O=b[0]+d ; for(i=1+d;i<=O;i++) if((a[i]-=b[i-d]*c)<0) a[i+1]+=(a[i]-DEPTH+1)/DEPTH,a[i]-=(a[i]-DEPTH+1)/DEPTH*DEPTH ; for(;a[i]<0;a[i+1]+=(a[i]-DEPTH+1)/DEPTH,a[i]-=(a[i]-DEPTH+1)/DEPTH*DEPTH,i++); for(;!a[a[0]]&&a[0]>1;a[0]--); } void mul(bignum_t c,const bignum_t a,const bignum_t b) { int i,j ; memset((void*)c,0,sizeof(bignum_t)); for(c[0]=a[0]+b[0]-1,i=1;i<=a[0];i++) for(j=1;j<=b[0];j++) if((c[i+j-1]+=a[i]*b[j])>=DEPTH) c[i+j]+=c[i+j-1]/DEPTH,c[i+j-1]%=DEPTH ; for(c[0]+=(c[c[0]+1]>0);!c[c[0]]&&c[0]>1;c[0]--); } void mul(bignum_t a,const int b) { int i ; for(a[1]*=b,i=2;i<=a[0];i++) { a[i]*=b ; if(a[i-1]>=DEPTH) a[i]+=a[i-1]/DEPTH,a[i-1]%=DEPTH ; } for(;a[a[0]]>=DEPTH;a[a[0]+1]=a[a[0]]/DEPTH,a[a[0]]%=DEPTH,a[0]++); for(;!a[a[0]]&&a[0]>1;a[0]--); } void mul(bignum_t b,const bignum_t a,const int c,const int d) { int i ; memset((void*)b,0,sizeof(bignum_t)); for(b[0]=a[0]+d,i=d+1;i<=b[0];i++) if((b[i]+=a[i-d]*c)>=DEPTH) b[i+1]+=b[i]/DEPTH,b[i]%=DEPTH ; for(;b[b[0]+1];b[0]++,b[b[0]+1]=b[b[0]]/DEPTH,b[b[0]]%=DEPTH); for(;!b[b[0]]&&b[0]>1;b[0]--); } void div(bignum_t c,bignum_t a,const bignum_t b) { int h,l,m,i ; memset((void*)c,0,sizeof(bignum_t)); c[0]=(b[0]<a[0]+1)?(a[0]-b[0]+2):1 ; for(i=c[0];i;sub(a,b,c[i]=m,i-1),i--) for(h=DEPTH-1,l=0,m=(h+l+1)>>1;h>l;m=(h+l+1)>>1) if(comp(b,m,i-1,a))h=m-1 ; else l=m ; for(;!c[c[0]]&&c[0]>1;c[0]--); c[0]=c[0]>1?c[0]:1 ; } void div(bignum_t a,const int b,int&c) { int i ; for(c=0,i=a[0];i;c=c*DEPTH+a[i],a[i]=c/b,c%=b,i--); for(;!a[a[0]]&&a[0]>1;a[0]--); } void sqrt(bignum_t b,bignum_t a) { int h,l,m,i ; memset((void*)b,0,sizeof(bignum_t)); for(i=b[0]=(a[0]+1)>>1;i;sub(a,b,m,i-1),b[i]+=m,i--) for(h=DEPTH-1,l=0,b[i]=m=(h+l+1)>>1;h>l;b[i]=m=(h+l+1)>>1) if(comp(b,m,i-1,a))h=m-1 ; else l=m ; for(;!b[b[0]]&&b[0]>1;b[0]--); for(i=1;i<=b[0];b[i++]>>=1); } int digit(const bignum_t a,const int b) { int i,ret ; for(ret=a[(b-1)/DIGIT+1],i=(b-1)%DIGIT;i;ret/=10,i--); return ret%10 ; } #define SGN(x) ((x)>0?1:((x)<0?-1:0)) #define ABS(x) ((x)>0?(x):-(x)) int read(bignum_t a,int&sgn,istream&is=cin) { char str[MAX*DIGIT+2],ch,*buf ; int i,j ; memset((void*)a,0,sizeof(bignum_t)); if(!(is>>str))return 0 ; buf=str,sgn=1 ; if(*buf=='-')sgn=-1,buf++; for(a[0]=strlen(buf),i=a[0]/2-1;i>=0;i--) ch=buf[i],buf[i]=buf[a[0]-1-i],buf[a[0]-1-i]=ch ; for(a[0]=(a[0]+DIGIT-1)/DIGIT,j=strlen(buf);j<a[0]*DIGIT;buf[j++]='0'); for(i=1;i<=a[0];i++) for(a[i]=0,j=0;j<DIGIT;j++) a[i]=a[i]*10+buf[i*DIGIT-1-j]-'0' ; for(;!a[a[0]]&&a[0]>1;a[0]--); if(a[0]==1&&!a[1])sgn=0 ; return 1 ; } struct bignum { bignum_t num ; int sgn ; public : inline bignum() { memset(num,0,sizeof(bignum_t)); num[0]=1 ; sgn=0 ; } inline int operator!() { return num[0]==1&&!num[1]; } inline bignum&operator=(const bignum&a) { memcpy(num,a.num,sizeof(bignum_t)); sgn=a.sgn ; return*this ; } inline bignum&operator=(const int a) { memset(num,0,sizeof(bignum_t)); num[0]=1 ; sgn=SGN (a); add(num,sgn*a); return*this ; } ; inline bignum&operator+=(const bignum&a) { if(sgn==a.sgn)add(num,a.num); else if (sgn&&a.sgn) { int ret=comp(num,a.num); if(ret>0)sub(num,a.num); else if(ret<0) { bignum_t t ; memcpy(t,num,sizeof(bignum_t)); memcpy(num,a.num,sizeof(bignum_t)); sub (num,t); sgn=a.sgn ; } else memset(num,0,sizeof(bignum_t)),num[0]=1,sgn=0 ; } else if(!sgn) memcpy(num,a.num,sizeof(bignum_t)),sgn=a.sgn ; return*this ; } inline bignum&operator+=(const int a) { if(sgn*a>0)add(num,ABS(a)); else if(sgn&&a) { int ret=comp(num,ABS(a)); if(ret>0)sub(num,ABS(a)); else if(ret<0) { bignum_t t ; memcpy(t,num,sizeof(bignum_t)); memset(num,0,sizeof(bignum_t)); num[0]=1 ; add(num,ABS (a)); sgn=-sgn ; sub(num,t); } else memset(num,0,sizeof(bignum_t)),num[0]=1,sgn=0 ; } else if (!sgn)sgn=SGN(a),add(num,ABS(a)); return*this ; } inline bignum operator+(const bignum&a) { bignum ret ; memcpy(ret.num,num,sizeof (bignum_t)); ret.sgn=sgn ; ret+=a ; return ret ; } inline bignum operator+(const int a) { bignum ret ; memcpy(ret.num,num,sizeof (bignum_t)); ret.sgn=sgn ; ret+=a ; return ret ; } inline bignum&operator-=(const bignum&a) { if(sgn*a.sgn<0)add(num,a.num); else if (sgn&&a.sgn) { int ret=comp(num,a.num); if(ret>0)sub(num,a.num); else if(ret<0) { bignum_t t ; memcpy(t,num,sizeof(bignum_t)); memcpy(num,a.num,sizeof(bignum_t)); sub(num,t); sgn=-sgn ; } else memset(num,0,sizeof(bignum_t)),num[0]=1,sgn=0 ; } else if(!sgn)add (num,a.num),sgn=-a.sgn ; return*this ; } inline bignum&operator-=(const int a) { if(sgn*a<0)add(num,ABS(a)); else if(sgn&&a) { int ret=comp(num,ABS(a)); if(ret>0)sub(num,ABS(a)); else if(ret<0) { bignum_t t ; memcpy(t,num,sizeof(bignum_t)); memset(num,0,sizeof(bignum_t)); num[0]=1 ; add(num,ABS(a)); sub(num,t); sgn=-sgn ; } else memset(num,0,sizeof(bignum_t)),num[0]=1,sgn=0 ; } else if (!sgn)sgn=-SGN(a),add(num,ABS(a)); return*this ; } inline bignum operator-(const bignum&a) { bignum ret ; memcpy(ret.num,num,sizeof(bignum_t)); ret.sgn=sgn ; ret-=a ; return ret ; } inline bignum operator-(const int a) { bignum ret ; memcpy(ret.num,num,sizeof(bignum_t)); ret.sgn=sgn ; ret-=a ; return ret ; } inline bignum&operator*=(const bignum&a) { bignum_t t ; mul(t,num,a.num); memcpy(num,t,sizeof(bignum_t)); sgn*=a.sgn ; return*this ; } inline bignum&operator*=(const int a) { mul(num,ABS(a)); sgn*=SGN(a); return*this ; } inline bignum operator*(const bignum&a) { bignum ret ; mul(ret.num,num,a.num); ret.sgn=sgn*a.sgn ; return ret ; } inline bignum operator*(const int a) { bignum ret ; memcpy(ret.num,num,sizeof (bignum_t)); mul(ret.num,ABS(a)); ret.sgn=sgn*SGN(a); return ret ; } inline bignum&operator/=(const bignum&a) { bignum_t t ; div(t,num,a.num); memcpy (num,t,sizeof(bignum_t)); sgn=(num[0]==1&&!num[1])?0:sgn*a.sgn ; return*this ; } inline bignum&operator/=(const int a) { int t ; div(num,ABS(a),t); sgn=(num[0]==1&&!num [1])?0:sgn*SGN(a); return*this ; } inline bignum operator/(const bignum&a) { bignum ret ; bignum_t t ; memcpy(t,num,sizeof(bignum_t)); div(ret.num,t,a.num); ret.sgn=(ret.num[0]==1&&!ret.num[1])?0:sgn*a.sgn ; return ret ; } inline bignum operator/(const int a) { bignum ret ; int t ; memcpy(ret.num,num,sizeof(bignum_t)); div(ret.num,ABS(a),t); ret.sgn=(ret.num[0]==1&&!ret.num[1])?0:sgn*SGN(a); return ret ; } inline bignum&operator%=(const bignum&a) { bignum_t t ; div(t,num,a.num); if(num[0]==1&&!num[1])sgn=0 ; return*this ; } inline int operator%=(const int a) { int t ; div(num,ABS(a),t); memset(num,0,sizeof (bignum_t)); num[0]=1 ; add(num,t); return t ; } inline bignum operator%(const bignum&a) { bignum ret ; bignum_t t ; memcpy(ret.num,num,sizeof(bignum_t)); div(t,ret.num,a.num); ret.sgn=(ret.num[0]==1&&!ret.num [1])?0:sgn ; return ret ; } inline int operator%(const int a) { bignum ret ; int t ; memcpy(ret.num,num,sizeof(bignum_t)); div(ret.num,ABS(a),t); memset(ret.num,0,sizeof(bignum_t)); ret.num[0]=1 ; add(ret.num,t); return t ; } inline bignum&operator++() { *this+=1 ; return*this ; } inline bignum&operator--() { *this-=1 ; return*this ; } ; inline int operator>(const bignum&a) { return sgn>0?(a.sgn>0?comp(num,a.num)>0:1):(sgn<0?(a.sgn<0?comp(num,a.num)<0:0):a.sgn<0); } inline int operator>(const int a) { return sgn>0?(a>0?comp(num,a)>0:1):(sgn<0?(a<0?comp(num,-a)<0:0):a<0); } inline int operator>=(const bignum&a) { return sgn>0?(a.sgn>0?comp(num,a.num)>=0:1):(sgn<0?(a.sgn<0?comp(num,a.num)<=0:0):a.sgn<=0); } inline int operator>=(const int a) { return sgn>0?(a>0?comp(num,a)>=0:1):(sgn<0?(a<0?comp(num,-a)<=0:0):a<=0); } inline int operator<(const bignum&a) { return sgn<0?(a.sgn<0?comp(num,a.num)>0:1):(sgn>0?(a.sgn>0?comp(num,a.num)<0:0):a.sgn>0); } inline int operator<(const int a) { return sgn<0?(a<0?comp(num,-a)>0:1):(sgn>0?(a>0?comp(num,a)<0:0):a>0); } inline int operator<=(const bignum&a) { return sgn<0?(a.sgn<0?comp(num,a.num)>=0:1):(sgn>0?(a.sgn>0?comp(num,a.num)<=0:0):a.sgn>=0); } inline int operator<=(const int a) { return sgn<0?(a<0?comp(num,-a)>=0:1): (sgn>0?(a>0?comp(num,a)<=0:0):a>=0); } inline int operator==(const bignum&a) { return(sgn==a.sgn)?!comp(num,a.num):0 ; } inline int operator==(const int a) { return(sgn*a>=0)?!comp(num,ABS(a)):0 ; } inline int operator!=(const bignum&a) { return(sgn==a.sgn)?comp(num,a.num):1 ; } inline int operator!=(const int a) { return(sgn*a>=0)?comp(num,ABS(a)):1 ; } inline int operator[](const int a) { return digit(num,a); } friend inline istream&operator>>(istream&is,bignum&a) { read(a.num,a.sgn,is); return is ; } friend inline ostream&operator<<(ostream&os,const bignum&a) { if(a.sgn<0) os<<'-' ; write(a.num,os); return os ; } }; int main() { int t; cin>>t; while(t--) { int n; cin>>n; bignum jie,ans; jie+=1; for(int i=1;i<=n;i++) { jie*=i; } for(int i=1;i<=n;i++) { ans=ans+jie/i; } cout<<ans<<".0"<<endl; } }
有些小错误，有空再慢慢整理消化！ORZ原创教主。
It is your time to fight！
