4833: [Lydsy2017年4月月赛]最小公倍佩尔数

4833: [Lydsy2017年4月月赛]最小公倍佩尔数

Time Limit: 8 Sec Memory Limit: 128 MB
Submit: 119 Solved: 55
[Submit][Status][Discuss]
Description

令(1+sqrt(2))^n=e(n)+f(n)*sqrt(2),其中e(n),f(n)都是整数，显然有(1-sqrt(2))^n=e(n)-f(n)*sqrt(2)。令g(

n)表示f(1),f(2)…f(n)的最小公倍数，给定两个正整数n和p,其中p是质数，并且保证f(1),f(2)…f(n)在模p意义
下均不为0,请计算sigma(i*g(i)),1<=i<=n.其在模p的值。
Input

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

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

Sample Input

5

1 233

2 233

3 233

4 233

5 233
Sample Output

1

5

35

42

121
HINT

Source

鸣谢Tangjz提供试题

[Submit][Status][Discuss]

这感觉是道非常非常难的数学题。。。
以下为了方便，记$(a,b) = gcd(a,b)，[a,b] = lcm(a,b)$
好像数学上也是这样表示的？
不管怎么样先打个表，
发现$f(1) = 1,f(2) = 2,f(n) = 2f(n-1) + f(n-2)$
既然$f$函数这么有规律，那么$g$呢。。。。好吧死活没有
于是就去翻译题解了。。。。。。题解真难！
先是把递推式的转移矩阵列出来，然后手算几项
发现有规律$f(n+m) = f(n-1)f(m) + f(n)f(m+1)$
然后大力数学归纳一下。。。发现是能证明的。。
于是$(f(n+m),f(n)) = (f(n-1)f(m),f(n))$
由打出来的表不难发现，$f(n-1)$和$f(n)$是互质的
所以$(f(n-1)f(m),f(n)) = (f(m),f(n))$，
显然把互质部分的因子去掉等式仍然成立
考虑辗转相除法的过程，可以证明$(f(n),f(m)) = f((n,m))$
考虑任意两个位置$a,b$，$[f(a),f(b)] = \frac {f(a)f(b)}{(f(a),f(b))} = \frac {f(a)f(b)}{f((a,b))}$
考虑这个式子$([a_1,a_2,\dots,a_n],b)$
可以只关心每类素因子对答案的贡献
于是原式可化为$[(a_1,b),(a_2,b),\dots,(a_n,b)]$
考虑任意三个位置$i,j,k$
$[f(i),f(j),f(k)] = \frac {[f(i),f(j)]f(k)}{([f(i),f(j)],f(k))} = \frac {\frac {f(i)f(j)}{f((i,j))}f(k)}{[f((i,k)),f((j,k))]}$
下面那个$lcm$暴力展开就得到$\frac {f(i)f(j)f(k)f((i,j,k))}{f((i,j))f((i,k))f((k,j))}$
于是乎大胆猜结论？？？(一辈子猜不出来系列)

---

上面是小清新与现实的分割线
令$S_n = \{1,2,\dots,n\}$
$g(n) = \prod_{T\subseteq S}f(gcd_{i \in T}(i))^{(-1)^{|T| + 1}}$
这东西。。。听说可以归纳(这个部分苟蒻亲身写了式子。。然后吐了)
那就假装是可以用归纳法证明的吧！(皆大欢喜:D)
再定义$f(n) = \prod_{d|n}h(d)$
于是原式又化为$\prod_{T \subseteq S}\prod_{d|gcd_{i \in T}(i)}h(d)^{(-1)^{|T| + 1}}$
这个式子。。。题解化到结果的操作姿势太玄学。。完全看不懂
结果竟是$g(n) = \prod_{i=1}^n h(i)$
又$h(n) = \prod_{d|n}f(d)^{u(\frac nd)}$
这个用类似证明莫比乌斯反演的证明方式证明
反正剩下就是大力调和级数筛筛筛。。。
复杂度$O(nlogn)$

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;

const int maxn = 1E6 + 5;

int n,p,T,tot,f[maxn],g[maxn],h[maxn],mu[maxn],pri[maxn];
bool not_pri[maxn];

#define Mul(a,b) (1LL * (a) * (b) % p)
#define Add(a,b) ((a) + (b) < p ? (a) + (b) : (a) + (b) - p)

inline int ksm(int x)
{
    int ret = 1;
    for (int y = p - 2; y; y >>= 1)
    {
        if (y & 1) ret = Mul(ret,x);
        x = Mul(x,x);
    }
    return ret;
}

void Solve()
{
    scanf("%d%d",&n,&p);
    f[1] = h[1] = g[1] = 1;
    for (int i = 2; i <= n; i++)
    {
        f[i] = Add(Mul(2,f[i - 1]),f[i - 2]);
        g[i] = ksm(f[i]); h[i] = 1;
    }
    for (int i = 1; i <= n; i++)
    {
        if (mu[i] == 1)
        {
            for (int j = i,now = 1; j <= n; j += i,now++)
                h[j] = Mul(h[j],f[now]);
        }
        else if (mu[i] == -1)
        {
            for (int j = i,now = 1; j <= n; j += i,now++)
                h[j] = Mul(h[j],g[now]);
        }
    }
    int Ans = 0,tmp = 1;
    for (int i = 1; i <= n; i++)
    {
        tmp = Mul(tmp,h[i]);
        Ans = Add(Ans,Mul(i,tmp));
    }
    cout << Add(Ans,p) << endl;
}

int main()
{
    #ifdef DMC
        freopen("DMC.txt","r",stdin);
    #endif

    mu[1] = 1;
    for (int i = 2; i < maxn; i++)
    {
        if (!not_pri[i])
            pri[++tot] = i,mu[i] = -1;
        for (int j = 1; j <= tot; j++)
        {
            int Nex = i * pri[j];
            if (Nex >= maxn) break;
            not_pri[Nex] = 1;
            if (i % pri[j] == 0)
            {
                mu[Nex] = 0; break;
            }
            mu[Nex] = -mu[i];
        }
    }
    cin >> T; while (T--) Solve();
    return 0;
}