您现在的位置是：首页 > 其它

当前栏目

Lucas定理

定理

2023-09-14 09:06:47 时间

Lucas定理

$p$ 是质数，则
$C_n^m\equiv C_{n\ mod\ p}^{m\ mod\ p}C_{\lfloor \frac{n} {p}\rfloor}^{\lfloor \frac{m}{p}\rfloor}(mod\ p)$
其中 $C_{n}^{m}=0(n<m)$
另一种形式是
表示 $m, n$ 成 $p$ 进制
$n=\sum_{i=0}^{k}n_k p^k$
$m=\sum_{i=0}^{k}m_k p^k$
则
$C_{n}^{m}\equiv (\prod_{i=0}^{k}C_{n_k}^{m_k})(\ mod\ p)$

引理

$p$ 是质数， $0 < f < p$ ,则
$C_{p}^{f}\equiv 0(\ mod\ p)$

证明：
因为 $p$ 是质数，所以 $f!$ 和 $(p - f)!$ 中没有 $p$ 的因子,所以显然成立

证明

设
$n = s p + q (0 < q < p)$
$m = tp + r (0 < r < p)$
$\begin{aligned} (1+x)^{n} &\equiv (1+x)^{sp+q}(\ mod\ p)\\ &\equiv ((1+x)^{p})^{s}(1+x)^{q}(\ mod\ p)\\ &\equiv (\sum_{i=0}^{p}C_{p}^{i}x^i)^{s}(1+x)^{q}(\ mod\ p)\\ &\equiv (1+x^{p})^{s}(1+x)^{q}(\ mod\ p)\\ &\equiv \sum_{i=1}^{s}C_{s}^{i}x^{pi}\sum_{j=1}^{q}C_{q}^{j}x^{j}(\ mod\ p) \end{aligned}$
由二项式定理
$1+x)^{sp+q}$ 中 $x^{tp+r}$ 的系数为 $C_{sp+q}^{tp+r}$
对照 $x$ 的系数
$\begin{cases} tp=pi\\ r=j\\ \end{cases} \Rightarrow \begin{cases} t=i\\ r=j\\ \end{cases}$
$\Rightarrow C_{sp+q}^{tp+r}\equiv C_{s}^{t}C_{q}^{r}(\ mod\ p)$

代码

其中有几个地方可以改的
比如，在lucas那里，先判断取摸后的会不会 $m > n$
还有可以直接预处理阶乘逆元，然后用 $C_{n}^{m}=\frac{n!}{m!(n-m)!}$
反正怎么爽怎么来

#include<iostream>
using namespace std;
typedef long long int ll;
int n, m, p;
//快速幂a^m
ll quickMod(const ll &a,ll m) {
    ll base = a, ans = 1;
    while (m) {
        if (m & 1)
            ans = ans*base%p;
        base = base*base%p;
        m >>= 1;
    }
    return ans;
}
//费马小定理求逆元
ll modReverse(const ll &a, const ll p) {
    return quickMod(a, p - 2);
}
ll c(ll n, ll m) {
    if (m > n)
        return 0;
    ll a = 1, b = 1;
    while (m) {
        a = a*n%p;
        b = b*m%p;
        --n;
        --m;
    }
    return a*modReverse(b, p) % p;
}
//Lucas定理
ll Lucas(ll n, ll m) {
    if (m == 0)
        return 1;
    return c(n%p, m%p)*Lucas(n / p, m / p) % p;
}

扩展lucas

$C_n^m \mathop{mod} p$
其中 $p$ 不一定是质数

首先质因数分解
$p_1^{\alpha_1}p_2^{\alpha_2}\cdots$
求出
$a_i = C_n^m \mathop{mod} p_i^{\alpha_i}$
最后用中国剩余定理合并

现在求 $a_i = C_n^m \mathop{mod} p^k = \frac{n!}{m! \left(n-m\right)!} \mathop{mod} p^k$ 其中 $p$ 是质数
由于 $n!, (n - m)!$ 不一定与 $p$ 互质，所以逆元不一定存在，所以转换一下

$\frac{n!}{m! \left(n-m\right)!} \mathop{mod} p^k = \frac{\frac{n!}{p^x}}{\frac{m!}{p^y} \frac{\left(n-m\right)!}{p^z}} p^{x-y-z} \mathop{mod} p^k$
其中 $x$ 为 $n!$ 中 $p$ 的因数个数， $y$ 为 $m!$ 中 $p$ 的因数个数， $z$ 为 $(n - m)!$ 中 $p$ 的因数个数
那么现在问题就是
$\frac{n!}{p^x}\mathop{mod} p^k$

