• branch=trunk inherited from [9165bd3629]
• sym-trunk inherited from [9165bd3629]

Added SRM/518-U/1C.cpp version [1b6b6f1bf178804b]

            1  +#include <iostream>
            2  +#include <sstream>
            3  +#include <iomanip>
            4  +#include <vector>
            5  +#include <string>
            6  +#include <map>
            7  +#include <set>
            8  +#include <algorithm>
            9  +#include <numeric>
           10  +#include <iterator>
           11  +#include <functional>
           12  +#include <complex>
           13  +#include <queue>
           14  +#include <valarray>
           15  +#include <stack>
           16  +#include <cmath>
           17  +#include <cassert>
           18  +#include <cstring>
           19  +using namespace std;
           20  +typedef long long LL;
           21  +typedef complex<double> CMP;
           22  +
           23  +static const int MODVAL = 1000000007;
           24  +struct mint
           25  +{
           26  + int val;
           27  + mint():val(0){}
           28  + mint(int    x):val(x%MODVAL) {}
           29  + mint(size_t x):val(x%MODVAL) {}
           30  + mint(LL     x):val(x%MODVAL) {}
           31  +};
           32  +mint& operator+=(mint& x, mint y) { return x = x.val+y.val; }
           33  +mint& operator-=(mint& x, mint y) { return x = x.val-y.val+MODVAL; }
           34  +mint& operator*=(mint& x, mint y) { return x = LL(x.val)*y.val; }
           35  +mint POW(mint x, LL e) { mint v=1; for(;e;x*=x,e>>=1) if(e&1) v*=x; return v; }
           36  +mint& operator/=(mint& x, mint y) { return x *= POW(y, MODVAL-2); }
           37  +mint operator+(mint x, mint y) { return x+=y; }
           38  +mint operator-(mint x, mint y) { return x-=y; }
           39  +mint operator*(mint x, mint y) { return x*=y; }
           40  +mint operator/(mint x, mint y) { return x/=y; }
           41  +mint inv2 = mint(1) / 2;
           42  +
           43  +class Nim { public:
           44  + int count(int K, int L)
           45  + {
           46  +  valarray<mint> v(65536);
           47  +  v[slice(2,L-1,1)] = 1;
           48  +  for(unsigned p=2; p<=L; ++p)
           49  +   if( v[p].val )
           50  +    for(unsigned q=p*p; q<=L; q+=p)
           51  +     v[q] = 0;
           52  +
           53  +  pre(v, 0, 65536);
           54  +  valarray<mint> a(1, 65536);
           55  +  for(int i=1; i<=K; i<<=1, v*=v)
           56  +   if( K & i )
           57  +    a *= v;
           58  +  post(a, 0, 65536);
           59  +  return a[0].val;
           60  + }
           61  +
           62  + void pre(valarray<mint>& v, int s, int e)
           63  + {
           64  +  if( s+1 == e )
           65  +   return;
           66  +  int m = (s+e)/2;
           67  +  pre(v, s, m);
           68  +  pre(v, m, e);
           69  +  for(int i=s,j=m; i<m; ++i,++j) {
           70  +   mint vi=v[i], vj=v[j];
           71  +   v[i] = vi - vj;
           72  +   v[j] = vi + vj;
           73  +  }
           74  + }
           75  +
           76  + void post(valarray<mint>& v, int s, int e)
           77  + {
           78  +  if( s+1 == e )
           79  +   return;
           80  +  int m = (s+e)/2;
           81  +  for(int i=s,j=m; i<m; ++i,++j) {
           82  +   mint dif=v[i], sum=v[j];
           83  +   v[i] = (sum + dif) * inv2;
           84  +   v[j] = (sum - dif) * inv2;
           85  +  }
           86  +  post(v, s, m);
           87  +  post(v, m, e);
           88  + }
           89  +};
           90  +
           91  +// BEGIN CUT HERE
           92  +#include <ctime>
           93  +double start_time; string timer()
           94  + { ostringstream os; os << " (" << int((clock()-start_time)/CLOCKS_PER_SEC*1000) << " msec)"; return os.str(); }
           95  +template<typename T> ostream& operator<<(ostream& os, const vector<T>& v)
           96  + { os << "{ ";
           97  +   for(typename vector<T>::const_iterator it=v.begin(); it!=v.end(); ++it)
           98  +   os << '\"' << *it << '\"' << (it+1==v.end() ? "" : ", "); os << " }"; return os; }
           99  +void verify_case(const int& Expected, const int& Received) {
          100  + bool ok = (Expected == Received);
          101  + if(ok) cerr << "PASSED" << timer() << endl;  else { cerr << "FAILED" << timer() << endl;
          102  + cerr << "\to: \"" << Expected << '\"' << endl << "\tx: \"" << Received << '\"' << endl; } }
          103  +#define CASE(N) {cerr << "Test Case #" << N << "..." << flush; start_time=clock();
          104  +#define END  verify_case(_, Nim().count(K, L));}
          105  +int main(){
          106  +
          107  +CASE(0)
          108  + int K = 3;
          109  + int L = 7;
          110  + int _ = 6;
          111  +END
          112  +CASE(1)
          113  + int K = 4;
          114  + int L = 13;
          115  + int _ = 120;
          116  +END
          117  +CASE(2)
          118  + int K = 10;
          119  + int L = 100;
          120  + int _ = 294844622;
          121  +END
          122  +CASE(3)
          123  + int K = 123456789;
          124  + int L = 12345;
          125  + int _ = 235511047;
          126  +END
          127  +CASE(4)
          128  + int K = 1000000000;
          129  + int L = 50000;
          130  + int _ = 428193537;
          131  +END
          132  +/*
          133  +CASE(5)
          134  + int K = ;
          135  + int L = ;
          136  + int _ = ;
          137  +END
          138  +*/
          139  +}
          140  +// END CUT HERE

