CSP10

T2
我們首先可以看出是線性的
矩陣加速

矩陣乘法不滿足乘法交換律，所以$a\times b $ 不等於 \(b\times a\)，也就是說你想讓\(a\)的一行乘上\(b\)的一列,就把\(a\)放左邊
本題中\(b\)應放左邊

點選檢視程式碼

#include <bits/stdc++.h>
#define speed() ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
// #define ll long long
#define ull unsigned long long
#define lid (rt<<1)
#define rid (rt<<1|1)
// #define endl '\n'
//#define int long long
#define pb push_back
#define int long long
// #pragma comment(linker, “/STACK:512000000,512000000”) 
using namespace std;
#define ll __int128
long long n,mod;
struct Martix
{
	int n,m;
	ll a[4][4];
	void resize(int a,int b)
	{
		n=a;m=b;
	}
	Martix ()
	{
		memset(a,0,sizeof(a));
	}
	void one()
	{
		for(int i=1;i<=3;i++) a[i][i]=1;
	}
	friend Martix operator *(Martix A,Martix B)
	{
		Martix res;
		memset(res.a,0,sizeof(res.a));
		//res.resize(min(B.n,A.n),min(B.m,A.m));
		for(int i=1;i<=3;i++)
		{
			for(int j=1;j<=3;j++)
				for(int k=1;k<=3;k++)
					res.a[i][j]=(res.a[i][j]+A.a[i][k]*B.a[k][j]%mod)%mod;
		}
		return res;
	}
	void T()
	{
		for(int i=1;i<=3;i++)
		{
			for(int j=1;j<=3;j++)
			{
				cout<<(long long)a[i][j]<<" ";
			}
			cout<<endl;
		}
	}
}a,b;
Martix power(Martix a,Martix b,ll c)
{
	// Martix res;
	// res.one();
	while(c)
	{
		if(c&1)a=b*a;
		c>>=1;
		b=b*b;
	}
	return a;
}
signed main()
{
	speed();
	// freopen("in.in","r",stdin);
	// freopen("out.out","w",stdout);
	cin>>n>>mod;
	a.resize(3,3);
	// a.T();
	// cout<<"&&&&&&&&"<<endl;
	b.resize(3,3);
	a.a[2][1]=a.a[3][1]=1;
	b.a[1][2]=b.a[2][2]=b.a[2][3]=b.a[3][3]=1;
	//a.T();
	for(ll k=10;;k*=10)
	{
		//cout<<k<<endl;
		b.a[1][1]=k%mod;
		//cout<<b.a[1][1]<<endl;
		// b.T();
		if(n<k)
		{
			a=power(a,b,n-k/10+1);
			break;
		}
		 //cout<<k-k/10<<endl;
		a=power(a,b,k-k/10);

		//cout<<"======="<<endl;
		//a.T();
	}
	cout<<(long long)a.a[1][1]<<endl;
	return 0;
}

T3簡單題

題目背景

zbw 遇上了一道題，由於他急著去找 qby，所以他想讓你來幫他做。

題目描述

給你 \(n,k\) 求下式的值：

\(\sum\limits_{i=1}^n\sum\limits_{j=1}^n(i+j)^kf(\gcd(i,j))\gcd(i,j)\)

其中 \(\gcd(i,j)\) 表示 \(i,j\) 的最大公約數。

\(f\) 函式定義如下：

如果 \(k\) 有平方因子 \(f(k)=0\)，否則 \(f(k)=1\)。

Update:平方因子定義 如果存在一個整數 \(k(k>1),k^2|n\) 則 \(k\) 是 \(n\) 的一個平方因子

請輸出答案對 \(998244353\) 取模的值。

輸入格式

一行兩個整數 \(n,k\)。

輸出格式

一行一個整數，表示答案對 \(998244353\) 取模後的值。

樣例 #1

樣例輸入 #1

3 3

樣例輸出 #1

1216

樣例 #2

樣例輸入 #2

2 6

樣例輸出 #2

9714

樣例 #3

樣例輸入 #3

18 2

樣例輸出 #3

260108

樣例 #4

樣例輸入 #4

143 1

樣例輸出 #4

7648044

提示

測試點	\(n\)	\(k\)
\(1,2\)	\(\leq10^3\)	\(\leq10^3\)
\(3,4\)	\(\leq2 \times 10^3\)	\(\leq10^{18}\)
\(5 \sim8\)	\(\leq5 \times 10^4\)	\(\leq10^{18}\)
\(9\)	\(\leq 5\times10^6\)	\(=1\)
\(10,11\)	\(\leq 5\times10^6\)	\(=2\)
\(12,13\)	\(\leq 5\times10^6\)	\(\leq10^3\)
\(14 \sim20\)	\(\leq 5\times10^6\)	\(\leq10^{18}\)

