x / y * y == x となる確率について

Case 1: $x = 2^{p-1}$.

$x = y$ のとき明らかに $(x\oslash y)\otimes y = x$ であるから、以降 $x\lt y$ とする。このとき $\tfrac12\lt x/y\lt 1$ より $\hfloor{x/y} = \tfrac12$ である。

$2^{p-1}\le x/y\cdot 2^p\lt 2^p$ に注意し、 $$ \begin{aligned} |(x\oslash y)-(x/y)| &= |\roundcirc{x/y\cdot 2^p}\cdot 2^{-p} - x/y| \\ &= |\roundcirc{x/y\cdot 2^p} - x/y\cdot 2^p| \cdot 2^{-p} \\ &\le \tfrac12\cdot 2^{-p} = 2^{-p-1} \end{aligned} $$ が成り立つ。よって、ある $|\delta|\le 2^{-p-1}$ を用いて $x\oslash y = x/y+\delta$ と表せ、 $$ \begin{aligned} (x\oslash y)\otimes y &= (x/y+\delta)\otimes y \\ &= \roundcirc{(x/y+\delta)\cdot y} \\ &= \roundcirc{x+\delta y} \end{aligned} $$ である。$y\lt 2^p$ より $|\delta y| \lt \tfrac12$ であるから、 $$ (x\oslash y)\otimes y = \begin{cases} x-\tfrac12, & \text{if~} -\tfrac12 \lt \delta y \lt -\tfrac14; \\ x, & \text{if~} -\tfrac14 \le \delta y \lt \tfrac12 \end{cases} $$ が成り立つ。

Case 2: $x\gt 2^{p-1} \wedge x\le y$.

Case 1 同様にして $x\lt y$ とする。このとき $\tfrac12\lt x/y\lt 1$ より $\hfloor{x/y} = \tfrac12$ である。

ある $|\delta|\le 2^{-p-1}$ に対して $(x\oslash y)\otimes y = \roundcirc{x+\delta y}$ と表せ、$y\lt 2^p$ より $|\delta y| \lt \tfrac12$ なので $(x\oslash y)\otimes y = x$ が従う。

Case 3: $x\gt 2^{p-1} \wedge x\gt y$.

このとき $1\lt x/y\lt 2$ より $\hfloor{x/y} = 1$ である。

$2^{p-1}\le x/y\cdot 2^{p-1}\lt 2^p$ に注意し、 $$ \begin{aligned} |(x\oslash y)-(x/y)| &= |\roundcirc{x/y\cdot 2^{p-1}}\cdot 2^{-p+1} - x/y| \\ &= |\roundcirc{x/y\cdot 2^{p-1}} - x/y\cdot 2^{p-1}| \cdot 2^{-p+1} \\ &\le \tfrac12\cdot 2^{-p+1} = 2^{-p} \end{aligned} $$ が成り立つ。よって、ある $|\delta|\le 2^{-p}$ に対して$(x\oslash y)\otimes y = \roundcirc{x+\delta y}$ と表せる。 $y\lt 2^p$ より $|\delta y| \lt 1$ であるから、 $$ (x\oslash y)\otimes y \in \{x-1, x, x+1\} $$ が従う。$\qed$

Case 1: $r\lt 2^{p-2}$.

このとき $(x\oslash y)\otimes y = x$ を示す。$r/y \lt 2^{p-2}/2^{p-1} = \tfrac12$ が成り立つ。

ある $q\in\Z$ が存在して $x\cdot 2^{p-1}=qy+r$ と表せる。 $2^{p-1}\le x/y\cdot 2^{p-1}\lt 2^p$ より $$ \begin{aligned} x\oslash y &= \roundcirc{x/y\cdot 2^{p-1}}\cdot 2^{-p+1} \\ &= \roundcirc{q+r/y}\cdot 2^{-p+1} \\ &= q\cdot 2^{-p+1} \end{aligned} $$ であり、$r\cdot 2^{-p+1} \lt \tfrac12$ に注意して $$ \begin{aligned} (x\oslash y)\otimes y &= \roundcirc{qy}\cdot 2^{-p+1} \\ &= \roundcirc{x\cdot 2^{p-1}-r}\cdot 2^{-p+1} \\ &= \roundcirc{x-r\cdot 2^{-p+1}} \\ &= x \end{aligned} $$ が成り立つ。

