Ticket #14567: trac_14567_addon1-fc.patch

File trac_14567_addon1-fc.patch, 15.8 KB (added by chapoton, 7 years ago)
  • sage/rings/continued_fractions.py

    # HG changeset patch
    # User Frederic Chapoton <chapoton at math.univ-lyon1.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: 
    187187- Vincent Delecroix (2013): cleaning, refactorisation, documentation (:trac:`14567`)
    188188"""
    189189from sage.structure.unique_representation import UniqueRepresentation
    190 from sage.structure.element import Element
    191 from sage.structure.element import FieldElement
     190from sage.structure.element import Element, FieldElement
    192191
    193192from field import Field
    194193from integer import Integer
    195 from rational import Rational
    196194from infinity import Infinity
    197195from integer_ring import ZZ
    198196from rational_field import QQ
    ZZ_0 = Integer(0) 
    203201ZZ_1 = Integer(1)
    204202ZZ_2 = Integer(2)
    205203
     204
    206205def two_last_convergents(x):
    207206    """
    208207    Given the list ``x`` that consists of numbers, return the two last
    def two_last_convergents(x): 
    221220        sage: two_last_convergents([-1,1,3,2])
    222221        (-1, 4, -2, 9)
    223222    """
    224     p0,p1 = 0,1
    225     q0,q1 = 1,0
     223    p0, p1 = 0, 1
     224    q0, q1 = 1, 0
    226225    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
    230230
    231231class ContinuedFraction_generic(FieldElement):
    232232    r"""
    class ContinuedFraction_generic(FieldEle 
    301301            b = other.quotient(i)
    302302            test = cmp(a,b)
    303303            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
    308308                return 0
    309309            i += 1
    310310
    class ContinuedFraction_generic(FieldEle 
    405405                return -(R(m-1) << e)
    406406
    407407        # 3. positive non integer
    408         if self.quotient(0) == 0: # 0 <= self < 1
     408        if self.quotient(0) == 0:  # 0 <= self < 1
    409409            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:
    411411                # if self is of the form [0; 2^N] then we need the following
    412412                N -= 1
    413         else: # self > 1
     413        else:  # self > 1
    414414            N = R.prec() - self.quotient(0).nbits()
    415415
    416416        # even/odd convergents are respectively below/above
    417417        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)
    422424        while (m_odd - m_even) > 1:
    423425            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
    428432
    429433        assert m_odd.nbits() == R.prec() or m_even.nbits() == R.prec()
    430434
    431435        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
    433437
    434438        # check ordering
    435439        # m_even/2^N <= p_even/q_even <= self <= p_odd/q_odd <= m_odd/2^N
    436440        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
    438442        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)
    440444
    441445        rnd = R.rounding_mode()
    442         if rnd == 'RNDN': # round to the nearest
     446        if rnd == 'RNDN':  # round to the nearest
    443447            # in order to find the nearest approximation we possibly need to
    444448            # augment our precision on convergents.
    445449            while True:
    446                 assert not(p_odd<<(N+1) <= (2*m_odd-1) * q_odd) or not(p_even<<(N+1) >= (2*m_even+1) * q_even)
    447                 if p_odd<<(N+1) <= (2*m_odd-1) * q_odd:
     450                assert not(p_odd << (N+1) <= (2*m_odd-1) * q_odd) or not(p_even << (N+1) >= (2*m_even+1) * q_even)
     451                if p_odd << (N+1) <= (2*m_odd-1) * q_odd:
    448452                    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
    451455                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)
    454460        elif rnd == 'RNDU':  # round up (toward +infinity)
    455             return R(m_odd)>>N
     461            return R(m_odd) >> N
    456462        elif rnd == 'RNDD':  # round down (toward -infinity)
    457             return R(m_even)>>N
     463            return R(m_even) >> N
    458464        elif rnd == 'RNDZ':  # round toward zero
    459465            if m_even.sign() == 1:
    460                 return R(m_even)>>N
     466                return R(m_even) >> N
    461467            else:
    462                 return R(m_odd)>>N
     468                return R(m_odd) >> N
    463469        else:
    464             raise ValueError("%s unknown rounding mode"%rounding)
     470            raise ValueError("%s unknown rounding mode" % rnd)
    465471
    466472    def __float__(self):
    467473        """
    class ContinuedFraction_generic(FieldEle 
    601607                else:
    602608                    try:
    603609                        stop = n.stop.__index__()
    604                     except (AttributeError,ValueError):
     610                    except (AttributeError, ValueError):
    605611                        raise ValueError("stop should be an integer")
    606612                return self.parent()([self.quotient(i) for i in xrange(start,stop,step)])
    607613            start, stop, step = n.indices(self.length())
    608614            return self.parent()([self.quotient(i) for i in xrange(start,stop,step)])
    609615        try:
    610616            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)
    613619        if n < 0:
    614             raise ValueError("n (=%s) should be positive"%n)
     620            raise ValueError("n (=%s) should be positive" % n)
    615621        q = self.quotient(n)
    616622        if n > 0 and q == 0:
    617623            raise IndexError("index out of range")
    class ContinuedFraction_generic(FieldEle 
    887893            return ZZ_2
    888894        return Infinity
    889895
     896
    890897class ContinuedFraction_periodic(ContinuedFraction_generic):
    891898    r"""
    892899    Continued fraction associated with rational or quadratic number.
    class ContinuedFraction_periodic(Continu 
    980987        """
    981988        n = int(n)
    982989        if n < 0:
    983             raise ValueError("n (=%d) should be positive"%n)
     990            raise ValueError("n (=%d) should be positive" % n)
    984991        if n < len(self._x1):
    985992            return self._x1[n]
    986         return self._x2[(n-len(self._x1))%len(self._x2)]
     993        return self._x2[(n-len(self._x1)) % len(self._x2)]
    987994
    988995    def length(self):
    989996        r"""
    class ContinuedFraction_periodic(Continu 
    10331040                b = other.quotient(i)
    10341041                test = cmp(a,b)
    10351042                if test == 1:
    1036                     return -1 if i%2 else 1
     1043                    return -1 if i % 2 else 1
    10371044                if test == -1:
    1038                     return 1 if i%2 else -1
     1045                    return 1 if i % 2 else -1
    10391046            return 0
    10401047
    10411048        return ContinuedFraction_generic.__cmp__(self, other)
    class ContinuedFraction_periodic(Continu 
    11201127        from sage.misc.functional import squarefree_part
    11211128        D = (q0-p1)**2 + 4*q1*p0
    11221129        DD = squarefree_part(D)
    1123         Q = QuadraticField(DD, 'sqrt%d'%DD)
     1130        Q = QuadraticField(DD, 'sqrt%d' % DD)
    11241131        x = ((p1 - q0) + (D/DD).sqrt() * Q.gen()) / (2*q1)
    11251132
    11261133        # we add the preperiod
    class ContinuedFraction_periodic(Continu 
    11421149            sage: CFF([(),(1,3)])
    11431150            [(1, 3)*]
    11441151        """
    1145         if len(self._x2) == 1 and not self._x2[0]: # rational
     1152        if len(self._x2) == 1 and not self._x2[0]:  # rational
    11461153            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:]) + ']'
    11491156
    11501157        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
    11521159            return '[' + period + ']'
    11531160
    11541161        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 + ']'
    11571164
    11581165    def __reduce__(self):
    11591166        r"""
    class ContinuedFraction_periodic(Continu 
    12281235            return '0'
    12291236        s = str(v[0])
    12301237        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]
    12321239        s += '}'*(len(v)-1)
    12331240        return s
    12341241
    class ContinuedFraction_periodic(Continu 
    13381345            return self.parent()([(-x1[0]-1, x1[2]+1) + x1[3:], self._x2])
    13391346        return self.parent()([(-x1[0]-1, ZZ_1, x1[1]-1) + x1[2:], self._x2])
    13401347
     1348
    13411349class ContinuedFraction_irrational(ContinuedFraction_generic):
    13421350    r"""
    13431351    Continued fraction of an exact irrational number.
    class ContinuedFraction_irrational(Conti 
    14381446            sage: continued_fraction(pi) # indirect doctest
    14391447            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
    14401448        """
    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)) + ", ...]"
    14421450
    14431451    def quotient(self, n):
    14441452        r"""
    class ContinuedFraction_irrational(Conti 
    15521560        """
    15531561        return self._x0
    15541562
     1563
    15551564def check_and_reduce_pair(x1,x2=None):
    15561565    r"""
    15571566    There are often two ways to represent a given continued fraction. This
    def check_and_reduce_pair(x1,x2=None): 
    16091618    # check that y2 is not a pure power (in a very naive way!!)
    16101619    n2 = len(y2)
    16111620    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:]:
    16131622            y2 = y2[:i]
    16141623            break
    16151624
    def check_and_reduce_pair(x1,x2=None): 
    16201629
    16211630    return tuple(y1),tuple(y2)
    16221631
     1632
    16231633class ContinuedFractionField(UniqueRepresentation, Field):
    16241634    """
    16251635    A common parent for continued fractions.
    class ContinuedFractionField(UniqueRepre 
    17361746        if isinstance(x, Element) and x.parent() is self:
    17371747            return x
    17381748
    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
    17421750            if len(x) == 2 and isinstance(x[0], (list,tuple)) and isinstance(x[1], (list,tuple)):
    17431751                x1 = tuple(ZZ(a) for a in x[0])
    17441752                x2 = tuple(ZZ(a) for a in x[1])
    class ContinuedFractionField(UniqueRepre 
    17471755                x1, x2 = check_and_reduce_pair(x)
    17481756            return self._element_class_periodic(self, x1, x2)
    17491757
    1750         else: # a number
    1751             if x in QQ: # rational
     1758        else:  # a number
     1759            if x in QQ:  # rational
    17521760                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
    17541762                from sage.symbolic.ring import SR
    17551763                if x.parent() is SR:
    17561764                    # we try a conversion to the algebraic real field (in order
    class ContinuedFractionField(UniqueRepre 
    17791787                    try:
    17801788                        RIF(x)
    17811789                    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)
    17831791                    return self._element_class_number(self, x)
    17841792
    17851793            # now work with RIF
    class ContinuedFractionField(UniqueRepre 
    19011909        if preperiod:
    19021910            if random() > 0.1:
    19031911                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
    19061914                    x1.append(ZZ.random_element(x,y,distribution))
    19071915        if period:
    19081916            x2.append(ZZ.random_element(x,y,distribution))
    class ContinuedFractionField(UniqueRepre 
    19141922
    19151923CFF = ContinuedFractionField()
    19161924
     1925
    19171926def continued_fraction_list(x, type="std", partial_convergents=False, bits=None, nterms=None):
    19181927    r"""
    19191928    Returns the (finite) continued fraction of ``x`` as a list.
    def continued_fraction_list(x, type="std 
    20662075        limit = cf.length()
    20672076    else:
    20682077        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)
    20702079        limit = 20
    20712080    if partial_convergents:
    20722081        return [cf.quotient(i) for i in xrange(limit)], [(cf.pn(i),cf.qn(i)) for i in xrange(limit)]
    20732082    return [cf.quotient(i) for i in xrange(limit)]
    20742083
     2084
    20752085def continued_fraction(x, bits=None, nterms=None):
    20762086    """
    20772087    Return the continued fraction of ``x``.
    def continued_fraction(x, bits=None, nte 
    21792189    """
    21802190    return CFF(x, bits=bits, nterms=nterms)
    21812191
     2192
    21822193def Hirzebruch_Jung_continued_fraction_list(x, bits=None, nterms=None):
    21832194    r"""
    21842195    Return the Hirzebruch-Jung continued fraction of ``x`` as a list.
    def Hirzebruch_Jung_continued_fraction_l 
    22212232# sage.rings.continued_fractions in #14567
    22222233
    22232234from sage.structure.sage_object import register_unpickle_override
    2224 register_unpickle_override('sage.rings.contfrac','ContinuedFractionField_class',ContinuedFractionField)
     2235register_unpickle_override('sage.rings.contfrac', 'ContinuedFractionField_class',ContinuedFractionField)