$\begin{aligned} n! &= 1\cdot 2 \cdot 3 \cdots n \\ &=\left(p\cdot 2p \cdot 3p\cdots\right)\left(1\cdot 2\cdot 3\cdots\right)\\ &=p^{\left\lfloor\frac{n}{p}\right\rfloor}\left(\left\lfloor\frac{n}{p}\right\rfloor!\right)\prod_{i=1 \atop i \not\equiv 0 \left(\mathop{mod} p\right)}^{n}i\\ \end{aligned}$
举个例子
$3^7\left(1\cdot 2\cdots 7\right)\left(1\cdot 2\cdot 4 \cdot 5 \cdot 7 \cdot 8 \cdot 11 \cdot 13 \cdot 14 \cdot 16 \cdot 17 \cdot 19 \cdot 20 \cdot 22\right)$
显然最后一个部分，是有循环节的
$1\cdot 2\cdot 4 \cdot 5 \cdot 7 \cdot 8 \equiv \cdot 11 \cdot 13 \cdot 14 \cdot 16 \cdot 17 \cdot 19 \left(\mathop{mod} 3^2 \right)$
所以
$p^{\left\lfloor\frac{n}{p}\right\rfloor}\left(\left\lfloor\frac{n}{p}\right\rfloor!\right)\left(\prod_{i=1 \atop i \not\equiv 0 \left(\mathop{mod} p\right)}^{p^k}i\right)^{\left\lfloor\frac{n}{p^k}\right\rfloor}\prod_{i=\left\lfloor\frac{n}{p^k}\right\rfloor p^k \atop i \not\equiv 0\left(\mathop{mod} p\right)}^{n}i$
然而 $\left\lfloor\frac{n}{p}\right\rfloor!$ 可能还有 $p^k$ 的因数，所以要继续递归
定义函数
$f\left(n\right) = \frac{n!}{p^k}$
显然
$f\left(n \right) = f\left(\left\lfloor\frac{n}{p}\right\rfloor\right)\left(\prod_{i=1 \atop i \not\equiv 0 \left(\mathop{mod} p\right)}^{p^k}i\right)^{\left\lfloor\frac{n}{p^k}\right\rfloor}\prod_{i=\left\lfloor\frac{n}{p^k}\right\rfloor p^k \atop i \not\equiv 0\left(\mathop{mod} p\right)}^{n}i$
边界条件 $f (0) = 1$

现在要求 $p^{x -y -z}$
定义函数 $g\left(n\right)$ 为 $n!$ 中 $p$ 的因数个数
类似于 $f$ ，可以得到
$g\left(n \right) = \left\lfloor\frac{n}{p}\right\rfloor + g\left(\left\lfloor\frac{n}{p}\right\rfloor\right)$
所以
$a_i = f\left(n\right) * inv\left(f\left(m\right)\right)inv\left(f\left(n-m\right)\right) g^{x-y-z} \left(\mathop{mod} p_i^{\alpha_i}\right)$
最后用中国剩余定理合并
$\equiv a_i\left(\mathop{mod} p_i^{\alpha_i}\right)$

代码

const int maxN = 1000000;
typedef long long LL;
LL a[maxN], b[maxN], cnt;

LL gcd(LL a, LL b) {
	LL c;
	while (b) {
		c = a % b;
		a = b;
		b = c;
	}
	return a;
}

LL extgcd(LL a, LL b, LL& x, LL& y) {
	if (b == 0) {
		x = 1;
		y = 0;
		return a;
	}
	LL ans = extgcd(b, a % b, y, x);
	y -= a / b * x;
	return ans;
}

LL inv(LL a, LL p) {
	LL x, y;
	extgcd(a, p, x, y);
	return (x % p + p) % p;
}

LL quick_mul(LL a, LL b, LL c) {
	a = (a % c + c) % c;
	b = (b % c + c) % c;
	LL res = 0;
	while (b) {
		if (b & 1)res = (res + a) % c;
		a = (a + a) % c;
		b >>= 1;
	}
	return res;
}

LL quick_pow(LL a, LL b, LL c) {
	LL res = 1;
	while (b) {
		if (b & 1)res = (res * a) % c;
		a = (a * a) % c;
		b >>= 1;
	}
	return res;
}

LL f(LL n, LL p, LL pk) {
	if (n == 0) {
		return 1;
	}
	LL ans = 1;
	//循环节
	for (LL i = 1; i < pk; ++i) {
		if (i % p != 0) {
			ans = ans * i % pk;
		}
	}
	//循环节的n/pk次方
	ans = quick_pow(ans, n / pk, pk);
	//余项
	for (LL i = n / pk * pk; i <= n; ++i) {
		if (i % p != 0) {
			ans = i % pk * ans % pk;
		}
	}
	return ans * (f(n / p, p, pk) % pk) % pk;
}