Case 2: $r = 2^{p-2}$.

このとき $x\bmod 2=0 \iff (x\oslash y)\otimes y = x$ を示す。

Case 1 同様にして $$ \begin{aligned} (x\oslash y)\otimes y &= \roundcirc{x-r\cdot 2^{-p+1}} \\ &= \roundcirc{x-\tfrac12} \end{aligned} $$ であり、Property 0.1 より従う。

Case 3: $r\gt 2^{p-2}$.

このとき $r/y\lt \tfrac12 \iff (x\oslash y)\otimes y = x-1$ を示す。

Case 3-1: $r/y\lt \tfrac12$.

$r\cdot 2^{-p+1} \gt \tfrac12$ に注意しつつ Case 1 同様にして $$ \begin{aligned} (x\oslash y)\otimes y &= \roundcirc{x-r\cdot 2^{-p+1}} \\ &= x-1 \end{aligned} $$ が成り立つ。

Case 3-2: $r/y = \tfrac12$.

これを満たす $(x, y)$ は存在しないことを背理法で示す。

$r = \tfrac y2$ とすると、$x\cdot 2^{p-1} = qy+\tfrac y2$、すなわち $x\cdot 2^p = (2q+1)\cdot y$ が成り立つ。 $x$ は整数であるから左辺は $2^p$ で割り切れるが、$y\lt 2^p$ なので右辺は $2^p$ で割り切れないため矛盾。

Case 3-3: $r/y \gt \tfrac12$.

$2^{p-1}\le x/y\cdot 2^{p-1}\lt 2^p$ より $$ \begin{aligned} x\oslash y &= \roundcirc{x/y\cdot 2^{p-1}}\cdot 2^{-p+1} \\ &= \roundcirc{q+r/y}\cdot 2^{-p+1} \\ &= (q+1)\cdot 2^{-p+1}. \end{aligned} $$ $r\lt y$ より $$ \begin{aligned} (x\oslash y)\otimes y &= ( (q+1)\cdot 2^{-p+1})\otimes y \\ &= \roundcirc{( (q+1)\cdot 2^{-p+1})\cdot y} \\ &= \roundcirc{qy + y}\cdot 2^{-p+1} \\ &= \roundcirc{x\cdot 2^{p-1}+y-r}\cdot 2^{-p+1} \\ &= \roundcirc{x + (y-r)\cdot 2^{-p+1})} \\ &\ge \roundcirc{x}\cdot 2^{-p+1} \\ &= x\cdot 2^{-p+1}. \quad\qed \end{aligned} $$

$0\lt y-2^{p-1}\lt 2^{p-1}\lt y$ が成り立つことに注意する。任意の $0\le i\lt k-1$ に対し、 $$ \begin{aligned} r_y(x+i+1) &= (x+i+1)\cdot 2^{p-1}\bmod y \\ &= ( (x+i)\cdot 2^{p-1} + 2^{p-1})\bmod y \\ &= ( (x+i)\cdot 2^{p-1} - (y - 2^{p-1}) )\bmod y \\ &= r_y(x+i) - (y-2^{p-1}) \end{aligned} $$ が成り立つ。よって、$\angled{r_y(x), r_y(x+1), \dots, r_y(x+k-1)}$ は交差 $-(y-2^{p-1})$ の等差数列である。

Case 1: $2^{p-2} \in \{r_y(x), r_y(x+1), \dots, r_y(x+k-1)\}$.

