Advertisement
Guest User

Untitled

a guest
Apr 20th, 2019
82
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.45 KB | None | 0 0
  1. #include <iostream>
  2. #include <vector>
  3. #include <algorithm>
  4. #include <complex>
  5. #include <math.h>
  6. using namespace std;
  7.  
  8. struct DNA {
  9. DNA(int M) : A(M, 0), C(M, 0), G(M, 0), T(M, 0) {}
  10.  
  11. vector<complex<long double>> A, C, G, T;
  12. };
  13.  
  14. void dfft(vector<complex<long double>> &vec, bool inv) {
  15. int sz = vec.size(), bit;
  16. for (int i = 1, j = 0; i < sz; i++) {
  17. bit = sz >> 1;
  18. for (; j >= bit; bit >>= 1)
  19. j -= bit;
  20. j += bit;
  21. if (i < j)
  22. swap(vec[i], vec[j]);
  23. }
  24. for (int len = 2; len <= sz; len <<= 1) {
  25. long double ang = 2 * M_PI / len * (inv ? -1 : 1);
  26. complex<long double> wlen(cos(ang), sin(ang));
  27. for (int i = 0; i < sz; i += len) {
  28. complex<long double> w(1);
  29. for (int j = 0; j < len / 2; j++) {
  30. complex<long double> u = vec[i + j], v = w * vec[i + j + (len >> 1)];
  31. vec[i + j] = u + v;
  32. vec[i + j + (len >> 1)] = u - v;
  33. w *= wlen;
  34. }
  35. }
  36. }
  37. if (inv)
  38. for (int i = 0; i < sz; i++)
  39. vec[i] /= sz;
  40. }
  41.  
  42. void parse(const string &s, DNA &d) {
  43. int t = 0;
  44.  
  45. for (auto it = s.begin(); it != s.end(); it++) {
  46. switch (*it) {
  47. case 'A':
  48. d.A[t] = 1;
  49. ++t;
  50. break;
  51. case 'C':
  52. d.C[t] = 1;
  53. ++t;
  54. break;
  55. case 'G':
  56. d.G[t] = 1;
  57. ++t;
  58. break;
  59. case 'T':
  60. d.T[t] = 1;
  61. ++t;
  62. break;
  63. }
  64. }
  65. }
  66.  
  67. int main() {
  68. int M;
  69. string s1, s2;
  70. cin >> M >> s1 >> s2;
  71. reverse(s1.begin(), s1.end());
  72. DNA first(2 * M), second(2 * M);
  73. parse(s1, first), parse(s2, second);
  74.  
  75. vector<complex<long double>> a_mul(2 * M, 0), c_mul(2 * M, 0), g_mul(2 * M, 0), t_mul(2 * M, 0);
  76. vector<unsigned long long> a_real(2 * M, 0), c_real(2 * M, 0), g_real(2 * M, 0), t_real(2 * M, 0);
  77.  
  78. vector<unsigned long long> result(M, 0);
  79.  
  80. for (int i = 0; i < M; i++) {
  81. second.A[M + i] = second.A[i];
  82. second.C[M + i] = second.C[i];
  83. second.G[M + i] = second.G[i];
  84. second.T[M + i] = second.T[i];
  85. }
  86.  
  87. dfft(first.A, false), dfft(second.A, false);
  88. dfft(first.C, false), dfft(second.C, false);
  89. dfft(first.G, false), dfft(second.G, false);
  90. dfft(first.T, false), dfft(second.T, false);
  91.  
  92. for (int i = 0; i < 2 * M; i++) {
  93. a_mul[i] = first.A[i] * second.A[i];
  94. c_mul[i] = first.C[i] * second.C[i];
  95. g_mul[i] = first.G[i] * second.G[i];
  96. t_mul[i] = first.T[i] * second.T[i];
  97. }
  98.  
  99. dfft(a_mul, true); dfft(c_mul, true);
  100. dfft(g_mul, true); dfft(t_mul, true);
  101.  
  102. for (int i = 0; i < 2 * M; i++) {
  103. a_real[i] = (unsigned long long int)(a_mul[i].real() + 0.5);
  104. c_real[i] = (unsigned long long int)(c_mul[i].real() + 0.5);
  105. g_real[i] = (unsigned long long int)(g_mul[i].real() + 0.5);
  106. t_real[i] = (unsigned long long int)(t_mul[i].real() + 0.5);
  107. }
  108.  
  109. for (int i = M; i < 2 * M; i++) {
  110. unsigned long long int t = a_real[i] + c_real[i] + g_real[i] + t_real[i];
  111. result[i - M] = t;
  112. }
  113.  
  114. cout << *max_element(result.begin(), result.end()) << ' ';
  115. cout << ((result.end() - max_element(result.begin(), result.end()) - 1) % M);
  116. return 0;
  117. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement