BZOJ 2813 奇妙的Fibonacci 线性筛

题目大意：给定i，求斐波那契数列中有多少F[j]是F[i]的约数，以及这些j的平方和

定理：Gcd(F[i],F[j])=F[Gcd(i,j)]

证明见 http://hi.baidu.com/wyl8899/item/b4ed30e6b9f404acce2d4f68

那么当F[j]|F[i]时，必有Gcd(F[j],F[i])=F[j]

则此时F[Gcd(j,i)]=F[j]

若Gcd(j,i)==j，则j|i

若Gcd(j,i)!=j，由于斐波那契数列中相等的两项只有F[1]=F[2]=1，故有i=2k+1,j=2

那么我们只要求出i的约数个数和约数平方和即可，如果i是奇数，约数个数要+1，平方和要+4

然后就是线性筛的问题了- - 熟悉约数和筛法的应该不难YY出这两个东西怎么筛吧- - 我不赘述了吧- -

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 10001000
#define MOD 1000000007
using namespace std;
int prime[1001001],tot;
int cnt[M],a[M],p_a[M];
long long sigma_p_2i[M],square_sum[M];
bool not_prime[M];
//a[n]表示n的最小质因数的次数
//p_a[n]表示n的最小质因数的a[n]次
//sigma_p_2i[n]表示Σ[0<=i<=a]p^2i
long long ans1,ans2;
void Linear_Shaker()
{
	int i,j;
	cnt[1]=1;square_sum[1]=1;
	for(i=2;i<=10000000;i++)
	{
		if(!not_prime[i])
		{
			prime[++tot]=i;
			cnt[i]=2;
			a[i]=1;
			p_a[i]=i;
			square_sum[i]=sigma_p_2i[i]=(long long)i*i+1;
		}
		for(j=1;prime[j]*i<=10000000;j++)
		{
			not_prime[prime[j]*i]=1;
			if(i%prime[j]==0)
			{
				cnt[prime[j]*i]=cnt[i]/(a[i]+1)*(a[i]+2);
				a[prime[j]*i]=a[i]+1;
				p_a[prime[j]*i]=p_a[i]*prime[j];
				sigma_p_2i[prime[j]*i]=sigma_p_2i[i]+(long long)p_a[prime[j]*i]*p_a[prime[j]*i];
				square_sum[prime[j]*i]=square_sum[i]/sigma_p_2i[i]*sigma_p_2i[prime[j]*i];
				break;
			}
			cnt[prime[j]*i]=cnt[i]<<1;
			a[prime[j]*i]=1;
			p_a[prime[j]*i]=prime[j];
			sigma_p_2i[prime[j]*i]=(long long)prime[j]*prime[j]+1;
			square_sum[prime[j]*i]=square_sum[i]*((long long)prime[j]*prime[j]+1);
		}
	}
}
void Check()
{
	int i,j;
	for(i=1;i<=10000000;i++)
	{
		int _cnt=0;long long _square_sum=0;
		for(j=1;j*j<=i;j++)
			if(i%j==0)
			{
				_cnt++,_square_sum+=(long long)j*j;
				if(j*j!=i)
					_cnt++,_square_sum+=(long long)(i/j)*(i/j);
			}
		if(_cnt!=cnt[i]||_square_sum!=square_sum[i])
			printf("%d\n",1/0);
	}
}
int main()
{
	Linear_Shaker();
	//Check();
	int T;
	long long n,a,b,c;
	for(cin>>T>>n>>a>>b>>c;T--;n=(n*a+b)%c+1)
	{
		if(~n&1)
			(ans1+=cnt[n])%=MOD,(ans2+=square_sum[n])%=MOD;
		else
			(ans1+=cnt[n]+1)%=MOD,(ans2+=square_sum[n]+4)%=MOD;
	}
	cout<<ans1<<endl<<ans2<<endl;
	return 0;
}