ある実数 $\alpha, \beta, \varepsilon > 0$ と整数 $N$ が与えられたとき、$$ \min_{0 \le i \le N} \{ \{\alpha i + \beta\} \}, \\ \mathop{\min \{i \mid i \in \mathbb{N}, \{\alpha i + \beta\} \lt \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