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

Description

令 $(1+\sqrt 2)^n=e(n)+\sqrt 2\cdot f(n)$ ，其中 $e(n),f(n)$ 都是整数，显然有 $(1-\sqrt 2)^n=e(n)-\sqrt 2\cdot 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}i\cdot g(i)$ 模 $p$ 的值。

Input

第一行包含一个正整数 $T$ ，表示有 $T$ 组数据，满足 $T\leq 2^{10}$ 。接下来是测试数据。每组测试数据只占一行，包含两个正整数 $n$ 和 $p$ ，满足 $1\leq n\leq 10^6,2\leq p\leq 10^9+7$ 。保证所有测试数据的 $n$ 之和不超过 $3\times 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)=a\cdot f(n-1)+b\cdot 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\neq \varnothing $）

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

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

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

得到式子：$$\begin{align*} g(n)&=\prod _{T\subset 2^{[n]}}\left ( \prod _{d|gcd_{i\in T}(i)}h(d) \right )^{(-1)^{|T|+1}}\\ &=\prod _{d=1}^{n}h(d)^{\sum _{T\subset 2^{~[\lfloor \frac{n}{d}\rfloor ]~}}(-1)^{|T|+1}} \end{align*}$$

又由二项式定理可证：$$\sum _{T\subset 2^{[\lfloor \frac{n}{d}\rfloor ]}}(-1)^{|T|+1}=-\sum _{i=1}^{\frac{n}{d}}(-1)^i\binom{\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\n",ans);
    }
    return 0;
}