Cod sursa(job #2418187)

Utilizator stoianmihailStoian Mihail stoianmihail Data 4 mai 2019 01:59:29
Problema Potrivirea sirurilor Scor 40
Compilator cpp-64 Status done
Runda Arhiva educationala Marime 11.8 kb
#include <fstream>
#include <cmath>
#include <cassert>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <vector>

using namespace std;

#define MAX_LEN 2000000

const int SIGMA = 26;
const int ALL = 1 << 8;
int alphaTable[ALL + 1];

void init() {
  for (char c = 'a'; c <= 'z'; c++)
    alphaTable[c] = c - 'a' + 1;
  for (char c = 'A'; c <= 'Z'; c++)
    alphaTable[c] = c - 'A' + SIGMA + 1;
  for (char c = '0'; c <= '9'; c++)
    alphaTable[c] = c - '0' + SIGMA * 2 + 1;
}

inline int encode(char c) {
  return alphaTable[c];
}

const double PI = acos(0) * 2;

class complex {
public:
	double a, b;
	complex() {a = 0.0; b = 0.0;}
	complex(double na, double nb) {a = na; b = nb;}
	const complex operator+(const complex &c) const
		{return complex(a + c.a, b + c.b);}
  void operator += (const complex &c) {
    *this = *this + c;
  }
	const complex operator-(const complex &c) const
		{return complex(a - c.a, b - c.b);}
	const complex operator*(const complex &c) const
		{return complex(a*c.a - b*c.b, a*c.b + b*c.a);}
  void operator *= (const complex &c) {
    *this = *this * c;
  }

	double magnitude() {return sqrt(a*a+b*b);}
	void print() {
    printf("(%.3f %.3f)\n", a, b);
	}
};

complex *roots;

class mirrorBits {
  #define MAX_SIZE (1 << 22)
private:
  unsigned size;
  unsigned *rev;
  char* lg;

  void compute() {
    int currLg = 0, mask = 0;
    rev[1] = 1;
    lg[1] = 0;
    for (unsigned i = 2; i <= size; i++) {
      if (!(i & (i - 1))) {
        ++currLg;
        mask = (mask << 1) | 1;
      }
      // nimm die Bits nach dem Bit von log und shifte danach mit wieviel Nullen dazwischen vorliegen.
      rev[i] = (rev[mask & i] << (currLg - lg[mask & i])) | 1;
      lg[i] = currLg;
    }
    /*for (unsigned i = 0; i < 20; ++i)
      cerr << i << " : " << rev[i] << " " << (int)lg[i] << endl;
    */
  }
public:
  mirrorBits() {
    this -> size = MAX_SIZE;
    rev = new unsigned[size + 1]();
    lg = new char[size + 1]();
    compute();
  }

  // TODO: speichere falls moeglich alles in rev.
  inline int getRev(int space, int x) {
    //cerr << space << " mit " << x << " aber rev[x] = " << rev[x] << " ret = " << (rev[x] << (space - lg[x] - 1)) << endl;
    return rev[x] << (space - lg[x] - 1);
  }

  inline int getLog(int x) {
    return lg[x];
  }

  inline int getPower(int x) {
    return (!(x & (x - 1))) ? x : (1 << (lg[x] + 1));
  }
} spiegel;

inline void swap(complex& a, complex& b) {
  complex tmp = a; a = b; b = tmp;
}

unsigned power(char val, int pow) {
  val = encode(val);
  unsigned ret = 1;
  while (pow--)
    ret *= val;
  return ret;
}

// reverse the positions as shown in emaxx.
void reversePositions(complex** transformation, unsigned log) {
  unsigned size = 1 << log;
  for (unsigned index = 0; index < size; ++index) {
    unsigned moveTo = spiegel.getRev(log, index);
    if (index < moveTo) {
      swap(transformation[index][0], transformation[moveTo][0]);
      swap(transformation[index][1], transformation[moveTo][1]);
    }
  }
}

void transform(complex** data, unsigned log, bool inverse = false) {
  int i, j, k, n = 1 << log;
  for (i = 1; i <= log; i++) {
    int m = (1 << i), md2 = m / 2;
    int start = 0, increment = (1 << (log - i));
    if (inverse) {
      start = n;
      increment *= -1;
    }
    complex t;
    for (k = 0; k < n; k += m) {
      int index = start;
      for (j = k; j < md2+k; j++) {
        // for normal polynomial
        t = roots[index] * data[j+md2][0];
        data[j+md2][0] = data[j][0] - t;
        data[j][0] = data[j][0] + t;

        // for polynomial ^ 2.
        t = roots[index] * data[j+md2][1];
        data[j+md2][1] = data[j][1] - t;
        data[j][1] = data[j][1] + t;
        index += increment;
      }
    }
  }
  if (inverse)
    for (i = 0; i < n; i++) {
      for (j = 0; j < 2; j++) {
        data[i][j].a /= n;
        data[i][j].b /= n;
      }
    }
}

unsigned logOfPow(unsigned x) {
  unsigned log = 0;
  while ((1 << log) < x)
    ++log;
  return log;
}

char pattern[MAX_LEN + 1], text[MAX_LEN + 1];
unsigned *partialSums, totalSum;

unsigned reverseSum(unsigned pos) {
  assert(pos);
  return totalSum - partialSums[pos - 1];
}

int rou(complex x) {
  return (int)(x.a + 0.5);
}

void print(complex **trans, int part, unsigned size) {
  cerr << "Print trans\n";
  for (unsigned i = 0; i < size; ++i) {
    fprintf(stderr, "(%.3f, %.3f) ", trans[i][part].a, trans[i][part].b);
  }
  cerr << endl;
}

int main(void) {
  ifstream cin("strmatch.in");
  ofstream cout("strmatch.out");

  init();

  cin >> pattern >> text;
  unsigned patternLen = strlen(pattern);
  unsigned textLen = strlen(text);

  cerr << patternLen << " mit " << textLen << endl;

  if (patternLen > textLen) {
    cout << 0 << endl;
    return 0;
  }
  unsigned logOfLength = logOfPow(3 * patternLen - 1);
  unsigned length = 1 << logOfLength;

  //cerr << "logOfLength = " << logOfLength << endl;

  roots = new complex[length + 1]();
  roots[0] = complex(1, 0);
  complex mult = complex(cos(2 * PI / length), sin(2 * PI / length));
  for (unsigned i = 1; i <= length; i++)
    roots[i] = roots[i - 1] * mult;

  complex **patternTransformation;
  patternTransformation = new complex*[1 << logOfLength];

  // positions of pattern should be mirrored.
  partialSums = new unsigned[patternLen];
  for (unsigned index = 0; index < patternLen; ++index) {
    //cerr << power(pattern[index], 3);
    partialSums[index] = power(pattern[index], 3) + ((!index) ? 0 : partialSums[index - 1]);

    patternTransformation[index] = new complex[2]();
    patternTransformation[index][0] = complex(power(pattern[patternLen - index - 1], 1), 0);
    patternTransformation[index][1] = complex(power(pattern[patternLen - index - 1], 2), 0);
  }
  totalSum = partialSums[patternLen - 1];

  //cerr << "totalSum = " << totalSum << endl;

  // Alloc for the rest of the vector.
  for (unsigned index = patternLen; index < length; ++index)
    patternTransformation[index] = new complex[2]();

  reversePositions(patternTransformation, logOfLength);

  //print(patternTransformation, 0, length);
  //print(patternTransformation, 1, length);

  transform(patternTransformation, logOfLength);

  cerr << "transform of pattern!" << endl;

  complex **textTransformation = new complex*[1 << logOfLength];
  for (unsigned index = 0; index < length; ++index)
    textTransformation[index] = new complex[2]();

  vector <int> matches;
  unsigned countOfMatches = 0;
  unsigned currIndex = 0;
  unsigned lastCheck = 0;
  unsigned blockSize = (1 << logOfLength) - patternLen + 1;
  unsigned *stack = new unsigned[patternLen];
  while (currIndex < textLen) {
    cerr << "lastCheck = " << lastCheck << " and currIndex = " << currIndex << endl;
    if (lastCheck == 1) {
      cerr << "enter the most dna" << endl;
      // check the string that should be continued!
      // At this point we should save all the pattern.
      int val = patternLen - stack[0] - 1;
      int k = val;
      unsigned start = 0;
      while (k < patternLen && text[currIndex + start++] == pattern[k])
        ++k;
      if (k == patternLen) {
        if (++countOfMatches <= 1000)
          matches.push_back(currIndex + stack[0]);
        //cerr << "found";
        //goto found;
      }
      lastCheck = 0;
    }
    // because we let the first if to enter this one, if lastCheck gets changed.
    if (!lastCheck) {
      cerr << "push the boundaries " << endl;
      unsigned index = currIndex;
      while (index < textLen && text[index] != pattern[0])
        index++;
      if (index == textLen) {
        goto terminate;
      } else {
        currIndex = index;
      }
      //cerr << "now is currIndex = " << currIndex << endl;
    }

    if (currIndex + blockSize > textLen)
      blockSize = textLen - currIndex;

    //cerr << "current blocksize = " << blockSize << endl;

    // if it is in free modus, and the rest of the text isn't enough.

    //cerr << "blockSize = " << blockSize << endl;

    // dies soll immer gecleared werden.
    for (unsigned index = 0; index < length; ++index) {
      textTransformation[index][0] = complex(0, 0);
      textTransformation[index][1] = complex(0, 0);
    }

    // Compute FFT from currIndex.
    for (unsigned index = 0; index < blockSize; ++index) {
      unsigned moveTo = spiegel.getRev(logOfLength, index);
      if (index < moveTo) {
        textTransformation[index][0] = complex((moveTo < blockSize) ? power(text[currIndex + moveTo], 1) : 0, 0);
        textTransformation[moveTo][0] = complex(power(text[currIndex + index], 1), 0);
        textTransformation[index][1] = complex((moveTo < blockSize) ? power(text[currIndex + moveTo], 2) : 0, 0);
        textTransformation[moveTo][1] = complex(power(text[currIndex + index], 2), 0);
      } else if (index == moveTo) {
        textTransformation[index][0] = complex(power(text[currIndex + index], 1), 0);
        textTransformation[index][1] = complex(power(text[currIndex + index], 2), 0);
      }
    }
    transform(textTransformation, logOfLength);

    cerr << "transofmat of text\n";

    // The order here is quite important.
    for (unsigned i = 0; i < length; ++i) {
      textTransformation[i][0] *= patternTransformation[i][1];
      textTransformation[i][1] *= patternTransformation[i][0];
    }

    // combine this with the convolution.
    reversePositions(textTransformation, logOfLength);
    transform(textTransformation, logOfLength, true);

    cerr << "reverese taernawe\n";

    //print(textTransformation, 0, length);
    //print(textTransformation, 1, length);

    // Don't forget to return the leftmost position!
    if (lastCheck) {
      for (unsigned i = 0; i < lastCheck; ++i) {
        unsigned val = stack[i];
        //cerr << "in lastCheck\n";
        //cerr << rou(textTransformation[val][1]) << "+" << reverseSum(patternLen - val - 1) << " vs " << (2 * rou(textTransformation[val][0])) << endl;
        if (rou(textTransformation[val][1]) + reverseSum(patternLen - val - 1) == 2 * rou(textTransformation[val][0])) {
          if (++countOfMatches <= 1000)
            matches.push_back(currIndex + val);
        }
      }
    }
    if (blockSize < patternLen) {
      goto terminate;
    }

    cerr << "hier arrives with blockSize = " << blockSize << endl;
    for (int i = 0; i < blockSize - patternLen + 1; ++i) {
      //cerr << (rou(textTransformation[i + patternLen - 1][1])) << " + " << totalSum << " vs " << (2 * rou(textTransformation[i + patternLen - 1][0])) << endl;
      if (rou(textTransformation[i + patternLen - 1][1]) + totalSum == 2 * rou(textTransformation[i + patternLen - 1][0])) {
        //cerr << "found normal!" << (currIndex + i) << endl;
        if (++countOfMatches <= 1000)
          matches.push_back(currIndex + i + patternLen - 1);
      }
    }
    cerr << "over the easy part\n";
    unsigned goBack = patternLen - 2;
    lastCheck = 0;
    for (int i = blockSize - patternLen + 1; i < blockSize; ++i) {
      if (rou(textTransformation[i + patternLen - 1][1]) + partialSums[goBack] == 2 * rou(textTransformation[i + patternLen - 1][0])) {
        stack[lastCheck++] = i - blockSize + patternLen - 1;
        //cerr << "push on stack " << (i - blockSize + patternLen - 1) << endl;
      }
      goBack--;
    }
    currIndex += blockSize;
  }
  terminate : {
    //cerr << "No of matches = " << matches << endl;
    cout << countOfMatches << endl;
    if (matches.size() < countOfMatches) {
      countOfMatches = matches.size();
    }
    for (unsigned i = 0; i < countOfMatches; ++i)
      cout << (matches[i] - patternLen + 1) << " ";
    return 0;
  }
  noFound : {
    cerr << "not found!";
    return 0;
  }
  return 0;
}