BerlekampMassey.d

/**
Berlekamp-Massey Algorithm の実装
東京大学理学部数学科 数学講究XB
2012/05/15 松本眞教授「互除法と Berlekamp-Massey 法」
の講義ノートによる
*/
module BerlekampMassey;

alias bool B; /// bool 配列を多用するためエイリアスにしたが、本当はこうするのはよくない。
import std.algorithm;
import std.stdio;
import std.typecons : Tuple, Proxy;
import std.format;
import std.conv;
//alias  QR;

version = demo;

/**
除算の結果: 商 (Polynomial) と余り (Laurent) とを保持する struct
*/
struct BPL
{
    alias Tuple!(BP, "quotient", BL, "remainder") QR;
    private QR value;
    mixin Proxy!value;
    this (B[] quotient, B[] remainder)
    {
        value = QR(BP(quotient), BL(remainder));
    }
    string toString() // 文字列化の方法を変更
    {
        return this.quotient.toString() ~ " ... " ~ this.remainder.toString();
    }
}

/// Binary Polynomial
struct BP
{
    B[] coefs; /// c[i] は x^i の係数
    this (B[] coefs) /// コンストラクタ: 余分な 0 を削除する
    {
        foreach_reverse (i, c; coefs)
        {
            if (c)
            {
                this.coefs = coefs[0..i+1];
                break;
            }
        }
    }

    string toString() // 文字列化
    {
        string ret = "";
        foreach_reverse (b; this.coefs)
        {
            ret ~= b ? "1" : "0";
        }
        if (ret == "")
        {
            return "0";
        }
        return ret;
    }

    BP opUnary(string s)() if (s == "+" || s == "-") // 単項演算子は identity
    {
        return BP(this.coefs);
    }

    BP opBinary(string s)(BP other) if (s == "+" || s == "-") // 2項演算: 加減は同じ
    {
        B[] coefs = this.coefs.dup;
        coefs.length = this.coefs.length.max(other.coefs.length);
        foreach (i; 0..other.coefs.length)
        {
            coefs[i] ^= other.coefs[i];
        }
        return BP(coefs);
    }

    BP opBinary(string s)(BP other) if (s == "*") // 乗算
    {
        if (this.degree < 0 || other.degree < 0)
        {
            return BP([]);
        }
        B[] coefs;
        coefs.length = this.degree + other.degree + 1;
        foreach (i, c; this.coefs)
        {
            if (!c)
            {
                continue;
            }
            foreach (j, d; other.coefs)
            {
                coefs[i + j] ^= d;
            }
        }
        return BP(coefs);
    }

    BP opBinary(string s)(int shift) if (s == ">>") // シフト演算
    {
        if (shift < 0)
        {
            return this << -shift;
        }
        return BP(this.coefs[shift..$]);
    }

    BP opBinary(string s)(int shift) if (s == "<<") // シフト演算
    {
        if (shift < 0)
        {
            return this >> -shift;
        }
        B[] zeros;
        zeros.length = shift;
        return BP(zeros + this.coefs);
    }

    @property int degree() // 次数
    {
        int ret = -1;
        foreach (c, b; this.coefs)
        {
            if (b)
            {
                ret = c;
            }
        }
        return ret;
    }
}

version (demo) unittest
{
    writeln("Polynomial test:");
    writefln("  x^4 + x^3 + x^2 + x^1 + x^0 = %s", BP([true, true, true, true, true]));
    writefln("  x^4 + x^3 + x^2 + x^1 = %s", BP([false, true, true, true, true]));
    writefln("  x^2 + x^1 + x^0 = %s", BP([true, true, true]));
    writefln("  0x^2 + x^1 + x^0 = %s", BP([true, true, false]));
    auto x = BP([true, true]);
    auto y = BP([true, false, true]);
    auto z = BP([false, true, false, false, true]);
    writefln("  %s + %s = %s", x, y, (x + y));
    writefln("  %s + %s = %s", x, z, (x + z));
    writefln("  %s + %s = %s", y, z, (y + z));
    writefln("  %s * %s = %s", x, y, (x * y));
    writefln("  %s * %s = %s", x, z, (x * z));
    writefln("  %s * %s = %s", y, z, (y * z));
    writeln();
}


