Advertisement
Guest User

Untitled

a guest
Aug 2nd, 2011
519
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 2.65 KB | None | 0 0
  1. #include <iostream>
  2. #include <stdio.h>
  3. #include <algorithm>
  4. #include <math.h>
  5. #include <stdlib.h>
  6. #include <string>
  7. #include <string.h>
  8. #include <vector>
  9. #include <map>
  10. #include <valarray>
  11. #include <set>
  12. #include <limits.h>
  13. #include <ctime>
  14. #include <stack>
  15. #include <queue>
  16. #include <ctype.h>
  17. #include <complex>
  18. using namespace std;
  19.  
  20. typedef complex<double> base;
  21. const double PI = acos(-1.0);
  22. int R[280000];
  23. base W[280000][2];
  24.  
  25. int rev(int x,const int& l) {
  26.     int sol = 0;
  27.     for(int i=0;i<=l;i++) {
  28.         sol <<= 1;
  29.         sol |= (x&1);
  30.         x >>= 1;
  31.     }
  32.     return sol;
  33. }
  34.  
  35. void fft(vector<base>& a,const int& invert) {
  36.     int n = a.size();
  37.     for(int i=0;i<n;i++) {
  38.         if(i<R[i]) {
  39.             swap(a[i],a[R[i]]);
  40.         }
  41.     }
  42.     for(int p=2;p<=n;(p<<=1)) {
  43.         int wn = n/p;
  44.         for(int i=0;i<n;i+=p) {
  45.             int ind = 0;
  46.             for(int j=0;j<(p>>1);j++) {
  47.                 base u = a[i+j], t = a[i+j+(p>>1)] * W[ind][invert];
  48.                 ind = ((ind+wn) & (n-1));
  49.                 a[i+j] = u + t;
  50.                 a[i+j+(p>>1)] = u - t;
  51.             }
  52.         }
  53.     }
  54.     if (invert)
  55.         for(int i=0;i<n;i++)
  56.             a[i] /= n;
  57. }
  58.  
  59. char s[131172];
  60.  
  61. int main() {
  62.     //freopen("input.txt","r",stdin);
  63.     //freopen("output.txt","w",stdout);
  64.     freopen("a2.in","r",stdin);
  65.     freopen("a2.out","w",stdout);
  66.     //int Time = clock();
  67.     int n;
  68.     scanf("%d\n",&n);
  69.     gets(s);
  70.     vector<int> a(n),b(n);
  71.     for(int i=0;i<n;i++) {
  72.         switch (s[n-i-1]) {
  73.             case ('A'):
  74.                 a[i] = 0;
  75.                 break;
  76.             case ('C'):
  77.                 a[i] = 1;
  78.                 break;
  79.             case ('G'):
  80.                 a[i] = 2;
  81.                 break;
  82.             default:
  83.                 a[i] = 3;
  84.                 break;
  85.         }
  86.     }
  87.     gets(s);
  88.     for(int i=0;i<n;i++) {
  89.         switch (s[n-i-1]) {
  90.             case ('A'):
  91.                 b[i] = 0;
  92.                 break;
  93.             case ('C'):
  94.                 b[i] = 1;
  95.                 break;
  96.             case ('G'):
  97.                 b[i] = 2;
  98.                 break;
  99.             default:
  100.                 b[i] = 3;
  101.                 break;
  102.         }
  103.     }
  104.     int o=0;
  105.     while((1 << o)<n) o++;
  106.     vector<int> m(n,0);
  107.     for(int i=0;i<(n<<1);i++) {
  108.             R[i] = rev (i,o);
  109.     }
  110.     W[0][0] = 1;
  111.     W[0][1] = 1;
  112.     double ang = 2*PI/(n<<1);
  113.     base w(cos(ang),sin(ang));
  114.     W[1][0] = w;
  115.     for(int i=2;i<=(n<<1);i++) {
  116.         W[i][0] = W[i-1][0] * w;
  117.     }
  118.     ang = -2*PI/(n<<1);
  119.     base wn(cos(ang),sin(ang));
  120.     W[1][1] = wn;
  121.     for(int i=2;i<=(n<<1);i++) {
  122.         W[i][1] = W[i-1][1] * wn;
  123.     }
  124.     for(int p=0;p<4;p++) {
  125.         vector<base> x(n<<1,0),y(n<<1,0);
  126.         for(int i=0;i<n;i++) {
  127.             x[n-i-1] = (a[i] == p) ? 1 : 0;
  128.             y[i] = (b[i] == p) ? 1 : 0;
  129.             y[i+n] = y[i];
  130.         }
  131.         fft(x,0);
  132.         fft(y,0);
  133.         for(int i=0;i<(n<<1);i++) {
  134.             x[i] *= y[i];
  135.         }
  136.         fft(x,1);
  137.         for(int i = 0; i<n; i++) {
  138.             m[i] += int(x[n-1+i].real()+0.5);
  139.         }
  140.     }
  141.     int ps = 0;
  142.     for(int i=1;i<n;i++) {
  143.         if(m[i] > m[ps]) ps = i;
  144.     }
  145.     printf("%d %d\n",m[ps],ps);
  146.     //printf("%.3f",float(clock()-Time)/1000);
  147.     return 0;
  148. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement