NEUOJ 1660 (容斥＋矩阵)

1660: Sequence

时间限制: 5 Sec 内存限制: 128 MB

提交: 13 解决: 3

[提交][状态][讨论版]

题目描述

Easy question once more!
define f(n) = f(n-1) + f(n-2) (n>=3), f(1) = f(2) = 1
you task is calculating g(n) = sigma(f(i)) (1<= i <=n and gcd(i, n) == 1)

输入

Multiple test cases (no more than 100). Each one contains an integer n(1 <= n <= 1e9).Process to end of file.

输出

For each case, output g(n) mod 1e9+7.

样例输入

2
3
4

样例输出

1
2
3

题目可以转化为求前f(1)+f(2)+..+f(n)-Σ f(i)
(gcd (i, n) != 1)．
对于利用容斥原理，就可以发现只需要搞定f(k)+f(2k)+f(3k)+..+f(tk) (t = n/k)．
对于前面部分的f(1)+f(2)+f(3)+...+f(n)也是一样．
这么一个下标等差的斐波那契数列的和可以二分求也可以纯矩阵求，答案就不难求出了．

#include <bits/stdc++.h>
#define pb push_back
using namespace std;
const long long mod = 1e9+7;

struct m {
    long long a[3][3];
    m operator + (m b) const {
        m ans;
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                ans.a[i][j] = a[i][j]+b.a[i][j];
                ans.a[i][j] %= mod;
            }
        }
        return ans;
    }
    m operator * (m b) const {
        m ans;
        memset (ans.a, 0, sizeof ans.a);
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                for (int k = 0; k < 3; k++) {
                    ans.a[i][j] += a[i][k]*b.a[k][j]%mod;
                    ans.a[i][j] %= mod;
                }
            }
        }
        return ans;
    }
    void show () {
        cout << ".........." << endl;
        for (int i = 0; i < 3; i++) {
            for (int  j= 0; j < 3; j++)
                cout << a[i][j] << " "; cout << endl;
        }
        cout << ".........." << endl;
    }
};
long long n;

m qpow (m a, long long b) {
    if (b == 1)
        return a;
    m ans = qpow (a, b>>1);
    ans = ans*ans;
    if (b&1)
        ans = ans*a;
    return ans;
}

long long cal (long long k) { //计算f(k)+f(2k)+f(3k)+..+f(tk),tk<=n
    long long xx[3][3] = {{0,1,0}, {1,1,0}, {0,0,1}}, yy[3][3] = {{1,0,1}, {0,1,0}, {0,0,1}};
    m x;
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++)
            x.a[i][j] = xx[i][j];
    }
    m y;
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++)
            y.a[i][j] = yy[i][j];
    }
    long long tot = n/k;
    m cur = qpow (x, k); cur = cur*y;
    m ans = qpow (cur, tot);
    return ans.a[1][2]%mod;
}

vector <long long> fac;
void get_fac () { //分解n
    fac.clear ();
    long long num = n;
    for (long long i = 2; i*i <= num; i++) {
        if (num%i == 0) {
            fac.pb (i);
            while (num%i == 0)
                num /= i;
        }
    }
    if (num > 1)
        fac.pb (num);
    return ;
}

int num_of_1 (long long num, long long &mul, long long Max) { //统计１的个数
    mul = 1;
    int ans = 0;
    for (long long bit = 1, i = 0; bit <= Max; bit <<= 1, i++) {
        if (num&bit) {
            ans++;
            mul *= fac[i];
        }
    }
    return ans;
}

void work () {
    long long ans = cal (1);
    long long bit = fac.size ();
    long long Max = (1<<bit)-1;
    for (long long i = 1; i <= Max; i++) {
        long long mul = 0; //这个数字对应的素因子乘积
        int cnt = num_of_1 (i, mul, Max); //需要加上还是减去
        if (cnt&1) { //加上
            ans -= cal (mul);
            while (ans < 0)
                ans += mod;
            ans %= mod;
        }
        else { //减去
            ans += cal (mul);
            ans %= mod;
        }
    }
    printf ("%lld\n", ans);
}

int main () {
    while (scanf ("%lld", &n) == 1) {
        get_fac ();
        work ();
    }
    return 0;
}