0192 最优有理逼近
* * * *
拉格朗日计划
* * * *
最优有理逼近

实数x的分母上限为$d\in\mathbb N$的最优有理逼近数,是一个最简分数$r/s$,其中$r,s\in\mathbb Z$,且$1\le s\le d$,使得只要$p,q\in\mathbb Z$满足 $$|\frac{p}{q}-x| < |\frac{r}{s}-x|$$ 就有$q>d$。

例如,$\sqrt{13}$的分母上限为20的最优有理逼近数是$18/5$,分母上限为$30$的最优有理逼近数是$101/28$。

对于所有满足$2\le n \le 100000$的非完全平方自然数$n$,找出其分母上限为$10^12$的最优有理逼近数,并求出所有这些最优有理逼近数的分母之和。

本题难度:



解答

本题乍看和第64题很相似,似乎只要求出$\sqrt n$的连分数表示就可以找出最优有理逼近数,但用题面中的数据尝试展开$\sqrt{13}$并观察分母可以发现5的下一项是33,因此第64题中的程序不能直接套用。

以下的解法使用法里序列的标准算法,这也是分母带上限的最优有理逼近问题的标准解法,最终结果是$57060635927998347$。

sq=set([i*i for i in range(2,317)])
bound=10**12
r=0
for n in range(2,100001):
    if n not in sq:
        a,b,c,d=0,1,1,0
        while b+d <= bound:
            if (a+c)*(a+c)-n*(b+d)*(b+d)>0:
                c,d=a+c,b+d
            else:
                a,b=a+c,b+d
        r+=b if 4*n*b*b*d*d < (a*d+b*c)*(a*d+b*c) else d
            
print r