実数関連の問題

$\min _{1 \le x \le N} \{\alpha x + \beta\}$ 問題

yukicoder No.219 巨大数の概算 など色々.

(2017/4/27: もっと効率のいい方法があるので直さないと!)

凡例

実数 $x$ に対して,その小数部を $\{ x \}$ で表す.

実数が絡んでくる問題

実数 $\alpha, \beta$ と正の整数 $N$ に対し,
\[
K^* := \mathop{\rm arg~min}\limits_{1 \le k \le N} \{ \{\alpha k + \beta\} \}
\]
を高速で求めることができると,誤差評価の厳しい例の生成に役立つ.

例えば,$a^k$ の先頭桁ができるだけ多く $9999...$ となっているような $k$ を見つけたい場合は $(\alpha, \beta) = (-\log_{10}{a}, 0)$ と置き,$700000....$ となっているような $k$ を見つけたい場合は $(\alpha, \beta) = (\log_{10}{a}, -\log_{10}{7})$ と置けばよい.

解き方

$N$ が固定されているならば,(特殊な場合を除き)$\alpha, \beta$ を十分な精度の分数で近似して解いても $K^*$ は変化しない.つまり,$\alpha \approx a/n, \beta \approx b/n$ とおいて,
\[
L^* := \mathop{\rm arg~min}\limits_{1 \le k \le N} \{ak + b \pmod{n}\}
\]
\[
= \mathop{\rm arg~min}\limits_{1 \le k \le N} \left\{ak + b - n\left\lfloor \frac{ak+b}{n} \right\rfloor \right\}
\]
を解けばよい.ここで,$b$ は定数なので無視すると,
\[
\mathop{\rm arg~min}\limits_{L \le k \le R} \left\{ak + b\left\lfloor \frac{ck+d}{e} \right\rfloor \right\}
\]
を解く問題となるが,これはユークリッドの互助法のような方法で $O(\log(\min\{c, e\}))$ 回ほどの算術演算を用いて解くことができる.

詳細は,kevinsogo さんの minion-of-the-year の Editorial を参照してください(要 login & unlock).

実装

def minimize_abcde(a, b, c, d, e, l, r):
    """
    Minimize {ax + b * floor( (cx + d) / e )} for l <= x <= r.
    Assume that e > 0 and l <= r.

    Returns:
        tuple(min, min(argmin))
    """
    if not (0 <= c < e):
        q, c = divmod(c, e) # 0 <= c < e
        return minimize_abcde(a + b * q, b, c, d, e, l, r)

    if not (0 <= d < e):
        q, d = divmod(d, e) # 0 <= d < e
        y, x = minimize_abcde(a, b, c, d, e, l, r)
        return (y + b * q, x)

    q1, q2 = (c * l + d) // e, (c * r + d) // e

    if q1 == q2:
        if a >= 0:
            return (a * l + b * q1, l)
        else:
            return (a * r + b * q2, r)

    if a > 0:
        y1, x1 = a * l + b * q1, l
        y2, x2 = minimize_abcde(b, a, e, c - d - 1, c, q1 + 1, q2)
        x2 = (e * x2 + c - d - 1) // c
        return (y1, x1) if y1 <= y2 else (y2, x2)
    else:
        y1, x1 = a * r + b * q2, r
        y2, x2 = minimize_abcde(b, a, e, e - d - 1, c, q1, q2 - 1)
        x2 = (e * x2 + e - d - 1) // c
        return (y1, x1) if y1 < y2 else (y2, x2) 

def minimize_Ax_plus_B_mod(a, b, mod, l, r):
    """
    Minimize {ax + b (modulo mod)} for l <= x <= r.

    Returns:
        tuple(min, min(argmin))
    """
    if mod == 1:
        return (0, l)
    
    a, b = a % mod, b % mod

    if a == 0:
        return (b, l)

    assert(0 <= l <= r < mod)

    y, x = minimize_abcde(a, -mod, a, b, mod, l, r)
    return (y + b, x)

def generate_hard_instances():
    from decimal import Decimal, getcontext

    N = 2 * 10 ** 9

    PREC = 50
    getcontext().prec = PREC

    denom = 10 ** PREC
    beta = Decimal(8).log10() - Decimal("1e-%d" % (PREC - 5))

    best = 0.0    
    for a in range(1, 2 * 10 ** 5):
        alpha = Decimal(a).log10()

        A = int(-alpha * denom) % denom
        B = int(beta * denom) % denom
        y, x = minimize_Ax_plus_B_mod(A, B, denom, 1, N)

        alpha *= x
        alpha -= int(alpha)
        prefix = pow(Decimal(10), alpha)
        if prefix > best:
            best = prefix
            print("%10d %10d %s" % (a, x, str(prefix)[:PREC - 10]))

generate_hard_instances()

$a^k$ が $7.99999... \times 10^e$ となるような例.

         1          1 1
         2 1923400333 7.99999999977710627498874172148782077364
        42  690093660 7.99999999977747818889741323883683423318
        43  425904478 7.99999999986089992830703537508895098551
       232 1951625800 7.99999999992265835214667774852499031572
       716  547661426 7.99999999997261291389957586642539855324
      1176  978611829 7.99999999998455115316335609401292081124
      1267  893290502 7.99999999998619163083093630649515095650
      2261  833222115 7.99999999998973694549821457004902549693
      2771 1472980651 7.99999999999480540389407415332527932469
      8169 1881050581 7.99999999999870878099499406511312160766
     10194  924928678 7.99999999999949471911802894284496606583
     22703 1133976529 7.99999999999986179625770083267966746172
    112898 1625668300 7.99999999999999968576074580814870538435

wolfram alpha: 112898^1625668300

その他

$\beta$ が 0 の場合は,$\alpha$ の最良近似分数を得た後に拡張ユークリッドの互除法などを使うことで, 良い $k$ が得られることもあります.