'''
Copyright (c) 2019 Serguei Zolotov

This code was built using following publication:
E. Charrier and L. Buzer, Approximating a real number by a rational number with a limited denominator: 
A geometric approach, Discrete Applied Mathematics, Volume 157 (2009) Pages 3473-3484.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as
published by the Free Software Foundation, either version 3 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

'''
#!/bin/python

# get rational close to pi
# use continued fraction as best approximation
# denominator shall be large enough (square of) compare to upper boundary
# use b001203.txt - A001203 sequence that contains cont fraction
def getpiratio():

....noms = [0, 1]
....denoms = [1, 0]
....fraction = []

....# load fraction coefficients
....# used first 5000 numbers - it gives 
....# denominator length of 2564 decimal digits
....# it is way more than we need: length of 2^1024 - 309 digits
....# open b-file or similar source
....fr = open("b001203.txt",'r')
....line = fr.readline()
....while line:
........cnt, value = line.split(" ")
........fraction.append(int(value))
........line = fr.readline()

....fr.close()

....# compute principal ratios
....for i in range(2, len(fraction)):
........a = fraction[i - 2]
........noms.append(a * noms[i - 1] + noms[i - 2])
........denoms.append(a * denoms[i - 1] + denoms[i - 2])

....return (noms[-1], denoms[-1])

# small utilities for computations
# we will be working with vectors

# dot product
def dotP(v1, v2):
....return v1[0]*v2[0] + v1[1]*v2[1]

# deternimant
def detP(v1, v2):
....return v1[0]*v2[1] - v1[1]*v2[0]

# square of length
def lenP2(v):
....return v[0]*v[0] + v[1]*v[1]

# addition with multiplier/sign
def addV(v1, v2, s):
....return (v1[0] + s*v2[0], v1[1] + s*v2[1])

# calculate distance to intersection of:
# line from point p0 - into direction of v
# with line L
# L  - boundary(line) if L == 0 - it is y = PI*x line
# e - number we are trying to approximate (will be pi in our case)
def intersect(p0, v, L, e):
....# simple calculation of intersection
....# a1x + b1y = c1
....# a2x + b2y = c2
....# (x, y) = (|c b|/|a b|, |a c|/|a b|)
....p1 = addV(p0, v, 1)

....a1, b1, c1 = p1[1] - p0[1], p0[0] - p1[0], detP(p0, p1)
....if L == 0:
........a2, b2, c2 = e[0], -e[1], 0 # intersection with y = pi*x
....else:
........a2, b2, c2 = 1, 0, L # intersection with right boundary

....# solve the system
....cxn = detP((c1, c2), (b1, b2))
....cyn = detP((a1, a2), (c1, c2))
....d = detP((a1, a2), (b1, b2))

....# intersection vector
....dvec = [cxn - p0[0]*d, cyn - p0[1]*d]
....# calculate projection
....n = dotP(dvec, v) 
....# return integer normalized coefficient according to original algorithm
....return n // (lenP2(v)*d)

# get next direction
def getUij(i, j, A, B):
....return addV(B[j], A[i], -1)

# initialize vectors
# t - left bound
# e - approximation target
def initV(t, e):
....# with Pi slope in mind checking only 2 options:
....if e[0] * t + e[0] < e[1]*(e[0] * t // e[1]) + e[1]*4:
........return (1, 3)
....else:
........return (1, 4)

# step update
def fupdate(p, v, u, g, e):
....p0 = addV(p, v, 1)
....k = intersect(p0, u, g, e)
....new_v = addV(v, u, k)
....return new_v

# main algorithm
def solve(l, r, vip):
....# initialize Pi
....pi_floor = vip[0] * l // vip[1]
....A = [(l, pi_floor + 1)]
....B = [(l, pi_floor)]
....# initialize steps
....U = {}
....U[(0,0)] = getUij(0, 0, A, B)
....V = {}
....V[(0,0)] = initV(l, vip)
....# set main direction
....N = (-vip[0], vip[1])
....i, j = 0, 0
....while True:
........vij = V[(i, j)]
........if dotP(vij, N) < 0:
............b = min(intersect(A[i], vij, 0, vip), intersect(A[i], vij, r, vip))
............if b >= 1:
................W = addV(A[i], vij, b)
................A.append(W)
................i = i+1
................# update
................uij = getUij(i, j, A, B)
................vij = fupdate(W, vij, uij, 0, vip)
................V[(i, j)] = vij
............else:
................break
........elif dotP(vij, N) > 0:
............b = min(intersect(B[j], vij, 0, vip), intersect(B[j], vij, r, vip))
............if b >= 1:
................W = addV(B[j], vij, b)
................B.append(W)
................j = j+1
................# update
................uij = getUij(i, j, A, B)
................# this step is required but was not illustarted in paper
................uij = addV((0, 0), uij, -1)
................vij = fupdate(W, vij, uij, 0, vip)
................V[(i, j)] = vij
............else:
................break
........else:
............break

....# here we need to check are we closer to pi from upper or lower curves
....# get upper and lower bounds numbers
....opt1 = A[-1]
....opt2 = B[-1]

....# it may seems confusing here
....# vip is a pi as {nom, denom} while optX as {denom, nom}
....# pick the right index for calculations
....d1 = abs(vip[0]*opt1[0] - vip[1]*opt1[1])
....d2 = abs(vip[0]*opt2[0] - vip[1]*opt2[1])
....if d1*opt2[0] < d2*opt1[0]:
........opt = opt1
....else:
........opt = opt2

....# return {nom, denom}
....return [opt[1], opt[0]]

if __name__ == '__main__':

....# get Pi
....vip = getpiratio()

....# open output files
....f8 = open("b325158_1.txt",'w')
....f9 = open("b325159_1.txt",'w')

....# loop through boundaries
....for n in range(1024):
........bl = 2**n
........br = 2**(n+1)-1
........ret = solve(bl, br, vip)
........print ("%d %d"%(n, ret[0]), file=f8)
........print ("%d %d"%(n, ret[1]), file=f9)

....# close output
....f8.close()
....f9.close()