$y - r_y(x+i) = 2^{p-2} \iff r_y(x+i+1) = 2^{p-2}$ であることと $(x+i)\bmod 2 = 1$ と $(x+i+1)\bmod 2 = 1$ のうちどちらか片方のみが成り立つことから従う。

Case 2: $2^{p-2} \notin \{r_y(x), r_y(x+1), \dots, r_y(x+k-1)\}$.

Lemmata 2–3 を踏まえると、$( (x+i)\oslash y)\otimes y\ne x+i$ であることは $r_y(x+i)$ が $\angled{2^{p-2}, \dots, y-2^{p-2}}$ に含まれることと同値である。 $\angled{2^{p-2}, \dots, y-2^{p-2}}$ は連続する $y-2^{p-1}+1$ 個の整数であり、交差 $-(y-2^{p-1})$ であり $2^{p-2}$ を含まない等差数列との共通部分はちょうど 1 つである。$\qed$

下記を満たす連続部分列の個数は、高々 1 つである。

• 各 $0\le i\lt k-1$ に対して $r_y(x+i) \gt r_y(x+i+1)$
• $x-1 = y$
• $r_y(x+k-1) \lt r_y(x+k)$

また、下記についても同様である。

• 各 $0\le i\lt k-1$ に対して $r_y(x+i) \gt r_y(x+i+1)$
• $r_y(x-1) \lt r_y(x)$
• $x+k-1 = 2^p$

これらの連続部分列に対応する $i$ のうち $( (x+i)\oslash y)\otimes y \ne x+i$ となるものは高々 1 つである。

各 $i$ に対して、$q_i$ が存在して $r_y(y+i) = r_y(y) - i(y-2^{p-1}) + q_iy$ が成り立つ。求める値は $q_{2^p-y}+O(1)$ となる。 $$ \begin{aligned} q_{2^p-y}y &= r_y(2^p) - r_y(y) + (2^p-y)(y-2^{p-1}) \\ &= (2^{2p-1}\bmod y) - 0 + (2^p-y)(y-2^{p-1}) \\ &= (2^p-y)(y-2^{p-1}) + O(y) \end{aligned} $$ であり、 $$ \begin{aligned} \frac{(y-2^{p-1})(2^p-y)}y &= \frac{2^p\cdot y - y^2 -2^{2p-1} + 2^{p-1}\cdot y}y \\ &= -y + 3\cdot 2^{p-1} - 2^{2p-1}/y \end{aligned} $$ であるから、 $$ q_{2^p-y} = -y + 3\cdot 2^{p-1} - 2^{2p-1}/y + O(1) $$ となる。$\qed$

余事象を考える。 まず、Lemma 1 より $( (x\oslash y)\otimes y)-x = -\tfrac12$ なる確率は $O(2^{-p})$ である。

$|( (x\oslash y)\otimes y)-x| = 1$ なる通り数は次の通りである。 $$ \begin{aligned} &\phantom{{}={}} \sum_{y=2^{p-1}+1}^{2^p-1} {-y + 3\cdot 2^{p-1} - 2^{2p-1}/y + O(1)} \\ &= -\tfrac12(2^{p-1}-1)(2^p+2^{p-1}) + 3\cdot 2^{p-1}\cdot (2^{p-1}-1) - 2^{2p-1}\cdot (H_{2^p-1}-H_{2^{p-1}}) + O(1) \\ &= 3\cdot 2^{p-1}\cdot (-\tfrac12(2^{p-1}-1) + (2^{p-1}-1) ) - 2^{2p-1}\cdot (H_{2^p-1}-H_{2^{p-1}}) + O(1) \\ &= 3\cdot 2^{p-1}\cdot \tfrac12(2^{p-1}-1) - 2^{2p-1}\cdot (H_{2^p-1}-H_{2^{p-1}}) + O(1) \\ &= 3\cdot 2^{p-2}\cdot (2^{p-1}-1) - 2^{2p-1}\cdot (H_{2^p-1}-H_{2^{p-1}}) + O(1) \\ &= 3\cdot (2^{2p-3}-2^{p-2}) - 2^{2p-1}\cdot (H_{2^p-1}-H_{2^{p-1}}) + O(1) \\ &= 3\cdot 2^{2p-3} - 2^{2p-1}\cdot (\ln(2^p)-\ln(2^{p-1}) ) + O(2^p) \\ &= 3\cdot 2^{2p-3} - 2^{2p-1}\ln(2) + O(2^p) \\ &= 2^{2(p-1)} \cdot \left(\tfrac32-\ln(4)\right) + O(2^p). \end{aligned} $$ これを $2^{2(p-1)}$ で割って、$\tfrac32-\ln(4)+O(2^{-p})$ を得る。これらより、求める確率は $$ 1-\left(\tfrac32-\ln(4)+O(2^{-p})\right) = \ln(4)-\tfrac12+O(2^{-p}) $$ となる。$\qed$

