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;
}
}
}