Hatena::Grouptopcoder

(iwi) { 反省します

TopCoder: [[iwi]] / Twitter: @iwiwi

 | 

2013-01-05

pow_mod, inverse の速度

14:22 | pow_mod, inverse の速度 - (iwi) { 反省します を含むブックマーク はてなブックマーク - pow_mod, inverse の速度 - (iwi) { 反省します pow_mod, inverse の速度 - (iwi) { 反省します のブックマークコメント

久々にコンテストの問題が解きたくなって AOJ で 問題 MinimumCostPath を解いた.個人的にはムチャクチャ面白いと思った.

ところで,別に TLE はしなかったんだが,mod の inverse がボトルネックだった.自分のライブラリの mod の inverse とからへん,コンテスト始めたころの大昔に作ってから変えてなくて,今やあまり気に入っていないと思いつつバグってはいないということで長らく放置していたが,この際ようやく弄ることにした.


ベンチマークコード

適当.g++ -O2.標準入力は適当に10とか入れてる.

const ll MOD = 1000000009;

const int N = 2000000;

int main() {
  int a;
  scanf("%d", &a);

  ll t = a;
  for (int k = 1; k <= N; ++k) {
    t = inverse(t, MOD) * k % MOD;
  }
  printf("%lld\n", t);

  return 0;
}

追記:手元の環境は

~/contests/tmp% cat /proc/cpuinfo [00:35:39]

processor : 0

vendor_id : GenuineIntel

cpu family : 6

model : 26

model name : Intel(R) Core(TM) i7 CPU 960 @ 3.20GHz

stepping : 5

cpu MHz : 1600.000

cache size : 8192 KB

~/contests/tmp% gcc -v [00:36:44]

Using built-in specs.

Target: x86_64-linux-gnu

Configured with: ../src/configure -v --with-pkgversion='Ubuntu 4.4.3-4ubuntu5.1' --with-bugurl=file:///usr/share/doc/gcc-4.4/README.Bugs --enable-languages=c,c++,fortran,objc,obj-c++ --prefix=/usr --enable-shared --enable-multiarch --enable-linker-build-id --with-system-zlib --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --with-gxx-include-dir=/usr/include/c++/4.4 --program-suffix=-4.4 --enable-nls --enable-clocale=gnu --enable-libstdcxx-debug --enable-plugin --enable-objc-gc --disable-werror --with-arch-32=i486 --with-tune=generic --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu

Thread model: posix

gcc version 4.4.3 (Ubuntu 4.4.3-4ubuntu5.1)

今までのやつ

駆け出しだった頃の僕が,適当に bmerry か誰かのやつを参考にして書いた.中で再帰している上にいちいち丁寧に正にするために % を 2 回もとってて遅そう.(でもマイナスが帰るとそれが原因のバグをつくりそうで恐ろしいという気持ちは理解できる)

inline ll mod(ll a, ll m) { return (a % m + m) % m; }

ll inverse(ll a, ll m) {
  if ((a = mod(a, m)) == 1) return 1;
  return mod((1 - m * inverse(m % a, a)) / a, m);
}

~/contests/tmp% time ./a.out <<< 10 [14:10:35]

249011238

./a.out <<< 10 2.84s user 0.01s system 99% cpu 2.870 total

~/contests/tmp% time ./a.out <<< 10 [14:10:39]

249011238

./a.out <<< 10 2.85s user 0.00s system 99% cpu 2.870 total

~/contests/tmp% time ./a.out <<< 10 [14:10:43]

249011238

./a.out <<< 10 2.83s user 0.02s system 99% cpu 2.868 total


inline

再帰関数だからあんま意味分からんと思うが一応 inverse に inline つけてみた.

~/contests/tmp% g++ -O2 inverse.cpp [14:11:39]

~/contests/tmp% time ./a.out <<< 10 [14:11:41]

249011238

./a.out <<< 10 2.75s user 0.00s system 99% cpu 2.763 total

~/contests/tmp% time ./a.out <<< 10 [14:11:44]

249011238

./a.out <<< 10 2.75s user 0.00s system 99% cpu 2.762 total

~/contests/tmp% time ./a.out <<< 10 [14:11:49]

249011238

./a.out <<< 10 2.75s user 0.00s system 100% cpu 2.748 total

ちょっぴりはやくなったw



pow_mod 1

今まで使っていた pow_mod を使う.これも再帰していて丁寧に mod を呼んでいて見るからに遅そうである.

inline ll mod(ll a, ll m) { return (a % m + m) % m; }

ll pow_mod(ll a, ll k, ll m) {
  if (k == 0) return 1;
  ll t = pow_mod(a, k / 2, m);
  return mod(mod(t * t, m) * (k % 2 ? a : 1), m);
}

inline ll inverse(ll a, ll m) {
  return pow_mod(a, m - 2, m);
}

~/contests/tmp% time ./a.out <<< 10 [14:13:19]

249011238

./a.out <<< 10 3.39s user 0.00s system 99% cpu 3.407 total

~/contests/tmp% time ./a.out <<< 10 [14:13:23]

249011238

./a.out <<< 10 3.39s user 0.00s system 99% cpu 3.407 total

~/contests/tmp% time ./a.out <<< 10 [14:13:28]

249011238

./a.out <<< 10 3.37s user 0.02s system 99% cpu 3.411 total

遅くなったw



pow_mod 2

pow_mod を高速そうな実装にしてみる.

inline ll pow_mod(ll a, ll k, ll m) {
  ll r = 1;
  for (; k > 0; k >>= 1) {
    if (k & 1) (r *= a) %= m;
    (a *= a) %= m;
  }
  return r;
}

~/contests/tmp% time ./a.out <<< 10 [14:16:27]

249011238

./a.out <<< 10 0.45s user 0.01s system 99% cpu 0.463 total

~/contests/tmp% time ./a.out <<< 10 [14:16:28]

249011238

./a.out <<< 10 0.47s user 0.00s system 100% cpu 0.468 total

~/contests/tmp% time ./a.out <<< 10 [14:16:31]

249011238

./a.out <<< 10 0.45s user 0.01s system 99% cpu 0.463 total

かなりはやくなったw

少し予想はしていたが,今までのはやはり爆遅だった



非再帰 exgcd

Mi_Sawa さんに教えてもらった(ありがとうございます).確かに速そう.

inline ll exgcd(ll a, ll b, ll &x, ll &y){
  ll u[] = {a, 1, 0}, v[] = {b, 0, 1};
  while(*v){
    ll t = *u / *v;
    rep(i, 3) swap(u[i] -= t * v[i], v[i]);
  }
  x = u[1]; y = u[2];
  return u[0];
}
inline ll inverse(ll a, ll m){
  ll p, q;
  exgcd(a, m, p, q);
  return (p%m+m)%m;
}
手元環境

~/contests/tmp% g++ -O2 inverse.cpp [23:26:03]

~/contests/tmp% time ./a.out <<< 10 [23:26:08]

249011238

./a.out <<< 10 0.65s user 0.00s system 100% cpu 0.649 total

~/contests/tmp% time ./a.out <<< 10 [23:26:10]

249011238

./a.out <<< 10 0.64s user 0.01s system 100% cpu 0.644 total

~/contests/tmp% time ./a.out <<< 10 [23:26:16]

249011238

./a.out <<< 10 0.64s user 0.01s system 100% cpu 0.644 total

自分の環境だと僅かに pow_mod より遅い.

ideone

しかし,ideone に貼りつけたら,exgcd のほうが 2 倍ぐらい速い.

Codeforces

ideone と似たっけかになった.

pow_mod 0: 1.015625 sec

pow_mod 1: 1.014648 sec

pow_mod 2: 1.014649 sec

exgcd 0: 0.695312 sec

exgcd 1: 0.694336 sec

exgcd 2: 0.694336 sec

TopCoder

やっぱしexgcd のほうが 2 倍ぐらい速い.

pow_mod 0: 0.369712 sec

pow_mod 1: 0.351320 sec

pow_mod 2: 0.351127 sec

exgcd 0: 0.171135 sec

exgcd 1: 0.170425 sec

exgcd 2: 0.170366 sec

結論

exgcd のほうが速い環境がどうやら多い,pow_mod のほうが速い時も差は僅か

exgcd のほうがよさそうかも



非再帰 exgcd ループアンローリング

ループアンローリングしたらちょっと速くなった

http://ideone.com/Dz1LLI



整理して

inline ll inverse(ll a, ll m){
  ll b = m, u = 1, v = 0;
  while (b) {
    ll t = a / b;
    swap(a -= t * b, b);
    swap(u -= t * v, v);
  }
  return (u % m + m) % m;
}

Mi_SawaMi_Sawa2013/01/05 20:38自分の環境ではexgcdの方が速かったです.
(pow_modだとコンスタントにlog Mくらいかかりますが, exgcdだと速い時がそれなりにあるので, 良さそうでもある)

inline ll exgcd(ll a, ll b, ll &x, ll &y){
ll u[] = {a, 1, 0}, v[] = {b, 0, 1};
while(*v){
ll t = *u / *v;
rep(i, 3) swap(u[i] -= t * v[i], v[i]);
}
x = u[1]; y = u[2];
return u[0];
}
inline ll invMod(ll a, ll m){
ll p, q;
exgcd(a, m, p, q);
return (p%m+m)%m;
}

iwiwiiwiwi2013/01/06 00:23コメントありがとうございます!!
試してみました,詳しくは記事に追記します.自分の環境だと僅かに pow_mod のほうが高速でしたが,ideone で実行したら確かに貼ってくださった exgcd のほうが 2 倍とか速かったです.

Mi_SawaMi_Sawa2013/01/06 14:36もう一つよいのは, Mが素数である必要がなく, 逆元が無い時も,
ax = gcd(a,M)
なるxを見つけてくれる所ですね.
負の数にちゃんと対応していなかったり(普通に正規化していれば大丈夫ですが, 一度これで引っかかった), 逆元の存在のチェックをしていないので, ロバストにするなら,
ll u[] = {llabs(a), 1, 0}, v[] = {llabs(b), 0, 1};
にしたり,
assert(exgcd(a, m, p, q) == 1)
にしておくとよいです.

Mi_SawaMi_Sawa2013/01/06 16:22あ, ll u[] = {llabs(a), 1, 0}, v[] = {llabs(b), 0, 1};だとダメな気がします, コメント汚して申し訳ないです…
ループ終わったあと, u[0]が負ならrep(i, 3) u[i] = -u[i];とかすればよい気がします…

iwiwiiwiwi2013/01/07 11:35ありがとうございます.張り付けで使うならちょっと長くなっても周りのちょっとしたミスも吸収できるようにできるだけロバストにしておくのが良さそうですよね.

utevaxisuutevaxisu2017/04/29 07:42http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

onecodijonecodij2017/04/29 07:47http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

uigiqikoleuigiqikole2017/04/29 08:04http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

uduluvouduluvo2017/04/29 08:36http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

aqoosotawopaqoosotawop2017/04/29 08:39http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

oqawaaahobooqawaaahobo2017/04/29 08:39http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

eqafifedozovpeqafifedozovp2017/04/29 08:59http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

epesameepesame2017/04/29 09:00http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

itejifuqaitejifuqa2017/04/29 09:16http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

edkuyerugowaedkuyerugowa2017/04/29 09:30http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

omococucomococuc2017/04/29 09:52http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

obefodaeqobefodaeq2017/04/29 14:35http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

iwohaxemiiwohaxemi2017/04/30 01:40http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

oyelujloyelujl2017/04/30 12:01http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

oafocibebupovoafocibebupov2017/04/30 12:25http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

hopuyqahopuyqa2017/05/01 03:27http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

ocusnasitewocusnasitew2017/05/01 03:46http://100mgcheapest-price-viagra.com/ - 100mgcheapest-price-viagra.com.ankor <a href="http://tadalafil-buy-5mg.com/">tadalafil-buy-5mg.com.ankor</a> http://20mgprednisone-order.com/

 |