0905 minimal distance to pi

원문: https://www.hackerrank.com/challenges/minimal-distance-to-pi
참조한 블로그 : www.libragold.com

임의의 무리수에 가장 근접한 유리수를 찾는 문제.
단, 유리수를 \(\frac{p}{q}\)로 나타냈을 때, \(q\)가 가질 수 있는 범위를 한정한다.

일단, Farey sequence[1]를 알면 좋다.
위키피디아에 들어가서 보면 정렬된 그림이 있는데, 이해하기 좋다.
Farey sequence는 수열을 얻는 과정이 쉽고, 근사하는 속도가 매우 빠르다.

본 문제는 \(\pi-3\)에 근사하는 유리수를 찾지만, 임의의 수(\(\in \Q\))에 적용 가능하다.
일단 \(r(q)\)를 ‘\(\pi-3\)에 가장 근접하는 유리수 중 분모가 q이하인 것’을 뜻한다고 하자. 아래는 예시.
\(r(4) = 1/4 \)
\(r(8) = 1/7\) (\(\because 1/7\)이 \(1/8\)보다 \(\pi-3\)에 더 가깝다. \(8/57\) 이전까지 \(1/7\)의 오차가 가장 작다. )
\(r\)의 argument로 (inclusive) range를 허용하면 위는 \(r(1,4), r(1,8)\)이 된다.
구하려고 하는 것은 임의의 \(a\lt b \in \N\)에서 \(r(a,b)\)이다.
\(b-a\)가 매우 크기 때문에 다 해 볼 수는 없고, a와 b사이의 어떤 값이 r(a,b)의 분모가 될 것인지 빨리 알아내는 것이 핵심이다.

#cont.

9/29 언제 정리를… 아는데까지만이라도 적고 싶은데, 그림이 몇개 필요한듯 해서 계속 미루기만 ㅋ

#!/bin/python

import sys

def nextfarey(before,x,C): # real value = x/C
    a,b,c,d = before
    assert a*C/b <= x
    assert c*C/d >= x
    mediant = (a+c)*C/(b+d)
    if mediant < x:
        return a+c, b+d, c, d
    else:
        return a, b, a+c, b+d
 
def fareyseq(x, C, denom_max):
    seq = ((0,1,1,1),)
    while True:
        nf = nextfarey(seq[-1], x, C)
        seq += (nf,)
        if nf[1] >= denom_max or nf[3] >= denom_max:
            return seq
        
def fareypairs(x, C, denom_max):
    seq = fareyseq(x, C, denom_max)
    seqs = ()
    for s in seq:
        s1 = (s[0], s[1])
        s2 = (s[2], s[3])
        if s1 not in seqs:
            seqs += (s1,)
        if s2 not in seqs:
            seqs += (s2,)
    return seqs

def err(p, q, r, s, C=10**100):
        return abs(p*C/q - r*C/s)
        
def sol(lower, upper):
    x= 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
    C= 10**100
    #closest numerator
    q = lower
    limit = upper
    p = int(x*q/C)

    #print seqs
    remainlen = limit - q
    movecnt = 0
    seqs = fareypairs(x, C, limit)
    while remainlen > 0:
        #print 'remainlen=', remainlen
        minerr = err(p, q, x, C)
        #print 'minerr = ', minerr
        #print 'minerr=',minerr
        p1 = 0
        q1 = 0
        for s in seqs:
            #print s
            if q+s[1] > limit :
                break
            newerr = err(p+s[0], q+s[1], x, C)
            #print 'newerr = ', newerr
            if minerr > newerr:
                p1 = s[0]
                q1 = s[1]
                minerr = newerr
                break
        if q1 == 0 :
            break
        else:
            #print p1, q1
            movecnt += 1
            p += p1
            q += q1
            #print 'add:',p1, q1, 'now:',p, q, 'error =', err(p, q, x, C)
            remainlen -= p1
    print '%d/%d'%(p+q*3,q)

        
min,max = raw_input().strip().split(' ')
min,max = [long(min),long(max)]
sol(min, max)