OI题解 - 牛客2019多校赛 Big Integer

多校赛 = 被神犇队友带飞 ⊙ 被更仙的神犇踩爆……
然后这是一道挺有意思的数学题。
Tab. 好久没更博客了……

---

· 题目

（手动翻译，如有偏差，请在邮箱内告知~勿喷）

定义函数 $f(x)=\sum_{i=0}^{x-1}10^i$ ，例如 $f(1)=1,f(3)=111,f(4)=1111….$ 。给定 $n,m,p$ ，保证 $p$ 为质数。求有多少对 $(a,b)$ 满足：

• $a,b$ 均为整数；
• $a\in[1,n],b\in[1,m]$；
• $p|f(a^b)$ （$x|y$ 表示 $y$ 被 $x$ 整除）；

---

· 解析

- Part 1. 转换同余方程

根据 $f(x)$ 的定义式，我们不难看出那是一个等比数列，于是利用等比数列的公式（或者手推一下也可以）可以得到 $f(x)=\frac{10^x-1}9$ 。

再把题目中的 “$p|f(a^b)$” 换一下 —— 即 “$f(a^b)\equiv 0\pmod p$”，把 $f(x)=\frac{10^x-1}9$ 代入该式可得：

$$\frac{10^{a^b}-1}9\equiv 0\pmod p$$

这个时候我们发现这个分数看起来很烦……就想到把分母 $9$ 消去，也就是两边同时乘 $9$；
但是不能直接乘，对于同余式，只有当 $c,p$ 互质时，存在 $a\equiv b\Leftrightarrow ac\equiv bc$ ；
所以说必须要让 $9$ 与 $p$ 互质……$p$ 是质数，则只有 $p=3$ 时不能直接将 $9$ 消去 —— 那么我们特判 $p=3$；【第一个特判的质数】

于是就可以愉快地消去 $9$ 了，所以式子变为：

$$\begin{align} 10^{a^b}-1&\equiv 0\pmod p\\ 10^{a^b}&\equiv 1\pmod p \end{align}$$

根据上面这个式子，显然得到当 $p=2$ 或 $p=5$ 时，不可能存在 $(a,b)$ 满足条件 —— 因为此时 $10^{a^b}\bmod p=0$ 。那么 $p=2$ 和 $p=5$ 特判一下，答案为 $0$；【第二个和第三个特判及其对应的答案】

- Part 2. 数论阶

排除掉以上的需要特判的情况，我们的目标就是求解 $10^{a^b}\equiv 1\pmod p$ 。那么这里引入一个叫作“数论阶”的东西 ——

对于 $\gcd(a,n)=1$ 的整数，满足 $a^r\equiv1\pmod n $ 的最小整数 $r$，称为 $a$ 模 $n$ 的阶。
（copy from 百度百科，有修改）

可以证明（但是这里就不展开了）：

对于同余方程 $a^x\equiv b\pmod p$ ；若 $a$ 与 $p$ 互质，且 $r$ 是 $a$ 模 $p$ 的阶；
则该方程的解 $x$ 一定是 $r$ 的倍数，且 $r$ 的倍数一定是方程的解 （充分必要性）。

所以对于之前我们求得的同余方程 $10^{a^b}\equiv 1\pmod p$ ，我们求出 $10$ 模 $p$ 的阶 $r$ ，那么 $a^b$ 一定是 $r$ 的倍数，且若 $a^b$ 是 $r$ 的倍数，它一定是方程的一个解。（Tab. 怎么求解 $r$ 一会再说）

- Part 3. 求解答案

这时候我们想到——把 $r,a$ 分解质因数！如果 $a^b$ 是 $r$ 的倍数：

