# HG changeset patch
# User Frederic Chapoton <chapoton at math.univlyon1.fr>
# Date 1370286082 7200
# Node ID ee5a28955de5246c1e03241e1ef59e6baaf4722f
# Parent 05cbe12e7bfd588fc1c746176dbea93cf4a4bd5c
trac 14567 minor details, review patch
diff git a/sage/rings/continued_fractions.py b/sage/rings/continued_fractions.py
a

b

AUTHORS: 
187  187   Vincent Delecroix (2013): cleaning, refactorisation, documentation (:trac:`14567`) 
188  188  """ 
189  189  from sage.structure.unique_representation import UniqueRepresentation 
190   from sage.structure.element import Element 
191   from sage.structure.element import FieldElement 
 190  from sage.structure.element import Element, FieldElement 
192  191  
193  192  from field import Field 
194  193  from integer import Integer 
195   from rational import Rational 
196  194  from infinity import Infinity 
197  195  from integer_ring import ZZ 
198  196  from rational_field import QQ 
… 
… 
ZZ_0 = Integer(0) 
203  201  ZZ_1 = Integer(1) 
204  202  ZZ_2 = Integer(2) 
205  203  
 204  
206  205  def two_last_convergents(x): 
207  206  """ 
208  207  Given the list ``x`` that consists of numbers, return the two last 
… 
… 
def two_last_convergents(x): 
221  220  sage: two_last_convergents([1,1,3,2]) 
222  221  (1, 4, 2, 9) 
223  222  """ 
224   p0,p1 = 0,1 
225   q0,q1 = 1,0 
 223  p0, p1 = 0, 1 
 224  q0, q1 = 1, 0 
226  225  for a in x: 
227   p0,p1 = p1,a*p1+p0 
228   q0,q1 = q1,a*q1+q0 
229   return p0,q0,p1,q1 
 226  p0, p1 = p1, a*p1+p0 
 227  q0, q1 = q1, a*q1+q0 
 228  return p0, q0, p1, q1 
 229  
230  230  
231  231  class ContinuedFraction_generic(FieldElement): 
232  232  r""" 
… 
… 
class ContinuedFraction_generic(FieldEle 
301  301  b = other.quotient(i) 
302  302  test = cmp(a,b) 
303  303  if test == 1: # a > b 
304   return 1 if i%2 else 1 
305   if test == 1: # b > a 
306   return 1 if i%2 else 1 
307   if a == 0 and b == 0 and i: # rational case 
 304  return 1 if i % 2 else 1 
 305  if test == 1: # b > a 
 306  return 1 if i % 2 else 1 
 307  if a == 0 and b == 0 and i: # rational case 
308  308  return 0 
309  309  i += 1 
310  310  
… 
… 
class ContinuedFraction_generic(FieldEle 
405  405  return (R(m1) << e) 
406  406  
407  407  # 3. positive non integer 
408   if self.quotient(0) == 0: # 0 <= self < 1 
 408  if self.quotient(0) == 0: # 0 <= self < 1 
409  409  N = R.prec() + self.quotient(1).nbits()  1 
410   if self.quotient(2) == 0 and self.quotient(1)%(1<<(self.quotient(1).nbits()1)) == 0: 
 410  if self.quotient(2) == 0 and self.quotient(1) % (1 << (self.quotient(1).nbits()1)) == 0: 
411  411  # if self is of the form [0; 2^N] then we need the following 
412  412  N = 1 
413   else: # self > 1 
 413  else: # self > 1 
414  414  N = R.prec()  self.quotient(0).nbits() 
415  415  
416  416  # even/odd convergents are respectively below/above 
417  417  k = 0 
418   p_even = self.pn(2*k); p_odd = self.pn(2*k+1) 
419   q_even = self.qn(2*k); q_odd = self.qn(2*k+1) 
420   m_even = (p_even<<N) // q_even # floor((2^N p_even) / q_even) 
421   m_odd = (p_odd<<N + q_odd  1) // q_odd # ceil((2^N p_odd) / q_odd) 
 418  p_even = self.pn(2*k) 
 419  p_odd = self.pn(2*k+1) 
 420  q_even = self.qn(2*k) 
 421  q_odd = self.qn(2*k+1) 
 422  m_even = (p_even << N) // q_even # floor((2^N p_even) / q_even) 
 423  m_odd = (p_odd << N + q_odd  1) // q_odd # ceil((2^N p_odd) / q_odd) 
422  424  while (m_odd  m_even) > 1: 
423  425  k += 1 
424   p_even = self.pn(2*k); p_odd = self.pn(2*k+1) 
425   q_even = self.qn(2*k); q_odd = self.qn(2*k+1) 
426   m_even = (p_even<<N) // q_even 
427   m_odd = ((p_odd<<N) + q_odd  1) // q_odd 
 426  p_even = self.pn(2*k) 
 427  p_odd = self.pn(2*k+1) 
 428  q_even = self.qn(2*k) 
 429  q_odd = self.qn(2*k+1) 
 430  m_even = (p_even << N) // q_even 
 431  m_odd = ((p_odd << N) + q_odd  1) // q_odd 
428  432  
429  433  assert m_odd.nbits() == R.prec() or m_even.nbits() == R.prec() 
430  434  
431  435  if m_even == m_odd: # no need to worry (we have a decimal number) 
432   return R(m_even)>>N 
 436  return R(m_even) >> N 
433  437  
434  438  # check ordering 
435  439  # m_even/2^N <= p_even/q_even <= self <= p_odd/q_odd <= m_odd/2^N 
436  440  assert m_odd == m_even + 1 
437   assert m_even / (ZZ_1<<N) <= p_even/q_even 
 441  assert m_even / (ZZ_1 << N) <= p_even/q_even 
438  442  assert p_even / q_even <= p_odd / q_odd 
439   assert p_odd / q_odd <= m_odd / (ZZ_1<<N) 
 443  assert p_odd / q_odd <= m_odd / (ZZ_1 << N) 
440  444  
441  445  rnd = R.rounding_mode() 
442   if rnd == 'RNDN': # round to the nearest 
 446  if rnd == 'RNDN': # round to the nearest 
443  447  # in order to find the nearest approximation we possibly need to 
444  448  # augment our precision on convergents. 
445  449  while True: 
446   assert not(p_odd<<(N+1) <= (2*m_odd1) * q_odd) or not(p_even<<(N+1) >= (2*m_even+1) * q_even) 
447   if p_odd<<(N+1) <= (2*m_odd1) * q_odd: 
 450  assert not(p_odd << (N+1) <= (2*m_odd1) * q_odd) or not(p_even << (N+1) >= (2*m_even+1) * q_even) 
 451  if p_odd << (N+1) <= (2*m_odd1) * q_odd: 
448  452  return R(m_even) >> N 
449   if p_even<<(N+1) >= (2*m_even+1) * q_even: 
450   return R(m_odd)>>N 
 453  if p_even << (N+1) >= (2*m_even+1) * q_even: 
 454  return R(m_odd) >> N 
451  455  k += 1 
452   p_even = self.pn(2*k); p_odd = self.pn(2*k+1) 
453   q_even = self.qn(2*k); q_odd = self.qn(2*k+1) 
 456  p_even = self.pn(2*k) 
 457  p_odd = self.pn(2*k+1) 
 458  q_even = self.qn(2*k) 
 459  q_odd = self.qn(2*k+1) 
454  460  elif rnd == 'RNDU': # round up (toward +infinity) 
455   return R(m_odd)>>N 
 461  return R(m_odd) >> N 
456  462  elif rnd == 'RNDD': # round down (toward infinity) 
457   return R(m_even)>>N 
 463  return R(m_even) >> N 
458  464  elif rnd == 'RNDZ': # round toward zero 
459  465  if m_even.sign() == 1: 
460   return R(m_even)>>N 
 466  return R(m_even) >> N 
461  467  else: 
462   return R(m_odd)>>N 
 468  return R(m_odd) >> N 
463  469  else: 
464   raise ValueError("%s unknown rounding mode"%rounding) 
 470  raise ValueError("%s unknown rounding mode" % rnd) 
465  471  
466  472  def __float__(self): 
467  473  """ 
… 
… 
class ContinuedFraction_generic(FieldEle 
601  607  else: 
602  608  try: 
603  609  stop = n.stop.__index__() 
604   except (AttributeError,ValueError): 
 610  except (AttributeError, ValueError): 
605  611  raise ValueError("stop should be an integer") 
606  612  return self.parent()([self.quotient(i) for i in xrange(start,stop,step)]) 
607  613  start, stop, step = n.indices(self.length()) 
608  614  return self.parent()([self.quotient(i) for i in xrange(start,stop,step)]) 
609  615  try: 
610  616  n = n.__index__() 
611   except (AttributeError,ValueError): 
612   raise ValueError("n (=%s) should be an integer"%n) 
 617  except (AttributeError, ValueError): 
 618  raise ValueError("n (=%s) should be an integer" % n) 
613  619  if n < 0: 
614   raise ValueError("n (=%s) should be positive"%n) 
 620  raise ValueError("n (=%s) should be positive" % n) 
615  621  q = self.quotient(n) 
616  622  if n > 0 and q == 0: 
617  623  raise IndexError("index out of range") 
… 
… 
class ContinuedFraction_generic(FieldEle 
887  893  return ZZ_2 
888  894  return Infinity 
889  895  
 896  
890  897  class ContinuedFraction_periodic(ContinuedFraction_generic): 
891  898  r""" 
892  899  Continued fraction associated with rational or quadratic number. 
… 
… 
class ContinuedFraction_periodic(Continu 
980  987  """ 
981  988  n = int(n) 
982  989  if n < 0: 
983   raise ValueError("n (=%d) should be positive"%n) 
 990  raise ValueError("n (=%d) should be positive" % n) 
984  991  if n < len(self._x1): 
985  992  return self._x1[n] 
986   return self._x2[(nlen(self._x1))%len(self._x2)] 
 993  return self._x2[(nlen(self._x1)) % len(self._x2)] 
987  994  
988  995  def length(self): 
989  996  r""" 
… 
… 
class ContinuedFraction_periodic(Continu 
1033  1040  b = other.quotient(i) 
1034  1041  test = cmp(a,b) 
1035  1042  if test == 1: 
1036   return 1 if i%2 else 1 
 1043  return 1 if i % 2 else 1 
1037  1044  if test == 1: 
1038   return 1 if i%2 else 1 
 1045  return 1 if i % 2 else 1 
1039  1046  return 0 
1040  1047  
1041  1048  return ContinuedFraction_generic.__cmp__(self, other) 
… 
… 
class ContinuedFraction_periodic(Continu 
1120  1127  from sage.misc.functional import squarefree_part 
1121  1128  D = (q0p1)**2 + 4*q1*p0 
1122  1129  DD = squarefree_part(D) 
1123   Q = QuadraticField(DD, 'sqrt%d'%DD) 
 1130  Q = QuadraticField(DD, 'sqrt%d' % DD) 
1124  1131  x = ((p1  q0) + (D/DD).sqrt() * Q.gen()) / (2*q1) 
1125  1132  
1126  1133  # we add the preperiod 
… 
… 
class ContinuedFraction_periodic(Continu 
1142  1149  sage: CFF([(),(1,3)]) 
1143  1150  [(1, 3)*] 
1144  1151  """ 
1145   if len(self._x2) == 1 and not self._x2[0]: # rational 
 1152  if len(self._x2) == 1 and not self._x2[0]: # rational 
1146  1153  if len(self._x1) == 1: 
1147   return '[%d]'%self._x1[0] 
1148   return '[%d; '%self._x1[0] + ', '.join(str(a) for a in self._x1[1:]) + ']' 
 1154  return '[%d]' % self._x1[0] 
 1155  return '[%d; ' % self._x1[0] + ', '.join(str(a) for a in self._x1[1:]) + ']' 
1149  1156  
1150  1157  period = '(' + ', '.join(str(a) for a in self._x2) + ')*' 
1151   if not self._x1: # purely periodic case 
 1158  if not self._x1: # purely periodic case 
1152  1159  return '[' + period + ']' 
1153  1160  
1154  1161  if len(self._x1) == 1: 
1155   return '[%d; '%self._x1[0] + period + ']' 
1156   return '[%d; '%self._x1[0] + ', '.join(str(a) for a in self._x1[1:]) + ', ' + period + ']' 
 1162  return '[%d; ' % self._x1[0] + period + ']' 
 1163  return '[%d; ' % self._x1[0] + ', '.join(str(a) for a in self._x1[1:]) + ', ' + period + ']' 
1157  1164  
1158  1165  def __reduce__(self): 
1159  1166  r""" 
… 
… 
class ContinuedFraction_periodic(Continu 
1228  1235  return '0' 
1229  1236  s = str(v[0]) 
1230  1237  for i in range(1,len(v)): 
1231   s += '+ \\frac{\\displaystyle 1}{\\displaystyle %s'%v[i] 
 1238  s += '+ \\frac{\\displaystyle 1}{\\displaystyle %s' % v[i] 
1232  1239  s += '}'*(len(v)1) 
1233  1240  return s 
1234  1241  
… 
… 
class ContinuedFraction_periodic(Continu 
1338  1345  return self.parent()([(x1[0]1, x1[2]+1) + x1[3:], self._x2]) 
1339  1346  return self.parent()([(x1[0]1, ZZ_1, x1[1]1) + x1[2:], self._x2]) 
1340  1347  
 1348  
1341  1349  class ContinuedFraction_irrational(ContinuedFraction_generic): 
1342  1350  r""" 
1343  1351  Continued fraction of an exact irrational number. 
… 
… 
class ContinuedFraction_irrational(Conti 
1438  1446  sage: continued_fraction(pi) # indirect doctest 
1439  1447  [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...] 
1440  1448  """ 
1441   return '[%d; '%self.quotient(0) + ', '.join(str(self.quotient(i)) for i in xrange(1,20)) + ", ...]" 
 1449  return '[%d; ' % self.quotient(0) + ', '.join(str(self.quotient(i)) for i in xrange(1,20)) + ", ...]" 
1442  1450  
1443  1451  def quotient(self, n): 
1444  1452  r""" 
… 
… 
class ContinuedFraction_irrational(Conti 
1552  1560  """ 
1553  1561  return self._x0 
1554  1562  
 1563  
1555  1564  def check_and_reduce_pair(x1,x2=None): 
1556  1565  r""" 
1557  1566  There are often two ways to represent a given continued fraction. This 
… 
… 
def check_and_reduce_pair(x1,x2=None): 
1609  1618  # check that y2 is not a pure power (in a very naive way!!) 
1610  1619  n2 = len(y2) 
1611  1620  for i in xrange(1,(n2+2)/2): 
1612   if n2%i == 0 and y2[:i] == y2[i:]: 
 1621  if n2 % i == 0 and y2[:i] == y2[i:]: 
1613  1622  y2 = y2[:i] 
1614  1623  break 
1615  1624  
… 
… 
def check_and_reduce_pair(x1,x2=None): 
1620  1629  
1621  1630  return tuple(y1),tuple(y2) 
1622  1631  
 1632  
1623  1633  class ContinuedFractionField(UniqueRepresentation, Field): 
1624  1634  """ 
1625  1635  A common parent for continued fractions. 
… 
… 
class ContinuedFractionField(UniqueRepre 
1736  1746  if isinstance(x, Element) and x.parent() is self: 
1737  1747  return x 
1738  1748  
1739   if isinstance(x, (list,tuple)): # a digit expansion 
1740   cls = self._element_class_periodic 
1741   
 1749  if isinstance(x, (list, tuple)): # a digit expansion 
1742  1750  if len(x) == 2 and isinstance(x[0], (list,tuple)) and isinstance(x[1], (list,tuple)): 
1743  1751  x1 = tuple(ZZ(a) for a in x[0]) 
1744  1752  x2 = tuple(ZZ(a) for a in x[1]) 
… 
… 
class ContinuedFractionField(UniqueRepre 
1747  1755  x1, x2 = check_and_reduce_pair(x) 
1748  1756  return self._element_class_periodic(self, x1, x2) 
1749  1757  
1750   else: # a number 
1751   if x in QQ: # rational 
 1758  else: # a number 
 1759  if x in QQ: # rational 
1752  1760  return QQ(x).continued_fraction() 
1753   if bits is None and nterms is None and isinstance(x, Element): # an irrational 
 1761  if bits is None and nterms is None and isinstance(x, Element): # an irrational 
1754  1762  from sage.symbolic.ring import SR 
1755  1763  if x.parent() is SR: 
1756  1764  # we try a conversion to the algebraic real field (in order 
… 
… 
class ContinuedFractionField(UniqueRepre 
1779  1787  try: 
1780  1788  RIF(x) 
1781  1789  except StandardError: 
1782   raise ValueError("the number %s does not seem to be a real number"%x) 
 1790  raise ValueError("the number %s does not seem to be a real number" % x) 
1783  1791  return self._element_class_number(self, x) 
1784  1792  
1785  1793  # now work with RIF 
… 
… 
class ContinuedFractionField(UniqueRepre 
1901  1909  if preperiod: 
1902  1910  if random() > 0.1: 
1903  1911  x1.append(ZZ.random_element(x,y,distribution)) 
1904   while random() > 0.5: # from here they may be arbitrarily many 
1905   # elements but the mean length is 2 
 1912  while random() > 0.5: # from here they may be arbitrarily many 
 1913  # elements but the mean length is 2 
1906  1914  x1.append(ZZ.random_element(x,y,distribution)) 
1907  1915  if period: 
1908  1916  x2.append(ZZ.random_element(x,y,distribution)) 
… 
… 
class ContinuedFractionField(UniqueRepre 
1914  1922  
1915  1923  CFF = ContinuedFractionField() 
1916  1924  
 1925  
1917  1926  def continued_fraction_list(x, type="std", partial_convergents=False, bits=None, nterms=None): 
1918  1927  r""" 
1919  1928  Returns the (finite) continued fraction of ``x`` as a list. 
… 
… 
def continued_fraction_list(x, type="std 
2066  2075  limit = cf.length() 
2067  2076  else: 
2068  2077  import warnings 
2069   warnings.warn("the continued fraction of %s seems infinite, return only the first 20 terms"%x) 
 2078  warnings.warn("the continued fraction of %s seems infinite, return only the first 20 terms" % x) 
2070  2079  limit = 20 
2071  2080  if partial_convergents: 
2072  2081  return [cf.quotient(i) for i in xrange(limit)], [(cf.pn(i),cf.qn(i)) for i in xrange(limit)] 
2073  2082  return [cf.quotient(i) for i in xrange(limit)] 
2074  2083  
 2084  
2075  2085  def continued_fraction(x, bits=None, nterms=None): 
2076  2086  """ 
2077  2087  Return the continued fraction of ``x``. 
… 
… 
def continued_fraction(x, bits=None, nte 
2179  2189  """ 
2180  2190  return CFF(x, bits=bits, nterms=nterms) 
2181  2191  
 2192  
2182  2193  def Hirzebruch_Jung_continued_fraction_list(x, bits=None, nterms=None): 
2183  2194  r""" 
2184  2195  Return the HirzebruchJung continued fraction of ``x`` as a list. 
… 
… 
def Hirzebruch_Jung_continued_fraction_l 
2221  2232  # sage.rings.continued_fractions in #14567 
2222  2233  
2223  2234  from sage.structure.sage_object import register_unpickle_override 
2224   register_unpickle_override('sage.rings.contfrac','ContinuedFractionField_class',ContinuedFractionField) 
 2235  register_unpickle_override('sage.rings.contfrac', 'ContinuedFractionField_class',ContinuedFractionField) 