ある実数 $\alpha, \beta, \varepsilon$ と整数 $N$ が与えられたとき、
$$ \min_{1 \le i \le N} \{\alpha i + \beta\}, \\ \mathop{\min \{i \in \mathbb{N} \mid \{\alpha i + \beta\} < \varepsilon\} } $$
を求める問題について考える。ここで、$\{x\}$ は $x$ の小数部分を表すものとする。
関連する問題
例えば、「$2$ の冪乗であって、上 $10$ 桁が $1234567890$ であるものを求めよ」といった問題は、上に示した問題に帰着できる。これを利用すると、(旧)yukicoder No.219 巨大数の概算 に意地の悪い例題を与えることができる。
高速な解き方
この問題は、貪欲的なアルゴリズムによって $O(\log{n})$ 回の演算で解くことができる。
まず、
$$ \begin{aligned} x_0 &:= 1 - \alpha & y_0 := \alpha \\ r_0 &:= \beta & n_0 := 0 \\ u_0 &:= 1 & v_0 := 1 \\ \end{aligned} $$
と置く。以後 $k = 0, 1, \ldots$ について、
$$ \begin{aligned} r_{k} &= q x_k + r_{k+1}\, (0 \le r_{k+1} < r_k) \\ n_{k+1} &:= n_{k} + qu_{k}\\ y_{k} &= q x_k + y_{k+1}\, (0 \le y_{k+1} < y_k) \\ v_{k+1} &:= v_k + qu_{k} \\ x_{k} &= q y_{k+1} + x_{k+1}\, (q := \min\left( \left\lceil{\frac{x_k-r_{k+1}}{y_{k+1}}}\right\rceil, \left\lfloor\frac{x_k}{y_{k+1}}\right\rfloor \right) ) \\ u_{k+1} &:= u_k + qv_{k+1} \\ \end{aligned} $$
を計算し、$n_k > N$ となった地点、あるいは $r_k < \epsilon$ となった地点、で計算を打ち切る。
実装コード
以下は、上で挙げた 2 種の問題の計算をするためのコードです。
$a^b$ の上位桁が $7.999999...$ になるような例を計算するためのコードもついています。
def solve1(N, a, b, mod): """ Return x (0 <= x <= N) which minimizes (a * x + b) % mod. """ x, y, r = -a % mod, a % mod, b % mod u, v = 1, 1 ret = 0 while x: if x <= r: q, r = divmod(r, x) ret += q * u if ret > N: q = (ret - N + u - 1) // u r += q * x ret -= q * u break q, y = divmod(y, x) if y == 0: break v += q * u q = min((x - r + y - 1) // y, x // y) x -= q * y u += q * v return (ret, r) def solve2(a, b, c, mod): """ Return the minimum x such that (a * x + b) % mod <= c. """ x, y, r = -a % mod, a % mod, b % mod u, v = 1, 1 ret = 0 if r > c: while x: if x <= r: q, r = divmod(r, x) ret += q * u if r <= c: q = (c - r) // x r += q * x ret -= q * u break q, y = divmod(y, x) if y == 0: break v += q * u q = min((x - r + y - 1) // y, x // y) x -= q * y u += q * v return (ret, r, r <= c) def generate_hard_instances(): from decimal import Decimal, getcontext N = 2 * 10 ** 9 PREC = 40 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 x = solve1(N, A, B, denom)[0] 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()
実行結果
1 0 1 2 1923400333 7.9999999997771062749887417214 42 690093660 7.9999999997774781888974132388 43 425904478 7.9999999998608999283070353750 232 1951625800 7.9999999999226583521466777485 716 547661426 7.9999999999726129138995758664 1176 978611829 7.9999999999845511531633560940 1267 893290502 7.9999999999861916308309363065 2261 833222115 7.9999999999897369454982145700 2771 1472980651 7.9999999999948054038940741533 8169 1881050581 7.9999999999987087809949940651 10194 924928678 7.9999999999994947191180289428 22703 1133976529 7.9999999999998617962577008326 112898 1625668300 7.9999999999999996857607458081