1. 程式人生 > 實用技巧 >BZOJ2813 - 奇妙的Fibonacci(篩法求積性函式)

BZOJ2813 - 奇妙的Fibonacci(篩法求積性函式)

題目

Fibonacci數列是這樣一個數列:F1 = 1, F2 = 1, F3 = 2 . . .Fi = Fi-1 + Fi-2 (當 i >= 3)

對於某一個Fibonacci數項Fi,有多少個Fj能夠整除Fi (i可以等於j);所有j的平方之和是多少。

題解

Fibonacci數列有許多性質,其中一個性質就是

\[f_{\gcd(i,j)}=\gcd(f_i, f_j) \]

因此求 \(f_j | f_i\) 可以推出 \(j | i\)。故問題變為求每個數i的約數個數 和 約數平方和。(但是要額外考慮\(f_2=1\)的情況,因為前面那個條件不是充要條件)

i的約數個數 和 約數平方和 都是積性函式,而積性函式基本都可以用線性篩篩出來。關鍵在於若列舉的質數p有i%p==0(即p為i*p的平方因子)時的處理。

求約數個數

\(f(i)\)為i約數個數,\(num(i)\)為i最小質因數對應的次方數。
\(i = p_1^{c_1}p_2^{c_2}...p_k^{c_k}\),那麼\(f(i) = \prod (c_i+1)\),為積性函式。

若i%p==0,p為i*p的最小質因數,也是i*p的平方因子

\(num(i*p)=num(i)+1\)

\(f(i*p)=\frac{f(i)}{num(i)}\cdot (num(i)+1)\)

求約束平方和

\(f(i)\)為i約數平方和,\(g(i)\)為i最小質因數p的(\(p^{2\cdot 0}+p^{2\cdot 1}+...+p^{2\cdot k}\)

)。
\(i = p_1^{c_1}p_2^{c_2}...p_k^{c_k}\),那麼\(f(i) = \prod (p^{2\cdot 0}+p^{2\cdot 1}+...+p^{2\cdot c_i})\),為積性函式。

若i%p==0,p為i*p的最小質因數,也是i*p的平方因子

\(g(i \cdot p)=g(i) \cdot p^2 + 1\)

\(f(i \cdot p)=\frac{f(i)}{g(i)}\cdot (g(i) \cdot p^2 +1)\)

#include <bits/stdc++.h>

#define endl '\n'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define FILE freopen(".//data_generator//in.txt","r",stdin),freopen("res.txt","w",stdout)
#define FI freopen(".//data_generator//in.txt","r",stdin)
#define FO freopen("res.txt","w",stdout)
#define mp make_pair
#define seteps(N) fixed << setprecision(N) 
typedef long long ll;

using namespace std;
/*-----------------------------------------------------------------*/

ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f

const ll N = 1e7 + 10;
const ll M = 1e9 + 7;
const double eps = 1e-5;

bool vis[N];
ll prime[N];
ll cnt;
ll num[N]; // 最小約數的k
ll g[N];   // 最小約數2k方和
ll f1[N];  // 所有約數的個數
ll f2[N];  // 所有約數的平方和

void init() {
    f1[1] = f2[1] = num[1] = 1;
    for(ll i = 2; i < N; i++) {
        if(!vis[i]) {
            num[i] = 1;
            g[i] = i * i + 1;
            f1[i] = 2;
            f2[i] = g[i];
            vis[i] = 1;
            prime[cnt++] = i;
        }
        for(ll j = 0; j < cnt && prime[j] * i < N; j++) {
            ll p = prime[j];
            vis[p * i] = 1;
            if(i % p == 0) {
                num[p * i] = num[i] + 1;
                f1[p * i] = f1[i] / num[p * i] * (num[p * i] + 1);
                g[p * i] = g[i] * p * p + 1;
                f2[p * i] = f2[i] / g[i] * g[p * i];
                break;
            } else {
                f1[p * i] = 2 * f1[i];
                num[p * i] = 1;
                f2[p * i] = f2[p] * f2[i];
                g[p * i] = p * p + 1;
            }
        }
    }
}

int main() {
    IOS;
    init();
    ll q;
    cin >> q;
    ll tq, a, b, c;
    cin >> tq >> a >> b >> c;
    ll ans1 = 0;
    ll ans2 = 0;
    for(ll i = 1; i <= q; i++) {
        ll p = tq;
        tq = (tq * a + b) % c + 1;
        if(p % 2) {
            ans1 = (ans1 + f1[p] + 1) % M;
            ans2 = (ans2 + f2[p] + 4) % M;
        } else {
            ans1 = (ans1 + f1[p]) % M;
            ans2 = (ans2 + f2[p]) % M;
        }
        ans1 %= M;
        ans2 %= M;
    }
    cout << ans1 << endl;
    cout << ans2 << endl;
}