半在线卷积的实现

rogeryoungh · 2021-12-15 22:41:52 · 个人记录

为了简明，只写二叉。

递归的方式已经被大家熟知，即 CDQ 分治。

struct PolyCDQ {
    int deg, lim;
    Poly ret, F, nf[32];
    std::function<void(int)> relax;

    PolyCDQ(Poly tf) : F(tf) {
        deg = F.size(), lim = get_lim(deg);
        ret.resize(lim);
    }

    void cdq(int l, int r) {
        int len = r - l, mid = (l + r) / 2;
        if (r - l <= 64) {
            for (int i = l; i < std::min(r, deg); ++i) {
                for (int j = l; j < i; ++j)
                    ret[i] += ret[j] * F[i - j];
                relax(i);
            }
            return;
        }
        cdq(l, mid);
        if (mid > deg)
            return;
        Poly a = Poly(ret.begin() + l, len / 2);
        Poly &b = nf[std::__lg(len)];
        if (b.empty()) {
            b = F.cut(len);
        }
        mul(a, b, len);
        for (int i = len / 2; i < len; i++)
            ret[l + i] += a[i];
        cdq(mid, r);
    }

    Poly exp() {
        pre_inv(lim);
        for (int i = 0; i < deg; i++)
            F[i] *= i;
        relax = [&](int i) {
            ret[i] = i == 0 ? 1 : ret[i] * Inv[i];
        };
        cdq(0, lim);
        return ret.cut(deg);
    }

    Poly inv() {
        m32 iv = F[0].inv();
        relax = [&](int i) {
            ret[i] = i == 0 ? iv : -ret[i] * iv;
        };
        cdq(0, lim);
        return ret.cut(deg);
    }

    Poly quo(Poly h) { // 注意是 h / f
        h.resize(lim);
        m32 iv = F[0].inv();
        relax = [&](int i) {
            ret[i] = i == 0 ? h[0] * iv : (h[i] - ret[i]) * iv;
        };
        cdq(0, lim);
        return ret.cut(deg);
    }
};

我想了想，好像迭代树和树状数组挺像的。被 hly1204 的 next 方式启发，于是有了以下的非递归写法。

struct PolyCDQ {
    int deg, lim, now = 1;
    const int M = 64;
    Poly ret, F, nf[32];

    PolyCDQ(Poly tf) : F(tf) {
        deg = F.size(), lim = get_lim(deg);
        ret.resize(lim);
    }

    m32 next() {
        int len = now & -now, l = now - len;
        if (len < M) {
            for (int j = now & -M; j < now; ++j)
                ret[now] += ret[j] * F[now - j];
        } else {
            Poly a = ret.cut(len, l), &b = nf[__lg(len)];
            if (b.empty())
                b = F.cut(len * 2);
            mul(a, b, len * 2);
            for (int i = len; i < len * 2; i++)
                ret[l + i] += a[i];
        }
        return ret[now++];
    }

    Poly inv() {
        m32 iv = F[0].inv();
        ret[0] = iv;
        for (int i = 1; i < deg; i++)
            ret[i] = -next() * iv;
        return ret.cut(deg);
    }

    Poly exp() {
        pre_inv(lim);
        for (int i = 0; i < deg; i++)
            F[i] *= i;
        ret[0] = 1;
        for (int i = 1; i < deg; i++)
            ret[i] = next() * Inv[i];
        return ret.cut(deg);
    }

    Poly quo(Poly h) { // 注意是 h / f
        h.resize(deg);
        m32 iv = F[0].inv();
        ret[0] = h[0] * iv;
        for (int i = 1; i < deg; i++)
            ret[i] = (h[i] - next()) * iv;
        return ret.cut(deg);
    }
};

在线的还在搞。

感觉 sqrt 可以半在线的实现，但是我不会。 弄出来了一个，但是跑的好慢。

struct PolyCDQS {
    int deg, lim, now = 1;
    static const int M = 64;
    Poly ret, &F = ret, nf[32], h;

    PolyCDQS(Poly tf) : h(tf) {
        deg = h.size(), lim = get_lim(deg);
        ret.resize(lim);
    }

    m32 next() {
        int len = now & -now, l = now - len;
        if (len < M) {
            if (now >= M) {
                for (int j = now & -M; j < now; ++j)
                    ret[now] += ret[j] * F[now - j] * 2;
            } else {
                for (int j = now & -M; j < now; ++j)
                    ret[now] += ret[j] * F[now - j];
            }
        } else {
            Poly a = ret.cut(len, l), &b = nf[__lg(len)];
            if (b.empty())
                b = F.cut(len * 2);
            if (l == len * 2)
                b = F.cut(len * 2) * Poly{2};
            mul(a, b, len * 2);
            for (int i = len; i < len * 2; i++) {
                ret[l + i] += a[i];
            }
        }
        return ret[now++];
    }

    Poly sqrt() {
        fill(F.begin(), F.end(), 0);
        F[0] = ret[0] = h[0].sqrt();
        m32 iv = (F[0] * 2).inv();
        for (int i = 1; i < deg; i++) {
            F[i] = ret[i] = (h[i] - next()) * iv;
        }
        return ret.cut(deg);
    }
};

在线的也造出来了。见 #150. 挑战多项式