關於演算法競賽中快速乘的一些優化

01-27

0x00 寫在前面

在一些演算法競賽題目中經常遇到一些需要取模的乘法運算，一般都是會爆 $C++$ 中 $long;long$ 範圍的。

以下的計算均只考慮正整數的情況。

所以有人按著快速冪的思想提出了 $O(logW)$ （ $W$ 是乘數）的快速乘。

inline long long mmul ( long long a, long long b, const long long& Mod ) {ntlong long rt ( 0 ) ;ntfor ( ; b ; b >>= 1, ( a <<= 1 ) %= Mod )nttif ( b & 1 ) {nttt( rt += a ) %= Mod ;ntt}ntreturn rt ;n}n

複雜度好高吖！！！

也有加以思考提出 $O(1)$ 的使用大範圍的 $long;double$ 的快速乘的。

inline long long mmul ( long long a, long long b, const long long& Mod ) {n a %= Mod, b %= Mod ;n return ( a * b - ( long long ) ( ( ( long double ) a * b + 0.5 ) / Mod ) * Mod ) % Mod ;n}n

對於模數大於 $sqrt{LLong_MAX}$ 的，它也會爆掉（計算 $acdot b$ 爆掉）。

但是上面的方法都有時間複雜度過高或者是常數大的缺點。

為什麼不使用小學知識的乘法分配律來優化呢？

0x01 O(1)的常數較小的快速乘

我們可以考慮提出一個把大數拆成兩個數之和的形式，使得可以在 $long;long$ 範圍內計算。

我們要計算 $a cdot b;mod;p$

設 $b=lf+rg$

那麼原式為

$acdot(lf+rg);mod;p=((acdot lf);mod;p+(acdot rg);mod;p);mod;p$

我們把 $lf$ 欽定為 $b$ 的二進位前 $x$ 位， $rg$ 為 $b$ 的後 $(64-x)$ 位。

$x$ 根據實際情況來定吧，實在不行把 $a,b$ 都拆了，可以應對模數大於 $sqrt{LLong_MAX}$ 的情況。

參考代碼：

inline long long mmul ( long long a, long long b, const long long& Mod ) {ntlong long lf = a * ( b >> 25LL ) % Mod * ( 1LL << 25 ) % Mod ;n long long rg = a * ( b & ( ( 1LL << 25 ) - 1 ) ) % Mod ;ntreturn ( lf + rg ) % Mod ;n}n

0x02 速度測試

我隨機了 $1e7$ 組數據，數據範圍均為 $1e10$ 以內（對於第二種方法可能小小有失公正）。

使用環境（機房電腦）為：

i5-4590 @3.30GHz
4GB 內存
Windows 7系統（這個機房沒有Linux系統......）
G++ 4.9.2

編譯命令：

-O2 -std=c++11 -Wl,--stack=80000000 -Wall -Wconversion -Wextran

這是生成數據的 $generator$ ：

# include <bits/stdc++.h>nnlong long r ( ) {ntreturn ( rand ( ) << 15LL | rand ( ) ) << 16LL | ( ( rand ( ) << 15LL | rand ( ) ) ) ;n}nnint main ( ) {ntsrand ( time ( 0 ) ) ;ntfreopen ( "in.txt", "w", stdout ) ;ntint T = 10000000 ;ntconst long long Mod = 10000000000LL ;nntwhile ( T -- ) {nttlong long a = r ( ) % Mod, b = r ( ) % Mod ;nttif ( a < 0 ) a += Mod ;nttif ( b < 0 ) b += Mod ;nttassert ( a >= 0 && b >= 0 ) ;nttprintf ( "%lld %lldn", a, b ) ;nt}ntreturn 0 ;n}n

這是測試用的 $cpp$ 文件：

# include <bits/stdc++.h>nninline long long mmul1 ( long long a, long long b, const long long& Mod ) {ntlong long lf = a * ( b >> 25LL ) % Mod * ( 1LL << 25 ) % Mod ;n long long rg = a * ( b & ( ( 1LL << 25 ) - 1 ) ) % Mod ;ntreturn ( lf + rg ) % Mod ;n}nninline long long mmul2 ( long long a, long long b, const long long& Mod ) {ntlong long rt ( 0 ) ;ntfor ( ; b ; b >>= 1LL, ( a <<= 1LL ) %= Mod )nttif ( b & 1 ) {nttt( rt += a ) %= Mod ;ntt}ntreturn rt ;n}nninline long long mmul3 ( long long a, long long b, const long long& Mod ) {n a %= Mod, b %= Mod ;n return ( a * b - ( long long ) ( ( ( long double ) a * b + 0.5 ) / Mod ) * Mod ) % Mod ;n}nnint main ( ) {ntfreopen ( "in.txt", "r", stdin ) ;ntfreopen ( "result.txt", "w", stdout ) ;ntint t1 ( 0 ), t2 ( 0 ), t3 ( 0 ) ;ntlong long a, b, ans1, ans2, ans3 ;nt# define MOD 1000000007LLntwhile ( ~ scanf ( "%lld%lld", & a, & b ) ) {nttint start = clock ( ) ;nttans1 = mmul1 ( a, b, MOD ) ;nttt1 += clock ( ) - start ;nttstart = clock ( ) ;nttans2 = mmul2 ( a, b, MOD ) ;nttt2 += clock ( ) - start ;nttstart = clock ( ) ;nttans3 = mmul3 ( a, b, MOD ) ;nttt3 += clock ( ) - start ;nt//tfprintf ( stderr, "%lld %lld %lldn", ans1, ans2, ans3 ) ;nttassert ( ans1 == ans2 && ans2 == ans3 ) ;nt}nt# undef MODntputs ( "test result :" ) ;ntprintf ( "mmult1 : %.3lfsn", 1.0 * t1 / CLOCKS_PER_SEC ) ;ntprintf ( "mmult2 : %.3lfsn", 1.0 * t2 / CLOCKS_PER_SEC ) ;ntprintf ( "mmult3 : %.3lfsn", 1.0 * t3 / CLOCKS_PER_SEC ) ;ntreturn 0 ;n}n