区間代入/区間積 O(logN)/query で償却計算量にしないアルゴリズム

changelog

2024-02-23 問題文中の積の説明を改善。実行時間測定用の例題の具体的な問題を説明。

概要
1. 実用性
問題
方法
実装
速度
（追記）その他の方法

概要

こういうものがあります。

[1] :

その記事によると、計算量を償却することで空間 $O(N)$ を達成しています。この記事では計算量を償却しない場合の方法を考えます。

実用性

そんなものはない

問題

考える問題は次の通りです。

長さ $N$ の列 $A=(a_0,a_1,\ldots ,a_{N-1})$ がある。列の要素には結合的な積が定義されている。以降で必要な積の計算過程のすべてにおいて、 $2$ 項の積 $1$ 回あたり $O(1)$ 時間で計算できる。この列に対して以下の操作をたくさん処理せよ。

$l,r,x$ が与えられる。 $a_l,a_{l+1},\ldots ,a_{r-1}$ の値をそれぞれ $x$ に変更する。
$l,r$ が与えられる。積 $a_la_{l+1}\cdots a_{r-1}$ を計算する。

[1] では、クエリが与えられるごとにすぐに（つまり、オンラインで）処理するとして、 $O(N)$ の空間で各クエリを償却 $O(\log N)$ 時間で処理できることが示されています。この記事では、空間 $O(N)$ は守りつつ各クエリを償却なしの $O(\log N)$ 時間で処理する方法を紹介します。

方法

以降では、 $N$ は $2$ 冪であるとします。 $n=\log_2 N$ とします。

$x$ で値を上書きするクエリが与えられるごとに、 $x$ を作用素 $k$ として $k=0,1,2,\ldots $ の順に番号を振ります。この番号は $k=2N$ を $k=0$ と同一視して周期 $2N$ で変化させます。また、作用素の $2^e$ 乗（ $0\leq e\leq n$ ）を管理するために、長さ $2N$ の配列 $B$ を使います。また、整数 $i$ で $B[i]$ とアクセスしたときは、代わりに $B$ の $i\bmod {2N}$ 番目の要素にアクセスすることにします。

作用素 $k$ の $2^e$ 乗は、 $B[k+2N-2^{n-e}]$ に上書きで保存することにします。

必要な作用素が上書きされて失われるのを防ぐために、作用素 $k$ の更新処理が終わった直後に、 $k\bmod N$ に $n$ ビットのビットリバーサルをした値を $k^{\prime}$ として、 $a_{k^{\prime}}$ に作用する遅延伝播をすべて解消します。

$2^{e}$ 乗を保存したものが上書きで失われるまでに $2^{n-e}$ 回の猶予があります。ビットリバーサルの順で遅延伝播を解消することで、 $2^{e}$ 乗が必要なノードの遅延伝播がすべて解消されます。

実装

実装例 (C++)

Elem どうしの掛け算が定義されている必要があります。

#include <vector>
#include <algorithm>