/// Binary Laurent Polynomial: 非正次数に限る
struct BL
{
    B[] coefs; /// c[i] は -i 次の係数
    this (B[] coefs, long n = -1) // 次数を指定すると足りなければ 0 で埋める。
    {
        this.coefs = coefs.dup;
        if (coefs.length < n)
        {
            this.coefs.length = cast(size_t)n;
        }
    }

    @property int degree() // 次数
    {
        foreach (i, c; this.coefs)
        {
            if (c)
            {
                return -i;
            }
        }
        return 1;
    }

    bool opCast(T)() if (is(T == bool)) /// non-zero check
    {
        return reduce!("a || b")(false, coefs);
    }

    string toString()
    {
        string ret = "";
        bool first = true;
        int leading_zeros;
        foreach (b; coefs)
        {
            if (first && !b)
            {
                leading_zeros += 1;
                continue;
            }
            ret ~= (b ? "1" : "0") ~ (first ? "." : "");
            first = false;
        }
        return first ? "0" : // 1 が最初であるというフラグが寝ている: 0 を示している
            ret ~ (leading_zeros ? "e-" ~ leading_zeros.to!string() : "");
    }

    BL opUnary(string s)() if (s == "+" || s == "-")
    {
        return BL(this.coefs);
    }

    BL opBinary(string s)(BL other) if (s == "+" || s == "-") // やることは一緒なので BP と BP の加算を流用
    {
        auto n = min(this.coefs.length, other.coefs.length); // 有効桁数は少ない方に合わせる
        return BL((BP(this.coefs) + BP(other.coefs)).coefs[0..($.max(n))], n); // 桁が縮んでいる可能性もあるので n に合わせる
    }

    BL opBinary(string s)(BP multiplier) if (s == "*")
    {
        B[] product;
        product.length = this.coefs.length;
        foreach (i, b; multiplier.coefs)
        {
            if (!b)
            {
                continue;
            }
            foreach_reverse (j, c; this.coefs)
            {
                try
                {
                    product[j - i] ^= c;
                }
                catch
                {
                    break;
                }
            }
        }
        return BL(product);
    }

    BL opBinary(string s)(int shift) if (s == ">>")
    {
        if (shift < 0)
        {
            return this << -shift;
        }
        B[] zeros;
        zeros.length = shift;
        return BL(zeros ~ this.coefs);
    }

    BL opBinary(string s)(int shift) if (s == "<<")
    {
        if (shift < 0)
        {
            return this >> -shift;
        }
        return BL(this.coefs[shift..$]);
    }

    BL opBinary(string s)(BL other) if (s == "*") // BP*BP を流用
    {
        return BL((BP(this.coefs) * BP(other.coefs)).coefs);
    }

    BPL opBinary(string s)(BL divisor) if (s == "/") // 除算
    {
        assert (this.degree >= divisor.degree); // 次数に関する制約: 0 による除算もdivisor.degree = 1 で弾くことができる。
        B[] remainder = this.coefs.dup;
        B[] quotient;
        sizediff_t remainder_deg, quotient_deg;

        while (true)
        {
            if (
                0 < (remainder_deg = -remainder.countUntil(true)) || // 割り切れた
                0 > (quotient_deg = remainder_deg - divisor.degree) // 余りが出た
            )
            {
                return BPL(quotient, remainder);
            }
            quotient.length = quotient.length.max(quotient_deg + 1);
            quotient[quotient_deg] = true;
            foreach_reverse (i, c; divisor.coefs)
            {
                if (i < quotient_deg)
                {
                    break;
                }
                remainder[i - quotient_deg] ^= divisor.coefs[i];
            }
        }
    }
}

