Day 159-Q4 — 2026-09-24 — HyperLogLog（基数推定・「先頭の 0 の長さの最大値」から相異なる要素数を 16 KB で当てる）

問題

長さ $N$ の整数列 $A_1,\dots,A_N$ と精度パラメータ $b$ が与えられる。以下に定める HyperLogLog 推定値 $E$ を計算し、$\lfloor E\rfloor$ を出力せよ（ハッシュと計算手順を固定するので答えは一意に決まる）。

ハッシュ: 64 ビット演算（$\bmod 2^{64}$）で

h(x): x += 0x9E3779B97F4A7C15
      x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9
      x = (x ^ (x >> 27)) * 0x94D049BB133111EB
      return x ^ (x >> 31)

（splitmix64 の finalizer。Vigna, http://xorshift.di.unimi.it/splitmix64.c ）

手順: $m=2^b$ 個のレジスタ $M[0..m-1]$ を 0 で初期化。各 $A_i$ について $h=h(A_i)$ とし、上位 $b$ ビットを添字 $j$、残り $64-b$ ビットを $w$ とする。$\rho(w)$ を「$w$ の最上位側から見て最初に 1 が現れる位置（1-indexed）」とし（$w=0$ なら $65-b$）、$M[j]\leftarrow\max(M[j],\rho(w))$。

$$E_{raw}=\alpha_m\, m^2 \Big/ \sum_{j} 2^{-M[j]},\qquad \alpha_{16}=0.673,\ \alpha_{32}=0.697,\ \alpha_{64}=0.709,\ \alpha_m=\frac{0.7213}{1+1.079/m}\ (m\ge128)$$

$E_{raw}\le 2.5m$ かつ $V=\#\{j: M[j]=0\}>0$ なら $E = m\ln(m/V)$（線形計数）、それ以外は $E=E_{raw}$。

入力形式

N b
A_1 A_2 ... A_N

制約

入出力例

10 4
1 2 3 2 1 4 5 5 5 6

4

（真の相異なる要素数は 6。$m=16$ では 6 要素が 4 個のレジスタにしか当たらず $V=12$、$E=16\ln(16/12)=4.60$。$m=16$ の理論誤差は $1.04/\sqrt{16}=26\%$ なので、この程度の外れは想定内——後述の実測表参照）

この問題は Flajolet, Fusy, Gandouet, Meunier「HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm」（AofA 2007、https://algo.inria.fr/flajolet/Publications/FlFuGaMe07.pdf 、2026-09-24 参照）の Fig.2 のアルゴリズムを、ハッシュを固定して決定的にしたものである。Redis の PFADD/PFCOUNT は $b=14$（16384 レジスタ、12 KB）でこれを実装している（https://redis.io/docs/latest/commands/pfcount/ 、同上参照。Redis は Google の HLL++ 由来のバイアス補正も併用）。

概念図：ハッシュの分割とレジスタ

ヒント（クリックで段階的に開示）

shift = 64 - b
for v in A:
    h = splitmix64(v)
    j = h >> shift                 # 上位 b ビット
    w = h & ((1 << shift) - 1)     # 残り
    r = shift - w.bit_length() + 1 # 先頭 0 の数 + 1（w == 0 なら shift + 1）
    if r > M[j]: M[j] = r
E = alpha(m) * m * m / sum(2.0 ** -x for x in M)
if E <= 2.5 * m and (V := M.count(0)) > 0: E = m * log(m / V)

模範解答 (Python)

📋 コピーimport sys, math

MASK = (1 << 64) - 1


def splitmix64(x):
    x = (x + 0x9E3779B97F4A7C15) & MASK
    x = ((x ^ (x >> 30)) * 0xBF58476D1CE4E5B9) & MASK
    x = ((x ^ (x >> 27)) * 0x94D049BB133111EB) & MASK
    return x ^ (x >> 31)


def alpha(m):
    if m == 16: return 0.673
    if m == 32: return 0.697
    if m == 64: return 0.709
    return 0.7213 / (1 + 1.079 / m)


def main():
    data = sys.stdin.buffer.read().split()
    n, b = int(data[0]), int(data[1])
    m = 1 << b
    shift = 64 - b
    low_mask = (1 << shift) - 1
    M = [0] * m
    for i in range(n):
        h = splitmix64(int(data[2 + i]))
        j = h >> shift
        w = h & low_mask
        r = shift - w.bit_length() + 1      # w == 0 のとき bit_length() == 0 で r = shift + 1
        if r > M[j]:
            M[j] = r
    Z = 0.0
    for x in M:
        Z += 2.0 ** (-x)
    E = alpha(m) * m * m / Z
    if E <= 2.5 * m:
        V = M.count(0)
        if V > 0:
            E = m * math.log(m / V)          # 線形計数（小さい基数の補正）
    print(math.floor(E))


main()

Step-by-Step 解説

迷った他案と採らなかった理由

よくあるミス

次のステップ

• 発展1: 2 つの HLL を レジスタごとの max でマージすると和集合の推定になる（分散カウントの根拠）。100 個に分割して個別に作りマージした推定と、1 本で作った推定が ビット単位で一致 することを確かめる
• 発展2: 包除 $|A\cap B|=|A|+|B|-|A\cup B|$ で積集合を推定し、誤差が和集合の誤差の差で膨らむことを実測する（積集合には MinHash の方が向く）
• 扱わなかった範囲: HLL++ の疎表現とバイアス補正表、レジスタを 6 ビットに詰めるメモリ最適化、削除（HLL は削除不可。必要なら Counting 版や Sliding HLL）

自己評価

自分の回答

気づき・メモ