template<class Elem>
struct RangeUpdateRangeProduct {
private:
    int N;
    int logN;
    std::vector<int> lazy;
    std::vector<int> nodeRank;
    std::vector<Elem> lazyTable;
    std::vector<Elem> prod;
    int tableIterator;
    int bitReverse(int x){
        int y = 0;
        for(int i=0; i<logN; i++){
            y = (y << 1) | (x & 1);
            x >>= 1;
        }
        return y;
    }
    void apply(int p, int k){
        if(p < N) lazy[p] = k;
        prod[p] = lazyTable[(k+N*2-nodeRank[p]) & (N*2-1)];
    }
    void propagate(int p){
        if(p >= N || lazy[p] == -1) return;
        apply(p*2, lazy[p]);
        apply(p*2+1, lazy[p]);
        lazy[p] = -1;
    }
    void propagateFor(int p){
        p += N;
        for(int d=logN; d>=1; d--) propagate(p>>d);
    }
    void rangeUpdateRecursion(int l, int r, int k, int i, int a, int b){
        if(b <= l || r <= a) return;
        if(l <= a && b <= r){ apply(i, k); return; }
        propagate(i);
        rangeUpdateRecursion(l, r, k, i*2, a, (a+b)/2);
        rangeUpdateRecursion(l, r, k, i*2+1, (a+b)/2, b);
        prod[i] = prod[i*2] * prod[i*2+1];
    }
    Elem rangeProductRecursion(int l, int r, int i, int a, int b){
        if(l <= a && b <= r) return prod[i];
        propagate(i);
        if(r <= (a+b)/2) return rangeProductRecursion(l, r, i*2, a, (a+b)/2);
        if((a+b)/2 <= l) return rangeProductRecursion(l, r, i*2+1, (a+b)/2, b);
        Elem lq = rangeProductRecursion(l, r, i*2, a, (a+b)/2);
        Elem rq = rangeProductRecursion(l, r, i*2+1, (a+b)/2, b);
        return lq * rq;
    }
public:
    // すべての要素を filler で初期化
    RangeUpdateRangeProduct(int n, Elem filler){
        N = 1; logN = 0; while(N < n){ N *= 2; logN++; }
        tableIterator = 0;
        lazy.assign(N, -1);
        nodeRank.assign(N*2, N);
        for(int i=N-1; i>=1; i--) nodeRank[i] = nodeRank[i*2] >> 1;
        lazyTable.assign(N*2, filler);
        prod.assign(N*2, filler);
        for(int i=N-1; i>=1; i--) prod[i] = prod[i*2] * prod[i*2+1];
    }
    void rangeUpdate(int l, int r, Elem x){
        int k = tableIterator++;
        if(tableIterator == N*2) tableIterator = 0;
        lazyTable[(k+N*2-N) & (N*2-1)] = x;
        for(int i=1; i<=logN; i++){
            x = x * x;
            lazyTable[(k+N*2-(N>>i)) & (N*2-1)] = x;
        }
        rangeUpdateRecursion(l, r, k, 1, 0, N);
        propagateFor(bitReverse(k&(N-1)));
    }
    // l < r
    Elem rangeProd(int l, int r){
        return rangeProductRecursion(l, r, 1, 0, N);
    }
};

速度

実行時間測定用の例題として、要素が $998\, 244\, 353$ 未満の非負整数で、「積」が通常の積を $998\, 244\, 353$ で割ったあまりである場合を使います。この問題について、 [1] からリンクされていた実装と、上で示した実装それぞれ AtCoder のコードテストで実行時間を計測しました。要素数 $N$ とクエリの個数 $Q$ の大きさを変えながらそれぞれ $3$ つのシード値を入力して計測して得た値の中央値を表にまとめました。

実験に用いたソースコード

クリックで展開

#include <vector>
#include <algorithm>

namespace nachia {

template<class Elem>
struct assign_segment_tree {
private:
    int N;
    int logN;
    std::vector<int> lazy;
    std::vector<int> nodeRank;
    std::vector<Elem> lazyTable;
    std::vector<Elem> prod;
    int tableIterator;
    int bitReverse(int x){
        int y = 0;
        for(int i=0; i<logN; i++){
            y = (y << 1) | (x & 1);
            x >>= 1;
        }
        return y;
    }
    void apply(int p, int k){
        if(p < N) lazy[p] = k;
        prod[p] = lazyTable[(k+N*2-nodeRank[p]) & (N*2-1)];
    }
    void propagate(int p){
        if(p >= N || lazy[p] == -1) return;
        apply(p*2, lazy[p]);
        apply(p*2+1, lazy[p]);
        lazy[p] = -1;
    }
    void propagateFor(int p){
        p += N;
        for(int d=logN; d>=1; d--) propagate(p>>d);
    }
    void rangeUpdateRecursion(int l, int r, int k, int i, int a, int b){
        if(b <= l || r <= a) return;
        if(l <= a && b <= r){ apply(i, k); return; }
        propagate(i);
        rangeUpdateRecursion(l, r, k, i*2, a, (a+b)/2);
        rangeUpdateRecursion(l, r, k, i*2+1, (a+b)/2, b);
        prod[i] = prod[i*2] * prod[i*2+1];
    }
    Elem rangeProductRecursion(int l, int r, int i, int a, int b){
        if(l <= a && b <= r) return prod[i];
        propagate(i);
        if(r <= (a+b)/2) return rangeProductRecursion(l, r, i*2, a, (a+b)/2);
        if((a+b)/2 <= l) return rangeProductRecursion(l, r, i*2+1, (a+b)/2, b);
        Elem lq = rangeProductRecursion(l, r, i*2, a, (a+b)/2);
        Elem rq = rangeProductRecursion(l, r, i*2+1, (a+b)/2, b);
        return lq * rq;
    }
public:
    assign_segment_tree(int n, Elem filler){
        N = 1; logN = 0; while(N < n){ N *= 2; logN++; }
        tableIterator = 0;
        lazy.assign(N, -1);
        nodeRank.assign(N*2, N);
        for(int i=N-1; i>=1; i--) nodeRank[i] = nodeRank[i*2] >> 1;
        lazyTable.assign(N*2, filler);
        prod.assign(N*2, filler);
        for(int i=N-1; i>=1; i--) prod[i] = prod[i*2] * prod[i*2+1];
    }
    void assign(int l, int r, Elem x){
        int k = tableIterator++;
        if(tableIterator == N*2) tableIterator = 0;
        lazyTable[(k+N*2-N) & (N*2-1)] = x;
        for(int i=1; i<=logN; i++){
            x = x * x;
            lazyTable[(k+N*2-(N>>i)) & (N*2-1)] = x;
        }
        rangeUpdateRecursion(l, r, k, 1, 0, N);
        propagateFor(bitReverse(k&(N-1)));
    }
    // l < r
    Elem fold(int l, int r){
        return rangeProductRecursion(l, r, 1, 0, N);
    }
};

}