相加相乗平均の関係に注意して、$|( (x\oslash y)\otimes y)-x| = 1$ なる通り数は $$ \begin{aligned} &\phantom{{}={}} \max_y {(-y + 3\cdot 2^{p-1} - 2^{2p-1}/y)} \\ &= 3\cdot 2^{p-1} - \min_y {(y + 2^{2p-1}/y)} \\ &= 3\cdot 2^{p-1} - 2\cdot (2^{2p-1})^{1/2} \\ &= 3\cdot 2^{p-1} - 2^{p+1/2} \\ &= 3\cdot 2^{p-1} - 2\sqrt 2\cdot 2^{p-1} \\ &= (3-2\sqrt 2)\cdot 2^{p-1} \end{aligned} $$ である。上界を達成するのは $y=2^{p-1/2}$ である。$( (x\oslash y)\otimes y)-x = -\tfrac12$ なる確率は $O(2^{-p})$ であるから、求める確率は $$ 1 - ( (3-2\sqrt 2) + O(2^{-p}) ) = 2\sqrt 2-2 + O(2^{-p}) $$ である。$\qed$

#include <iomanip>
#include <iostream>
#include <map>
#include <quadmath.h>
#include <random>

template <typename F> F random_1_2(std::mt19937& mt) { return F(); }

template <> float random_1_2(std::mt19937& mt) {
  return (mt() & ~(~0u << 23)) | 1u << 23;
}

template <> double random_1_2(std::mt19937& mt) {
  unsigned long res = mt();
  res <<= 32;
  res |= mt();
  res &= ~(~0uL << 52);
  res |= 1uL << 52;
  return res;
}

template <> long double random_1_2(std::mt19937& mt) {
  unsigned long res = mt();
  res <<= 32;
  res |= mt();
  res |= 1uL << 63;
  return res;
}

template <> __float128 random_1_2(std::mt19937& mt) {
  unsigned __int128 res = mt();
  res <<= 32;
  res |= mt();
  res <<= 32;
  res |= mt();
  res <<= 32;
  res |= mt();
  res &= ~(~(unsigned __int128)0u << 112);
  res |= (unsigned __int128)(1u) << 112;
  return res;
}

template <typename F> void count_f(long n) {
  std::map<F, unsigned long> c;
  std::cerr << std::fixed << std::setprecision(1);
  std::random_device rd;
  std::mt19937 mt(rd());
  for (int i = 0; i < n; ++i) {
    volatile F x = random_1_2<F>(mt);
    volatile F y = random_1_2<F>(mt);
    volatile F z = F(x / y) * y;
    ++c[F(z - x)];
  }
  std::cout << "---\n";
  for (auto [x, y]: c) {
    std::cout << double(x) << '\t' << y << '\n';
  }
}

int main() {
  long const N = 1e9;
  count_f<float>(N);
  count_f<double>(N);
  count_f<long double>(N);
  count_f<__float128>(N);
}