実数関連の問題

ある実数 $\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$ の小数部分を表すものとする。

高速な解き方

この問題は、貪欲的なアルゴリズムによって $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$ となった地点、で計算を打ち切る。

証明

😇

（Handbook of Floating-Point Arithmetic, p436- の Lefèvre algorithm を $O(\log n)$ にしたものです。）

実装コード

以下は、上で挙げた 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

参考: wolfram alpha: 112898^1625668300

memo

実数関連の問題

関連する問題

高速な解き方

証明

実装コード

実行結果