Advertisement
Guest User

Untitled

a guest
Dec 10th, 2019
114
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 4.60 KB | None | 0 0
  1. namespace fft
  2. {
  3. struct num
  4. {
  5. double x,y;
  6. num() {x=y=0;}
  7. num(double x,double y):x(x),y(y){}
  8. };
  9. inline num operator+(num a,num b) {return num(a.x+b.x,a.y+b.y);}
  10. inline num operator-(num a,num b) {return num(a.x-b.x,a.y-b.y);}
  11. inline num operator*(num a,num b) {return num(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
  12. inline num conj(num a) {return num(a.x,-a.y);}
  13.  
  14. int base=1;
  15. vector<num> roots={{0,0},{1,0}};
  16. vector<int> rev={0,1};
  17. const double PI=acosl(-1.0);
  18.  
  19. void ensure_base(int nbase)
  20. {
  21. if(nbase<=base) return;
  22. rev.resize(1<<nbase);
  23. for(int i=0;i<(1<<nbase);i++)
  24. rev[i]=(rev[i>>1]>>1)+((i&1)<<(nbase-1));
  25. roots.resize(1<<nbase);
  26. while(base<nbase)
  27. {
  28. double angle=2*PI/(1<<(base+1));
  29. for(int i=1<<(base-1);i<(1<<base);i++)
  30. {
  31. roots[i<<1]=roots[i];
  32. double angle_i=angle*(2*i+1-(1<<base));
  33. roots[(i<<1)+1]=num(cos(angle_i),sin(angle_i));
  34. }
  35. base++;
  36. }
  37. }
  38.  
  39. void fft(vector<num> &a,int n=-1)
  40. {
  41. if(n==-1) n=a.size();
  42. assert((n&(n-1))==0);
  43. int zeros=__builtin_ctz(n);
  44. ensure_base(zeros);
  45. int shift=base-zeros;
  46. for(int i=0;i<n;i++)
  47. if(i<(rev[i]>>shift))
  48. swap(a[i],a[rev[i]>>shift]);
  49. for(int k=1;k<n;k<<=1)
  50. {
  51. for(int i=0;i<n;i+=2*k)
  52. {
  53. for(int j=0;j<k;j++)
  54. {
  55. num z=a[i+j+k]*roots[j+k];
  56. a[i+j+k]=a[i+j]-z;
  57. a[i+j]=a[i+j]+z;
  58. }
  59. }
  60. }
  61. }
  62.  
  63. vector<num> fa,fb;
  64.  
  65. vector<int> multiply(vector<int> &a, vector<int> &b)
  66. {
  67. int need=a.size()+b.size()-1;
  68. int nbase=0;
  69. while((1<<nbase)<need) nbase++;
  70. ensure_base(nbase);
  71. int sz=1<<nbase;
  72. if(sz>(int)fa.size()) fa.resize(sz);
  73. for(int i=0;i<sz;i++)
  74. {
  75. int x=(i<(int)a.size()?a[i]:0);
  76. int y=(i<(int)b.size()?b[i]:0);
  77. fa[i]=num(x,y);
  78. }
  79. fft(fa,sz);
  80. num r(0,-0.25/sz);
  81. for(int i=0;i<=(sz>>1);i++)
  82. {
  83. int j=(sz-i)&(sz-1);
  84. num z=(fa[j]*fa[j]-conj(fa[i]*fa[i]))*r;
  85. if(i!=j) fa[j]=(fa[i]*fa[i]-conj(fa[j]*fa[j]))*r;
  86. fa[i]=z;
  87. }
  88. fft(fa,sz);
  89. vector<int> res(need);
  90. for(int i=0;i<need;i++) res[i]=fa[i].x+0.5;
  91. return res;
  92. }
  93.  
  94. vector<int> multiply_mod(vector<int> &a,vector<int> &b,int m,int eq=0)
  95. {
  96. int need=a.size()+b.size()-1;
  97. int nbase=0;
  98. while((1<<nbase)<need) nbase++;
  99. ensure_base(nbase);
  100. int sz=1<<nbase;
  101. if(sz>(int)fa.size()) fa.resize(sz);
  102. for(int i=0;i<(int)a.size();i++)
  103. {
  104. int x=(a[i]%m+m)%m;
  105. fa[i]=num(x&((1<<15)-1),x>>15);
  106. }
  107. fill(fa.begin()+a.size(),fa.begin()+sz,num{0,0});
  108. fft(fa,sz);
  109. if(sz>(int)fb.size()) fb.resize(sz);
  110. if(eq) copy(fa.begin(),fa.begin()+sz,fb.begin());
  111. else
  112. {
  113. for(int i=0;i<(int)b.size();i++)
  114. {
  115. int x=(b[i]%m+m)%m;
  116. fb[i]=num(x&((1<<15)-1),x>>15);
  117. }
  118. fill(fb.begin()+b.size(),fb.begin()+sz,num{0,0});
  119. fft(fb,sz);
  120. }
  121. double ratio=0.25/sz;
  122. num r2(0,-1),r3(ratio,0),r4(0,-ratio),r5(0,1);
  123. for(int i=0;i<=(sz>>1);i++)
  124. {
  125. int j=(sz-i)&(sz-1);
  126. num a1=(fa[i]+conj(fa[j]));
  127. num a2=(fa[i]-conj(fa[j]))*r2;
  128. num b1=(fb[i]+conj(fb[j]))*r3;
  129. num b2=(fb[i]-conj(fb[j]))*r4;
  130. if(i!=j)
  131. {
  132. num c1=(fa[j]+conj(fa[i]));
  133. num c2=(fa[j]-conj(fa[i]))*r2;
  134. num d1=(fb[j]+conj(fb[i]))*r3;
  135. num d2=(fb[j]-conj(fb[i]))*r4;
  136. fa[i]=c1*d1+c2*d2*r5;
  137. fb[i]=c1*d2+c2*d1;
  138. }
  139. fa[j]=a1*b1+a2*b2*r5;
  140. fb[j]=a1*b2+a2*b1;
  141. }
  142. fft(fa,sz);fft(fb,sz);
  143. vector<int> res(need);
  144. for(int i=0;i<need;i++)
  145. {
  146. ll aa=fa[i].x+0.5;
  147. ll bb=fb[i].x+0.5;
  148. ll cc=fa[i].y+0.5;
  149. res[i]=(aa+((bb%m)<<15)+((cc%m)<<30))%m;
  150. }
  151. return res;
  152. }
  153. vector<int> square_mod(vector<int> &a,int m)
  154. {
  155. return multiply_mod(a,a,m,1);
  156. }
  157. };
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement