Description：

对于某一个 $Fibonacci$ 数 $F_i$ ，求有多少个 $F_j$ 能够整除 $F_i$ 和所有 $j$ 的平方之和。

$1\le n \le 10^7$

Solution：

对于每个斐波那契数，有 $\gcd(f[n],f[n+1])=\gcd(f[n],f[n]+f[n-1])=\gcd(f[n],f[n-1])(辗转相减)=\cdots=\gcd(f[2],f[1])=1$ 。

还有 $f[n+m]=f[n]\times f[m-1]+f[m]\times f[n+1]$ 。

$\gcd(f[n],f[n+m])=\gcd(f[n],f[n]\times f[m-1]+f[m]\times f[n+1])=\gcd(f[n],f[m]\times f[n+1])=\gcd(f[n],f[m])$

$\gcd(f[n],f[m])=\gcd(f[n],f[m-n])$ ，和辗转相减非常像，最后就是 $\gcd(f[0],f[gcd(n,m)])=f[\gcd(n,m)]$

所以 $f[i]|f[j]\Longleftrightarrow i|j(i\ne 2)$ ，因为 $f[2]=1$ 所以要特判，这题实际求的就是 $\sigma_0$ 和 $\sigma_2$ ，线性筛一下就好了。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cstring>

using namespace std;

#define MOD 1000000007

typedef long long ll;

int q;

#define MAXN 10000001

bool isprime[MAXN];

int prime[MAXN],tot = 0;

ll sigma[MAXN],sigma2[MAXN];

int mindiv[MAXN],mindivs[MAXN];

int main()

{

	for(register int i = 2;i < MAXN;++i)isprime[i] = true;

	sigma[1] = sigma2[1] = 1;

	mindiv[1] = mindivs[1] = 1;

	for(register int i = 2;i < MAXN;++i)

	{

		if(isprime[i])

		{

			prime[++tot] = i;

			sigma[i] = 2;

			sigma2[i] = (1 + (ll)i * i) % MOD;

			mindiv[i] = i;mindivs[i] = i;

		}

		for(register int j = 1;j <= tot && i * prime[j] < MAXN;++j)

		{

			register int k = i * prime[j];

			isprime[k] = false;

			mindiv[k] = prime[j];

			if(i % prime[j] != 0)

			{

				sigma[k] = sigma[i] * 2 % MOD;

				mindivs[k] = prime[j];

				sigma2[k] = (ll)sigma2[i] * sigma2[prime[j]] % MOD; 

			}

			else

			{

				sigma[k] = (ll)sigma[i / mindivs[i]] * (sigma[mindivs[i]] + 1) % MOD;

				sigma2[k] = (ll)sigma2[i / mindivs[i]] * (sigma2[mindivs[i]] * prime[j] % MOD * prime[j] % MOD + 1) % MOD;

				mindivs[k] = mindivs[i] * prime[j];

				break;

			}

		}

	}

	scanf("%d",&q);

	register int s,a,b,c;

	scanf("%d%d%d%d",&s,&a,&b,&c);

	register int ans1 = 0,ans2 = 0;

	for(register int i = 1;i <= q;++i)

	{

		ans1 = (ans1 + sigma[s] + (s & 1)) % MOD;ans2 = (ans2 + sigma2[s] + (s & 1) * 4) % MOD;

		s = ((ll)s * a + b) % c + 1;

	}

	cout << ans1 << endl << ans2 << endl;

	return 0;

}