18寒假的数论

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("
");
}



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
", 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
", 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
", ans);
}

 exgcd:http://www.cnblogs.com/frog112111/archive/2012/08/19/2646012.html

            https://blog.csdn.net/shengtao96/article/details/51234355

CRT :http://yzmduncan.iteye.com/blog/1323599/

catalan:https://blog.csdn.net/doc_sgl/article/details/8880468

欧拉函数：https://blog.csdn.net/sentimental_dog/article/details/52002608

//直接求解欧拉函数  
int euler(int n){ //返回euler(n)   
     int res=n,a=n;  
     for(int i=2;i*i<=a;i++){  
         if(a%i==0){  
             res=res/i*(i-1);//先进行除法是为了防止中间数据的溢出   
             while(a%i==0) a/=i;  
         }  
     }  
     if(a>1) res=res/a*(a-1);  
     return res;  
} 

euler[1]=1;    
     for(int i=2;i<Max;i++)    
       euler[i]=i;    
     for(int i=2;i<Max;i++)    
        if(euler[i]==i)    
           for(int j=i;j<Max;j+=i)    
              euler[j]=euler[j]/i*(i-1);//先进行除法是为了防止中间数据的溢出