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);
}
原文地址：https://www.cnblogs.com/EdSheeran/p/8453013.html
时间： 2024-10-05 22:41:03