結果

問題 No.313 π
ユーザー Min_25
提出日時 2015-12-06 00:04:17
言語 PyPy2
(7.3.15)
結果
AC  
実行時間 726 ms / 5,000 ms
コード長 4,696 bytes
コンパイル時間 843 ms
コンパイル使用メモリ 77,056 KB
実行使用メモリ 114,688 KB
最終ジャッジ日時 2024-09-14 14:33:33
合計ジャッジ時間 27,880 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge4
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 2
other AC * 32
権限があれば一括ダウンロードができます

ソースコード

diff # 

import sys
import math
from math import sqrt

def isqrt(n):
  if n <= 0:
    return 0

  x = int(sqrt(n) * (1 + 1e-12))
  while True:
    y = (x + n // x) >> 1
    if y >= x:
      return x
    x = y

def fast_str(N):
  def rec(N, e, zerofill=False):
    if e <= 10:
      if not zerofill:
        strs.append(str(N))
      else:
        strs.append("%0*d" % (1 << (e + 1), N))
    else:
      shamt = 1 << e
      q, r = fast_divmod(N >> shamt, fives[e])
      r = (r << shamt) | (N & masks[e])
      if zerofill or q:
        rec(q, e - 1, zerofill)
        rec(r, e - 1, True)
      else:
        rec(r, e - 1, zerofill)

  strs = []

  if N < 0:
    N = -N
    strs.append("-")

  fives = []
  masks = []

  five = 5
  mask = 1
  
  e = 0

  while 1:
    shamt = 1 << e
    if (five << shamt) > N:
      break
    fives.append(five)
    masks.append(mask)
    if (five.bit_length() + shamt) * 2 - 1 > N.bit_length():
      break

    e += 1
    five *= five
    mask |= mask << shamt

  rec(N, e)
  return "".join(strs)

def fast_divmod(a, b):
  """
  Burnikel & Ziegler
  """

  DIVMOD_THRESHOLD = 8000

  def _pack(pack, shamt):
    size = len(pack)
    while size > 1:
      npack = []
      for i in range(0, size - 1, 2):
        npack.append(pack[i] | (pack[i+1] << shamt))
      if size & 1:
        npack.append(pack[-1])
      pack = npack
      size = (size + 1) >> 1
      shamt <<= 1
    return pack[0]

  def _unpack(M, size, shamt):
    needed_sizes = []
    s = size
    while s > 1:
      needed_sizes.append(s)
      s = (s + 1) >> 1
    ret = [M]
    for needed_size in needed_sizes[::-1]:
      mask = (1 << shamt) - 1
      nret = []
      for c in ret:
        nret.append(c & mask)
        nret.append(c >> shamt)
      ret = nret[:needed_size]
      shamt >>= 1
    return ret

  def _div32(a21, a0, b, b1, b0, n):
    if (a21 >> n) >= b1:
      q, r = (1 << n) - 1, a21 - (b1 << n) + b1
    else:
      q, r = _div21(a21, b1, n)
    r = ((r << n) | a0) - q * b0
    while r < 0:
      q -= 1
      r += b
    return q, r

  def _div21(a, b, n):
    if a < b:
      return (0, a)
    if n <= DIVMOD_THRESHOLD:
      return divmod(a, b)
    odd = n & 1
    if odd:
      a, b, n = a << 1, b << 1, n + 1
    nh = n >> 1
    mask = (1 << nh) - 1
    b1, b0 = b >> nh, b & mask
    q1, r = _div32(a >> n, (a >> nh) & mask, b, b1, b0, nh)
    q0, r = _div32(r, a & mask, b, b1, b0, nh)
    if odd:
      r >>= 1
    return (q1 << nh) | q0, r

  if a < b:
    return (0, a)

  m = a.bit_length()
  n = b.bit_length()

  if n <= DIVMOD_THRESHOLD:
    return divmod(a, b)

  nqs = m // n
  a_blocks = _unpack(a, nqs + 1, (n << (nqs.bit_length() - 1)))

  qs = []
  r = a_blocks[-1]
  for i in range(nqs - 1, -1, -1):
    q, r = _div21((r << n) | a_blocks[i], b, n)
    qs.append(q)
  q = _pack(qs[::-1], n)
  return q, r

def fast_isqrt(a, threshold=1024):
  """
  GMP Karatsuba Square Root
  """

  def _isqrt(a):
    if a < threshold:
      r = isqrt(a)
      return r, a - r * r

    n = a.bit_length()
    adjusted = ((n - 1) & 2) == 0
    if adjusted:
      a <<= 2
      n += 2

    nq = (n + 1) >> 2
    nh = nq << 1
    mask = (1 << nh) - 1
    ah, al = a >> nh, a & mask
    mask >>= nq

    a1, a0 = al >> nq, al & mask
    s1, r1 = _isqrt(ah)

    q, u = fast_divmod((r1 << nq) + a1, s1 << 1)
    s = (s1 << nq) + q
    r = (u << nq) + a0 - q * q
    if r < 0:
      r = r + (s << 1) - 1
      s -= 1

    if adjusted:
      r >>= 2
      if s & 1:
        s >>= 1
        r += s + 1
      else:
        s >>= 1

    return s, r

  return _isqrt(a)[0]

def pi_chudnovsky_bs(digits):
  # (c) Nick Craig-Wood
  def bs(a, b):
    if b - a == 1:
      if a == 0:
        Pab = Qab = 1
      else:
        Pab = (6 * a - 5) * (2 * a - 1) * (6 * a - 1)
        Qab = a * a * a * C3_OVER_24
      Tab = Pab * (13591409 + 545140134 * a)
      if a & 1:
        Tab = -Tab
    else:
      m = (a + b) // 2
      Pam, Qam, Tam = bs(a, m)
      Pmb, Qmb, Tmb = bs(m, b)
      Pab = Pam * Pmb
      Qab = Qam * Qmb
      Tab = Qmb * Tam + Pam * Tmb
    return Pab, Qab, Tab

  C = 640320
  C3_OVER_24 = C ** 3 // 24
  DIGITS_PER_TERM = math.log10(C3_OVER_24 / 72.)
  N = int(digits / DIGITS_PER_TERM + 1)
  _, Q, T = bs(0, N)

  shamt = 0
  # prec = int(digits * 3.3219280948873626)
  # shamt = Q.bit_length() - prec
  numer = (Q >> shamt) * 426880 * fast_isqrt(10005 * 100 ** digits)
  denom = T >> shamt
  return fast_divmod(numer, denom)[0]

def prob313():
  rl = sys.stdin.readline
  st = pi_chudnovsky_bs(200000)
  s = fast_str(st)
  s = s[0:1] + "." + s[1:]

  A = rl().rstrip()
  for a, b in zip(A, s):
    if a != b:
      print("%s %s" % (a, b))
      break

prob313()
0