$$\begin{align} &r=\prod_{i=1}^{k}r_i^{p_i}\\ &a=\left(\prod_{i=1}^{k}r_i^{q_i}\right)\cdot a'&&(\gcd\{a',r\}=1)\\ &a^b=\left[\left(\prod_{i=1}^{k}r_i^{q_i}\right)\cdot a'\right]^b=\prod_{i=1}^{k}r_i^{q_ib}\cdot a'^b \end{align}$$

则有： $\forall_i q_ib\ge p_i$ ，即 $\forall_i q_i\ge \left\lceil\frac{p_i}b\right\rceil$；换句话说，令：

$$G=\prod_{i=1}^kr_i^{\left\lceil\frac{p_i}b\right\rceil}$$

则 $G|a$ ，那么对于给定的 $b$，$a$ 的个数就是 $\left\lfloor\frac nG\right\rfloor$ （这个 $n$ 是题目里的 “$a\in[1,n]$”）；所以我们枚举 $b$ —— ？？？$b$ 的范围好像很大没法枚举？？？
emm……其实我们会发现，当 $b\ge \max\{p_i\}$ 的时候，$\left\lceil\frac{p_i}b\right\rceil$ 就都为 $1$ 了——而且我们高兴地发现：$p_i$ 最大就是 $(\log_2n)$ ~
所以我们只需要枚举 $b$ 从 $1$ 到 $\max\{p_i\}$ 就可以了，当 $b$ 超过 $\max\{p_i\}$ 时，$G$ 恒等于 $\prod_{i=1}^kr_i$ ，可以直接计算答案。

- Part 4. 求解数论阶

我用的是 BSGS 直接求解 $10^r\equiv1\pmod p$ 的最小 $r$ ……（感谢 Tiw 神犇）

还有一种方法是利用了数论阶的性质：

$r|\phi(p)$ ；
因为 $p$ 是质数，则 $\phi(p)=p-1$；
所以 $r|p-1$

对 $(p-1)$ 分解因数，从小到大检验是否是数论阶即可。时间复杂度 $O(\sqrt p)$ 。

---

· 源代码

/*Lucky_Glass*/
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
 
typedef long long ll;
const int LOGN=50;
 
int cas;
int n,m,p;
int xp[LOGN],yp[LOGN],np; //p-1= xp[i]^yp[i] (i=1~np)
int d[LOGN],xd[LOGN][LOGN];

/*BSGS*/
//write by Tiw_Air_OAO
typedef pair<int, int> pii;
const int HASHSIZE = 1000037;
vector<pii>h[HASHSIZE];
int hash_find(int k) {
	int x = k%HASHSIZE;
	for(int i=0;i<h[x].size();i++)
		if( h[x][i].first == k )
			return h[x][i].second;
	return -1;
}
void hash_insert(int k, int n) {
	h[k%HASHSIZE].push_back(make_pair(k, n));
}
int  P,SQ;
void hash_clear(int A, int B) {
	for(int i=1,j=1LL*A*B%P;i<=SQ;i++,j=1LL*j*A%P)
		h[j%HASHSIZE].clear();
}
int BSGS(int A, int B) {
	int tmp; SQ = sqrt(P);
	for(int i=1,j=A;i<=SQ;i++,j=1LL*j*A%P) {
		hash_insert(1LL*j*B%P, i);
		if( i == SQ ) tmp = j;
	}
	for(int i=0,j=tmp;i<=P/SQ;i++,j=1LL*j*tmp%P) {
		int x = hash_find(j);
		if( x != -1 && 1LL*(i+1)*SQ - x != 0 ) {
			hash_clear(A, B);
			return 1LL*(i+1)*SQ - x;
		}
	}
	hash_clear(A, B);
	return -1;
}
/*end-BSGS*/ 
int main() {
	scanf("%d",&cas);
	while(cas--) {
		np=0;
		scanf("%d%d%d",&p,&n,&m);
		if(p==2 || p==5) {
			printf("0\n");
			continue;
		}
		if(p==3) {
			printf("%lld\n",1ll*(n/3)*m);
			continue;
		}
		P=p;
		int tmp=BSGS(10,1);
		for(ll i=2; i*i<=tmp; i++)
			if(tmp%i==0) {
				np++;
				xp[np]=i;
				xd[np][0]=1;yp[np]=0;
				while(tmp%i==0)
					tmp/=i,yp[np]++,
					xd[np][yp[np]]=xd[np][yp[np]-1]*i;
			}
		if(tmp) {
			np++;
			xp[np]=tmp;xd[np][0]=1;
			yp[np]=1;xd[np][1]=tmp;
		}
		int b=1;
		ll ans=0;
		int cur;
		while(b<=m) {
			cur=1;
			bool Break=false;
			for(int i=1; i<=np; i++){
				d[i]=(yp[i]+b-1)/b,
				cur*=xd[i][d[i]];
				if(d[i]!=1) Break=true;
			}
			if(!Break) break;
			ans+=n/cur;
			b++;
		}
		ans+=(m-b+1ll)*(n/cur);
		printf("%lld\n",ans);
	}
	return 0;
}

---

The End

Thanks for reading!

Email: lucky_glass@foxmail.com ，欢迎提问~