fork download
copy 
  1. import math
  2.  
  3. def expand(n, numerator, denominator):
  4. # 1回連分数展開して整数部分を返します．
  5. # numerator / (sqrt(n) + denominator)
  6. # = 1 / (sqrt(23) - 4)
  7. # = 1 + 1 / (7 / (sqrt(23) - 3))
  8. # = ret + 1 / (numerator / (sqrt(N) + denominator))
  9. # ret, numerator, denominator : 返り値
  10.  
  11.     ret = 1
  12.     _numerator = (n - math.pow(denominator,2)) // numerator
  13.     _denominator = -1 * denominator
  14.     while True:
  15.         _denominator = _denominator - _numerator
  16.         if n < pow(_numerator - _denominator, 2):
  17.             return ret, _numerator, _denominator
  18.         ret = ret + 1
  19.  
  20. def circulation_list(D):
  21. # [c0;c1,c2,...,ck]を返します．(k:周期)
  22.  
  23.     list = []
  24.     list.append(int(math.floor(math.sqrt(D))))
  25.     numerator = first_numerator = 1
  26.     denominator = first_denominator = -1 * list[0]
  27.     while True:
  28.         c, numerator, denominator = expand(D, numerator, denominator)
  29.         list.append(c)
  30.         if numerator == first_numerator and denominator == first_denominator:
  31.             return list
  32.  
  33. def is_square_num(n):
  34.     if n == 0:
  35.         return 1
  36.  
  37.     i = 1
  38.     while i * i <= n:
  39.         if i * i == n:
  40.             return True
  41.         i = i + 1
  42.  
  43.     return False
  44.  
  45. max = 0
  46. maxD = 2
  47. for D in range(2, 1001):
  48.     if is_square_num(D):
  49.         continue
  50.  
  51.     list = circulation_list(D)
  52.     period = len(list) - 1
  53.     if period % 2 == 0:
  54.         m = 1
  55.     else:
  56.         m = 2
  57.         list = list[:] + list[1:]
  58.     # p,qと添字をそろえる
  59.     list = [0] + list[:]
  60.  
  61.     p = [1, list[1]]
  62.     q = [0, 1]
  63.     for i in range(2, m * period + 1):
  64.         p.append(list[i] * p[i - 1] + p[i - 2])
  65.         q.append(list[i] * q[i - 1] + q[i - 2])
  66.  
  67.     #print('D = {0} : (x, y) = ({1}, {2})'.format(D, p[-1], q[-1]))
  68.     if max < p[-1]:
  69.         max = p[-1]
  70.         maxD = D
  71.  
  72. print(maxD)# your code goes here
Success #stdin #stdout 0.05s 27720KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
661
https://ideone.com/WF6K5w
language:
Python 3 (python  3.9.5)
created:
7 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!