欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

【小技巧】O(1)快速乘

程序员文章站 2022-04-28 15:36:59
...
  • 问题:求 a×b mod pa\times b\ mod\ pa,b,pa,b,plong long 范围内。
  • CRT 等算法中应用广泛。
  • 为了处理模数在 int 范围外的情况,就是两数相乘可能会爆 long long 时,我们不能直接用整型的乘法来计算。
  • 首先我们可以进行二进制拆分,化乘法为加法,类似快速幂那样,写出一个 O(logn)O(\log n) 的快速乘
typedef long long s64; 

inline void add(s64 &a, const s64 &b)
{
	a += b; 
	if (a >= mod)
		a -= mod; 
}

inline s64 qmul(s64 a, s64 b, const s64 &mod)
{
	a = (a % mod + mod) % mod; 
	b = (b % mod + mod) % mod; //这两行依据情况不写
	s64 res = 0; 
	for (; b; b >>= 1, add(a, a, mod))
		if (b & 1)
			add(res, a, mod); 
	return res; 
}
  • 多次使用时,为了避免毒瘤出题人卡时间(或是为了优化常数),我们可以利用 long double 写出一个优秀的 O(1)O(1) 快速乘。
  • 简单原理:a×b mod p=a×ba×bp×pa\times b\ mod\ p=a\times b-\lfloor \frac{a\times b}{p}\rfloor\times p
  • 利用 long double 来处理这个 a×bp\lfloor \frac{a\times b}{p}\rfloor
  • 然后处理一下浮点误差就可以了。
  • 模数较大时可能会出锅。
  • 不过出锅概率很小
typedef long long s64; 
typedef long double ld; 

inline s64 qmul(s64 a, s64 b, s64 mod)
{
	a = (a % mod + mod) % mod; 
	b = (b % mod + mod) % mod; //这两行依据情况不写
	s64 res = a * b - (s64)((ld)a / mod * b + 1e-8) * mod; 
	return res < 0 ? res + mod : res;  
}
相关标签: 小技巧