Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <vector>
- #include <algorithm>
- #include <complex>
- #include <math.h>
- using namespace std;
- struct DNA {
- DNA(int M) : A(M, 0), C(M, 0), G(M, 0), T(M, 0) {}
- vector<complex<long double>> A, C, G, T;
- };
- void dfft(vector<complex<long double>> &vec, bool inv) {
- int sz = vec.size(), bit;
- for (int i = 1, j = 0; i < sz; i++) {
- bit = sz >> 1;
- for (; j >= bit; bit >>= 1)
- j -= bit;
- j += bit;
- if (i < j)
- swap(vec[i], vec[j]);
- }
- for (int len = 2; len <= sz; len <<= 1) {
- long double ang = 2 * M_PI / len * (inv ? -1 : 1);
- complex<long double> wlen(cos(ang), sin(ang));
- for (int i = 0; i < sz; i += len) {
- complex<long double> w(1);
- for (int j = 0; j < len / 2; j++) {
- complex<long double> u = vec[i + j], v = w * vec[i + j + (len >> 1)];
- vec[i + j] = u + v;
- vec[i + j + (len >> 1)] = u - v;
- w *= wlen;
- }
- }
- }
- if (inv)
- for (int i = 0; i < sz; i++)
- vec[i] /= sz;
- }
- void parse(const string &s, DNA &d) {
- int t = 0;
- for (auto it = s.begin(); it != s.end(); it++) {
- switch (*it) {
- case 'A':
- d.A[t] = 1;
- ++t;
- break;
- case 'C':
- d.C[t] = 1;
- ++t;
- break;
- case 'G':
- d.G[t] = 1;
- ++t;
- break;
- case 'T':
- d.T[t] = 1;
- ++t;
- break;
- }
- }
- }
- int main() {
- int M;
- string s1, s2;
- cin >> M >> s1 >> s2;
- reverse(s1.begin(), s1.end());
- DNA first(2 * M), second(2 * M);
- parse(s1, first), parse(s2, second);
- vector<complex<long double>> a_mul(2 * M, 0), c_mul(2 * M, 0), g_mul(2 * M, 0), t_mul(2 * M, 0);
- vector<unsigned long long> a_real(2 * M, 0), c_real(2 * M, 0), g_real(2 * M, 0), t_real(2 * M, 0);
- vector<unsigned long long> result(M, 0);
- for (int i = 0; i < M; i++) {
- second.A[M + i] = second.A[i];
- second.C[M + i] = second.C[i];
- second.G[M + i] = second.G[i];
- second.T[M + i] = second.T[i];
- }
- dfft(first.A, false), dfft(second.A, false);
- dfft(first.C, false), dfft(second.C, false);
- dfft(first.G, false), dfft(second.G, false);
- dfft(first.T, false), dfft(second.T, false);
- for (int i = 0; i < 2 * M; i++) {
- a_mul[i] = first.A[i] * second.A[i];
- c_mul[i] = first.C[i] * second.C[i];
- g_mul[i] = first.G[i] * second.G[i];
- t_mul[i] = first.T[i] * second.T[i];
- }
- dfft(a_mul, true); dfft(c_mul, true);
- dfft(g_mul, true); dfft(t_mul, true);
- for (int i = 0; i < 2 * M; i++) {
- a_real[i] = (unsigned long long int)(a_mul[i].real() + 0.5);
- c_real[i] = (unsigned long long int)(c_mul[i].real() + 0.5);
- g_real[i] = (unsigned long long int)(g_mul[i].real() + 0.5);
- t_real[i] = (unsigned long long int)(t_mul[i].real() + 0.5);
- }
- for (int i = M; i < 2 * M; i++) {
- unsigned long long int t = a_real[i] + c_real[i] + g_real[i] + t_real[i];
- result[i - M] = t;
- }
- cout << *max_element(result.begin(), result.end()) << ' ';
- cout << ((result.end() - max_element(result.begin(), result.end()) - 1) % M);
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement