Elias-Fano Encodingで遊ぶ

はじめに

読んでる論文に出てきてたElias-Fano Encodingをちょっと書いて遊んでみた。

Elias-Fano Encodingとは

  • 単調増加整数列の表現方法のひとつ
  • 厳密にはsuccinctではないが、succinctに近い表現
    • quasi-succinct representationと言っている
  • 整数列の各値を「上位bit列」と「下位bit列」に分割し、整数列全体では、上位bit列を「(negated) unary code表現」、下位bit列を「各下位bitを連結したbit列」で表現したもの
    • 上位、下位のbit数は等分ではなく、上位bitがceil(lg(n))個
    • 例: {3,4,5}という整数列の場合
    • ceil(lg(3))=2として、上位2bitと下位1bitに分ける
    • 3は011、4は100、5は101なので、上位bitは01と10と10、下位bitは1と0と1
    • 上位bit列はnegated unary codeで表現するので、00が0回、01が1回、10が2回、11が0回なので、「0101100」と表現する
      • 1^{00の個数}0 1^{01の個数}0 1^{10の個数}0 ...と表現
    • 下位bit列は、連結するので「101」と表現

コード

http://pages.di.unipi.it/pibiri/slides/seminar_ef.pdf
を参考に「正の整数の単調増加列(1以上の整数、同じ整数はなし、ソート済み)」を想定した場合のElias-Fano Encodingするコードを書いてみた。
操作は「access: S[i]の値」と「successor(x): min{S[i] | S[i] >= x}」の2つ。
いくつかのランダムケースで試して問題なさそうなことまでしか確認していない。

#include <vector>
#include <map>
#include <cstdint>
#include <algorithm>
#include <iostream>
#include <cmath>

class FID {
  static const int BIT_SIZE = 64;
  using BLOCK_TYPE = uint64_t;
  int size;
  int block_size;
  std::vector<BLOCK_TYPE> blocks;
  std::vector<int> s_rank, s0_rank;
public:
  BLOCK_TYPE popcount(BLOCK_TYPE x){
    x = ((x & 0xaaaaaaaaaaaaaaaaULL) >> 1) + (x & 0x5555555555555555ULL);
    x = ((x & 0xccccccccccccccccULL) >> 2) + (x & 0x3333333333333333ULL);
    x = ((x & 0xf0f0f0f0f0f0f0f0ULL) >> 4) + (x & 0x0f0f0f0f0f0f0f0fULL);
    x = ((x & 0xff00ff00ff00ff00ULL) >> 8) + (x & 0x00ff00ff00ff00ffULL);
    x = ((x & 0xffff0000ffff0000ULL) >> 16) + (x & 0x0000ffff0000ffffULL);
    x = ((x & 0xffffffff00000000ULL) >> 32) + (x & 0x00000000ffffffffULL);
    return x;
  }
public:
  FID(int size):
  size(size),
  block_size(((size + BIT_SIZE - 1) / BIT_SIZE) + 1),
  blocks(block_size, 0),
  s_rank(block_size, 0),
  s0_rank(block_size, 0)
  {}

  void init(int sz){
    blocks.clear();
    s_rank.clear();
    size = sz;
    block_size = ((size + BIT_SIZE - 1) / BIT_SIZE) + 1;
    blocks.resize(block_size, 0);
    s_rank.resize(block_size, 0);
    s0_rank.resize(block_size, 0);
  }
  
  void set(int i){
    blocks[i/BIT_SIZE] |= 1ULL << (i%BIT_SIZE);
  }
  
  void finalize(){
    s_rank[0] = 0;
    for(int i=1; i<block_size; i++){
      s_rank[i] = s_rank[i-1] + popcount(blocks[i-1]);
    }
    s0_rank[0] = 0;
    for(int i=1; i<block_size; i++){
      s0_rank[i] = s0_rank[i-1] + (BIT_SIZE - popcount(blocks[i-1]));
    }
  }
  
  bool access(int i){
    return (blocks[i/BIT_SIZE] >> (i%BIT_SIZE)) & 1ULL;
  }

  int rank(int i){
    BLOCK_TYPE mask = (1ULL << (i%BIT_SIZE)) - 1;
    return s_rank[i/BIT_SIZE] + popcount(mask & blocks[i/BIT_SIZE]);
  }

  int select(int x){
    if(rank((block_size-1) * BIT_SIZE) <= x) return -1;
    int lb = 0, ub = block_size-1;
    while(ub-lb>1){
      int m = (lb+ub)/2;
      if(s_rank[m]<=x) lb = m;
      else ub = m;
    }
    int lbb = lb*BIT_SIZE, ubb = (lb+1)*BIT_SIZE;
    while(ubb-lbb>1){
      int m = (lbb+ubb)/2;
      if(rank(m)<=x) lbb = m;
      else ubb = m;
    }
    return lbb;
  }

  int select0(int x){
    if((block_size-1) * BIT_SIZE - rank((block_size-1) * BIT_SIZE) <= x) return -1;
    int lb = 0, ub = block_size-1;
    while(ub-lb>1){
      int m = (lb+ub)/2;
      if(s0_rank[m]<=x) lb = m;
      else ub = m;
    }
    int lbb = lb*BIT_SIZE, ubb = (lb+1)*BIT_SIZE;
    while(ubb-lbb>1){
      int m = (lbb+ubb)/2;
      if(m-rank(m)<=x) lbb = m;
      else ubb = m;
    }
    return lbb;
  }
};

class EliasFano {
  uint32_t n, u, lgn, lgun;
  FID low, high;

  //ceil(lg(x))
  uint32_t ceillg(double x){
    return ceil(log(x)/log(2));
  }

  uint32_t low_get(size_t i){
    uint32_t ret = 0;
    for(uint32_t p=0; p<lgun; p++){
      ret |= (low.access(lgun*i+p)) << (lgun-1-p);
    }
    return ret;
  }
  
public:
  EliasFano():low(0),high(0){}

  size_t size(){ return n; }
  
  void build(const std::vector<uint32_t>& v){
    if(v.size() == 0) return;
    n = v.size();
    u = v.back();
    lgun = ceillg(u/(double)n);
    lgn = ceillg(u+1) - lgun;

    {//low bits
      low.init(lgun * n);
      int pos = lgun * n - 1;
      for(auto it=v.rbegin(); it!=v.rend(); ++it){
        for(uint32_t b=0; b<lgun; b++){
          if(((*it)>>b)&1){
            low.set(pos);
          }
          pos--;
        }
      }
      low.finalize();
    }
    {//high bits
      high.init((1<<lgn) + n);
      std::vector<uint32_t> cnt(1<<lgn, 0);
      for(uint32_t x : v){
        cnt[(x>>lgun)]++;
      }
      int pos = 0;
      for(uint32_t x : cnt){
        for(uint32_t i=0; i<x; i++){
          high.set(pos);
          pos++;
        }
        pos++;
      }
      high.finalize();
    }
  }
  
  uint32_t access(int i){
    if(n == 0) return 0;
    uint32_t ret = 0;
    ret |= (high.select(i) - i) << lgun;
    ret |= low_get(i);
    return ret;
  }

  uint32_t successor(uint32_t x){
    if(n == 0) return 0;
    if(u < x) return 0;
    uint32_t h = (x>>lgun);
    uint32_t l = x&((1<<lgun)-1);
    uint32_t p1 = (h==0)?0:(high.select0(h-1)-h+1);
    uint32_t p2 = high.select0(h)-h;

    if(x <= access(p1)) return access(p1);
    
    while(p2-p1>1){
      uint32_t m = (p2+p1)/2;
      if(low_get(m) < l) p1 = m;
      else p2 = m;
    }

    return access(p2);
  }

  void dump(){
    std::cout << "n = " << n << ", ";
    std::cout << "u = " << u << ", ";
    std::cout << "ceil(lg(n)) = " << lgn << ", ";
    std::cout << "ceil(lg(u/n)) = " << lgun << std::endl;

    std::cout << "L = ";
    for(uint32_t i=0; i<lgun*n; i++){
      if(low.access(i)) std::cout << 1;
      else std::cout << 0;
    }
    std::cout << std::endl;
    std::cout << "H = ";
    for(uint32_t i=0; i<(1<<lgn)+n; i++){
      if(high.access(i)) std::cout << 1;
      else std::cout << 0;
    }
    std::cout << std::endl;
  }
};

int main(){
  //正の整数の単調増加列(1以上の整数, 同じ整数はなし, ソート済み)
  std::vector<uint32_t> v{3,4,7,13,14,15,21,43};
    
  EliasFano ef;
  ef.build(v);

  ef.dump();

  std::cout << "[access()]" << std::endl;
  for(size_t i=0; i<ef.size(); i++){
    std::cout << i << "\t" << v[i] << "\t" << ef.access(i) << std::endl;
  }

  std::cout << "[successor()]" << std::endl;
  for(uint32_t x=0; x<50; x++){
    std::cout << x << "\t" << ef.successor(x) << std::endl;
  }
  
  return 0;
}

結果

プレゼンの例と同じになっていることを確認。

n = 8, u = 43, ceil(lg(n)) = 3, ceil(lg(u/n)) = 3
L = 011100111101110111101011
H = 1110111010001000
[access()]
0       3       3
1       4       4
2       7       7
3       13      13
4       14      14
5       15      15
6       21      21
7       43      43
[successor()]
0       3
1       3
2       3
3       3
4       4
5       7
6       7
7       7
8       13
9       13
10      13
11      13
12      13
13      13
14      14
15      15
16      21
17      21
18      21
19      21
20      21
21      21
22      43
23      43
24      43
25      43
26      43
27      43
28      43
29      43
30      43
31      43
32      43
33      43
34      43
35      43
36      43
37      43
38      43
39      43
40      43
41      43
42      43
43      43
44      0
45      0
46      0
47      0
48      0
49      0