Newcoder 18 F.Course（数论+矩阵快速幂）

Description

$Aria$正面临算设课程的考试。

设$F(n)=\prod _{i=1}^{n}i!$

对于给定的$n,m$（其中$n$为质数），求$F(n^m)$中质因子$n$的出现次数，即求一个最大的非负整数$e$满足$n^e$整除$F(n^m)$。

由于是算设课的考试，答案当然是对$10^9+7$取模的。

Input

一行，两个正整数$n,m$，其中$n$是质数

$(n,m\le 10^9)$

Output

一行，一个整数，质因子$n$的出现次数，答案对$10^9+7$取模。

Sample Input

233 2333

Sample Output

861127073

Solution

记所给素数为$p$，$S(n)$表示$n!$中素数$p$的幂指数，那么有$S(n)=\lfloor \frac{n}{p}\rfloor +S(\lfloor \frac{n}{p}\rfloor )$

考虑$F(m)=\sum _{i=0}^{p^m-1}S(i)$，答案即为$F(m)+S(p^m)$，而$S(p^m)=p^{m-1}+...+1=\frac{p^m-1}{p-1}$
$$F(m)=\sum _{i=0}^{p^m-1}(\lfloor \frac{i}{p}\rfloor +S(\lfloor \frac{i}{p}\rfloor ))=p\cdot \sum _{i=0}^{p^{m-1}-1}S(i)+G(m)$$
其中
$$G(m)=\sum _{i=0}^{p^m-1}\lfloor \frac{i}{p}\rfloor =\sum _{i=1}^{p^{m-1}-1}p\cdot i=\frac{p^m(p^{m-1}-1)}{2}$$
故有$F(m)=p\cdot F(m-1)+\frac{p^m(p^{m-1}-1)}{2}=p\cdot F(m-1)+\frac{1}{2}\cdot p^m+\frac{1}{2}\cdot p^{2m-1}$，用矩阵快速幂维护$F(m),p^m,p^{2m-1}$的值即可，时间复杂度$O(logm)$

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int INF=0x3f3f3f3f,maxn=100001;
#define mod 1000000007
int add(int x,int y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x; 
}
int mul(int x,int y)
{
	ll z=1ll*x*y;
	return z-z/mod*mod;
}
int Pow(int x,int y)
{
	ll z=1;
	while(y)
	{
		if(y&1)z=mul(z,x);
		x=mul(x,x);
		y>>=1;
	}
	return z;
}
typedef int M[3][3];
void MUL(M &A,M B)
{
	M C;
	memset(C,0,sizeof(C));
	for(int i=0;i<3;i++)
		for(int j=0;j<3;j++)
			for(int k=0;k<3;k++)
				C[i][j]=add(C[i][j],mul(A[i][k],B[k][j]));
	for(int i=0;i<3;i++)
		for(int j=0;j<3;j++)
			A[i][j]=C[i][j];
}
void POW(M &A,int n)
{
	M B;
	for(int i=0;i<3;i++)
		for(int j=0;j<3;j++)
			B[i][j]=(i==j);
	while(n)
	{
		if(n&1)MUL(B,A);
		MUL(A,A);
		n>>=1;
	}
	for(int i=0;i<3;i++)
		for(int j=0;j<3;j++)
			A[i][j]=B[i][j];
}
int main()
{
	int p,m;
	M A;
	while(~scanf("%d%d",&p,&m))
	{
		A[0][0]=p,A[0][1]=(mod+1)/2,A[0][2]=mod-(mod+1)/2;
		A[1][0]=0,A[1][1]=mul(p,p),A[1][2]=0;
		A[2][0]=0,A[2][1]=0,A[2][2]=p;
		POW(A,m);
		int ans=add(mul(A[0][1],p),mul(A[0][2],p));
		ans=add(ans,mul(add(Pow(p,m),mod-1),Pow(p-1,mod-2)));
		printf("%d\n",ans);
	}
	return 0;
}