BZOJ 2813 奇妙的Fibonacci

题目：http://www.lydsy.com/JudgeOnline/problem.php?id=2813

题意：
定义Fibonacci数列为$F_1=1,F_2=1,F_i=F_{i-1}+F_{i-2}(i>=3)$。
对于某一个数$F_i$，求有多少个$F_j$能够整除$F_i$ (i可以等于j)，以及所有满足条件的j的平方之和。询问次数$\leq3*10^6$，$i\leq10^7$。

题解：
由定义不难得到

$$gcd(fib_i, fib_{i + 1}) = gcd(fib_{i - 1}, fib_i) = ... = gcd(fib_1, fib_2) = 1$$ 以及 $$fib_{n + m} = fib_{n - 1} * fib_m + fib_n * fib_{m + 1}$$
所以 $$gcd(fib_{n + m}, fib_n) = gcd(fib_{n - 1} * fib_m, fib_n) \\= gcd(fib_{n - 1}, fib_n) * gcd(fib_m, fib_n) \\= gcd(fib_m, fib_n)$$
所以可以归纳得出 $$gcd(fib_i, fib_j) = fib_{gcd(i, j)}$$
那么$fib_j | fib_i$等价于$j | i$。

对于每个询问直接算因子明显会超时，考虑用线性筛法来做，埃式筛法复杂度为$O(nloglogn)$，欧式筛法复杂度为$O(n)$，这里使用的是欧拉筛法。
记$e[i]$为$i$的最小质因数次数，$d[i]$为$i$除去最小质因数的数，$g[i]$为$i$的因数个数，$f[i]$为$i$的因数平方和。
对于欧式筛中：

• $i \; mod \; prime[j] > 0$，$i*prime[j]$第一次遇见最小质因数，因数个数乘2，$f[i*prime[j]]$为$f[i] * prime[j] ^ 2 + f[i]$（每个因子乘或不乘$prime[j]$）。
• $i \; mod \; prime[j] = 0$，$prime[j]$为$i*prime[j]$的最小质因数但已经遇见过，最小质因数加1，因数个数关于最小质因数部分重算，$d[i*prime[j]]$不变，$f[i*prime[j]]$为$f[i] * prime[j] ^ 2 + f[d[i*prime[j]]]$。
• 注意$fib_1 = fib_2 = 1$，这种情况并未算入某些奇数的答案中。

代码：

#include <cstdio>
typedef long long LL;
const int maxc = 664580, maxv = 10000001, mod = 1000000007;
int cnt, prime[maxc], e[maxv], d[maxv], g[maxv], f[maxv], sumg, sumf;
bool vis[maxv];
int main()
{
    int Q, q, A, B, C;
    scanf("%d%d%d%d%d", &Q, &q, &A, &B, &C);
    A %= C, B %= C;
    g[1] = f[1] = 1;
    for(int i = 2; i <= C; ++i)
    {
        if(!vis[i])
        {
            prime[cnt++] = i;
            e[i] = d[i] = 1;
            g[i] = 2;
            f[i] = ((LL)i * i + 1) % mod;
        }
        for(int j = 0, k; j < cnt && (k = i * prime[j]) <= C; ++j)
        {
            vis[k] = 1;
            if(i % prime[j] == 0)
            {
                e[k] = e[i] + 1;
                g[k] = (g[i] / (e[i] + 1)) * (e[k] + 1);
                d[k] = d[i];
                f[k] = (f[i] * ((LL)prime[j] * prime[j] % mod) + f[d[i]]) % mod;
                break;
            }
            else
            {
                e[k] = 1;
                d[k] = i;
                g[k] = g[i] << 1;
                f[k] = f[i] * (((LL)prime[j] * prime[j] + 1) % mod) % mod;
            }
        }
    }
    while(Q--)
    {
        sumg += g[q] + (q & 1);
        sumf += f[q] + 4 * (q & 1);
        if(sumg >= mod) sumg -= mod;
        if(sumf >= mod) sumf -= mod;
        q = (q * (LL)A + B) % C + 1;
    }
    printf("%d\n%d\n", sumg, sumf);
    return 0;
}