[Lydsy1704月赛] 最小公倍佩尔数

4833: [Lydsy1704月赛]最小公倍佩尔数

Time Limit: 8 Sec Memory Limit: 128 MB
Submit: 202 Solved: 99
[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提供试题

一眼看去就是一道毒瘤数论题233333

首先要把 f() 这个数列搞出来，把f(n)的表达式写出来，发现是 ∑ [i%2==1] * C(n,i) * (sqrt(2)^(i-1))。

(打个表就可以发现 f(n) = 2 * f(n-1) + f(n-2) 了嘛 2333)

显然类菲波那切数列(满足 f(n+m) = f(n+1) * f(m) + f(n) * f(m-1)) 都是满足 gcd(f(x) , f(y)) = f( gcd(x,y) )的啊，之后对指数做min_max容斥，就可以得到一坨子关于子集gcd乘乘除除的玩意，看起来还有点麻烦。。。

于是可以再设 f(n) = π h(d) * [d|n] ，然后样把所有 f(gcd) 替换成 π h(),化简一下就可以发现： g(i) = π h(j) * [j<=i] ，于是就可以直接做了QWQ

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=1e6+5;

inline int add(int x,int y,const int ha){ x+=y; return x>=ha?x-ha:x;}
inline void ADD(int &x,int y,const int ha){ x+=y; if(x>=ha) x-=ha;}
inline int mul(int x,int y,const int ha){ return x*(ll)y%ha;}

inline int ksm(int x,int y,const int ha){
	int an=1;
	for(;y;y>>=1,x=mul(x,x,ha)) if(y&1) an=mul(an,x,ha);
	return an;
}

int T,f[maxn],h[maxn],n,ans=0,p,now;

inline void solve(){
	f[1]=1; for(int i=2;i<=n;i++) f[i]=add(add(f[i-1],f[i-1],p),f[i-2],p);
	
	for(int i=1,inv;i<=n;i++){
		h[i]=f[i],inv=ksm(h[i],p-2,p);
		
		for(int j=i*2;j<=n;j+=i) f[j]=mul(f[j],inv,p);
		
		now=mul(now,h[i],p);
		ADD(ans,mul(now,i,p),p);
	}
}

int main(){
	scanf("%d",&T);
	while(T--){
		scanf("%d%d",&n,&p);
	    
		ans=0,now=1,solve();
		
		printf("%d
",ans);
	}
	
	return 0;
}