#include <atcoder/modint>
using Modint = atcoder::static_modint<998244353>;

namespace noshi91 {
    
using T = Modint;
T op(T l, T r) { return l * r; }
T id() { return 1; }

class assign_segment_tree {
private:
  class node {
  public:
    T sum;
    T *lazy;
    node() : sum(id()), lazy(nullptr) {}
    T get() {
      if (lazy) {
        return *lazy;
      } else {
        return sum;
      }
    }
  };

  int height;
  std::vector<node> tree;
  std::vector<T> table;

  void push(int index) {
    if (tree[index].lazy) {
      tree[index * 2].lazy = tree[index].lazy - 1;
      tree[index * 2 + 1].lazy = tree[index].lazy - 1;
      tree[index].sum = *tree[index].lazy;
      tree[index].lazy = nullptr;
    }
  }

  T fold(int index, int n_left, int n_right, int q_left, int q_right) {
    if (q_left <= n_left && n_right <= q_right) {
      return tree[index].get();
    }
    if (n_right <= q_left || q_right <= n_left) {
      return id();
    }
    push(index);
    int n_mid = (n_left + n_right) / 2;
    return op(fold(index * 2, n_left, n_mid, q_left, q_right),
              fold(index * 2 + 1, n_mid, n_right, q_left, q_right));
  }

  void assign(int index, int n_left, int n_right, int q_left, int q_right,
              T *lazy) {
    if (q_left <= n_left && n_right <= q_right) {
      tree[index].lazy = lazy;
      return;
    }
    if (n_right <= q_left || q_right <= n_left) {
      return;
    }
    push(index);
    int n_mid = (n_left + n_right) / 2;
    assign(index * 2, n_left, n_mid, q_left, q_right, lazy - 1);
    assign(index * 2 + 1, n_mid, n_right, q_left, q_right, lazy - 1);
    tree[index].sum = op(tree[index * 2].get(), tree[index * 2 + 1].get());
  }

public:
  assign_segment_tree(int n) {
    height = 1;
    int s = 1;
    while (s < n) {
      s *= 2;
      ++height;
    }
    tree.assign(s * 2, node());
    table.reserve(s * 2);
  }

  int size() { return tree.size() / 2; }

  T fold(int first, int last) { return fold(1, 0, size(), first, last); }

  void assign(int first, int last, T value) {
    for (int i = 0; i < height; ++i) {
      table.push_back(value);
      value = op(value, value);
    }
    assign(1, 0, size(), first, last, &table.back());
    if (table.capacity() - table.size() < height) {
      for (int i = 1; i < size(); ++i) {
        push(i);
      }
      for (int i = size(); i != tree.size(); ++i) {
        if (tree[i].lazy) {
          tree[i].sum = *tree[i].lazy;
          tree[i].lazy = nullptr;
        }
      }
      table.clear();
    }
  }
};

}

