51nod 1147 连分数

二次无理数的连分数是循环的,循环节从\(a_1\)开始然后到某一个\(a_i=2 a_1\)的时候结束,后面以此循环,即\(\left[a_0;a_1,a_2,a_3,\text{...},2 a_1\right]\),wiki百科里有一个针对\(\sqrt{n}\)的连分数求法,不涉及浮点数求倒数。

\[\begin{align} m_0&=0, & d_0&=1, & a_0&=\lfloor\sqrt{S}\rfloor\\ m_{n+1}&=d_na_n-m_n, & d_{n+1}&=\frac{S-m_{n+1}^2}{d_n}, & a_{n+1}&=\lfloor\frac{a_0+m_{n+1}}{d_{n+1}}\rfloor \end{align} \]

记连分数的循环节长度为\(L\),计算的时候会计算\(\left\lfloor \frac{k+1}{L}\right\rfloor\)个完整循环和从\(a_1\)\(a_ {(k + 1) \% L}\)剩余部分,用矩阵乘法和矩阵快速幂分别算好再合并就好了,注意乘的顺序,别忘了最后的\(a_0\)

def contfrac(n):
    cf=[]
    s=int(n**0.5)
    m0,d0,a0=0,1,s
    cf.append(a0)
    m1=d0*a0-m0
    d1=(n-m1*m1)//d0
    a1=(s+m1)//d1
    cf.append(a1)
    while a1!=2*s:
        m0,d0,a0=m1,d1,a1
        m1=d0*a0-m0
        d1=(n-m1*m1)//d0
        a1=(s+m1)//d1
        cf.append(a1)
    return cf

def mat_mult(A, B, p):
    C = [[0,0],[0,0]]
    for i in range(2):
        for j in range(2):
            for k in range(2):
                C[i][j] += A[i][k]*B[k][j]
            C[i][j] %= p
    return C

def mat_pow(A, p, m):
    if p==0:return [[1,0],[0,1]]
    if p == 1:
        return A
    if p % 2:
        return mat_mult(A, mat_pow(A, p-1, m), m)
    X = mat_pow(A, p//2, m)
    return mat_mult(X, X, m)

n,k=map(int,raw_input().split())
mod=10**9+7
k=k+1
frac=contfrac(n)
Len=len(frac)-1
remain,loop=k%Len,k//Len
mp,mr=[[1,0],[0,1]],[[1,0],[0,1]]
for i in range(Len,0,-1):
    A=[[frac[i],1],[1,0]]
    mp=mat_mult(mp,A,mod)
mp=mat_pow(mp,loop,mod)
for i in range(remain,0,-1):
    A=[[frac[i],1],[1,0]]
    mr=mat_mult(mr,A,mod)
ma=mat_mult(mr,mp,mod)
p,q=ma[1][0]%mod,ma[1][1]%mod
t,q=q,p
p=(frac[0]*p+t)%mod
print "%d/%d"%(p,q)
posted @ 2017-02-22 01:03  czp001  阅读(213)  评论(0编辑  收藏  举报