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;
    }