struct Problem {
    int N;
    int Q;
    Modint init;
    std::vector<int> T;
    std::vector<int> L;
    std::vector<int> R;
    std::vector<Modint> val;
};

#include <unordered_map>
#include <cassert>
#include <cstdint>


namespace nachia{

class Xoshiro256pp{
public:

    using i32 = int32_t;
    using u32 = uint32_t;
    using i64 = int64_t;
    using u64 = uint64_t;


private:
    uint64_t s[4];

    // https://prng.di.unimi.it/xoshiro256plusplus.c
    static inline uint64_t rotl(const uint64_t x, int k) noexcept {
        return (x << k) | (x >> (64 - k));
    }
    inline uint64_t gen(void) noexcept {
        const uint64_t result = rotl(s[0] + s[3], 23) + s[0];
        const uint64_t t = s[1] << 17;
        s[2] ^= s[0];
        s[3] ^= s[1];
        s[1] ^= s[2];
        s[0] ^= s[3];
        s[2] ^= t;
        s[3] = rotl(s[3], 45);
        return result;
    }

    // https://xoshiro.di.unimi.it/splitmix64.c
    u64 splitmix64(u64& x) {
        u64 z = (x += 0x9e3779b97f4a7c15);
        z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
        z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
        return z ^ (z >> 31);
    }

public:

    void seed(u64 x = 7001){
        assert(x != 0);
        s[0] = x;
        for(int i=1; i<4; i++) s[i] = splitmix64(x);
    }

    Xoshiro256pp(){ seed(); }
    
    u64 rng64() { return gen(); }
    u64 operator()(){ return gen(); }

    // generate x : l <= x <= r
    u64 random_unsigned(u64 l,u64 r){
        assert(l<=r);
        r-=l;
        auto res = rng64();
        if(res<=r) return res+l;
        u64 d = r+1;
        u64 max_valid = 0xffffffffffffffff/d*d;
        while(true){
            auto res = rng64();
            if(res<=max_valid) break;
        }
        return res%d+l;
    }

    // generate x : l <= x <= r
    i64 random_signed(i64 l,i64 r){
        assert(l<=r);
        u64 unsigned_l = (u64)l ^ (1ull<<63);
        u64 unsigned_r = (u64)r ^ (1ull<<63);
        u64 unsigned_res = random_unsigned(unsigned_l,unsigned_r) ^ (1ull<<63);
        return (i64)unsigned_res;
    }


    // permute x : n_left <= x <= n_right
    // output r from the front
    template<class Int>
    std::vector<Int> random_nPr(Int n_left, Int n_right, Int r){
        Int n = n_right-n_left;

        assert(n>=0);
        assert(r<=(1ll<<27));
        if(r==0) return {};  
        assert(n>=r-1);

        std::vector<Int> V;
        std::unordered_map<Int,Int> G;
        for(int i=0; i<r; i++){
            Int p = random_signed(i,n);
            Int x = p - G[p];
            V.push_back(x);
            G[p] = p - (i - G[i]);
        }

        for(Int& v : V) v+=n_left;
        return V;
    }



    // V[i] := V[perm[i]]
    // using swap
    template<class E,class PermInt_t>
    void permute_inplace(std::vector<E>& V,std::vector<PermInt_t> perm){
        assert(V.size() == perm.size());
        int N=V.size();
        for(int i=0; i<N; i++){
            int p=i;
            while(perm[p]!=i){
                assert(0 <= perm[p] && perm[p] < N);
                assert(perm[p] != perm[perm[p]]);
                std::swap(V[p],V[perm[p]]);
                int pbuf = perm[p]; perm[p] = p; p = pbuf;
            }
            perm[p] = p;
        }
    }

    template<class E>
    std::vector<E> shuffle(const std::vector<E>& V){
        int N=V.size();
        auto P = random_nPr(0,N-1,N);
        std::vector<E> res;
        res.reserve(N);
        for(int i=0; i<N; i++) res.push_back(V[P[i]]);
        return res;
    }

    // shuffle using swap
    template<class E>
    void shuffle_inplace(std::vector<E>& V){
        int N=V.size();
        permute_inplace(V,random_nPr(0,N-1,N));
    }


};

} // namespace nachia

