【bzoj 4833】[Lydsy1704月赛]最小公倍佩尔数

Description

令 $(1+sqrt 2)^n=e(n)+sqrt 2cdot f(n)$ ，其中 $e(n),f(n)$ 都是整数，显然有 $(1-sqrt 2)^n=e(n)-sqrt 2cdot f(n)$。令 $g(n)$ 表示 $f(1),f(2),cdots ,f(n)$ 的最小公倍数，给定两个正整数 $n$ 和 $p$ ，其中 $p$ 是质数，并且保证 $f(1),f(2),cdots ,f(n)$ 在模 $p$ 意义下均不为 $0$，请计算 $sum _{i=1}^{n}icdot g(i)$ 模 $p$ 的值。

Input

第一行包含一个正整数 $T$ ，表示有 $T$ 组数据，满足 $Tleq 2^{10}$ 。接下来是测试数据。每组测试数据只占一行，包含两个正整数 $n$ 和 $p$ ，满足 $1leq nleq 10^6,2leq pleq 10^9+7$ 。保证所有测试数据的 $n$ 之和不超过 $3 imes 10^6$ 。

Output

对于每组测试数据，输出一行一个非负整数，表示这组数据的答案。

在开始推导前先观察两个式子：

$$gcd(fib(a),fib(b))=fib(gcd(a,b))$$

$$gcd(x^a-1,x^b-1)=x^{gcd(a,b)}-1$$

形如 $f(n)=acdot f(n-1)+bcdot f(n-2)$ 的式子具有性质 $gcd(f(x),f(y))=f(gcd(x,y))$ 。

而题目中的式子等价于： $f(0)=0,f(1)=1,f(n)=2f(n-1)+f(n-2)$，同样满足这个性质。

（以下集合 $T$ 均满足 $T eq varnothing $）

再由式子：$$lcm(S)=prod_{Tsubset S}gcd(T)^{(-1)^{|T|+1}}$$

可以得到：$$g(n)=prod _{Tsubset 2^{[n]}}f(gcd_{iin T}(i))^{(-1)^{|T|+1}}$$

构造出 $h$ 满足 $f(n)=prod _{d|n}h(d)$

得到式子：$$egin{align*} g(n)&=prod _{Tsubset 2^{[n]}}left ( prod _{d|gcd_{iin T}(i)}h(d) ight )^{(-1)^{|T|+1}}\ &=prod _{d=1}^{n}h(d)^{sum _{Tsubset 2^{~[lfloor frac{n}{d} floor ]~}}(-1)^{|T|+1}} end{align*}$$

又由二项式定理可证：$$sum _{Tsubset 2^{[lfloor frac{n}{d} floor ]}}(-1)^{|T|+1}=-sum _{i=1}^{frac{n}{d}}(-1)^iinom{frac{n}{d}}{i}=1$$

所以 $g(n)=prod _{d=1}^{n}h(d)$

问题解决，时间复杂度 $O(nlogn)$。

#include<cstdio>
#include<algorithm>
#include<cstring>
#define LL long long
using namespace std;
const int N=1e6+5;
int T,n,mod,inv,sum,ans,f[N],h[N];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
int main()
{
    T=read();
    while(T--)
    {
        n=read();mod=read();
        f[0]=0;h[1]=f[1]=1;
        for(int i=2;i<=n;i++)
            h[i]=f[i]=(1ll*f[i-1]*2+f[i-2])%mod;
        for(int i=1;i<=n;i++)
        {
            inv=power(h[i],mod-2);
            for(int j=i+i;j<=n;j+=i)h[j]=1ll*h[j]*inv%mod;
        }
        sum=1;ans=0;
        for(int i=1;i<=n;i++)
            sum=1ll*sum*h[i]%mod,ans=(ans+1ll*sum*i)%mod;
        printf("%d
",ans);
    }
    return 0;
}