Added lib/doc/fft.txt version [893844b7e5ed0726]

            1  +
            2  +http://apps.topcoder.com/wiki/display/tc/SRM+518
            3  +
            4  +
            5  +以下のような関数を考える。
            6  +
            7  +int[] xor_count(int[] a, int[] b, int N)
            8  +{
            9  + int[] c = {0} * N;
           10  + for(int i=0; i<N; ++i)
           11  +  for(int k=0; k<N; ++k)
           12  +   c[i ^ k] += a[i] * b[k];
           13  + return c;
           14  +}
           15  +
           16  +高速化。
           17  +
           18  +int[] xor_count_ultra(int[] a, int[] b, int N)
           19  +{
           20  + a = xor_pre(a);
           21  + b = xor_pre(b);
           22  + c = a * b; // element-wise product
           23  + return xor_post(c);
           24  +}
           25  +
           26  +このようになるマジカルな変換 xor_pre, xor_post が存在する。
           27  +
           28  +N=1 の場合、pre : {x} --> {x}
           29  +N=2 の場合、pre : {x,y} --> {x-y, x+y}
           30  +  xor_count( {x,y}, {u,v} ) = {xu+yv, yu+xv}
           31  +  と {x-y, x+y} * {u-v, u+v} = {xu+yv-xu-xv, xu+yv+yu+xv}
           32  +  から確認できる。
           33  +一般にNが2のベキの場合、xorの作用するビットのパターンが再帰的なので
           34  +  pre( [X,Y] ) = [pre(X)-pre(Y), pre(X)+pre(Y)]
           35  +と再帰的にすれば成り立つ
           36  +
           37  +
           38  +
           39  +
           40  +
           41  +
           42  +xorでなくて足し算の場合
           43  +
           44  +int[] plus_count(int[] a, int[] b, int N)
           45  +{
           46  + int[] c = {0} * 2N;
           47  + for(int i=0; i<N; ++i)
           48  +  for(int k=0; k<N; ++k)
           49  +   c[i + k] += a[i] * b[k];
           50  + return c;
           51  +}
           52  +
           53  +も同様に。
           54  +
           55  +int[] plus_count_ultra(int[] a, int[] b, int N)
           56  +{
           57  + a = plus_pre(a);
           58  + b = plus_pre(b);
           59  + c = a * b; // element-wise product
           60  + return plus_post(c);
           61  +}
           62  +
           63  +とできるのが FFT