version (demo) unittest
{
    writeln("Laurent test:");
    writefln("  x^-4 + x^-3 + x^-2 + x^-1 + x^-0 = %s", BL([true, true, true, true, true]));
    writefln("  x^-4 + x^-3 + x^-2 + x^-1 = %s", BL([false, true, true, true, true]));
    writefln("  x^-2 + x^-1 + x^-0 = %s", BL([true, true, true]));
    writefln("  0x^-2 + x^-1 + x^-0 = %s", BL([true, true, false]));
    auto x = BL([true, true]);
    auto y = BL([true, false, true]);
    auto z = BL([false, true, false, false, true]);
    writefln("  %s + %s = %s", x, y, (x + y));
    writefln("  %s + %s = %s", x, z, (x + z));
    writefln("  %s + %s = %s", y, z, (y + z));
    writefln("  %s * %s = %s", x, y, (x * y));
    writefln("  %s * %s = %s", x, z, (x * z));
    writefln("  %s * %s = %s", y, z, (y * z));
    writeln();
}

version (demo) unittest
{
    writeln("Multiplication test:");
    auto xp = BP([true, true]);
    auto yp = BP([true, false, true]);
    auto zp = BP([false, true, false, false, true]);
    auto xl = BL([true, true]);
    auto yl = BL([true, false, true]);
    auto zl = BL([false, true, false, false, true]);
    foreach (l; [xl, yl, zl])
    {
        foreach (p; [xp, yp, zp])
        {
            writefln("  %s * %s = %s", l, p, (l*p));
        }
    }
    writeln();
}
 
version (demo) unittest
{
    writeln("Division test:");
    void test_division(BL dividend, BL divisor)
    {
        writefln("  %s / %s = %s", dividend, divisor, (dividend / divisor));
    }

    auto dividend = BL([true], 6);
    foreach (divisor; [[true, true], [false, true, true], [false, false, true, false, true]])
    {
        test_division(dividend, BL(divisor, 6));
    }
    writeln();
}


/**
BP を要素とする 2*2 行列
*/
struct matrix
{
    BP[2][2] elem;

    this (BP q) /// 左上に q を、左下と右上に 1 を、右下に 0 を、それぞれ入れる
    {
        this.elem[0][0] = -q;
        this.elem[0][1] = BP([true]);
        this.elem[1][0] = BP([true]);
    }

    this (
        BP a, BP b,
        BP c, BP d
    ) /// 4成分を指定して行列を作る
    {
        this.elem[0][0] = a; this.elem[0][1] = b;
        this.elem[1][0] = c; this.elem[1][1] = d;
    }

    string toString()
    {
        return "matrix(\n" ~ 
            this.elem[0][0].toString() ~ " " ~ this.elem[0][1].toString() ~  "\n" ~ 
            this.elem[1][0].toString() ~ " " ~ this.elem[1][1].toString() ~ "\n)";
    }

    @property BP[2] lowerRow()
    {
        return this.elem[1];
    }

    @property BP[2] upperRow()
    {
        return this.elem[0];
    }

    matrix opBinary(string s)(matrix other) if (s == "*") // 乗算
    {
        return matrix(
            this.elem[0][0] * other.elem[0][0] + this.elem[0][1] * other.elem[1][0],
            this.elem[0][0] * other.elem[0][1] + this.elem[0][1] * other.elem[1][1],
            this.elem[1][0] * other.elem[0][0] + this.elem[1][1] * other.elem[1][0],
            this.elem[1][0] * other.elem[0][1] + this.elem[1][1] * other.elem[1][1]
        );
    }

    static matrix unit() // 単位行列
    {
        return matrix(BP([true]), BP([]), BP([]), BP([true]));
    }
}


