bzoj 2813 奇妙的Fibonacci （线性筛求约束平方和）

2813: 奇妙的Fibonacci

Time Limit: 20 Sec   Memory Limit: 512 MB
Submit: 369   Solved: 87
[ Submit][ Status][ Discuss]

Description

Fibonacci数列是这样一个数列：

F₁ = 1， F₂ = 1， F₃ = 2 . . .

F_i = F_i-1 + F_i-2 (当 i >= 3)

pty忽然对这个古老的数列产生了浓厚的兴趣，他想知道：对于某一个Fibonacci数F_i，有多少个F_j能够整除F_i (i可以等于j)，他还想知道所有j的平方之和是多少。

Input

第一行一个整数Q,表示Q个询问。

第二行四个整数:Q₁, A, B, C

第i个询问Q_i = (Q_i-1 * A + B) mod C + 1(当i >= 2)

Output

Ai代表第i个询问有多少个F_j能够整除F_Qi。

Bi代表第i个询问所有j的平方之和。

输出包括两行：

第一行是所有的Ai之和。

第二行是所有的Bi之和。

由于答案过大，只需要输出除以1000000007得到的余数即可。

Sample Input

2

2 2 1 8

Sample Output

6

55

HINT

对于100%的数据保证：Q <= 3*10^6，C <= 10^7，A <= 10^7，B <= 10^7，

1 <= Q1<= C

Source

题解：用线性筛在线性时间内求解。之前已经分析过线性筛求约数个数了，那么如何通过线性筛求解约数平方和呢？

1. i%prime[j]!=0
这时prime[j]是A的最小质因数且为1次

分析A的因数的组成，可以知道他由两部分组成，
一部分是i的所有因数 另一部分是i的所有因数各乘上prime[j]得到的 
因此g[A]=g[i]*2  f[A]=f[i]*(prime[j]^2+1)
2. i%prime[j]==0
此时prime[j]是A的最小质因数 为e[i]+1次 首先更新t[i*prime[j]]的值
令B=i/(prime[j]^t[i]) 即B是i除去最小质因数后得到的数
f[i]=f[B]* ( sigma(k=0..t[i])(prime[j]^(2*k) )
f[A]=f[B]*( sigma(k=0..t[i]+1)(prime[j]^(2*k) )
第二个式子可以分解成f[A]=f[B]+f[B]*(sigma(k=1..t[i]+1)  prime[j]^(2*k))

                        =f[B]+f[B]*sigma(k=0..t[i]) prime[j]^(2*k))*prime[j]^2

                       代入f[i]=f[B]* ( sigma(k=0..t[i])(prime[j]^(2*k) )

                      所以f[A]=f[B]+f[i]*prime[j]^2  这样就可以利用线性筛求解了，另外需要一个数组来记录当前数/（最小质因数^它对应的指数）  

这道题刚开始全开的long long ，结果TLE了，把一些不必要的数组改成int 就过了T_T。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long 
using namespace std;
int  d[10000100],t[10000100],p[10000100],g[10000003],prime[10000100],num,n,m,a,b,c;
long long sum,n1,f[10000003],ans;
const int  inf=1e9+7;
void get_d()
{
  p[1]=1; d[1]=1; t[1]=1;  f[1]=1; num=0;
  for (long long i=2;i<=10000000;i++)
   {
     if (!p[i])
      {
        d[i]=2;t[i]=1;  g[i]=1;
        f[i]=(i*i+1)%inf;
        num++; prime[num]=i;
      }
     for (long long j=1;j<=num;j++)
      {
        if (prime[j]*i>10000000) break;
        p[i*prime[j]]=1;
        if (i%prime[j]==0)
         {
            t[i*prime[j]]=t[i]+1;
            d[i*prime[j]]=d[i]/(t[i]+1)*(t[i*prime[j]]+1);
            g[i*prime[j]]=g[i];
            f[i*prime[j]]=(long long)(f[i]*prime[j]*prime[j]%inf+f[g[i]])%inf;
            break;
         }
        else
         {
            d[i*prime[j]]=d[i]*2;
            t[i*prime[j]]=1;
            g[i*prime[j]]=i;
            f[i*prime[j]]=(long long)(f[i]*prime[j]*prime[j]%inf+f[i])%inf;
         }
      }
   }
 
}
int main()
{
  scanf("%d%d%d%d%d",&m,&n,&a,&b,&c);
  get_d();
  n1=n;
  ans=0;
  for (long long i=1;i<=m;i++)
   {
     sum=(sum+d[n1])%inf; ans=(ans+f[n1])%inf;
     if (n1%2!=0)  sum=(sum+1)%inf,ans=(ans+4)%inf;
     n1=(long long)(n1%c*a%c+b%c)%c+1;
   }
  printf("%lld\n",sum);
  printf("%lld\n",ans);
}