fixed.py

from dataclasses import dataclass
from math import gcd
from typing import Self, Tuple


@dataclass
class Fraction:
  p: int
  q: int

  def assert_coprime(self):
    assert gcd(self.p, self.q) == 1

  def __lt__(self, other: Self):
    return self.p * other.q < self.q * other.p

  def __sub__(self, other: Self):
    return Fraction(self.p * other.q - self.q * other.p, self.q * other.q)

  def __mul__(self, other: Self):
    return Fraction(self.p * other.p, self.q * other.q)


def make_fraction(a: int, b: int):
  g = gcd(a, b)
  return Fraction(a // g, b // g)


def make_power(a: int, b: int):
  return Fraction(a ** b, 1) if b >= 0 else Fraction(1, a ** -b)


def extended_gcd(a: int, b: int) -> Tuple[int, int, int]:
  if a == 0:
    return (b, 0, 1)
  g, x, y = extended_gcd(b % a, a)
  return (g, y - (b // a) * x, x)


def get_bounds(x: Fraction, nmax: int) -> Tuple[Fraction, Fraction]:
  x.assert_coprime()

  if x.q <= nmax:
    g, a, b = extended_gcd(x.p, x.q)
    assert g == 1 and a * x.p + b * x.q == 1
    v = (-a) % x.q
    v += (nmax - v) // x.q * x.q
    assert v <= nmax
    return (x, make_fraction(x.p * v + 1, x.q * v))

  cfrac = []
  p, q = x.p, x.q
  while 1:
    d, r = p // q, p % q
    cfrac.append(d)
    if r != 0:
      p, q = q, r
    else:
      break

  prev = Fraction(0, 1)
  curr = Fraction(1, 0)

  for i, a in enumerate(cfrac):
    nxt = Fraction(prev.p + a * curr.p, prev.q + a * curr.q)
    if nxt.q > nmax:
      s = (nmax - prev.q) // curr.q
      f = Fraction(prev.p + s * curr.p, prev.q + s * curr.q)
      if i % 2:
        return (curr, f)
      else:
        return (f, curr)
    prev, curr = curr, nxt
  assert False


# from random import randint as ri

# def get_bounds_stupid(x: Fraction, nmax: int) -> Tuple[Fraction, Fraction]:
#   lower = Fraction(0, 1)
#   upper = Fraction(1, 0)
#   for n in range(1, nmax + 1):
#     lower = max(lower, make_fraction(n * x.p // x.q, n))
#     upper = min(upper, make_fraction(n * x.p // x.q + 1, n))
#   return (lower, upper)

# while 1:
#   x, y, z = ri(1, 1000), ri(1, 1000), ri(1, 1000)
#   print(x, y, z)
#   assert get_bounds_stupid(make_fraction(x, y), z) == get_bounds(make_fraction(x, y), z)

q = 64
p = 52

eta = 1

nmax = 2 ** (p + 1) - 1
bias = 2 ** (q - p - 2) - 1
emin = 1 - bias - p
emax = 2 ** (q - p - 1) - 2 - bias - p

cnt = 0
maxq = 0

for e in range(emin, emax + 1):
  digits_min = 0 if e >= 0 else e
  digits_max = int(f"{float(2 ** e * (2 ** (p + 1) - 1)):e}".split("e")[1])
  print(e, digits_min, digits_max, flush=True)

  for digit in range(digits_min, digits_max + 1):
    # last digit in segment is `digit`
    left, right = get_bounds(make_power(2, e - digit) * make_power(5, -digit), nmax)
    distance = right - left
    l = distance.p
    r = 10 ** eta * distance.q
    q = max(0, r.bit_length() - l.bit_length())
    if (l << q) < r:
      q += 1
      assert l << q >= r
    maxq = max(maxq, q)
    cnt += 1

print(cnt, maxq)