/**
Berlekamp-Massey Algorithm
Params:
    arg = 最初の2N項    
Returns:
    [分母, 分子]: arg = 分子/分母
*/
BP[2] Berlekamp_Massey(B[] arg)
{
    auto N = arg.length >> 1; // 最大次数
    BL x = BL([false] ~ arg);
    BL y = BL([true], arg.length+1);
    /*
    講義ノートでは (1, y) からスタートしており、最終的に y = a[n] / b[n] を得ている。
    ここでは y を x で割ることを繰り返したいので、(x, 1) で始めることにする。
    */
    matrix qm = matrix.unit();
    while (true)
    {
        auto qr = y / x;
        debug writefln("%s / %s = %s", y, x, qr);
        qm = matrix(qr.quotient) * qm; // 行列を更新
        debug writeln(qm);
        assert (qm.lowerRow[0].degree < qm.upperRow[0].degree);
        assert (qm.lowerRow[1].degree < qm.upperRow[1].degree);
        /*
        次数が増加していることをチェック: レポート問題4 (1) に対応。
        もちろん Assersion failure が出ないことで証明になるわけではないが、
        プログラムの誤りで無限ループになってしまったとき検出の役に立つかもしれない。
    証明:
        まず、すべてのステップで商の次数は正であることを示す。
        最初のステップでは次数の大きい 1 を次数の小さい 0.***... で割るため、商の次数は正である。
        それ以降のステップでは前のステップの divisor を remainder で割ることになる。
        前者の方が大きい次数をもつため、商の次数は正である。
        
        次に、次数が単調増加になっていることを示す。
        n (非負整数) 回の更新を終えた時点の qm の値を A[n] とする。
        ここでは行列の要素を 0-origin で A[n](0, 0) のように参照する。
            a[n] := A[n](0, 0)
            b[n] := A[n](0, 1)
        とすると、matrix(q) をかけることによって
            A[n+1](1, 0) = a[n]
            A[n+1](1, 1) = b[n]
        となる。これから数学的帰納法を用いる。
        initial step:
            A[1] を見ると a[0] = 1, a[1] = q, b[0] = 0, b[1] = 1 である。
            したがって、a[0].deg < a[1].deg および b[0].deg < b[1].deg が成り立つ。
        inductive step:
            n >0 とし、a[n-1].deg < a[n].deg, b[n-1].deg < b[n].deg を仮定する。
            a[n] および b[n] は同じ形の漸化式 p[n+1] = q[n]p[n] + p[n-1] をみたす。
            ここで q[n] たちは次数が正であるから、右辺第1項の次数は p[n] の次数より真に大きく、
            かつ右辺第2項の次数よりも大きいから、それは左辺の次数となる。
            したがって a[n].deg < a[n+1].deg, b[n].deg < b[n+1].deg を得る。
        */
        if (N < qm.upperRow[0].degree) // a[n] の次数が N を超えた
        {
            debug writefln("N < b[n+1].degree");
            return qm.lowerRow; // 欲しい (a, b)[n] は下段
        }
        if (!qr.remainder) // 余り 0: 割り切れたときは上段が欲しい。
        {
            debug writefln("x | y");
            return qm.upperRow;
        }
        y = x;
        x = qr.remainder;
    }
}

/// 入力パース: 最初に現れるスペースまでの '1' を true へ、それ以外の文字 (普通は '0') を false へ変換
auto toBs(string s)
{
    B[] ret;
    foreach (c; s)
    {
        if (c == ' ')
        {
            break;
        }
        ret ~= c == '1';
    }
    return ret;
}


auto BerlekampMassey(string s)
{
    return s.toBs().Berlekamp_Massey();
}


version (demo) unittest
{
    writeln("BMA test");
    foreach (seq; [
        "1",
        "11111111", // 1; 11
        "10101010", // 10; 101
        "11011011", // 11; 111
        "1111010110010001111010110010001111010110010001111010110010001111", // 1111; 11001
        "1101011001000111101011001000111101011001000111101011001000111101", // 1101; 11001
        "0010001111010110010001111010110010001111010110010001111010110010", // 0010; 11001
        "1111111101011001000111101011001000111101011001000111101011001000", // 11111111; 11001
        "0000000000001101101101101101101101101101101101101101101101101101"
    ])
    {
        auto res = seq.BerlekampMassey();
        "input: %s".writefln(seq);
        "output: %s".writefln(res);
        writeln();
    }
}

main.d

module main;

import std.stdio;
import std.string;
import std.algorithm;
import BerlekampMassey;

alias bool B;

void main()
{
    while (true)
    {
        auto buf = readln().strip();
        if (buf.length == 0)
        {
            break;
        }
        writefln("%s", buf.toBs().Berlekamp_Massey());
    }
}