nachia::Xoshiro256pp rng;

Problem generate(int n, int q){
    Problem res;
    res.N = n;
    res.Q = q;
    res.init = 1;
    for(int i=0; i<q; i++){
        int l = 0, r = 0;
        while(l == r){
            l = rng.random_signed(0, n);
            r = rng.random_signed(0, n);
            if(l > r) std::swap(l, r);
        }
        res.T.push_back(rng.random_signed(0, 1));
        res.L.push_back(l);
        res.R.push_back(r);
        res.val.push_back(rng.random_signed(1, 992844352));
    }
    return res;
}

std::vector<int> solve_nachia(Problem problem){
    std::vector<int> ans;
    nachia::assign_segment_tree<Modint> A(problem.N, problem.init);
    for(int qi=0; qi<problem.Q; qi++){
        int t = problem.T[qi];
        int l = problem.L[qi];
        int r = problem.R[qi];
        Modint e = problem.val[qi];
        if(t == 0){
            A.assign(l, r, e);
        } else {
            ans.push_back((int)A.fold(l, r).val());
        }
    }
    return ans;
}

std::vector<int> solve_noshi91(Problem problem){
    std::vector<int> ans;
    noshi91::assign_segment_tree A(problem.N);
    for(int qi=0; qi<problem.Q; qi++){
        int t = problem.T[qi];
        int l = problem.L[qi];
        int r = problem.R[qi];
        Modint e = problem.val[qi];
        if(t == 0){
            A.assign(l, r, e);
        } else {
            ans.push_back((int)A.fold(l, r).val());
        }
    }
    return ans;
}

#include <iostream>
#include <chrono>

int main(){
    std::ios::sync_with_stdio(false); std::cin.tie(nullptr);
    unsigned long long seed; std::cin >> seed;
    int nq; std::cin >> nq;
    // case 1 : seed = 144262995209904656
    // case 2 : seed = 7577536898627496
    // case 3 : seed = 795742045562236292
    rng.seed(seed);
    auto problem = generate(nq, nq);
    auto t0 = std::chrono::high_resolution_clock::now();
    auto ans = solve_nachia(std::move(problem));
    // auto ans = solve_noshi91(std::move(problem));
    auto t1 = std::chrono::high_resolution_clock::now();
    std::cout << "T = " << (t1-t0).count() << std::endl;
    for(auto a : ans) std::cout << a << '\n';
    return 0;
}

結果を次の表に示します。

$N,Q$	nachia [ms]	noshi91 [ms]
$2\times 10^5$	$105$	$85$
$5\times 10^5$	$298$	$242$
$10\times 10^5$	$723$	$746$
$50\times 10^5$	$6081$	$4580$

多くの場合で noshi91 によるもの [1] のほうが高速でした。表には出ていませんが $\pm 15\%$ くらいのぶれはあったのでたくさん実験すれば $10\times 10^5$ の場合も noshi91 のもののほうが良い成績になると思います。

一旦おわり

（追記）その他の方法

↓ noshi91 さんに情報を頂きました。

kopricky さんに指摘されたものと記憶しているのですが、単方向連結リストでメモを管理して、どこからも参照されなくなったメモを適宜破棄していくと、log(高さ) の和になって空間が O(n) になるという方法もあります。
— 熨斗袋 (@noshi91) December 28, 2023

本編で紹介した方法は、古い作用素による遅延伝播を強制的に押し下げるアプローチによります。しかしこれを行わなくても保存される作用素を柔軟に切り替えれば必要なメモは全体で $O(N)$ のサイズに収まることが示されます。作用素の前計算を単方向連結リストで実装して参照カウントで管理し、不要になったものは直ちにメモリを解放してゆけば、空間 $O(N)$ が保たれます。

実はこのアプローチも検討には入れていましたが、この方法で実験するためには（ C++ の new , delete が遅いため）メモリ管理を実装する必要があり、また参照カウンタも必要になります。このため実行時高速な実装ができないと思い、実験を省きました。

一方で、本編の方法と比べてこの方法には、特殊なアドレスの計算・ビットリバーサルの計算・追加の遅延伝播が要らないという利点があります。

追記おわり