﻿ Deorowicz and Danek Algorithm_MLCS

# Deorowicz and Danek Algorithm_MLCS(Bit-Parallel, 2013)

Main features
• searching phase in O(⌈r / wmn logw) time complexity
Abstract of the paper [Deorowicz 2013]

It is often a necessity to compare some sequences to ﬁnd 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.

C++ source code
```int MLCS_Bitpar(vector< int> T, vector< int> A, vector< int> B, int Tlen, int Alen, int Blen, int size, int WordSize) {
vector< vector< uint64_t> > W1(Blen + 1, vector< uint64_t>(ceil((Tlen + 1) / WordSize + 1)));
vector< vector< uint64_t> > W2(Blen + 1, vector< uint64_t>(ceil((Tlen + 1) / WordSize + 1)));
vector< vector< uint64_t> > Y(size, vector< uint64_t>(ceil((Tlen + 1) / WordSize + 1), 0));
uint64_t carrybitW = 0;
uint64_t MaxNUM64 = 0XFFFFFFFFFFFFFFFF;//2^64 - 1
uint64_t modLastWord = 0X0000020000000000;//2^42
//uint64_t modLastWord = 0X0000000000000010;
int WordNum = ceil((Tlen + 1) / WordSize + 1) - 1;

//Constructing array of Y

uint64_t shiftNum = 1;
for (int i = 1; i <= min(Tlen, WordSize - 1); ++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;
for (int k = 1; k <= Blen; ++k) {
carrybitW = 0;
for (int i = 0; i <= WordNum - 1; ++i)
{
tmpa = W2[k - 1][i];
U = tmpa & Y[B[k] - 1][i];
W2[k][i] = (tmpa + U + carrybitW) | (W2[k - 1][i] - 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);
//detect overflow
if (U > MaxNUM64 - tmpa)
carrybitW = 1;
else
carrybitW = 0;
}

int layer = 1;
while (1) {
if (layer % 2 == 1) {
carrybitW = 0;
for (int i = 0; i <= WordNum - 1; ++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);
//detect overflow
if (U > MaxNUM64 - tmpa)
carrybitW = 1;
else
carrybitW = 0;

//Main calculations
uint64_t Up = 0, Upp = 0, Wp = 0, Wpp = 0, Ut = 0, Wt = 0, Vt = 0, carrybitWp = 0, carrybitWpp = 0;
int tmpv = ceil(log2(WordSize)) - 1;
for (int k = 1; k <= Blen; ++k) {
carrybitWp = 0, carrybitWpp = 0;
bool f = false;
for (int i = 0; i <= WordNum - 1; ++i)
{
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;

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;

Ut = Wp | Wpp;
Wt = Wp ^ Wpp;
Vt = Wt;

for (int ipp = 0; ipp <= tmpv; ++ipp) {
Vt = Vt ^ (Vt << (1 << ipp));
}
if (f == true) Vt = ~Vt;
if ((Vt & 0X8000000000000000) > 0) f = true;
else f = false;
W1[k][i] = ~(Wt & Vt) & Ut;
}
tmpa = W2[k][WordNum];
Up = tmpa & Y[A[layer] - 1][WordNum];
Wp = ((tmpa + Up + carrybitWp) % modLastWord) | (tmpa - Up);
if (Up > MaxNUM64 - tmpa)
carrybitWp = 1;
else
carrybitWp = 0;

tmpa = W1[k - 1][WordNum];
Upp = tmpa & Y[B[k] - 1][WordNum];
Wpp = ((tmpa + Upp + carrybitWpp) % modLastWord) | (tmpa - Upp);
if (Upp > MaxNUM64 - tmpa)
carrybitWpp = 1;
else
carrybitWpp = 0;

Ut = Wp | Wpp;
Wt = Wp ^ Wpp;
Vt = Wt;

for (int ipp = 0; ipp <= tmpv; ++ipp) {
Vt = Vt ^ (Vt << (1 << ipp));
}
if (f == true) Vt = ~Vt;
if ((Vt & 0X8000000000000000) > 0) f = true;
else f = false;
W1[k][WordNum] = ~(Wt & Vt) & Ut;
}

if (layer == Alen)
{
int z = 0;
uint64_t Vz = 0;
for (int i = 0; i <= WordNum - 1; ++i) {
Vz = ~W1[Blen][i];
while (Vz != 0) {
Vz = Vz & (Vz - 1);
++z;
}
}
Vz = (~W1[Blen][WordNum]) % modLastWord;
while (Vz != 0) {
Vz = Vz & (Vz - 1);
++z;
}
return z;
}
++layer;
}
else {

carrybitW = 0;
for (int i = 0; i <= WordNum - 1; ++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);
//detect overflow
if (U > MaxNUM64 - tmpa)
carrybitW = 1;
else
carrybitW = 0;

//Main calculations
uint64_t Up = 0, Upp = 0, Wp = 0, Wpp = 0, Ut = 0, Wt = 0, Vt = 0, carrybitWp = 0, carrybitWpp = 0;
int tmpv = ceil(log2(WordSize)) - 1;
for (int k = 1; k <= Blen; ++k) {
carrybitWp = 0, carrybitWpp = 0;
bool f = false;
for (int i = 0; i <= WordNum - 1; ++i)
{
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;

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;

Ut = Wp | Wpp;
Wt = Wp ^ Wpp;
Vt = Wt;

for (int ipp = 0; ipp <= tmpv; ++ipp) {
Vt = Vt ^ (Vt << (1 << ipp));
}
if (f == true) Vt = ~Vt;
if ((Vt & 0X8000000000000000) > 0) f = true;
else f = false;
W2[k][i] = ~(Wt & Vt) & Ut;
}
tmpa = W1[k][WordNum];
Up = tmpa & Y[A[layer] - 1][WordNum];
Wp = ((tmpa + Up + carrybitWp) % modLastWord) | (tmpa - Up);
if (Up > MaxNUM64 - tmpa)
carrybitWp = 1;
else
carrybitWp = 0;

tmpa = W2[k - 1][WordNum];
Upp = tmpa & Y[B[k] - 1][WordNum];
Wpp = ((tmpa + Upp + carrybitWpp) % modLastWord) | (tmpa - Upp);
if (Upp > MaxNUM64 - tmpa)
carrybitWpp = 1;
else
carrybitWpp = 0;

Ut = Wp | Wpp;
Wt = Wp ^ Wpp;
Vt = Wt;

for (int ipp = 0; ipp <= tmpv; ++ipp) {
Vt = Vt ^ (Vt << (1 << ipp));
}
if (f == true) Vt = ~Vt;
if ((Vt & 0X8000000000000000) > 0) f = true;
else f = false;
W2[k][WordNum] = ~(Wt & Vt) & Ut;
}
if (layer == Alen)
{
int z = 0;
uint64_t Vz = 0;
for (int i = 0; i <= WordNum - 1; ++i) {
Vz = ~W2[Blen][i];
while (Vz != 0) {
Vz = Vz & (Vz - 1);
++z;
}
}
Vz = (~W2[Blen][WordNum]) % modLastWord;
while (Vz != 0) {
Vz = Vz & (Vz - 1);
++z;
}
return z;
}
++layer;
}
}
}

```
The files
All_MLCS.cpp

Reference
Deorowicz, S., Danek, A., 2013. Bit-parallel algorithms for the merged longest common sub- sequence problem. International Journal of Foundations of Computer Science 24, 1281–1298.