Smiley face

Deorowicz and Danek Algorithm_BMLCS(Bit-Parallel, 2013)


Main features
Abstract of the paper [Deorowicz 2014]

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.

C++ source code
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;
        }

    }
}

The files
  All_BMLCS.cpp

Reference
Danek, A., Deorowicz, S., 2014. Bit-parallel algorithm for the block variant of the merged longest common subsequence problem. Advances in Intelligent Systems and Computing 242, 173–181.