1.容斥原理找1到M中与N互质的数
int primes[33], ptot; vector<pair<int,int> > split(int n) { vector<pair<int,int> > stk; for(int i = 2; i * i <= n; i++) { if(n % i == 0) { stk.push_back(make_pair(i, 0));//一维表示因子，二维表示个数 while(n % i == 0) { n /= i; stk.back().second++; } } } if(n != 1) { stk.push_back(make_pair(n, 1));//特判 } return stk; } void print(int s) { for(int i = 7; i >= 0; i--) printf("%d", (s>>i)&1); printf("\n"); } void subset(int u) { // int u = 0x5; //0x5 = 5 0000 0101 for(int s = u; s; s = (s - 1) & u) { print(s); } print(0); } */ #include <vector> #include <cstdio> #include <algorithm> using namespace std; int countbits(int s) { int rt = 0; while(s) { rt++; s -= s & -s; } return rt; } int count(int m, int n) { vector<int> stk; for(int i = 2; i * i <= n; i++) { if(n % i == 0) { stk.push_back(i); while(n % i == 0) n /= i; } } if(n != 1) stk.push_back(n); if(stk.size() == 0) return 1; int rt = 0; int u = (1<<stk.size()) - 1;//造全集 for(int s = u; s; s = (s - 1) & u) {//造子集 int mul = 1; for(int t = 0; t < (int)stk.size(); t++) if((s>>t) & 1) mul *= stk[t]; if(countbits(s) & 1) //奇加偶减 rt += m / mul; else rt -= m / mul; } return m - rt; } int main() { while(1) { int m, n; scanf("%d%d", &m, &n); printf("%d\n", count(m,n)); } }
2.矩阵快速幂
主要分清单位矩阵和初始矩阵以及向量之间的关系
#include <bits/stdc++.h> using namespace std; #define ll long long ll a0,a1,a2,b,c,d,e,n; const ll M = 1000000000000000000LL; struct Maxtri{ ll w[4][4]; void unit() {//单位矩阵 for(int i = 0; i < 4; i++) for(int j = 0; j < 4; j++) w[i][j] = (i == j); } ll *operator[](int i){ return w[i]; } }; ll add(ll a,ll b){ ll c = a + b; if(c >= M)c -= M; return c; } ll quick(ll a,ll b){ ll rt = 0; for(;b;b>>=1,a=add(a,a)) if(b & 1)rt=add(a,rt); return rt; } Maxtri operator*(Maxtri l, Maxtri r){ Maxtri c; for(int i = 0; i < 4; i++) for(int j = 0; j < 4; j++){ c[i][j] = 0; for(int k = 0; k < 4; k++) c[i][j] = add( c[i][j] , quick(l[i][k], r[k][j]) ); } return c; } Maxtri mul(ll b, Maxtri a){ Maxtri rt; for(rt.unit(); b; b>>=1, a=a*a) if(b & 1) rt = rt*a; return rt; } int main(){ //freopen("seq.in","r",stdin); //freopen("seq.out","w",stdout); scanf("%I64d%I64d%I64d%I64d%I64d%I64d%I64d%I64d",&a0,&a1,&a2,&b,&c,&d,&e,&n); Maxtri a; for(int i = 0; i < 4; i++) for(int j = 0; j < 4; j++) a[i][j] = 0; a[2][0] = d, a[2][1] = c, a[2][2] = b, a[2][3] = e; a[0][1] = a[1][2] = a[3][3] = 1;//初始矩阵 Maxtri q = mul(n, a); Maxtri p; p[0][0] = a0, p[0][1] = a1, p[0][2] = a2,p[0][3] = 1;//向量 ll ans = 0; for(int i = 0; i < 4; i++) ans = add(ans, quick( q[0][i] , p[i][0] ));//最后点乘向量 stringstream ss; string sans; ss << ans; ss >> sans; for(int i = 1; i < 18 - (int)sans.size(); i++) printf("0"); cout<<sans<<endl; return 0; }
再贴一个板子
const int Mod = 1e9 + 7; struct Matrix { int w[2][2]; void zero() { for(int i = 0; i < 2; i++) for(int j = 0; j < 2; j++) w[i][j] = 0; } void unit() {//构造单位矩阵，不是初始矩阵 for(int i = 0; i < 2; i++) for(int j = 0; j < 2; j++) w[i][j] = (i == j); } int *operator[](int i) {//重载【】 return w[i]; } }; Matrix operator*(const Matrix &r, const Matrix &s) {//重载乘号 Matrix c; for(int i = 0; i < 2; i++) for(int j = 0; j < 2; j++) { c[i][j] = 0; for(int k = 0; k < 2; k++) c[i][j] = (c[i][j] + 1LL * r[i][k] * s[k][j])% Mod; } return c; } Matrix mpow(Matrix a, int b) {//矩阵快速幂 Matrix rt; for(rt.unit(); b; b>>=1,a=a*a) if(b & 1) rt=rt*a; return rt; }
然后是数论代码板子了→
#include <bits/stdc++.h> using namespace std; const int N = 1000100; namespace NT { bool isnot[N]; int mu[N], phi[N]; int primes[N], ptot; //筛素数 void normal_sieve( int n ) { isnot[1] = true; for( register int i = 2; i <= n; i++ ) { if( !isnot[i] ) { for( register int j = i + i; j <= n; j += i ) { isnot[j] = true; } } } } void linear_sieve( int n ) { isnot[1] = true; for( int i = 2; i <= n; i++ ) { if( !isnot[i] ) //isnot是合数 primes[ptot++] = i; for( int t = 0; t < ptot; t++ ) { int j = primes[t] * i; if( j > n ) break; isnot[j] = true; if( i % primes[t] == 0 ) break; } } } void linear_sieve_more( int n ) { isnot[1] = true; mu[1] = 1; phi[1] = 1; for( int i = 2; i <= n; i++ ) { if( !isnot[i] ) { primes[ptot++] = i; mu[i] = -1; phi[i] = i - 1; } for( int t = 0; t < ptot; t++ ) { int j = primes[t] * i; if( j > n ) break; isnot[j] = true; mu[j] = mu[primes[t]] * mu[i]; phi[j] = phi[primes[t]] * phi[i]; if( i % primes[t] == 0 ) { mu[j] = 0; phi[j] = primes[t] * phi[i]; break; } } } } // greatest common divisor long long gcd( long long a, long long b ) { return b == 0 ? a : gcd( b, a % b ); } long long lcm( long long a, long long b ) { return a / gcd(a,b) * b;//先除后乘 } // gcd(a,b) = a * x + b * y long long exgcd( long long a, long long b, long long &x, long long &y ) { if( b == 0 ) { x = 1; y = 0; return a; } else { long long x0, y0; long long cd = exgcd( b, a % b, x0, y0 ); x = y0; y = x0 - (a/b) * y0; return cd; } } int main() { int a, b; while( scanf("%d%d",&a,&b) == 2 ) { long long x, y; long long cd = exgcd(a,b,x,y); printf( "%lld = %d * %lld + %d * %lld\n", cd, a, x, b, y ); } } // ax + by = c bool linear_equation( long long a, long long b, long long c, long long &x, long long &y, long long &xinc, long long &yinc ) { long long d = gcd(a,b); if( c % d != 0 ) return false; a /= d, b /= d, c /= d; exgcd( a, b, x, y ); x *= c; y *= c; xinc = b; yinc = -a; return true; } // lucas‘t theorem // C(n,m) = C(n / p, m / p) * C(n % p, m % p) ( mod p ) when 0 <= m <= n // C(n,m) = 0 ( mod p ) when m > n long long fac[N], vfac[N]; // 求组合数 /* 暴力 long long comb[N][N]; void init(int n) { for(int i = 0; i <= n; i++) { for(int j = 0; j <= i; j++) { if(j == 0 || j == i) comb[i][j] = 1; else comb[i][j] = (comb[i-1][j-1] + comb[i-1][j]) % Mod; } } } */ //公式 long long comb( long long n, long long m, long long p ) { if( m > n ) return 0; return fac[n] * vfac[m] % p * vfac[n-m] % p; } long long lucas( long long n, long long m, long long p ) { if( m > n ) return 0; if( n / p == 0 ) return 1;//key point! return lucas( n / p, m / p, p ) * comb( n % p, m % p, p ) % p; } } long long exgcd( long long a, long long b, long long &x, long long &y ) { if( b == 0 ) { x = 1; y = 0; return a; } else { long long x0, y0; long long cd = exgcd( b, a % b, x0, y0 ); x = y0; y = x0 - (a/b) * y0; return cd; } } // a^(-1) /* way 1 long long inverse( long long a, long long mod ) { // require mod is a prime return mpow(a, mod-2, mod); } */ // way 2 long long inverse( long long a, long long mod ) { long long x, y; exgcd( a, mod, x, y ); return (x % mod + mod) % mod; } // Chinese Remainder Theorem // x = a ( mod m ) long long crt( vector<long long> va, vector<long long> vm ) { int n = (int)va.size(); long long M = 1; for( int t = 0; t < n; t++ ) M *= vm[t]; long long ans = 0; for( int t = 0; t < n; t++ ) { long long ci = M / vm[t]; long long invci = inverse(ci,vm[t]); ans = (ans + ci * invci % M * va[t])% M; } return ans; } int main() { vector<long long> va, vm; va.push_back(2);vm.push_back(3); va.push_back(3);vm.push_back(5); va.push_back(0);vm.push_back(4); long long ans = crt(va, vm); printf("%I64d\n", ans); }