LL g(LL n, LL p) {
	if (n < p) {
		return 0;
	}
	return n / p + g(n / p, p);
}
// Cn_m % pk
LL c_pk(LL n, LL m, LL p, LL pk) {
	// (n!/ p^x)*inv(m!/p^y) * inv((n-m)!/p^z) * p^{z-y-x} % pk
	return f(n, p, pk) * inv(f(m, p, pk), pk) % pk * inv(f(n - m, p, pk), pk) % pk * quick_pow(p, g(n, p) - g(m, p) - g(n - m, p), pk);
}

//中国剩余定理, x=a(mod b) ,n是总共有几个同余方程，M用来接收最后的公倍数
LL china(LL a[], LL b[], LL n, LL& M) {
	LL ans = 0, lcm = 1, x = 0, y = 0;
	for (LL i = 0; i < n; ++i)lcm *= b[i];
	for (LL i = 0; i < n; ++i) {
		LL Mi = lcm / b[i];
		extgcd(Mi, b[i], x, y);
		LL ti = inv(Mi, b[i]);
		ans = (ans + a[i] % lcm * ti % lcm * Mi % lcm) % lcm;
	}
	M = lcm;
	//返回最小正整数解
	return (ans % lcm + lcm) % lcm;
}

LL exLucas(LL n, LL m, LL p) {
	LL temp = p;
	//质因数分解
	for (LL i = 2; i * i <= temp; ++i) {
		if (temp % i == 0) {
			LL pk = 1;
			while (temp % i == 0) {
				pk *= i;
				temp /= i;
			}
			b[cnt] = pk;
			a[cnt] = c_pk(n, m, i, pk);
			++cnt;
		}
	}
	if (temp > 1) {
		b[cnt] = temp;
		a[cnt] = c_pk(n, m, temp, temp);
		++cnt;
	}
	LL M;
	return china(a, b, cnt, M);
}

例题

hdu3037

#include<iostream>
using namespace std;
typedef long long int ll;
int n, m, p;
//快速幂a^m
ll quickMod(const ll &a,ll m) {
    ll base = a, ans = 1;
    while (m) {
        if (m & 1)
            ans = ans*base%p;
        base = base*base%p;
        m >>= 1;
    }
    return ans;
}
//费马小定理求逆元
ll modReverse(const ll &a, const ll p) {
    return quickMod(a, p - 2);
}
ll c(ll n, ll m) {
    if (m > n)
        return 0;
    ll a = 1, b = 1;
    while (m) {
        a = a*n%p;
        b = b*m%p;
        --n;
        --m;
    }
    return a*modReverse(b, p) % p;
}
//Lucas定理
ll Lucas(ll n, ll m) {
    if (m == 0)
        return 1;
    return c(n%p, m%p)*Lucas(n / p, m / p) % p;
}
int main() {
    int T;
    cin >> T;
    while (T--) {
        cin >> n >> m >> p;
        /*
        题目相当于不超过m个球放入n个盒子
        m个球放入n个盒子如果盒子不为空，则利用插板法，m-1个空插上n-1个挡板c(m-1,n-1)
        如果允许空，则相当于每个盒子先放一个球，然后不允许空，n+m个球放入n个盒子，c(n+m-1,n-1)=c(n+m-1,m)
        所以m个球放入n个盒子且允许空的方案个数为c(n+m-1,m)
        因为是题目是0-m个，则ans=c(n+0-1,0)+c(n+1-1,1)+...+c(n+m-1,m)
        由公式C(n,k) = C(n-1,k)+C(n-1,k-1)得
        ans=c(n-1,0)+c(n,1)+c(n+1,2)...c(n+m-1,m)=c(n,0)+c(n,1)+c(n+1,2)+...+c(n+m-1,m)
          =c(n+1,1)+c(n+1,2)+...+c(n+m-1,m)
          =c(n+m,m)
        */
        cout << Lucas(n + m, m) << endl;
    }
    return 0;
}

洛谷P4720

#include<cstdio>
using namespace std;
const int maxN = 1000000;
typedef long long LL;
LL a[maxN], b[maxN], cnt;

LL gcd(LL a, LL b) {
	LL c;
	while (b) {
		c = a % b;
		a = b;
		b = c;
	}
	return a;
}

LL extgcd(LL a, LL b, LL& x, LL& y) {
	if (b == 0) {
		x = 1;
		y = 0;
		return a;
	}
	LL ans = extgcd(b, a % b, y, x);
	y -= a / b * x;
	return ans;
}

LL inv(LL a, LL p) {
	LL x, y;
	extgcd(a, p, x, y);
	return (x % p + p) % p;
}

LL quick_mul(LL a, LL b, LL c) {
	a = (a % c + c) % c;
	b = (b % c + c) % c;
	LL res = 0;
	while (b) {
		if (b & 1)res = (res + a) % c;
		a = (a + a) % c;
		b >>= 1;
	}
	return res;
}

LL quick_pow(LL a, LL b, LL c) {
	LL res = 1;
	while (b) {
		if (b & 1)res = (res * a) % c;
		a = (a * a) % c;
		b >>= 1;
	}
	return res;
}

LL f(LL n, LL p, LL pk) {
	if (n == 0) {
		return 1;
	}
	LL ans = 1;
	//循环节
	for (LL i = 1; i < pk; ++i) {
		if (i % p != 0) {
			ans = ans * i % pk;
		}
	}
	//循环节的n/pk次方
	ans = quick_pow(ans, n / pk, pk);
	//余项
	for (LL i = n / pk * pk; i <= n; ++i) {
		if (i % p != 0) {
			ans = i % pk * ans % pk;
		}
	}
	return ans * (f(n / p, p, pk) % pk) % pk;
}

LL g(LL n, LL p) {
	if (n < p) {
		return 0;
	}
	return n / p + g(n / p, p);
}
// Cn_m % pk
LL c_pk(LL n, LL m, LL p, LL pk) {
	// (n!/ p^x)*inv(m!/p^y) * inv((n-m)!/p^z) * p^{z-y-x} % pk
	return f(n, p, pk) * inv(f(m, p, pk), pk) % pk * inv(f(n - m, p, pk), pk) % pk * quick_pow(p, g(n, p) - g(m, p) - g(n - m, p), pk);
}

//中国剩余定理, x=a(mod b) ,n是总共有几个同余方程，M用来接收最后的公倍数
LL china(LL a[], LL b[], LL n, LL& M) {
	LL ans = 0, lcm = 1, x = 0, y = 0;
	for (LL i = 0; i < n; ++i)lcm *= b[i];
	for (LL i = 0; i < n; ++i) {
		LL Mi = lcm / b[i];
		extgcd(Mi, b[i], x, y);
		LL ti = inv(Mi, b[i]);
		ans = (ans + a[i] % lcm * ti % lcm * Mi % lcm) % lcm;
	}
	M = lcm;
	//返回最小正整数解
	return (ans % lcm + lcm) % lcm;
}

LL exLucas(LL n, LL m, LL p) {
	LL temp = p;
	//质因数分解
	for (LL i = 2; i * i <= temp; ++i) {
		if (temp % i == 0) {
			LL pk = 1;
			while (temp % i == 0) {
				pk *= i;
				temp /= i;
			}
			b[cnt] = pk;
			a[cnt] = c_pk(n, m, i, pk);
			++cnt;
		}
	}
	if (temp > 1) {
		b[cnt] = temp;
		a[cnt] = c_pk(n, m, temp, temp);
		++cnt;
	}
	LL M;
	return china(a, b, cnt, M);
}
int main() {
	LL n, m, p;
	scanf("%lld%lld%lld", &n, &m, &p);
	printf("%lld\n", exLucas(n, m, p));
	return 0;
}

猜你喜欢

python 程序员进阶之路：从新手到高手的100个模块
【云栖大会】人工智能：智，在云端
CSDN如何获得2020技术圈认证（新徽章哦）
iOS开发之网络数据解析--GDataXML解析框架的使用
Fate1.6 支持的机器学习算法
HDFS中JAVA API的使用
[ES6] Spread Operator
ES6的新特性
[React] Setup 'beforeunload' listener
Java实现 LeetCode 75 颜色分类
zabbix导入数据出现ERROR 1046 (3D000) at line 1: No database selected
【Python3网络爬虫开发实战】 3.1.3-解析链接
PoolFormer实战：使用PoolFormer实现图像分类任务（一）
第七届蓝桥杯JavaC组国(决)赛真题
java实现罗马数字转十进制
解决办法：错误异常The superclass "javax.servlet.http.HttpServlet" was not found on the Java Build Path

相关主题

CAP定理
CAP 定理
中国剩余定理
2.3 中国剩余定理
中心极限定理
Lucas定理
矩阵树定理
卢卡斯定理
夹逼定理
采样定理
中值定理
中值定理们
裴蜀定理
欧拉定理

zl程序教程

当前栏目

Lucas定理

Lucas定理

引理

证明

代码

扩展lucas

代码

例题

hdu3037

洛谷P4720

相关文章