It is often a necessity to compare some sequences to find out how similar they are. There are many similarity measures that can be used, e.g., longest common subsequence, edit distance, sequence alignment. Recently a merged longest common subsequence (MergedLCS) problem was formulated with applications in bioinformatics. We propose the bit-parallel algorithms for the MergedLCS problem and evaluate them in practice showing that they are usually tens times faster than the already published methods.
int BMLCS_Bitpar(vector< int> T, vector< int> A, vector< int> B, int Tlen, int Alen, int Blen, int size, int WordSize, vector< int> EOBA, vector< int> EOBB) { int space = ceil((Tlen + 1) / WordSize + 1); vector< vector< uint64_t>> W1(Blen + 1, vector< uint64_t>(space)); vector< vector< uint64_t>> W2(Blen + 1, vector< uint64_t>(space)); vector< vector< uint64_t>> Y(size, vector< uint64_t>(space, 0)); uint64_t MaxNUM64 = 0XFFFFFFFFFFFFFFFF;//2^64 - 1 uint64_t modLastWord = 0X0000020000000000;//2^42 int WordNum = space - 1; int WordNumMinus = WordNum - 1; vector< int> FEOBA(Alen + 1, 0); vector< int> FEOBB(Blen + 1, 0); //construct FEOBA and FEOBB int sizeEOBA = EOBA.size(); int sizeEOBB = EOBB.size(); for (int i = 0; i < sizeEOBA; ++i) { FEOBA[EOBA[i]] = 1; } FEOBA[0] = 1; FEOBA[Alen] = 1; for (int i = 0; i < sizeEOBB; ++i) { FEOBB[EOBB[i]] = 1; } FEOBB[0] = 1; FEOBB[Blen] = 1; //Constructing array of Y uint64_t shiftNum = 1; for (int i = 1; i < WordSize; ++i) { shiftNum = shiftNum << 1; Y[T[i] - 1][0] |= shiftNum; } if (Tlen > 63) for (int i = WordSize; i <= Tlen; ++i) { shiftNum = shiftNum << 1; if (i % WordSize == 0) shiftNum = 1; Y[T[i] - 1][floor(i / WordSize)] |= shiftNum; } //Initialisation for (int i = 0; i <= WordNum; ++i) { W2[0][i] = MaxNUM64; } W2[0][WordNum] %= modLastWord; //Calculating boundaries uint64_t U = 0, tmpa = 0; uint64_t carrybitW = 0; for (int k = 1; k <= Blen; ++k) { carrybitW = 0; for (int i = 0; i <= WordNumMinus; ++i) { tmpa = W2[k - 1][i]; U = tmpa & Y[B[k] - 1][i]; W2[k][i] = (tmpa + U + carrybitW) | (tmpa - U); //detect overflow if (U > MaxNUM64 - tmpa) carrybitW = 1; else carrybitW = 0; } tmpa = W2[k - 1][WordNum]; U = tmpa& Y[B[k] - 1][WordNum]; W2[k][WordNum] = ((tmpa + U + carrybitW) % modLastWord) | (tmpa - U); } int layer = 1; while (1) { if (layer % 2 == 1) { carrybitW = 0; for (int i = 0; i <= WordNumMinus; ++i) { tmpa = W2[0][i]; U = tmpa & Y[A[layer] - 1][i]; W1[0][i] = (tmpa + U + carrybitW) | (tmpa - U); //detect overflow if (U > MaxNUM64 - tmpa) carrybitW = 1; else carrybitW = 0; } tmpa = W2[0][WordNum]; U = tmpa & Y[A[layer] - 1][WordNum]; W1[0][WordNum] = ((tmpa + U + carrybitW) % modLastWord) | (tmpa - U); //Main calculations uint64_t Up = 0, Upp = 0, Wp = 0, Wpp = 0, Ut = 0, Wt = 0, Vt = 0, carrybitWp = 0, carrybitWpp = 0, Ws = 0; int za = 0, zb = 0; int ta = 0, tb = 0; uint64_t one = 1; for (int k = 1; k <= Blen; ++k) { carrybitWp = 0, carrybitWpp = 0; za = 0; zb = 0; ta = FEOBA[layer]; tb = FEOBB[k]; if (ta + tb == 0) continue; for (int i = 0; i <= WordNumMinus; ++i) { if (tb == 1) { tmpa = W2[k][i]; Up = tmpa & Y[A[layer] - 1][i]; Wp = (tmpa + Up + carrybitWp) | (tmpa - Up); if (Up > MaxNUM64 - tmpa) carrybitWp = 1; else carrybitWp = 0; } if (ta == 1) { tmpa = W1[k - 1][i]; Upp = tmpa & Y[B[k] - 1][i]; Wpp = (tmpa + Upp + carrybitWpp) | (tmpa - Upp); if (Upp > MaxNUM64 - tmpa) carrybitWpp = 1; else carrybitWpp = 0; } if (ta == 1 && tb == 0) { W1[k][i] = Wpp; } else if (ta == 0 && tb == 1) { W1[k][i] = Wp; } else if (ta == 1 && tb == 1) { Ut = Wp & Wpp; Vt = Wp ^ Wpp; Ws = MaxNUM64; if (i == 0) for (int tmpi = 1; tmpi < WordSize; ++tmpi) { if ((Vt >> tmpi) % 2 == 1) { if ((Wp >> tmpi) % 2 == 0) { ++za; if (za > zb) Ws &= (~(one << tmpi)); } else { ++zb; if (zb > za) Ws &= (~(one << tmpi)); } } } else { for (int tmpi = 0; tmpi < WordSize; ++tmpi) { if ((Vt >> tmpi) % 2 == 1) { if ((Wp >> tmpi) % 2 == 0) { ++za; if (za > zb) Ws &= (~(one << tmpi)); } else { ++zb; if (zb > za) Ws &= (~(one << tmpi)); } } } } W1[k][i] = ((Wp | Wpp) & Ws) | Ut; } } if (tb == 1) { tmpa = W2[k][WordNum]; Up = tmpa & Y[A[layer] - 1][WordNum]; Wp = ((tmpa + Up + carrybitWp) % modLastWord) | (tmpa - Up); } if (ta == 1) { tmpa = W1[k - 1][WordNum]; Upp = tmpa & Y[B[k] - 1][WordNum]; Wpp = ((tmpa + Upp + carrybitWpp) % modLastWord) | (tmpa - Upp); } if (ta == 1 && tb == 0) { W1[k][WordNum] = Wpp; } else if (ta == 0 && tb == 1) { W1[k][WordNum] = Wp; } else if (ta == 1 && tb == 1) { Ut = Wp & Wpp; Vt = Wp ^ Wpp; Ws = MaxNUM64 % modLastWord; for (int tmpi = 0; tmpi < 41; ++tmpi) { if ((Vt >> tmpi) % 2 == 1) { if ((Wp >> tmpi) % 2 == 0) { ++za; if (za > zb) Ws &= (~(one << tmpi)); } else { ++zb; if (zb > za) Ws &= (~(one << tmpi)); } } } W1[k][WordNum] = (((Wp | Wpp) & Ws) | Ut) % modLastWord; } } if (layer == Alen) { int l = 0; uint64_t Vz = 0; for (int i = 0; i <= WordNum - 1; ++i) { Vz = ~W1[Blen][i]; while (Vz != 0) { Vz = Vz & (Vz - 1); ++l; } } Vz = (~W1[Blen][WordNum]) % modLastWord; while (Vz != 0) { Vz = Vz & (Vz - 1); ++l; } vector< vector< uint64_t>>().swap(W1); vector< vector< uint64_t>>().swap(Y); vector< int>().swap(FEOBA); vector< int>().swap(FEOBB); return l; } ++layer; } else { carrybitW = 0; for (int i = 0; i <= WordNumMinus; ++i) { tmpa = W1[0][i]; U = tmpa & Y[A[layer] - 1][i]; W2[0][i] = (tmpa + U + carrybitW) | (tmpa - U); //detect overflow if (U > MaxNUM64 - tmpa) carrybitW = 1; else carrybitW = 0; } tmpa = W1[0][WordNum]; U = tmpa & Y[A[layer] - 1][WordNum]; W2[0][WordNum] = ((tmpa + U + carrybitW) % modLastWord) | (tmpa - U); //Main calculations uint64_t Up = 0, Upp = 0, Wp = 0, Wpp = 0, Ut = 0, Wt = 0, Vt = 0, carrybitWp = 0, carrybitWpp = 0, Ws = 0; int za = 0, zb = 0; int ta = 0, tb = 0; uint64_t one = 1; for (int k = 1; k <= Blen; ++k) { carrybitWp = 0, carrybitWpp = 0; za = 0; zb = 0; ta = FEOBA[layer]; tb = FEOBB[k]; if (ta + tb == 0) continue; for (int i = 0; i <= WordNumMinus; ++i) { if (tb == 1) { tmpa = W1[k][i]; Up = tmpa & Y[A[layer] - 1][i]; Wp = (tmpa + Up + carrybitWp) | (tmpa - Up); if (Up > MaxNUM64 - tmpa) carrybitWp = 1; else carrybitWp = 0; } if (ta == 1) { tmpa = W2[k - 1][i]; Upp = tmpa & Y[B[k] - 1][i]; Wpp = (tmpa + Upp + carrybitWpp) | (tmpa - Upp); if (Upp > MaxNUM64 - tmpa) carrybitWpp = 1; else carrybitWpp = 0; } if (ta == 1 && tb == 0) { W2[k][i] = Wpp; } else if (ta == 0 && tb == 1) { W2[k][i] = Wp; } else if (ta == 1 && tb == 1) { Ut = Wp & Wpp; Vt = Wp ^ Wpp; Ws = MaxNUM64; if (i == 0) for (int tmpi = 1; tmpi < WordSize; ++tmpi) { if ((Vt >> tmpi) % 2 == 1) { if ((Wp >> tmpi) % 2 == 0) { ++za; if (za > zb) Ws &= (~(one << tmpi)); } else { ++zb; if (zb > za) Ws &= (~(one << tmpi)); } } } else { for (int tmpi = 0; tmpi < WordSize; ++tmpi) { if ((Vt >> tmpi) % 2 == 1) { if ((Wp >> tmpi) % 2 == 0) { ++za; if (za > zb) Ws &= (~(one << tmpi)); } else { ++zb; if (zb > za) Ws &= (~(one << tmpi)); } } } } W2[k][i] = ((Wp | Wpp) & Ws) | Ut; } } if (tb == 1) { tmpa = W1[k][WordNum]; Up = tmpa & Y[A[layer] - 1][WordNum]; Wp = ((tmpa + Up + carrybitWp) % modLastWord) | (tmpa - Up); } if (ta == 1) { tmpa = W2[k - 1][WordNum]; Upp = tmpa & Y[B[k] - 1][WordNum]; Wpp = ((tmpa + Upp + carrybitWpp) % modLastWord) | (tmpa - Upp); } if (ta == 1 && tb == 0) { W2[k][WordNum] = Wpp; } else if (ta == 0 && tb == 1) { W2[k][WordNum] = Wp; } else if (ta == 1 && tb == 1) { Ut = Wp & Wpp; Vt = Wp ^ Wpp; Ws = MaxNUM64 % modLastWord; for (int tmpi = 0; tmpi < 41; ++tmpi) { if ((Vt >> tmpi) % 2 == 1) { if ((Wp >> tmpi) % 2 == 0) { ++za; if (za > zb) Ws &= (~(one << tmpi)); } else { ++zb; if (zb > za) Ws &= (~(one << tmpi)); } } } W2[k][WordNum] = (((Wp | Wpp) & Ws) | Ut) % modLastWord; } } if (layer == Alen) { int l = 0; uint64_t Vz = 0; for (int i = 0; i <= WordNum - 1; ++i) { Vz = ~W2[Blen][i]; while (Vz != 0) { Vz = Vz & (Vz - 1); ++l; } } Vz = (~W2[Blen][WordNum]) % modLastWord; while (Vz != 0) { Vz = Vz & (Vz - 1); ++l; } vector< vector< uint64_t>>().swap(W2); vector< vector< uint64_t>>().swap(Y); vector< int>().swap(FEOBA); vector< int>().swap(FEOBB); return l; } ++layer; } } }