對於 \(100\%\) 的資料，\(1 \leq n \leq 5 \times 10^6\)，\(1 \leq k \leq 10^{18}\)。

Update on 2020/3/16:

時間調至 \(1\)s,卡掉了 \(O(n\log k)\),\(O(n\log mod)\) 的做法。

首先發現一個性質:我們設\(gcd(i,j)=d\),函式\(f(x)\)即為\(\mu^2(x)\),即為\(\sum_{i=1}^n\sum_{j=1}^n\mu^2(d)d(i+j)^k\)

\[\sum_{d=1}^n\mu^2(d)d\sum_{i=1}^n\sum_{j=1}^n(i+j)^k[gcd(i,j)=d] \]

\[\sum_{d=1}^n\mu^2(d)d^{k+1}\sum_{i=1}^{\lfloor {\frac{n}{d}} \rfloor}\sum_{j=1}^{\lfloor {\frac{n}{d}} \rfloor}(i+j)^k[gcd(i,j)=1] \]

根據\(\sum_{d|n}\mu (d)=[n=1]\)

\[\sum_{d=1}^n\mu^2(d)d^{k+1}\sum_{i=1}^{\lfloor {\frac{n}{d}} \rfloor}\sum_{j=1}^{\lfloor {\frac{n}{d}} \rfloor}(i+j)^k\sum_{t|gcd(i,j)}\mu (t) \]

\[\sum_{d=1}^n\mu^2(d)d^{k+1}\sum_{t=1}^{\lfloor {\frac{n}{d}} \rfloor}\mu(t)t^k\sum_{i=1}^{\lfloor {\frac{n}{dt}} \rfloor}\sum_{j=1}^{\lfloor {\frac{n}{dt}} \rfloor}(i+j)^k \]

然後我們肯定不能\(O(n^3)\)的去處理，最左邊\(\sum_{d=1}^n\mu^2(d)d^{k+1}\)和\(\sum_{t=1}^{\lfloor {\frac{n}{d}} \rfloor}\mu(t)t^k\sum_{i=1}^{\lfloor {\frac{n}{dt}} \rfloor}\)我們可以\(O(n)\)預處理出來,最後一部分如何做到\(O(n)\)預處理呢?
畫一下圖

^
7|
6|
5|
4|5 6 7 8
3|4 5 6 7
2|4 5 5 6
1|2 3 4 5
.------------------->
  1 2 3 4 5 6 7 8 9

設\(t(x)=x^k,sum(x)=\sum_{i=1}^x i^k\)發現\(t(x)=t(x-1)+2\times(sum(2i-1)-sum(i))+sum(2i)-sum(2i-1)\)
好了剩下的怎麼辦，用數列分塊就可以啦,計算次方可以用擴充套件尤拉定理!

點選檢視程式碼

#include <bits/stdc++.h>
#define speed() ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
#define ll long long
#define ull unsigned long long
#define lid (rt<<1)
#define rid (rt<<1|1)
// #define endl '\n'
//#define int long long
#define pb push_back
// #pragma comment(linker, “/STACK:512000000,512000000”) 
using namespace std;
const int N = 5e6+5,mod=998244353;
ll n,k;int mu[N],pr[N],tot;bool is[N];
inline ll qpow(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}
void getMu()
{
	mu[1]=1;
	for(register int i=2;i<=n;i=-~i)
	{
		if(!is[i])
		{
			pr[++tot]=i;mu[i]=-1;
		}
		for(register int j=1;j<=tot&&i*pr[j]<=n;j=-~j)
		{
			is[i*pr[j]]=1;
			if(i%pr[j]==0)
			{
				mu[i*pr[j]]=0;break;
			}
			mu[i*pr[j]]=-mu[i];
		}
	}
}
ll ud[N],td[N],tt[N];int sp[N*2],p[N*2];
inline ll slfk(ll x)
{
	ll res=0;
	register int l=1,r;
	while(l<=x)
	{
		r=x/(x/l);
		res=(1ll*(td[r]-td[l-1]+mod)%mod*tt[x/l]+res)%mod;
		l=r+1;
	}
	return res;
}
inline ll slfk2(ll x)
{
	ll res=0;
	register ll l=1,r;
	while(l<=x)
	{
		r=x/(x/l);
		ll rres=slfk(x/l);
		res=(1ll*(ud[r]-ud[l-1]+mod)%mod*rres+res)%mod;
		l=r+1;
	}
	return res;
}
signed main()
{
	speed();
	// freopen("in.in","r",stdin);
	// freopen("out.out","w",stdout);
	cin>>n>>k;
	// cout<<"**************8"<<endl;
	getMu();
	p[1]=1;
	for(int i=2;i<=2*n;i=-~i)//O(2e8)
	{
		p[i]=qpow(i,k%(mod-1));
		sp[i]=(sp[i-1]+p[i])%mod;
		// cout<<sp[i]<<endl;
	}
	for(ll d=1;d<=n;d=-~d)
	{
		ud[d]=(ud[d-1]+mu[d]*mu[d]*p[d]*d%mod)%mod;
	}
	for(ll t=1;t<=n;t=-~t)
	{
		td[t]=(td[t-1]+mu[t]*p[t])%mod;
	}
	tt[1]=p[2];
	for(int i=2;i<=n;i=-~i)
	{
		tt[i]=(tt[i-1]+2*(sp[2*i-1]-sp[i])%mod+p[2*i])%mod;
	}
	ll ans=0;
	ans=slfk2(n)%mod;
	cout<<ans;
	return 0;
}
/*
^
7|
6|
5|
4|5 6 7 8
3|4 5 6 7
2|4 5 5 6
1|2 3 4 5
.------------------->
  1 2 3 4 5 6 7 8 9
*/

當然，還沒完,學校題庫怎麼可能過得去呢（跑那麼慢）
預處理\(i^k\)的複雜度為\(O(nlogn)\)我們可以放進線性篩中

點選檢視程式碼

#include <bits/stdc++.h>
#define speed() ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
#define ll long long
#define ull unsigned long long
#define lid (rt<<1)
#define rid (rt<<1|1)
// #define endl '\n'
//#define int long long
#define pb push_back
// #pragma comment(linker, “/STACK:512000000,512000000”) 
using namespace std;
const int N = 5e6+5,mod=998244353;
ll n,k;int mu[N],pr[N],tot;bool is[N*2];
int sp[N*2];
inline ll qpow(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}
void getMu()
{
	mu[1]=1;
	for(register int i=2;i<=2*n;i=-~i)
	{
		if(!is[i])
		{
			pr[++tot]=i;
			if(i<=n)mu[i]=-1;
			sp[i]=qpow(i,k);
		}
		for(register int j=1;j<=tot&&i*pr[j]<=2*n;j=-~j)
		{
			is[i*pr[j]]=1;
			sp[i*pr[j]]=1ll*sp[i]*sp[pr[j]]%mod;
			if(i%pr[j]==0)
			{
				break;
			}
			if(i*pr[j]<=n)mu[i*pr[j]]=-mu[i];
		}
	}
}
ll ud[N],td[N],tt[N];
inline ll slfk(ll x)
{
	ll res=0;
	register int l=1,r;
	while(l<=x)
	{
		r=x/(x/l);
		res=(1ll*(td[r]-td[l-1]+mod)%mod*tt[x/l]+res)%mod;
		l=r+1;
	}
	return res;
}
inline ll slfk2(ll x)
{
	ll res=0;
	register ll l=1,r;
	while(l<=x)
	{
		r=x/(x/l);
		ll rres=slfk(x/l);
		res=(1ll*(ud[r]-ud[l-1]+mod)%mod*rres+res)%mod;
		l=r+1;
	}
	return res;
}
signed main()
{
	speed();
	// freopen("in.in","r",stdin);
	// freopen("out.out","w",stdout);
	cin>>n>>k;
	// cout<<"**************8"<<endl;
	getMu();
	sp[1]=1;
	for(int i=2;i<=2*n;i=-~i)//O(2e8)
	{
		// p[i]=qpow(i,k%(mod-1));
		sp[i]=1ll*(sp[i-1]+sp[i])%mod;
		// cout<<sp[i]<<endl;
	}
	for(ll d=1;d<=n;d=-~d)
	{
		ud[d]=1ll*(ud[d-1]+1ll*mu[d]*mu[d]*(sp[d]-sp[d-1])*d%mod)%mod;
	}
	for(ll t=1;t<=n;t=-~t)
	{
		td[t]=1ll*(td[t-1]+1ll*mu[t]*(sp[t]-sp[t-1]))%mod;
	}
	tt[1]=sp[2]-sp[1];
	for(int i=2;i<=n;i=-~i)
	{
		tt[i]=1ll*(tt[i-1]+2*(sp[2*i-1]-sp[i])%mod+(sp[2*i]-sp[2*i-1])%mod+mod)%mod;
	}
	ll ans=0;
	ans=slfk2(n)%mod;
	cout<<(ans+mod)%mod;
	return 0;
}
/*
^
7|
6|
5|
4|5 6 7 8
3|4 5 6 7
2|4 5 5 6
1|2 3 4 5
.------------------->
  1 2 3 4 5 6 7 8 9
*/

加強版看這裡