Ticket #14567: trac_14567-continued_fractions.patch

File trac_14567-continued_fractions.patch, 127.5 KB (added by vdelecroix, 6 years ago)
  • new file doc/en/reference/diophantine_approximation/conf.py

    # HG changeset patch
    # User Vincent Delecroix <20100.delecroix at gmail.com>
    # Date 1369482941 -7200
    # Node ID 3c6fae0bbd8ba5f0aaab0b0be572f40fa7fa9328
    # Parent  7c97a70b162921d91dca8c6c9cd792ea98dcfaec
    ticket 14567
    
    many improvements for continued fractions
     * properly initialized category
     * infinite continued fractions to deal with real numbers
     * ultimately periodic continued fractions and quadratic numbers
     * clean method for numerical approximation
     * link with sage.combinat.words
    
    diff --git a/doc/en/reference/diophantine_approximation/conf.py b/doc/en/reference/diophantine_approximation/conf.py
    new file mode 120000
    - +  
     1../conf_sub.py
     2 No newline at end of file
  • new file doc/en/reference/diophantine_approximation/index.rst

    diff --git a/doc/en/reference/diophantine_approximation/index.rst b/doc/en/reference/diophantine_approximation/index.rst
    new file mode 100644
    - +  
     1Diophantine approximation
     2=========================
     3
     4The diophantine approximation deals with the approximation of real numbers
     5(or real vectors) with rational numbers (or rational vectors).
     6See the article :wikipedia:`Diophantine_approximation` for more information.
     7
     8.. toctree::
     9   :maxdepth: 2
     10
     11   sage/rings/continued_fractions
     12
     13.. include:: ../footer.txt
  • doc/en/reference/index.rst

    diff --git a/doc/en/reference/index.rst b/doc/en/reference/index.rst
    a b Geometry and Topology 
    9595Number Theory, Algebraic Geometry
    9696---------------------------------
    9797
     98* :doc:`Diophantine approximation <diophantine_approximation/index>`
    9899* :doc:`Quadratic Forms <quadratic_forms/index>`
    99100* :doc:`L-Functions <lfunctions/index>`
    100101* :doc:`Schemes <schemes/index>`
  • sage/calculus/wester.py

    diff --git a/sage/calculus/wester.py b/sage/calculus/wester.py
    a b explicit calls to Maxima or other system 
    5353    sage: # (YES) Continued fraction of 3.1415926535
    5454    sage: a = 3.1415926535
    5555    sage: continued_fraction(a)
    56     [3, 7, 15, 1, 292, 1, 1, 6, 2, 13, 4]
     56    [3; 7, 15, 1, 292, 1, 1, 6, 2, 13, 4]
    5757
    5858::
    5959
  • sage/combinat/words/word_generators.py

    diff --git a/sage/combinat/words/word_generators.py b/sage/combinat/words/word_generators.py
    a b class LowerChristoffelWord(FiniteWord_li 
    198198            elif (p, q) == (1, 0):
    199199                w = [alphabet[1]]
    200200            else:
    201                 from sage.rings.all import QQ, CFF
    202                 cf = CFF(QQ((p, q)))
     201                from sage.rings.rational_field import QQ
     202                cf = QQ((p, q)).continued_fraction_list()
    203203                u = [alphabet[0]]
    204204                v = [alphabet[1]]
    205205                #do not consider the first zero if p < q
    class WordGenerator(object): 
    888888                msg = "The argument slope (=%s) must be in ]0,1[."%slope
    889889                raise ValueError, msg
    890890            from sage.rings.all import CFF
    891             cf = iter(CFF(slope, bits=bits))
    892             length = 'finite'
     891            cf = CFF(slope, bits=bits)
     892            if cf.length() == Infinity:
     893                length = Infinity
     894            else:
     895                length = 'finite'
     896            cf = iter(cf)
    893897        elif hasattr(slope, '__iter__'):
    894898            cf = iter(slope)
    895899            length = Infinity
    class WordGenerator(object): 
    927931        EXAMPLES::
    928932
    929933            sage: CFF(1/golden_ratio^2)[:8]
    930             [0, 2, 1, 1, 1, 1, 1, 1]
     934            [0; 2, 1, 1, 1, 1, 2]
    931935            sage: cf = iter(_)
    932936            sage: Word(words._CharacteristicSturmianWord_LetterIterator(cf))
    933             word: 0100101001001010010100100101001001
     937            word: 0100101001001010010100100101001010
    934938
    935939        ::
    936940
    937941            sage: alpha = (sqrt(3)-1)/2
    938942            sage: CFF(alpha)[:10]
    939             [0, 2, 1, 2, 1, 2, 1, 2, 1, 2]
     943            [0; 2, 1, 2, 1, 2, 1, 2, 1, 2]
    940944            sage: cf = iter(_)
    941945            sage: Word(words._CharacteristicSturmianWord_LetterIterator(cf))
    942946            word: 0100100101001001001010010010010100100101...
  • sage/geometry/cone.py

    diff --git a/sage/geometry/cone.py b/sage/geometry/cone.py
    a b def classify_cone_2d(ray0, ray1, check=T 
    12181218        ...       d, k = classify_cone_2d(ray0, ray1, check=True)
    12191219        ...       assert (d,k) == classify_cone_2d(ray1, ray0)
    12201220        ...       if d == 0: continue
    1221         ...       frac = Hirzebruch_Jung_continued_fraction_list(k/d)
     1221        ...       frac = (k/d).continued_fraction_list("hj")
    12221222        ...       if len(frac)>100: continue   # avoid expensive computation
    12231223        ...       hilb = Cone([ray0, ray1]).Hilbert_basis()
    12241224        ...       assert len(hilb) == len(frac) + 1
  • sage/modular/modsym/ambient.py

    diff --git a/sage/modular/modsym/ambient.py b/sage/modular/modsym/ambient.py
    a b class ModularSymbolsAmbient(space.Modula 
    640640        """
    641641        if alpha.is_infinity():
    642642            return self.manin_symbol((i,0,1), check=False)
    643         v, c = arith.continued_fraction_list(alpha._rational_(), partial_convergents=True)
     643        # v, c = arith.continued_fraction_list(alpha._rational_(), partial_convergents=True)
     644        cf = alpha._rational_().continued_fraction()
     645        v = list(cf)
     646        c = [(cf.pn(k),cf.qn(k)) for k in xrange(len(cf))]
    644647        a = self(0)
    645648        zero = rings.ZZ(0)
    646649        one = rings.ZZ(1)
  • sage/modular/modsym/modular_symbols.py

    diff --git a/sage/modular/modsym/modular_symbols.py b/sage/modular/modsym/modular_symbols.py
    a b import sage.modular.cusps as cusps 
    3636import sage.modular.modsym.manin_symbols
    3737from sage.structure.sage_object import SageObject
    3838import sage.structure.formal_sum as formal_sum
    39 import sage.rings.arith as arith
    4039from sage.rings.integer_ring import ZZ
    4140from sage.misc.latex import latex
    4241
    class ModularSymbol(SageObject): 
    333332        k = space.weight()
    334333        v = [(0,1), (1,0)]
    335334        if not alpha.is_infinity():
    336             v += [(x.numerator(), x.denominator()) for x in arith.convergents(alpha._rational_())]
     335            cf = alpha._rational_().continued_fraction()
     336            v.extend((cf.pn(k),cf.qn(k)) for k in xrange(len(cf)))
    337337        sign = 1
    338338        apply = sage.modular.modsym.manin_symbols.apply_to_monomial
    339339        mansym = sage.modular.modsym.manin_symbols.ManinSymbol
  • sage/plot/misc.py

    diff --git a/sage/plot/misc.py b/sage/plot/misc.py
    a b def _multiple_of_constant(n,pos,const): 
    233233        sage: plot(x^2, (x,0,10), ticks=[sqrt(2),8], tick_formatter=sqrt(2))
    234234    """
    235235    from sage.misc.latex import latex
    236     from sage.rings.arith import convergents
    237     c=[i for i in convergents(n/const.n()) if i.denominator()<12]
    238     return '$%s$'%latex(c[-1]*const)
     236    from sage.rings.continued_fractions import continued_fraction
     237    cf = continued_fraction(n/const)
     238    k = 1
     239    while cf.quotient(k) and cf.qn(k) < 12:
     240        k += 1
     241    return '$%s$'%latex(cf.convergent(k-1)*const)
    239242
    240243
    241244def get_matplotlib_linestyle(linestyle, return_type):
  • sage/rings/all.py

    diff --git a/sage/rings/all.py b/sage/rings/all.py
    a b Frac = FractionField 
    133133from fraction_field_element import is_FractionFieldElement
    134134
    135135# continued fractions
    136 from contfrac import continued_fraction, CFF, ContinuedFractionField
     136from continued_fractions import (
     137        continued_fraction, continued_fraction_list,
     138        CFF, ContinuedFractionField,
     139        Hirzebruch_Jung_continued_fraction_list)
     140
    137141
    138142# Arithmetic
    139143from arith import *
  • sage/rings/arith.py

    diff --git a/sage/rings/arith.py b/sage/rings/arith.py
    a b def farey(v, lim): 
    39753975            upper = mediant
    39763976
    39773977
    3978 ## def convergents_pnqn(x):
    3979 ##     """
    3980 ##     Return the pairs (pn,qn) that are the numerators and denominators
    3981 ##     of the partial convergents of the continued fraction of x.  We
    3982 ##     include (0,1) and (1,0) at the beginning of the list (these are
    3983 ##     the -2 and -1 th convergents).
    3984 ##     """
    3985 ##     v = pari(x).contfrac()
    3986 ##     w = [(0,1), (1,0)]
    3987 ##     for n in range(len(v)):
    3988 ##         pn = w[n+1][0]*v[n] + w[n][0]
    3989 ##         qn = w[n+1][1]*v[n] + w[n][1]
    3990 ##         w.append(int(pn), int(qn))
    3991 ##     return w
    3992    
    3993 def continued_fraction_list(x, partial_convergents=False, bits=None, nterms=None):
    3994     r"""
    3995     Returns the continued fraction of x as a list.
    3996 
    3997     The continued fraction expansion of `x` are the coefficients `a_i` in
    3998 
    3999     .. math::
    4000 
    4001         x = a_1 + 1/(a_2+1/(...) ... )
    4002 
    4003     with `a_1` integer and `a_2`, `...` positive integers.
    4004 
    4005     .. note::
    4006    
    4007        This may be slow for real number input, since it's implemented in pure
    4008        Python. For rational number input the PARI C library is used.
    4009    
    4010     .. SEEALSO::
    4011 
    4012          :func:`Hirzebruch_Jung_continued_fraction_list` for
    4013          Hirzebruch-Jung continued fractions.
    4014 
    4015     INPUT:
    4016 
    4017     - ``x`` -- exact rational or floating-point number. The number to
    4018       compute the continued fraction of.
    4019 
    4020     - ``partial_convergents`` -- boolean. Whether to return the partial convergents.
    4021 
    4022     - ``bits`` -- integer. the precision of the real interval field
    4023       that is used internally.
    4024 
    4025     - ``nterms`` -- integer. The upper bound on the number of terms in
    4026       the continued fraction expansion to return.
    4027 
    4028     OUTPUT:
    4029 
    4030     A lits of integers, the coefficients in the continued fraction
    4031     expansion of ``x``. If ``partial_convergents=True`` is passed, a
    4032     pair containing the coefficient list and the partial convergents
    4033     list is returned.
    4034 
    4035     EXAMPLES::
    4036    
    4037         sage: continued_fraction_list(45/17)
    4038         [2, 1, 1, 1, 5]
    4039         sage: continued_fraction_list(e, bits=20)
    4040         [2, 1, 2, 1, 1, 4, 1, 1]
    4041         sage: continued_fraction_list(e, bits=30)
    4042         [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8]
    4043         sage: continued_fraction_list(sqrt(2))
    4044         [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
    4045         sage: continued_fraction_list(sqrt(4/19))
    4046         [0, 2, 5, 1, 1, 2, 1, 16, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1]
    4047         sage: continued_fraction_list(RR(pi), partial_convergents=True)
    4048         ([3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3],
    4049          [(3, 1),
    4050           (22, 7),
    4051           (333, 106),
    4052           (355, 113),
    4053           (103993, 33102),
    4054           (104348, 33215),
    4055           (208341, 66317),
    4056           (312689, 99532),
    4057           (833719, 265381),
    4058           (1146408, 364913),
    4059           (4272943, 1360120),
    4060           (5419351, 1725033),
    4061           (80143857, 25510582),
    4062           (245850922, 78256779)])
    4063         sage: continued_fraction_list(e)
    4064         [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1]
    4065         sage: continued_fraction_list(RR(e))
    4066         [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1]
    4067         sage: continued_fraction_list(RealField(200)(e))
    4068         [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1,
    4069          14, 1, 1, 16, 1, 1, 18, 1, 1, 20, 1, 1, 22, 1, 1, 24, 1, 1,
    4070          26, 1, 1, 28, 1, 1, 30, 1, 1, 32, 1, 1, 34, 1, 1, 36, 1, 1, 38, 1, 1]
    4071 
    4072     TESTS::
    4073    
    4074         sage: continued_fraction_list(1 + 10^-10, nterms=3)
    4075         [1, 10000000000]
    4076         sage: continued_fraction_list(1 + 10^-20 - e^-100, bits=10, nterms=3)
    4077         [1, 100000000000000000000, 2688]
    4078         sage: continued_fraction_list(1 + 10^-20 - e^-100, bits=10, nterms=5)
    4079         [1, 100000000000000000000, 2688, 8, 1]
    4080         sage: continued_fraction_list(1 + 10^-20 - e^-100, bits=1000, nterms=5)
    4081         [1, 100000000000000000000, 2688, 8, 1]
    4082     """
    4083     if isinstance(x, (integer.Integer, int, long)):
    4084         if partial_convergents:
    4085             return [x], [(x,1)]
    4086         else:
    4087             return [x]
    4088 
    4089     if isinstance(x, sage.rings.rational.Rational):
    4090         if bits is not None and nterms is None:
    4091             x = RealIntervalField(bits)(x)
    4092         else:
    4093             # PARI is faster than the pure Python below, but doesn't give us the convergents.
    4094             v = pari(x).contfrac().python()
    4095             if nterms is not None:
    4096                 v = v[:nterms]
    4097             if partial_convergents:
    4098                 w = [(0,1), (1,0)]
    4099                 for a in v:
    4100                     pn = a*w[-1][0] + w[-2][0]
    4101                     qn = a*w[-1][1] + w[-2][1]
    4102                     w.append((pn, qn))
    4103                 return v, w[2:]
    4104             else:
    4105                 return v
    4106 
    4107     # Work in interval field, increasing precision as needed.
    4108     if bits is None:
    4109         try:
    4110             bits = x.prec()
    4111         except AttributeError:
    4112             bits = 53
    4113     RIF = RealIntervalField(bits)
    4114     v = []
    4115     w = [(0,1), (1,0)]
    4116     orig, x = x, RIF(x)
    4117    
    4118     while True:
    4119         try:
    4120             a = x.unique_floor()
    4121         except ValueError:
    4122             # Either we're done or we need more precision.
    4123             if nterms is None:
    4124                 break
    4125             else:
    4126                 RIF = RIF.to_prec(2*RIF.prec())
    4127                 x = RIF(orig)
    4128                 for a in v: x = ~(x-a)
    4129                 continue
    4130         if partial_convergents:
    4131             pn = a*w[-1][0] + w[-2][0]
    4132             qn = a*w[-1][1] + w[-2][1]
    4133             w.append((pn, qn))
    4134         v.append(a)
    4135         if x == a or nterms is not None and len(v) >= nterms:
    4136             break
    4137         x = ~(x-a)
    4138    
    4139     if partial_convergents:
    4140         return v, w[2:]
    4141     else:
    4142         return v
    4143 
    4144 
    4145 def Hirzebruch_Jung_continued_fraction_list(x, bits=None, nterms=None):
    4146     r"""
    4147     Return the Hirzebruch-Jung continued fraction of ``x`` as a list.
    4148 
    4149     The Hirzebruch-Jung continued fraction of `x` is similar to the
    4150     ordinary continued fraction expansion, but with minus signs. That
    4151     is, the coefficients `a_i` in
    4152 
    4153     .. math::
    4154 
    4155         x = a_1 - 1/(a_2-1/(...) ... )
    4156 
    4157     with `a_1` integer and `a_2`, `...` positive integers.
    4158 
    4159     .. SEEALSO::
    4160 
    4161          :func:`continued_fraction_list` for ordinary continued fractions.
    4162 
    4163     INPUT:
    4164 
    4165     - ``x`` -- exact rational or something that can be numerically
    4166       evaluated. The number to compute the continued fraction of.
    4167 
    4168     - ``bits`` -- integer (default: the precision of ``x``). the
    4169       precision of the real interval field that is used
    4170       internally. This is only used if ``x`` is not an exact fraction.
    4171 
    4172     - ``nterms`` -- integer (default: None). The upper bound on the
    4173       number of terms in the continued fraction expansion to return.
    4174 
    4175     OUTPUT:
    4176 
    4177     A lits of integers, the coefficients in the Hirzebruch-Jung continued
    4178     fraction expansion of ``x``.
    4179 
    4180     EXAMPLES::
    4181 
    4182         sage: Hirzebruch_Jung_continued_fraction_list(17/11)
    4183         [2, 3, 2, 2, 2, 2]
    4184         sage: Hirzebruch_Jung_continued_fraction_list(45/17)
    4185         [3, 3, 6]
    4186         sage: Hirzebruch_Jung_continued_fraction_list(e, bits=20)
    4187         [3, 4, 3, 2, 2, 2, 3, 7]
    4188         sage: Hirzebruch_Jung_continued_fraction_list(e, bits=30)
    4189         [3, 4, 3, 2, 2, 2, 3, 8, 3, 2, 2, 2, 2, 2, 2, 2, 3]
    4190         sage: Hirzebruch_Jung_continued_fraction_list(sqrt(2), bits=100)
    4191         [2, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4,
    4192          2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 2]
    4193         sage: Hirzebruch_Jung_continued_fraction_list(sqrt(4/19))
    4194         [1, 2, 7, 3, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, 7,
    4195          2, 2, 2, 7, 3, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
    4196         sage: Hirzebruch_Jung_continued_fraction_list(pi)
    4197         [4, 2, 2, 2, 2, 2, 2, 17, 294, 3, 4, 5, 16, 2, 2]
    4198         sage: Hirzebruch_Jung_continued_fraction_list(e)
    4199         [3, 4, 3, 2, 2, 2, 3, 8, 3, 2, 2, 2, 2, 2, 2, 2,
    4200          3, 12, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 10]
    4201         sage: Hirzebruch_Jung_continued_fraction_list(e, nterms=20)
    4202         [3, 4, 3, 2, 2, 2, 3, 8, 3, 2, 2, 2, 2, 2, 2, 2, 3, 12, 3, 2]
    4203         sage: len(_) == 20
    4204         True
    4205 
    4206     TESTS::
    4207 
    4208         sage: Hirzebruch_Jung_continued_fraction_list(1 - 10^-10, nterms=3)
    4209         [1, 10000000000]
    4210         sage: Hirzebruch_Jung_continued_fraction_list(1 - 10^-10 - e^-100, bits=100, nterms=5)
    4211         [1, 10000000000]
    4212         sage: Hirzebruch_Jung_continued_fraction_list(1 - 10^-20 - e^-100, bits=1000, nterms=5)
    4213         [1, 100000000000000000000, 2689, 2, 2]
    4214    """
    4215     if not isinstance(x, sage.rings.rational.Rational):
    4216         try:
    4217             x = QQ(x)
    4218         except TypeError:
    4219             # Numerically evaluate x
    4220             if bits is None:
    4221                 try:
    4222                     bits = x.prec()
    4223                 except AttributeError:
    4224                     bits = 53
    4225             x = QQ(x.n(bits))
    4226     v = []
    4227     while True:
    4228         div, mod = divmod(x.numerator(), x.denominator())
    4229         if mod == 0:
    4230             v.append(div)
    4231             break
    4232         v.append(div+1)
    4233         if nterms is not None and len(v) >= nterms:
    4234             break
    4235         x = 1/(div+1-x)
    4236     return v
    4237 
    4238 
    4239 def convergent(v, n):
    4240     r"""
    4241     Return the n-th continued fraction convergent of the continued
    4242     fraction defined by the sequence of integers v. We assume
    4243     `n \geq 0`.
    4244    
    4245     INPUT:
    4246    
    4247    
    4248     -  ``v`` - list of integers
    4249    
    4250     -  ``n`` - integer
    4251    
    4252    
    4253     OUTPUT: a rational number
    4254    
    4255     If the continued fraction integers are
    4256    
    4257     .. math::
    4258    
    4259        v = [a_0, a_1, a_2, \ldots, a_k]     
    4260    
    4261    
    4262     then ``convergent(v,2)`` is the rational number
    4263    
    4264     .. math::
    4265    
    4266        a_0 + 1/a_1     
    4267    
    4268     and ``convergent(v,k)`` is the rational number
    4269    
    4270     .. math::
    4271    
    4272        a1 + 1/(a2+1/(...) ... )         
    4273    
    4274     represented by the continued fraction.
    4275    
    4276     EXAMPLES::
    4277    
    4278         sage: convergent([2, 1, 2, 1, 1, 4, 1, 1], 7)
    4279         193/71
    4280     """
    4281     if hasattr(v, 'convergent'):
    4282         return v.convergent(n)
    4283     i = int(n)
    4284     x = QQ(v[i])
    4285     i -= 1
    4286     while i >= 0:
    4287         x = QQ(v[i]) + 1/x
    4288         i -= 1
    4289     return x
    4290 
    4291 
    4292 def convergents(v):
    4293     """
    4294     Return all the partial convergents of a continued fraction defined
    4295     by the sequence of integers v.
    4296    
    4297     If v is not a list, compute the continued fraction of v and return
    4298     its convergents (this is potentially much faster than calling
    4299     continued_fraction first, since continued fractions are
    4300     implemented using PARI and there is overhead moving the answer back
    4301     from PARI).
    4302    
    4303     INPUT:
    4304    
    4305    
    4306     -  ``v`` - list of integers or a rational number
    4307    
    4308    
    4309     OUTPUT:
    4310    
    4311    
    4312     -  ``list`` - of partial convergents, as rational
    4313        numbers
    4314    
    4315    
    4316     EXAMPLES::
    4317    
    4318         sage: convergents([2, 1, 2, 1, 1, 4, 1, 1])
    4319         [2, 3, 8/3, 11/4, 19/7, 87/32, 106/39, 193/71]
    4320     """
    4321     if hasattr(v, 'convergents'):
    4322         return v.convergents()
    4323     if not isinstance(v, list):
    4324         v = pari(v).contfrac()
    4325     w = [(0,1), (1,0)]
    4326     for n in range(len(v)):
    4327         pn = w[n+1][0]*v[n] + w[n][0]
    4328         qn = w[n+1][1]*v[n] + w[n][1]
    4329         w.append((pn, qn))
    4330     return [QQ(x) for x in w[2:]]
    4331 
     3978## Note: convergent, continued_fraction_list and convergents have been moved to
     3979## sage.rings.continued_fractions
    43323980
    43333981## def continuant(v, n=None):
    43343982##     """
    def continuant(v, n=None): 
    43754023        sage: q = continuant([1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10])
    43764024        sage: p/q
    43774025        517656/190435
    4378         sage: convergent([2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10],14)
     4026        sage: continued_fraction([2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10]).convergent(14)
    43794027        517656/190435
    43804028        sage: x = PolynomialRing(RationalField(),'x',5).gens()
    43814029        sage: continuant(x)
  • .py

    diff --git a/sage/rings/contfrac.py b/sage/rings/continued_fractions.py
    rename from sage/rings/contfrac.py
    rename to sage/rings/continued_fractions.py
    old new  
    11r"""
    2 Continued Fractions
     2Continued fractions
    33
    4 Sage implements the field ``ContinuedFractionField`` (or ``CFF``
    5 for short) of finite simple continued fractions.  This is really
    6 isomorphic to the field `\QQ` of rational numbers, but with different
    7 printing and semantics.  It should be possible to use this field in
    8 most cases where one could use `\QQ`, except arithmetic is slower.
     4A continued fraction is a representation of a real number in terms of a sequence
     5of integers denoted `[a_0; a_1, a_2, \ldots]`. The well known decimal expansion
     6is another way of representing a real number by a sequence of integers. The
     7value of a continued fraction is defined recursively as:
    98
    10 The ``continued_fraction(x)`` command returns an element of
    11 ``CFF`` that defines a continued fraction expansion to `x`.  The
    12 command ``continued_fraction(x,bits)`` computes the continued
    13 fraction expansion of an approximation to `x` with given bits of
    14 precision.  Use ``show(c)`` to see a continued fraction nicely
    15 typeset, and ``latex(c)`` to obtain the typeset version, e.g., for
    16 inclusion in a paper.
     9.. MATH::
     10
     11    [a_0; a_1, a_2, \ldots] = a_0 + \frac{1}{[a_1; a_2, \ldots]}
     12
     13In this expansion, all coefficients `a_n` are integers and only the value `a_0`
     14may be non positive. Note that `a_0` is nothing else but the floor (this remark
     15provides a way to build the continued fraction expansion from a given real
     16number).
     17
     18It is quite remarkable that
     19
     20- finite expansions correspond to rationals
     21- ultimately periodic expansions correspond to quadratic numbers (ie numbers of
     22  the form `a + b \sqrt{D}` with `a` and `b` rationals and `D` square free
     23  integer)
     24- two real numbers `x` and `y` have the same tail (up to a shift) if and only if
     25  there are integers `a,b,c,d` with `|ad - bc| = 1` and such that
     26  `y = (ax + b) / (cx + d)`.
     27
     28Moreover, the rational numbers obtained by truncation of the expansion of a real
     29number gives its so-called best approximations. For more informations on
     30continued fractions, you may have a look at :wikipedia:`Continued_fraction`.
    1731
    1832EXAMPLES:
    19 We create some example elements of the continued fraction field::
    2033
    21     sage: c = continued_fraction([1,2]); c
    22     [1, 2]
    23     sage: c = continued_fraction([3,7,15,1,292]); c
    24     [3, 7, 15, 1, 292]
    25     sage: c = continued_fraction(pi); c
    26     [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    27     sage: c.value()
    28     80143857/25510582
    29     sage: QQ(c)
    30     80143857/25510582
    31     sage: RealField(200)(QQ(c) - pi)
    32     -5.7908701643756732744264903067012149647564522968979302505514e-16
     34If you want to create the continued fraction of some real number you may either
     35use its method continued_fraction (if it exists) or call ``CFF`` (which
     36stands for real continued fraction) or the function
     37:func:`continued_fraction`::
    3338
    34 We can also create matrices, polynomials, vectors, etc., over the continued
    35 fraction field.
     39    sage: cf = CFF(13/27); cf
     40    [0; 2, 13]
     41    sage: 0 + 1/(2 + 1/13)
     42    13/27
     43    sage: cf.value()
     44    13/27
    3645
    37 ::
     46    sage: (22/45).continued_fraction()
     47    [0; 2, 22]
     48    sage: 0 + 1/(2 + 1/22)
     49    22/45
     50
     51    sage: continued_fraction(-132/17)
     52    [-8; 4, 4]
     53    sage: -8 + 1/(4 + 1/4)
     54    -132/17
     55
     56It is also possible to create a continued fraction from a list of digits::
     57
     58    sage: cf = CFF([-3,1,2,3,4,1,2,1])
     59    sage: cf.value()
     60    -465/202
     61
     62For quadratic numbers, the syntax is quite similar and the digits are
     63represented as a pair of tuples, the preperiod and the period::
     64
     65    sage: K.<sqrt2> = QuadraticField(2)
     66    sage: cf = CFF(sqrt2); cf
     67    [1; (2)*]
     68    sage: cf.value()
     69    sqrt2
     70    sage: cf.preperiod()
     71    (1,)
     72    sage: cf.period()
     73    (2,)
     74
     75    sage: (3*sqrt2 + 1/2).continued_fraction()
     76    [4; (1, 2, 1, 7)*]
     77
     78    sage: cf = CFF([(1,2,3),(1,4)]); cf
     79    [1; 2, 3, (1, 4)*]
     80    sage: cf.value()
     81    -2/23*sqrt2 + 36/23
     82
     83For an irrational and non quadratic number the continued fraction has, in general, no
     84particular structure. It is still possible to make computations::
     85
     86    sage: cf1 = continued_fraction(pi); cf1
     87    [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     88    sage: cf2 = continued_fraction(2*pi/7 + tan(2/5)); cf2
     89    [1; 3, 8, 3, 1, 33, 3, 1, 2, 1, 1, 1, 5, 2, 17, 3, 2, 1, 33, 1, ...]
     90    sage: cf1 + cf2
     91    [4; 2, 6, 13, 7, 2, 8, 2, 5, 1, 1, 5, 3, 9, 2, 14, 1, 1, 36, 2, ...]
     92
     93On the following we can remark how the tail may change even in the same
     94quadratic field::
     95
     96    sage: for i in xrange(20): print CFF(i*sqrt2)
     97    [0]
     98    [1; (2)*]
     99    [2; (1, 4)*]
     100    [4; (4, 8)*]
     101    [5; (1, 1, 1, 10)*]
     102    [7; (14)*]
     103    ...
     104    [24; (24, 48)*]
     105    [25; (2, 5, 6, 5, 2, 50)*]
     106    [26; (1, 6, 1, 2, 3, 2, 26, 2, 3, 2, 1, 6, 1, 52)*]
     107
     108Nevertheless, the tail is preserved under invertible integer homographies::
     109
     110    sage: apply_homography =  lambda m,z: (m[0,0]*z + m[0,1]) / (m[1,0]*z+m[1,1])
     111    sage: m1 = SL2Z.random_element()
     112    sage: m2 = SL2Z.random_element()
     113    sage: a = sqrt2/3
     114    sage: a.continued_fraction()
     115    [0; 2, (8, 4)*]
     116    sage: b = apply_homography(m1, a)
     117    sage: b.continued_fraction()
     118    [0; 1, 2, 1, 1, 1, 1, 6, (8, 4)*]
     119    sage: c = apply_homography(m2, a)
     120    sage: c.continued_fraction()
     121    [0; 1, 26, 1, 2, 2, (8, 4)*]
     122    sage: d = apply_homography(m1**2*m2**3, a)
     123    sage: d.continued_fraction()
     124    [0; 1, 2, 1, 1, 1, 1, 5, 2, 1, 1, 1, 1, 5, 26, 1, 2, 1, 26, 1, 2, 1, 26, 1, 2, 2, (8, 4)*]
     125
     126It is possible to make arithmetic operations on continued fractions::
     127
     128    sage: c1 = CFF([0,3,3,2,1,4,2]); c1
     129    [0; 3, 3, 2, 1, 4, 2]
     130    sage: c2 = CFF([-4,2,1,3]); c2
     131    [-4; 2, 1, 3]
     132    sage: c3 = -c1; c3
     133    [-1; 1, 2, 3, 2, 1, 4, 2]
     134    sage: c1+c2
     135    [-4; 1, 2, 628, 2]
     136    sage: c1+c3
     137    [0]
     138
     139    sage: c1/c2
     140    [-1; 1, 10, 1, 142]
     141    sage: ~c3
     142    [-4; 1, 2, 2, 1, 4, 2]
     143
     144And they can be used to create matrices, polynomials, vectors, etc::
    38145
    39146    sage: a = random_matrix(CFF, 4)
    40147    sage: a
    41     [    [-1, 2] [-1, 1, 94]      [0, 2]       [-12]]
    42     [       [-1]      [0, 2]  [-1, 1, 3]   [0, 1, 2]]
    43     [    [-3, 2]         [0]   [0, 1, 2]        [-1]]
    44     [        [1]        [-1]      [0, 3]         [1]]
     148    [      [0; 7]    [0; 4, 3]          [0]       [0; 9]]
     149    [[0; 1, 4, 4]          [1]          [0]    [0; 7, 5]]
     150    [      [0; 9] [0; 2, 8, 6]       [0; 9]    [0; 3, 3]]
     151    [   [0; 4, 5]       [0; 8]          [0]          [1]]
    45152    sage: f = a.charpoly()
    46153    sage: f
    47     [1]*x^4 + ([-2, 3])*x^3 + [14, 1, 1, 1, 9, 1, 8]*x^2 + ([-13, 4, 1, 2, 1, 1, 1, 1, 1, 2, 2])*x + [-6, 1, 5, 9, 1, 5]
     154    [1]*x^4 + ([-3; 1, 2, 1, 15])*x^3 + ... + [-1; 1, 165, 1, 1, 1, 1, 1, 1, 53, 1, 5]
    48155    sage: f(a)
    49156    [[0] [0] [0] [0]]
    50157    [[0] [0] [0] [0]]
    51158    [[0] [0] [0] [0]]
    52159    [[0] [0] [0] [0]]
    53160    sage: vector(CFF, [1/2, 2/3, 3/4, 4/5])
    54     ([0, 2], [0, 1, 2], [0, 1, 3], [0, 1, 4])
     161    ([0; 2], [0; 1, 2], [0; 1, 3], [0; 1, 4])
     162
     163The unary operations (negative and inversion) are quite fast but binary
     164operations are quite slow. It is then not adviced, if speed is needed, to use
     165them as the base class in a computation.
     166
     167.. TODO::
     168
     169    - Improve numerical approximation (the method
     170      :meth:`~ContinuedFraction_generic._mpfr_` is quite slow compared to the
     171      same method for an element of a number field)
     172
     173    - Make a class for infinite precision real number built from an infinite
     174      list (ie an infinite word)
     175
     176    - Make a class for non standard continued fractions of the form `a_0 +
     177      b_0/(a_1 + b_1/(...))` (the standard continued fractions are when all
     178      `b_n= 1` while the Hirzebruch-Jung continued fractions are the one for
     179      which `b_n = -1` for all `n`). See
     180      :wikipedia:`Generalized_continued_fraction`.
    55181
    56182AUTHORS:
    57183
    58 - Niles Johnson (2010-08): Trac #3893: ``random_element()`` should pass on ``*args`` and ``**kwds``.
     184- Niles Johnson (2010-08): ``random_element()`` should pass on ``*args`` and
     185  ``**kwds`` (:trac:`3893`).
    59186
     187- Vincent Delecroix (2013): cleaning, refactorisation, documentation (:trac:`14567`)
    60188"""
    61 from sage.structure.element         import FieldElement
    62 from sage.structure.parent_gens     import ParentWithGens
    63 from sage.libs.pari.all             import pari
     189from sage.structure.unique_representation import UniqueRepresentation
     190from sage.structure.element import Element
     191from sage.structure.element import FieldElement
    64192
    65 from field                          import Field
    66 from rational_field                 import QQ
    67 from integer_ring                   import ZZ
    68 from infinity                       import infinity
    69 from real_mpfr                      import is_RealNumber, RealField
    70 from real_double                    import RDF
    71 from arith                          import (continued_fraction_list,
    72                                             convergent, convergents)
     193from field import Field
     194from integer import Integer
     195from rational import Rational
     196from infinity import Infinity
     197from integer_ring import ZZ
     198from rational_field import QQ
     199from real_mpfr import RR
     200from real_mpfi import RealIntervalField
    73201
     202ZZ_0 = Integer(0)
     203ZZ_1 = Integer(1)
     204ZZ_2 = Integer(2)
    74205
    75 class ContinuedFractionField_class(Field):
     206def two_last_convergents(x):
    76207    """
    77     The field of all finite continued fraction of real numbers.
     208    Given the list ``x`` that consists of numbers, return the two last
     209    convergents `p_{n-1}, q_{n-1}, p_n, q_n`.
     210
     211    This function is principally used to compute the value of a ultimately
     212    periodic continued fraction.
    78213
    79214    EXAMPLES::
    80    
     215
     216        sage: from sage.rings.continued_fractions import two_last_convergents
     217        sage: two_last_convergents([])
     218        (0, 1, 1, 0)
     219        sage: two_last_convergents([0])
     220        (1, 0, 0, 1)
     221        sage: two_last_convergents([-1,1,3,2])
     222        (-1, 4, -2, 9)
     223    """
     224    p0,p1 = 0,1
     225    q0,q1 = 1,0
     226    for a in x:
     227        p0,p1 = p1,a*p1+p0
     228        q0,q1 = q1,a*q1+q0
     229    return p0,q0,p1,q1
     230
     231class ContinuedFraction_generic(FieldElement):
     232    r"""
     233    Generic class for (standard) continued fractions.
     234
     235    A derived class should implements at least:
     236
     237    - ``quotient(self, n)``: return the ``n``-th quotient of ``self``
     238
     239    - ``value(self)``: return the value of ``self`` (a real number)
     240
     241    This generic class provides:
     242
     243    - computation of convergents :meth:`convergent`, :meth:`pn` and :meth:`qn`
     244
     245    - comparison :meth:`__cmp__`
     246
     247    - elementary arithmetic function :meth:`floor`, :meth:`ceil`, :meth:`sign`
     248
     249    - numerical approximations :meth:`_mpfr_`
     250
     251    All other methods rely on :meth:`value` (and not on convergents) and may
     252    fail at execution.
     253    """
     254    def __init__(self, parent):
     255        r"""
     256        INPUT:
     257
     258        - ``parent`` -- the parent of ``self``
     259
     260        TESTS::
     261
     262            sage: TestSuite(CFF(3)).run()
     263        """
     264        FieldElement.__init__(self, parent)
     265        self._pn = [ZZ_0, ZZ_1]
     266        self._qn = [ZZ_1, ZZ_0]
     267
     268    def __abs__(self):
     269        """
     270        Return absolute value of self.
     271
     272        EXAMPLES::
     273
     274            sage: a = CFF(-17/389); a
     275            [-1; 1, 21, 1, 7, 2]
     276            sage: abs(a)
     277            [0; 22, 1, 7, 2]
     278            sage: QQ(abs(a))
     279            17/389
     280        """
     281        if self.quotient(0) >= 0:
     282            return self
     283        return self.__neg__()
     284
     285    def __cmp__(self, other):
     286        """
     287        EXAMPLES::
     288
     289            sage: a = CFF(-17/389)
     290            sage: b = CFF(1/389)
     291            sage: c = CFF([(),(1,)])     # the golden ratio
     292            sage: d = CFF([(-1,),(1,)])
     293            sage: d < a and a < b and b < c
     294            True
     295            sage: d >= a
     296            False
     297        """
     298        i = 0
     299        while True:
     300            a = self.quotient(i)
     301            b = other.quotient(i)
     302            test = cmp(a,b)
     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
     308                return 0
     309            i += 1
     310
     311    def _mpfr_(self, R):
     312        r"""
     313        Return a numerical approximation of ``self`` in the real mpfr ring ``R``.
     314
     315        It is guaranteed that the result is exact: when the rounding mode of
     316        ``R`` is 'RNDN' then the result is the nearest binary number of ``R`` to
     317        ``self``. The other rounding mode are 'RNDD' (toward +infinity), 'RNDU'
     318        (toward -infinity) and 'RNDZ' (toward zero).
     319
     320        EXAMPLES::
     321
     322            sage: continued_fraction(1/2).n()
     323            0.500000000000000
     324            sage: continued_fraction([0,4]).n()
     325            0.250000000000000
     326            sage: continued_fraction([12,1,3,4,2,2,3,1,2]).n(digits=4)
     327            12.76
     328
     329            sage: continued_fraction(12/7).n(digits=13) == (12/7).n(digits=13)
     330            True
     331            sage: continued_fraction(-14/333).n(digits=21) == (-14/333).n(digits=21)
     332            True
     333
     334            sage: a = (106*pi - 333) / (355 - 113*pi) - 292
     335            sage: cf = continued_fraction(a); cf
     336            [0; 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, ...]
     337            sage: cf.n(digits=3)
     338            0.635
     339            sage: cf.n(digits=4)
     340            0.6346
     341            sage: cf.n(digits=5)
     342            0.63459
     343            sage: cf.n(digits=6)
     344            0.634591
     345            sage: cf.n(digits=7)
     346            0.6345910
     347            sage: cf.n(digits=8)
     348            0.63459101
     349
     350            sage: K.<a> = NumberField(x^3-2, 'a', embedding=1.25)
     351            sage: b = 504/253*a^2 + 635/253*a + 661/253
     352            sage: cf = continued_fraction(b); cf
     353            [8; 1, 14, 1, 10, 2, 1, 4, 12, 2, 3, 2, 1, 3, 4, 1, 1, 2, 14, 3, ...]
     354            sage: cf.n(digits=3)
     355            8.94
     356            sage: cf.n(digits=6)
     357            8.93715
     358            sage: cf.n(digits=7)
     359            8.937154
     360            sage: cf.n(digits=8)
     361            8.9371541
     362            sage: cf.n(digits=9)
     363            8.93715414
     364            sage: cf.n(digits=10)
     365            8.937154138
     366            sage: cf.n(digits=11)
     367            8.9371541378
     368
     369        TESTS::
     370
     371        We check that the rounding works as expected, at least in the rational
     372        case::
     373
     374            sage: for _ in xrange(100):
     375            ....:     a = QQ.random_element(num_bound=1<<64)
     376            ....:     cf = continued_fraction(a)
     377            ....:     for prec in 17,24,53,128,256:
     378            ....:         for rnd in 'RNDN','RNDD','RNDU','RNDZ':
     379            ....:             R = RealField(prec=prec, rnd=rnd)
     380            ....:             assert R(cf) == R(a)
     381        """
     382        # 1. integer case
     383        if self.quotient(1) == 0:
     384            return R(self.quotient(0))
     385
     386        # 2. negative numbers
     387        # TODO: it is possible to deal with negative values. The only problem is
     388        # that we need to find the good value for N (which involves
     389        # self.quotient(k) for k=0,1,2)
     390        if self.quotient(0) < 0:
     391            rnd = R.rounding_mode()
     392            if rnd == 'RNDN' or rnd == 'RNDZ':
     393                return -R(-self)
     394            elif rnd == 'RNDD':
     395                r = R(-self)
     396                s,m,e = r.sign_mantissa_exponent()
     397                if e < 0:
     398                    return -(R(m+1) >> (-e))
     399                return -(R(m+1) << e)
     400            else:
     401                r = R(-self)
     402                s,m,e = r.sign_mantissa_exponent()
     403                if e < 0:
     404                    return -(R(m-1) >> (-e))
     405                return -(R(m-1) << e)
     406
     407        # 3. positive non integer
     408        if self.quotient(0) == 0: # 0 <= self < 1
     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:
     411                # if self is of the form [0; 2^N] then we need the following
     412                N -= 1
     413        else: # self > 1
     414            N = R.prec() - self.quotient(0).nbits()
     415
     416        # even/odd convergents are respectively below/above
     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)
     422        while (m_odd - m_even) > 1:
     423            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
     428
     429        assert m_odd.nbits() == R.prec() or m_even.nbits() == R.prec()
     430
     431        if m_even == m_odd:  # no need to worry (we have a decimal number)
     432            return R(m_even)>>N
     433
     434        # check ordering
     435        # m_even/2^N <= p_even/q_even <= self <= p_odd/q_odd <= m_odd/2^N
     436        assert m_odd == m_even + 1
     437        assert m_even / (ZZ_1<<N) <= p_even/q_even
     438        assert p_even / q_even <= p_odd / q_odd
     439        assert p_odd / q_odd <= m_odd / (ZZ_1<<N)
     440
     441        rnd = R.rounding_mode()
     442        if rnd == 'RNDN': # round to the nearest
     443            # in order to find the nearest approximation we possibly need to
     444            # augment our precision on convergents.
     445            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:
     448                    return R(m_even) >> N
     449                if p_even<<(N+1) >= (2*m_even+1) * q_even:
     450                    return R(m_odd)>>N
     451                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)
     454        elif rnd == 'RNDU':  # round up (toward +infinity)
     455            return R(m_odd)>>N
     456        elif rnd == 'RNDD':  # round down (toward -infinity)
     457            return R(m_even)>>N
     458        elif rnd == 'RNDZ':  # round toward zero
     459            if m_even.sign() == 1:
     460                return R(m_even)>>N
     461            else:
     462                return R(m_odd)>>N
     463        else:
     464            raise ValueError("%s unknown rounding mode"%rounding)
     465
     466    def __float__(self):
     467        """
     468        EXAMPLES::
     469
     470            sage: a = CFF(-17/389); a
     471            [-1; 1, 21, 1, 7, 2]
     472            sage: float(a)
     473            -0.043701799485861184
     474            sage: float(-17/389)
     475            -0.043701799485861184
     476        """
     477        return float(self._mpfr_(RR))
     478
     479    def pn(self, n):
     480        """
     481        Return the numerator of the `n`-th partial convergent of ``self``.
     482
     483        EXAMPLES::
     484
     485            sage: c = continued_fraction(pi); c
     486            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     487            sage: c.pn(0)
     488            3
     489            sage: c.pn(12)
     490            80143857
     491            sage: c.pn(152)
     492            3943771611212266962743738812600748213157266596588744951727393497446921245353005283
     493        """
     494        n = Integer(n)
     495
     496        p = self._pn
     497        q = self._qn
     498
     499        if n < -2:
     500            raise ValueError("n must be at least -2")
     501
     502        for k in xrange(len(p), n+3):
     503            x = self.quotient(k-2)
     504            if x == 0 and k != 2:
     505                return p[-1]
     506            p.append(x*p[k-1] + p[k-2])
     507            q.append(x*q[k-1] + q[k-2])
     508
     509        return p[n+2]
     510
     511    def qn(self, n):
     512        """
     513        Return the denominator of the ``n``-th partial convergent of ``self``.
     514
     515        EXAMPLES::
     516
     517            sage: c = continued_fraction(pi); c
     518            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     519            sage: c.qn(0)
     520            1
     521            sage: c.qn(12)
     522            25510582
     523            sage: c.qn(152)
     524            1255341492699841451528811722575401081588363886480089431843026103930863337221076748
     525        """
     526        self.pn(n)   # ! silent computation of qn
     527        if len(self._qn) < n+3:
     528            return self._qn[-1]
     529        return self._qn[n+2]
     530
     531    def convergent(self, n):
     532        """
     533        Return the ``n``-th partial convergent to self.
     534
     535        EXAMPLES::
     536
     537            sage: a = CFF(pi); a
     538            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     539            sage: a.convergent(3)
     540            355/113
     541            sage: a.convergent(15)
     542            411557987/131002976
     543        """
     544        return self.pn(n) / self.qn(n)
     545
     546    def __getitem__(self, n):
     547        r"""
     548        Return the ``n``-th partial quotient of ``self`` or a continued fraction
     549        associated to a sublist of the partial quotients of ``self``.
     550
     551        TESTS::
     552
     553            sage: cf1 = continued_fraction(pi); cf1
     554            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     555            sage: cf2 = continued_fraction(QuadraticField(2).gen()); cf2
     556            [1; (2)*]
     557            sage: cf3 = continued_fraction(4/17); cf3
     558            [0; 4, 4]
     559            sage: cf1[3:17]
     560            [1; 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2]
     561            sage: cf2[:10]
     562            [1; 2, 2, 2, 2, 2, 2, 2, 2, 2]
     563            sage: cf3[1:16]
     564            [4; 4]
     565
     566        Be careful that the truncation of an infinite continued fraction might
     567        be shorter by one::
     568
     569            sage: len(CFF(golden_ratio)[:8])
     570            7
     571        """
     572        if isinstance(n, slice):
     573            if self.length() == Infinity:
     574                if n.step is None:
     575                    step = 1
     576                else:
     577                    try:
     578                        step = n.step.__index__()
     579                    except (AttributeError,ValueError):
     580                        raise ValueError("step should be a non zero integer")
     581                if n.start is None:
     582                    if step < 0:
     583                        raise ValueError("if step is negative, start can not be None")
     584                    start = 0
     585                else:
     586                    try:
     587                        start = n.start.__index__()
     588                    except (AttributeError,ValueError):
     589                        raise ValueError("start should be an integer")
     590                if n.stop is None:
     591                    if step > 0:
     592                        raise ValueError("infinite list!")
     593                    stop = -1
     594                else:
     595                    try:
     596                        stop = n.stop.__index__()
     597                    except (AttributeError,ValueError):
     598                        raise ValueError("stop should be an integer")
     599                return self.parent()([self.quotient(i) for i in xrange(start,stop,step)])
     600            start, stop, step = n.indices(self.length())
     601            return self.parent()([self.quotient(i) for i in xrange(start,stop,step)])
     602        try:
     603            n = n.__index__()
     604        except (AttributeError,ValueError):
     605            raise ValueError("n (=%s) should be an integer"%n)
     606        if n < 0:
     607            raise ValueError("n (=%s) should be positive"%n)
     608        q = self.quotient(n)
     609        if n > 0 and q == 0:
     610            raise IndexError("index out of range")
     611
     612    def __iter__(self):
     613        r"""
     614        Iterate over the partial quotient of self.
     615
     616        EXAMPLES::
     617
     618            sage: cf = continued_fraction(pi)
     619            sage: i = iter(cf)
     620            sage: [i.next() for _ in xrange(10)]
     621            [3, 7, 15, 1, 292, 1, 1, 1, 2, 1]
     622            sage: [i.next() for _ in xrange(10)]
     623            [3, 1, 14, 2, 1, 1, 2, 2, 2, 2]
     624            sage: [i.next() for _ in xrange(10)]
     625            [1, 84, 2, 1, 1, 15, 3, 13, 1, 4]
     626        """
     627        yield self.quotient(0)
     628        i = 1
     629        while True:
     630            q = self.quotient(i)
     631            if not q:
     632                break
     633            yield q
     634            i += 1
     635
     636    def __int__(self):
     637        """
     638        EXAMPLES::
     639
     640            sage: a = CFF(-17/389); a
     641            [-1; 1, 21, 1, 7, 2]
     642            sage: int(a)
     643            -1
     644        """
     645        return int(self.quotient(0))
     646
     647    def __long__(self):
     648        """
     649        EXAMPLES::
     650
     651            sage: a = CFF(-17/389); a
     652            [-1; 1, 21, 1, 7, 2]
     653            sage: long(a)
     654            -1L
     655        """
     656        return long(self.quotient(0))
     657
     658    def sign(self):
     659        r"""
     660        Returns the sign of self as an Integer.
     661
     662        The sign is defined to be ``0`` if ``self`` is ``0``, ``1`` if ``self``
     663        is positive and ``-1`` if ``self`` is negative.
     664
     665        EXAMPLES::
     666
     667            sage: continued_fraction(tan(pi/7)).sign()
     668            1
     669            sage: continued_fraction(-34/2115).sign()
     670            -1
     671            sage: continued_fraction([0]).sign()
     672            0
     673        """
     674        if self.quotient(0) == 0:
     675            if self.quotient(1) == 0:
     676                return ZZ_0
     677            return ZZ_1
     678        return self.quotient(0).sign()
     679
     680    def floor(self):
     681        r"""
     682        Return the floor of ``self``.
     683
     684        EXAMPLES::
     685
     686            sage: cf = CFF([2,1,2,3])
     687            sage: cf.floor()
     688            2
     689        """
     690        return self.quotient(0)
     691
     692    def ceil(self):
     693        r"""
     694        Return the ceil of ``self``.
     695
     696        EXAMPLES::
     697
     698            sage: cf = CFF([2,1,3,4])
     699            sage: cf.ceil()
     700            3
     701        """
     702        if self.length() == 1:
     703            return self.quotient(0)
     704        return self.quotient(0)+1
     705
     706    def __neg__(self):
     707        r"""
     708        Return the opposite of ``self``.
     709
     710        EXAMPLES::
     711
     712            sage: -CFF(pi) == CFF(-pi)
     713            True
     714        """
     715        return self.parent()(-self.value())
     716
     717    def __invert__(self):
     718        r"""
     719        Return the inverse of ``self``.
     720
     721        EXAMPLES::
     722
     723            sage: ~CFF(pi) == CFF(~pi)
     724            True
     725        """
     726        return self.parent()(~self.value())
     727
     728    def _add_(self, other):
     729        """
     730        EXAMPLES::
     731
     732            sage: a = CFF(-17/389)
     733            sage: b = CFF(1/389)
     734            sage: c = a+b; c
     735            [-1; 1, 23, 3, 5]
     736            sage: c.value()
     737            -16/389
     738
     739        Note:
     740
     741        The algorithm is very naive !
     742        """
     743        try:
     744            return self.parent()(self.value() + other.value())
     745        except (ValueError,TypeError):
     746            raise NotImplementedError
     747
     748    def _sub_(self, other):
     749        """
     750        EXAMPLES::
     751
     752            sage: a = CFF(-17/389)
     753            sage: b = CFF(1/389)
     754            sage: c = a - b; c
     755            [-1; 1, 20, 1, 1, 1, 1, 3]
     756            sage: c.value()
     757            -18/389
     758        """
     759        try:
     760            return self.parent()(self.value() - other.value())
     761        except (ValueError,TypeError):
     762            raise NotImplementedError
     763
     764    def _mul_(self, other):
     765        """
     766        EXAMPLES::
     767
     768            sage: a = CFF(-17/389)
     769            sage: b = CFF(1/389)
     770            sage: c = a * b; c
     771            [-1; 1, 8900, 4, 4]
     772            sage: c.value(), (-1/389)*(17/389)
     773            (-17/151321, -17/151321)
     774        """
     775        try:
     776            return self.parent()(self.value() * other.value())
     777        except (ValueError,TypeError):
     778            raise NotImplementedError
     779
     780    def _div_(self, other):
     781        """
     782        EXAMPLES::
     783
     784            sage: a = CFF(-17/389)
     785            sage: b = CFF(1/389)
     786            sage: c = a / b; c
     787            [-17]
     788            sage: c.value(), (17/389) / (-1/389)
     789            (-17, -17)
     790        """
     791        try:
     792            return self.parent()(self.value() / other.value())
     793        except (ValueError,TypeError):
     794            raise NotImplementedError
     795
     796    def __nonzero__(self):
     797        """
     798        Return False if self is zero.
     799
     800        EXAMPLES::
     801
     802            sage: CFF(0).is_zero()    # indirect doctest
     803            True
     804            sage: CFF(1).is_zero()    # indirect doctest
     805            False
     806            sage: CFF([(),(1,)]).is_zero()     # indirect doctest
     807            False
     808            sage: CFF([(0,),(1,2)]).is_zero()  # indirect doctest
     809            False
     810        """
     811        return bool(self.quotient(0)) or bool(self.quotient(1))
     812
     813    def is_zero(self):
     814        r"""
     815        Test whether ``self`` is zero.
     816
     817        EXAMPLES::
     818
     819            sage: CFF(0).is_zero()
     820            True
     821            sage: CFF((0,1)).is_zero()
     822            False
     823            sage: CFF(-1/2).is_zero()
     824            False
     825            sage: CFF(pi).is_zero()
     826            False
     827        """
     828        return self.quotient(0) == ZZ_0 and self.quotient(1) == ZZ_0
     829
     830    def is_one(self):
     831        r"""
     832        Test whether ``self`` is one.
     833
     834        EXAMPLES::
     835
     836            sage: CFF(1).is_one()
     837            True
     838            sage: CFF(5/4).is_one()
     839            False
     840            sage: CFF(0).is_one()
     841            False
     842            sage: CFF(pi).is_one()
     843            False
     844        """
     845        return self.quotient(0) == ZZ_1 and self.quotient(1) == ZZ_0
     846
     847    def additive_order(self):
     848        """
     849        Return the additive order of this continued fraction,
     850        which we defined to be the additive order of its value.
     851
     852        EXAMPLES::
     853
     854            sage: CFF(-1).additive_order()
     855            +Infinity
     856            sage: CFF(0).additive_order()
     857            1
     858        """
     859        return Infinity if self else ZZ_1
     860
     861    def multiplicative_order(self):
     862        """
     863        Return the multiplicative order of this continued fraction,
     864        which we defined to be the multiplicative order of its value.
     865
     866        EXAMPLES::
     867
     868            sage: CFF(-1).multiplicative_order()
     869            2
     870            sage: CFF(1).multiplicative_order()
     871            1
     872            sage: CFF(pi).multiplicative_order()
     873            +Infinity
     874        """
     875        if self.is_zero():
     876            return Infinity
     877        if self.is_one():
     878            return ZZ_1
     879        if (-self).is_one():
     880            return ZZ_2
     881        return Infinity
     882
     883class ContinuedFraction_periodic(ContinuedFraction_generic):
     884    r"""
     885    Continued fraction associated with rational or quadratic number.
     886
     887    A rational number has a finite continued fraction expansion (or ultimately
     888    0). The one of a quadratic number, ie a number of the form `a + b \sqrt{D}`
     889    with `a` and `b` rational, is ultimately periodic.
     890
     891    .. NOTE::
     892
     893        This class stores a tuple ``_x1`` for the preperiod and a tuple ``_x2``
     894        for the period. In the purely periodic case ``_x1`` is empty while in
     895        the rational case ``_x2`` is the tuple ``(0,)``.
     896    """
     897    def __init__(self, parent, x1, x2, check=True):
     898        r"""
     899        TESTS::
     900
     901            sage: TestSuite(CFF([(1,2),(3,4,5)])).run()
     902        """
     903        ContinuedFraction_generic.__init__(self, parent)
     904        self._x1 = tuple(x1)
     905        self._x2 = tuple(x2)
     906
     907    def period(self):
     908        r"""
     909        Return the periodic part of ``self``.
     910
     911        EXAMPLES::
     912
     913            sage: K.<sqrt3> = QuadraticField(3)
     914            sage: cf = continued_fraction(sqrt3); cf
     915            [1; (1, 2)*]
     916            sage: cf.period()
     917            (1, 2)
     918
     919            sage: for k in xsrange(2,40):
     920            ....:     if not k.is_square():
     921            ....:         s = QuadraticField(k).gen()
     922            ....:         cf = continued_fraction(s)
     923            ....:         print '%2d %d %s'%(k, len(cf.period()), cf)
     924             2 1 [1; (2)*]
     925             3 2 [1; (1, 2)*]
     926             5 1 [2; (4)*]
     927             6 2 [2; (2, 4)*]
     928             7 4 [2; (1, 1, 1, 4)*]
     929             8 2 [2; (1, 4)*]
     930            10 1 [3; (6)*]
     931            11 2 [3; (3, 6)*]
     932            12 2 [3; (2, 6)*]
     933            13 5 [3; (1, 1, 1, 1, 6)*]
     934            14 4 [3; (1, 2, 1, 6)*]
     935            ...
     936            35 2 [5; (1, 10)*]
     937            37 1 [6; (12)*]
     938            38 2 [6; (6, 12)*]
     939            39 2 [6; (4, 12)*]
     940        """
     941        if self._x2[0] == ZZ_0:
     942            return ()
     943        return self._x2
     944
     945    def preperiod(self):
     946        r"""
     947        Return the preperiodic part of ``self``.
     948
     949        EXAMPLES::
     950
     951            sage: K.<sqrt3> = QuadraticField(3)
     952            sage: cf = continued_fraction(sqrt3); cf
     953            [1; (1, 2)*]
     954            sage: cf.preperiod()
     955            (1,)
     956
     957            sage: cf = continued_fraction(sqrt3/7); cf
     958            [0; 4, (24, 8)*]
     959            sage: cf.preperiod()
     960            (0, 4)
     961        """
     962        return self._x1
     963
     964    def quotient(self, n):
     965        r"""
     966        Return the ``n``-th quotient of ``self``.
     967
     968        EXAMPLES::
     969
     970            sage: cf = CFF([(12,5),(1,3)])
     971            sage: [cf.quotient(i) for i in xrange(10)]
     972            [12, 5, 1, 3, 1, 3, 1, 3, 1, 3]
     973        """
     974        n = int(n)
     975        if n < 0:
     976            raise ValueError("n (=%d) should be positive"%n)
     977        if n < len(self._x1):
     978            return self._x1[n]
     979        return self._x2[(n-len(self._x1))%len(self._x2)]
     980
     981    def length(self):
     982        r"""
     983        Returns the number of partial quotients of ``self``.
     984
     985        EXAMPLES::
     986
     987            sage: CFF(2/5).length()
     988            3
     989            sage: cf = CFF([(0,1),(2,)]); cf
     990            [0; 1, (2)*]
     991            sage: cf.length()
     992            +Infinity
     993        """
     994        if len(self._x2) > 1 or self._x2[0] != ZZ_0:
     995            return Infinity
     996        return ZZ(len(self._x1))
     997
     998    def __cmp__(self, other):
     999        r"""
     1000        EXAMPLES::
     1001
     1002            sage: a = CFF([(0,),(1,2,3,1,2,3,1)]); a.n()
     1003            0.694249167819459
     1004            sage: b = CFF([(0,),(1,2,3)]); b.n()
     1005            0.694254176766073
     1006            sage: c = CFF([(0,1),(2,3)]); c.n()
     1007            0.696140478029631
     1008            sage: d = CFF([(0,1,2),(3,)]); d.n()
     1009            0.697224362268005
     1010            sage: a < b and a < c and a < d
     1011            True
     1012            sage: b < c and b < d and c < d
     1013            True
     1014            sage: b == c
     1015            False
     1016            sage: c > c
     1017            False
     1018            sage: b >= d
     1019            False
     1020        """
     1021        if isinstance(other, ContinuedFraction_periodic):
     1022            n = max(len(self._x1) + 2*len(self._x2),
     1023                    len(other._x1) + 2*len(other._x2))
     1024            for i in xrange(n):
     1025                a = self.quotient(i)
     1026                b = other.quotient(i)
     1027                test = cmp(a,b)
     1028                if test == 1:
     1029                    return -1 if i%2 else 1
     1030                if test == -1:
     1031                    return 1 if i%2 else -1
     1032            return 0
     1033
     1034        return ContinuedFraction_generic.__cmp__(self, other)
     1035
     1036    def value(self):
     1037        r"""
     1038        Return the value of ``self`` as a quadratic number (with square free
     1039        discriminant).
     1040
     1041        EXAMPLES::
     1042
     1043        Some purely periodic examples::
     1044
     1045            sage: cf = CFF([(),(2,)]); cf
     1046            [(2)*]
     1047            sage: v = cf.value(); v
     1048            sqrt2 + 1
     1049            sage: v.continued_fraction()
     1050            [(2)*]
     1051
     1052            sage: cf = CFF([(),(1,2)]); cf
     1053            [(1, 2)*]
     1054            sage: v = cf.value(); v
     1055            1/2*sqrt3 + 1/2
     1056            sage: v.continued_fraction()
     1057            [(1, 2)*]
     1058
     1059            sage: cf = CFF([(),(3,2,1)])
     1060            sage: cf.value().continued_fraction() == cf
     1061            True
     1062            sage: cf = CFF([(),(1,3,1,5)])
     1063            sage: cf.value().continued_fraction() == cf
     1064            True
     1065
     1066        Some ultimately periodic but non periodic examples::
     1067
     1068            sage: cf = CFF([(1,),(2,)]); cf
     1069            [1; (2)*]
     1070            sage: v = cf.value(); v
     1071            sqrt2
     1072            sage: v.continued_fraction()
     1073            [1; (2)*]
     1074
     1075            sage: cf = CFF([(1,3),(1,2)]); cf
     1076            [1; 3, (1, 2)*]
     1077            sage: v = cf.value(); v
     1078            -sqrt3 + 3
     1079            sage: v.continued_fraction()
     1080            [1; 3, (1, 2)*]
     1081
     1082            sage: cf = CFF([(-5,18), (1,3,1,5)])
     1083            sage: cf.value().continued_fraction() == cf
     1084            True
     1085            sage: cf = CFF([(-1,),(1,)])
     1086            sage: cf.value().continued_fraction() == cf
     1087            True
     1088
     1089        TESTS::
     1090
     1091            sage: a1 = ((0,1),(2,3))
     1092            sage: a2 = ((-12,1,1),(2,3,2,4))
     1093            sage: a3 = ((1,),(1,2))
     1094            sage: a4 = ((-2,2),(1,124,13))
     1095            sage: a5 = ((0,),(1,))
     1096            sage: for a in a1,a2,a3,a4,a5:
     1097            ....:     cf = CFF(a)
     1098            ....:     assert cf.value().continued_fraction() == cf
     1099        """
     1100        if self._x1 and self._x1[0] < 0:
     1101            return -(-self).value()
     1102
     1103        if len(self._x2) == 1 and not self._x2[0]:
     1104            return self._rational_()
     1105
     1106        # determine the equation for the purely periodic cont. frac. determined
     1107        # by self._x2
     1108        p0,q0,p1,q1 = two_last_convergents(self._x2)
     1109
     1110        # now x is one of the root of the equation
     1111        #   q1 x^2 + (q0 - p1) x - p0 = 0
     1112        from sage.rings.number_field.number_field import QuadraticField
     1113        from sage.misc.functional import squarefree_part
     1114        D = (q0-p1)**2 + 4*q1*p0
     1115        DD = squarefree_part(D)
     1116        Q = QuadraticField(DD, 'sqrt%d'%DD)
     1117        x = ((p1 - q0) + (D/DD).sqrt() * Q.gen()) / (2*q1)
     1118
     1119        # we add the preperiod
     1120        p0,q0,p1,q1 = two_last_convergents(self._x1)
     1121        return (p1*x + p0) / (q1*x + q0)
     1122
     1123    def _repr_(self):
     1124        r"""
     1125        TESTS::
     1126
     1127            sage: a = continued_fraction(pi.n()); a
     1128            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3]
     1129            sage: a.rename('continued fraction of pi')
     1130            sage: a
     1131            continued fraction of pi
     1132
     1133            sage: CFF([(0,1),(2,)])
     1134            [0; 1, (2)*]
     1135            sage: CFF([(),(1,3)])
     1136            [(1, 3)*]
     1137        """
     1138        if len(self._x2) == 1 and not self._x2[0]: # rational
     1139            if len(self._x1) == 1:
     1140                return '[%d]'%self._x1[0]
     1141            return '[%d; '%self._x1[0] + ', '.join(str(a) for a in self._x1[1:]) + ']'
     1142
     1143        period = '(' + ', '.join(str(a) for a in self._x2) + ')*'
     1144        if not self._x1: # purely periodic case
     1145            return '[' + period + ']'
     1146
     1147        if len(self._x1) == 1:
     1148            return '[%d; '%self._x1[0] + period + ']'
     1149        return '[%d; '%self._x1[0] + ', '.join(str(a) for a in self._x1[1:]) + ', ' + period + ']'
     1150
     1151    def __reduce__(self):
     1152        r"""
     1153        Pickling support.
     1154
     1155        EXAMPLES::
     1156
     1157            sage: a = CFF([1,2,3])
     1158            sage: loads(dumps(a)) == a
     1159            True
     1160
     1161            sage: a = CFF([(-1,2),(3,1,4)])
     1162            sage: loads(dumps(a)) == a
     1163            True
     1164        """
     1165        return self.parent(),((self._x1,self._x2),)
     1166
     1167    def convergents(self):
     1168        """
     1169        Return the list of partial convergents of ``self``.
     1170
     1171        EXAMPLES::
     1172
     1173            sage: a = CFF(23/157); a
     1174            [0; 6, 1, 4, 1, 3]
     1175            sage: a.convergents()
     1176            [0, 1/6, 1/7, 5/34, 6/41, 23/157]
     1177        """
     1178        return [self.pn(n) / self.qn(n) for n in xrange(len(self))]
     1179
     1180    def __len__(self):
     1181        """
     1182        Return the number of terms in this continued fraction.
     1183
     1184        EXAMPLES::
     1185
     1186            sage: len(continued_fraction([1,2,3,4,5]) )
     1187            5
     1188        """
     1189        if len(self._x2) == 1 and self._x2[0] == ZZ_0:
     1190            return len(self._x1)
     1191        raise ValueError("the length is infinite")
     1192
     1193    def _rational_(self):
     1194        """
     1195        EXAMPLES::
     1196
     1197            sage: a = CFF(-17/389); a
     1198            [-1; 1, 21, 1, 7, 2]
     1199            sage: a._rational_()
     1200            -17/389
     1201            sage: QQ(a)
     1202            -17/389
     1203        """
     1204        if len(self._x2) > 1 or self._x2[0] != 0:
     1205            raise ValueError("this is not a rational!")
     1206        n = len(self)
     1207        return self.pn(n-1) / self.qn(n-1)
     1208
     1209    def _latex_(self):
     1210        """
     1211        EXAMPLES::
     1212
     1213            sage: a = CFF(-17/389)
     1214            sage: latex(a)
     1215            -1+ \frac{\displaystyle 1}{\displaystyle 1+ \frac{\displaystyle 1}{\displaystyle 21+ \frac{\displaystyle 1}{\displaystyle 1+ \frac{\displaystyle 1}{\displaystyle 7+ \frac{\displaystyle 1}{\displaystyle 2}}}}}
     1216        """
     1217        if self._x2 != (0,):
     1218            raise NotImplementedError("latex not implemented for non rational continued fractions")
     1219        v = self._x1
     1220        if len(v) == 0:
     1221            return '0'
     1222        s = str(v[0])
     1223        for i in range(1,len(v)):
     1224            s += '+ \\frac{\\displaystyle 1}{\\displaystyle %s'%v[i]
     1225        s += '}'*(len(v)-1)
     1226        return s
     1227
     1228    def __invert__(self):
     1229        """
     1230        Return the multiplicative inverse of self.
     1231
     1232        EXAMPLES::
     1233
     1234            sage: a = CFF(13/25)
     1235            sage: ~a == CFF(25/13)
     1236            True
     1237            sage: a * ~a
     1238            [1]
     1239
     1240            sage: a = CFF(-17/253)
     1241            sage: ~a == CFF(-253/17)
     1242            True
     1243            sage: a * ~a
     1244            [1]
     1245
     1246            sage: K.<sqrt5> = QuadraticField(5)
     1247            sage: a1 = (sqrt5+1)/2
     1248            sage: c1 = a1.continued_fraction(); c1
     1249            [(1)*]
     1250            sage: ~c1
     1251            [0; (1)*]
     1252            sage: ~c1 * c1
     1253            [1]
     1254            sage: c2 = (sqrt5/3 + 1/7).continued_fraction(); c2
     1255            [0; 1, (7, 1, 17, ..., 1, 2)*]
     1256            sage: ~c2
     1257            [1; (7, 1, 17, ..., 1, 2)*]
     1258            sage: c2 * ~c2
     1259            [1]
     1260        """
     1261        if not self:
     1262            raise ZeroDivisionError("Rational division by 0")
     1263        if self._x1:
     1264            if self._x1[0] < 0:
     1265                return -(-self).__invert__()
     1266            if self._x1[0] == 0:
     1267                return self.parent()([self._x1[1:],self._x2])
     1268        return self.parent()([(0,) + self._x1, self._x2])
     1269
     1270    def __neg__(self):
     1271        """
     1272        Return additive inverse of self.
     1273
     1274        EXAMPLES::
     1275
     1276            sage: a = CFF(-17/389); a
     1277            [-1; 1, 21, 1, 7, 2]
     1278            sage: -a
     1279            [0; 22, 1, 7, 2]
     1280            sage: -(-a)
     1281            [-1; 1, 21, 1, 7, 2]
     1282            sage: QQ(-a)
     1283            17/389
     1284
     1285            sage: -CFF([2])
     1286            [-2]
     1287            sage: -CFF([1,2])
     1288            [-2; 2]
     1289
     1290            sage: a = CFF([(1,2),(1,)]); a
     1291            [1; 2, (1)*]
     1292            sage: -a
     1293            [-2; (1)*]
     1294            sage: -(-a)
     1295            [1; 2, (1)*]
     1296
     1297            sage: a = CFF([(),(1,2,3)]); a
     1298            [(1, 2, 3)*]
     1299            sage: -a
     1300            [-2; 1, 1, (3, 1, 2)*]
     1301            sage: -(-a)
     1302            [(1, 2, 3)*]
     1303
     1304            sage: -CFF([0])
     1305            [0]
     1306            sage: -CFF([1])
     1307            [-1]
     1308            sage: -CFF([-1])
     1309            [1]
     1310        """
     1311        x1 = self._x1
     1312
     1313        # integer case
     1314        if len(self._x2) == 1 and self._x2[0] == 0 and len(self._x1) == 1:
     1315            return self.parent()(-self._x1[0])
     1316
     1317        # to make all other cases work, we need 3 elements in x1
     1318        if len(x1) < 3:
     1319            x1 = self._x1 + self._x2
     1320        if len(x1) < 3:
     1321            x1 += self._x2
     1322
     1323        # we are yet sure that there are at least two elements...
     1324        if x1[0] == 0 and x1[1] == 0:
     1325            return self.parent().zero()
     1326
     1327        if len(x1) < 3:
     1328            x1 += self._x2
     1329
     1330        if x1[1] == 1:
     1331            return self.parent()([(-x1[0]-1, x1[2]+1) + x1[3:], self._x2])
     1332        return self.parent()([(-x1[0]-1, ZZ_1, x1[1]-1) + x1[2:], self._x2])
     1333
     1334class ContinuedFraction_irrational(ContinuedFraction_generic):
     1335    r"""
     1336    Continued fraction of an exact irrational number.
     1337
     1338    EXAMPLES::
     1339
     1340        sage: continued_fraction(pi)
     1341        [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     1342        sage: continued_fraction(e)
     1343        [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1, ...]
     1344    """
     1345    def __init__(self, parent, x):
     1346        r"""
     1347        INPUT:
     1348
     1349        - ``parent`` -- the parent of that continued fraction
     1350
     1351        - ``x`` -- the real number from which we want the continued fraction
     1352
     1353        TESTS::
     1354
     1355            sage: TestSuite(CFF(pi)).run()
     1356        """
     1357        ContinuedFraction_generic.__init__(self, parent)
     1358        self._x0 = x
     1359
     1360        self._xa = RealIntervalField(53)(self._x0)   # an approximation of the
     1361                                                     # last element of the orbit
     1362                                                     # under the Gauss map
     1363        self._quotients = []
     1364
     1365    def length(self):
     1366        r"""
     1367        Return infinity
     1368
     1369        EXAMPLES::
     1370
     1371            sage: CFF(pi).length()
     1372            +Infinity
     1373        """
     1374        return Infinity
     1375
     1376    def __len__(self):
     1377        r"""
     1378        TESTS::
     1379
     1380            sage: len(CFF(pi))
     1381            Traceback (most recent call last):
     1382            ...
     1383            ValueError: the length is infinite!
     1384        """
     1385        raise ValueError("the length is infinite!")
     1386
     1387    def __reduce__(self):
     1388        r"""
     1389        Pickling support.
     1390
     1391        TESTS::
     1392
     1393            sage: cf = continued_fraction(pi)
     1394            sage: loads(dumps(cf)) == cf
     1395            True
     1396        """
     1397        return self.parent(),(self.value(),)
     1398
     1399    def __cmp__(self, other):
     1400        r"""
     1401        Comparison.
     1402
     1403        EXAMPLES::
     1404
     1405            sage: CFF(pi) > CFF(e)
     1406            True
     1407            sage: CFF(pi) > CFF(e+4)
     1408            False
     1409        """
     1410        try:
     1411            # The following is crazy and prevent us from using cmp(self.value(),
     1412            # other.value()). On sage-5.10.beta2:
     1413            #     sage: cmp(pi, 4)
     1414            #     -1
     1415            #     sage: cmp(pi, pi+4)
     1416            #     1
     1417            if self.value() == other.value():
     1418                return 0
     1419            if self.value() - other.value() > 0:
     1420                return 1
     1421            return -1
     1422        except StandardError:
     1423            return ContinuedFraction_generic.__cmp__(self, other)
     1424
     1425    def _repr_(self):
     1426        r"""
     1427        String representation.
     1428
     1429        EXAMPLES::
     1430
     1431            sage: continued_fraction(pi) # indirect doctest
     1432            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     1433        """
     1434        return '[%d; '%self.quotient(0) + ', '.join(str(self.quotient(i)) for i in xrange(1,20)) + ", ...]"
     1435
     1436    def quotient(self, n):
     1437        r"""
     1438        Returns the ``n``-th quotient of ``self``.
     1439
     1440        EXAMPLES::
     1441
     1442            sage: cf = continued_fraction(pi)
     1443            sage: cf.quotient(27)
     1444            13
     1445            sage: cf.quotient(2552)
     1446            152
     1447            sage: cf.quotient(10000)   # long time
     1448            5
     1449
     1450        The algorithm is not efficient with element of the symbolic ring and,
     1451        if possible, one can always prefer number fields elements. The reason is
     1452        that, given a symbolic element ``x``, there is no automatic way to
     1453        evaluate in ``RIF`` an expression of the form ``(a*x+b)/(c*x+d)`` where
     1454        both the numerator and the denominator are extremely small::
     1455
     1456            sage: a1 = pi
     1457            sage: c1 = continued_fraction(a1)
     1458            sage: p0 = c1.pn(12); q0 = c1.qn(12)
     1459            sage: p1 = c1.pn(13); q1 = c1.qn(13)
     1460            sage: num = (q0*a1 - p0); num.n()
     1461            1.49011611938477e-8
     1462            sage: den = (q1*a1 - p1); den.n()
     1463            -2.98023223876953e-8
     1464            sage: a1 = -num/den
     1465            sage: RIF(a1)
     1466            [-infinity .. +infinity]
     1467
     1468        The same computation with an element of a number field instead of
     1469        ``pi`` gives a very satisfactory answer::
     1470
     1471            sage: K.<a2> = NumberField(x^3 - 2, embedding=1.25)
     1472            sage: c2 = continued_fraction(a2)
     1473            sage: p0 = c2.pn(111); q0 = c2.qn(111)
     1474            sage: p1 = c2.pn(112); q1 = c2.qn(112)
     1475            sage: num = (q0*a2 - p0); num.n()
     1476            -4.56719261665907e46
     1477            sage: den = (q1*a2 - p1); den.n()
     1478            -3.65375409332726e47
     1479            sage: a2 = -num/den
     1480            sage: b2 = RIF(a2); b2
     1481            1.002685823312715?
     1482            sage: b2.absolute_diameter()
     1483            8.88178419700125e-16
     1484
     1485        The consequence is that the precision needed with ``c1`` grows when we
     1486        compute larger and larger partial quotients::
     1487
     1488            sage: c1.quotient(100)
     1489            2
     1490            sage: c1._xa.parent()
     1491            Real Interval Field with 353 bits of precision
     1492            sage: c1.quotient(200)
     1493            3
     1494            sage: c1._xa.parent()
     1495            Real Interval Field with 753 bits of precision
     1496            sage: c1.quotient(300)
     1497            5
     1498            sage: c1._xa.parent()
     1499            Real Interval Field with 1053 bits of precision
     1500
     1501            sage: c2.quotient(200)
     1502            6
     1503            sage: c2._xa.parent()
     1504            Real Interval Field with 53 bits of precision
     1505            sage: c2.quotient(500)
     1506            1
     1507            sage: c2._xa.parent()
     1508            Real Interval Field with 53 bits of precision
     1509            sage: c2.quotient(1000)
     1510            1
     1511            sage: c2._xa.parent()
     1512            Real Interval Field with 53 bits of precision
     1513        """
     1514        x = self._xa
     1515        for k in xrange(len(self._quotients), n+1):
     1516            if x.lower().is_infinity() or x.upper().is_infinity() or x.lower().floor() != x.upper().floor():
     1517                orbit = lambda z: -(self.qn(k-2)*z-self.pn(k-2))/(self.qn(k-1)*z-self.pn(k-1))
     1518                x = x.parent()(orbit(self._x0))
     1519
     1520                # It may happen that the above line fails to give an
     1521                # approximation with the expected number of digits (see in the
     1522                # examples). In that case, we augment the precision.
     1523                while x.lower().is_infinity() or x.upper().is_infinity() or x.lower().floor() != x.upper().floor():
     1524                    self._prec = x.parent().prec() + 100
     1525                    x = RealIntervalField(self._prec)(orbit(self._x0))
     1526
     1527            self._quotients.append(x.unique_floor())
     1528            x = (x-x.unique_floor())
     1529            if not x:
     1530                raise ValueError("the number is rational")
     1531            x = ~x
     1532
     1533        self._xa = x
     1534        return self._quotients[n]
     1535
     1536    def value(self):
     1537        r"""
     1538        Return the value of ``self`` (the number from which it was built).
     1539
     1540        EXAMPLES::
     1541
     1542            sage: cf = continued_fraction(e)
     1543            sage: cf.value()
     1544            e
     1545        """
     1546        return self._x0
     1547
     1548def check_and_reduce_pair(x1,x2=None):
     1549    r"""
     1550    There are often two ways to represent a given continued fraction. This
     1551    function makes it canonical.
     1552
     1553    In the very special case of the number `0` we return the pair
     1554    ``((0,),(0,))``.
     1555
     1556    TESTS::
     1557
     1558        sage: from sage.rings.continued_fractions import check_and_reduce_pair
     1559        sage: check_and_reduce_pair([])
     1560        ((0,), (0,))
     1561        sage: check_and_reduce_pair([-1,1])
     1562        ((0,), (0,))
     1563        sage: check_and_reduce_pair([1,1,1])
     1564        ((1, 2), (0,))
     1565        sage: check_and_reduce_pair([1,3],[2,3])
     1566        ((1,), (3, 2))
     1567        sage: check_and_reduce_pair([1,2,3],[2,3,2,3,2,3])
     1568        ((1,), (2, 3))
     1569        sage: check_and_reduce_pair([0,0],[0,0,0])
     1570        ((0,), (0,))
     1571        sage: check_and_reduce_pair([1,2,0],[0])
     1572        ((1, 2), (0,))
     1573    """
     1574    y1 = map(ZZ,x1)
     1575
     1576    if x2 is None or (not x2) or all(a == 0 for a in x2):
     1577        y2 = [0]
     1578        if not y1:
     1579            y1 = [ZZ_0]
     1580        elif len(y1) > 1 and y1[-1] == 1:
     1581            y1.pop(-1)
     1582            y1[-1] += 1
     1583
     1584    else:
     1585        y2 = map(ZZ,x2)
     1586        if any(b <= ZZ_0 for b in y2):
     1587            raise ValueError("the elements of the period can not be negative")
     1588
     1589    # add possibly some element of x1 into the period
     1590    while y1 and y1[-1] == y2[-1]:
     1591        y1.pop(-1)
     1592        y2.insert(0,y2.pop(-1))
     1593
     1594    # some special cases to treat
     1595    if len(y2) == 1 and y2[0] == 0:
     1596        if not y1:
     1597            y1 = [ZZ_0]
     1598        elif len(y1) > 1 and y1[-1] == 1:
     1599            y1.pop(-1)
     1600            y1[-1] += 1
     1601
     1602    # check that y2 is not a pure power (in a very naive way!!)
     1603    n2 = len(y2)
     1604    for i in xrange(1,(n2+2)/2):
     1605        if n2%i == 0 and y2[:-i] == y2[i:]:
     1606            y2 = y2[:i]
     1607            break
     1608
     1609    # check that at then end y1 has no zeros in it
     1610    for i in xrange(1,len(y1)):
     1611        if y1[i] <= 0:
     1612            raise ValueError("all quotient except the first must be positive")
     1613
     1614    return tuple(y1),tuple(y2)
     1615
     1616class ContinuedFractionField(UniqueRepresentation, Field):
     1617    """
     1618    A common parent for continued fractions.
     1619
     1620    The continued fraction is a field isomorphic to `\RR`. Nevertheless, it is
     1621    not a good idea to use continued fraction to make calculus: operations are
     1622    very slow.
     1623
     1624    .. SEEALSO::
     1625
     1626        :func:`continued_fraction`
     1627
     1628    EXAMPLES::
     1629
    811630        sage: CFF
    82         Field of all continued fractions
    83 
    84     The continued fraction field inherits from the base class
    85     :class:`sage.rings.ring.Field`. However it was initialised
    86     as such only since trac ticket #11900::
     1631        Real field with infinite precision (continued fractions)
     1632        sage: CFF([0,1,3,2])
     1633        [0; 1, 3, 2]
     1634        sage: CFF(133/25)
     1635        [5; 3, 8]
     1636        sage: CFF(cos(pi/12))
     1637        [0; 1, 28, 2, 1, 7, 21, 1, 8, 1, 3, 10, 1, 2, 3, 1, 8, 2, 1, 2, ...]
    871638
    881639        sage: CFF.category()
    891640        Category of fields
    90 
    911641    """
    921642    def __init__(self):
    93         """
    94         EXAMPLES::
    95        
    96             sage: ContinuedFractionField()
    97             Field of all continued fractions
     1643        r"""
     1644        TESTS::
    981645
    99         TESTS::
    100        
    101             sage: CFF._repr_option('element_is_atomic')
    102             False
     1646            sage: TestSuite(CFF).run()
    1031647        """
    1041648        Field.__init__(self, self)
    105         self._assign_names(('x'),normalize=False)
    1061649
    107     def __cmp__(self, right):
    108         """
    109         EXAMPLES::
    110        
    111             sage: CFF == ContinuedFractionField()
    112             True
    113             sage: CFF == CDF
    114             False
    115             sage: loads(dumps(CFF)) == CFF
    116             True
    117         """
    118         return cmp(type(self), type(right))
     1650        # build the main element classes
     1651        self.element_class = self._element_class_periodic = self.__make_element_class__(ContinuedFraction_periodic)
     1652        self._element_class_number = self.__make_element_class__(ContinuedFraction_irrational)
    1191653
    120     def __iter__(self):
    121         """
    122         EXAMPLES::
    123        
    124             sage: i = 0
    125             sage: for a in CFF:
    126             ...    print a
    127             ...    i += 1
    128             ...    if i > 5: break
    129             ...
    130             [0]
    131             [1]
    132             [-1]
    133             [0, 2]
    134             [-1, 2]
    135             [2]       
    136         """
    137         for n in QQ:
    138             yield self(n)
    139 
     1654        self._populate_coercion_lists_(convert_method_name='continued_fraction')
    1401655
    1411656    def _latex_(self):
    1421657        r"""
    1431658        EXAMPLES::
    144        
    145             sage: latex(CFF) 
     1659
     1660            sage: latex(CFF)
    1461661            \Bold{CFF}
    1471662        """
    1481663        return "\\Bold{CFF}"
    149    
    150     def _is_valid_homomorphism_(self, codomain, im_gens):
     1664
     1665    def _repr_(self):
    1511666        """
    152         Return whether or not the map to codomain by sending the
    153         continued fraction [1] of self to im_gens[0] is a
    154         homomorphism.
     1667        EXAMPLES::
     1668
     1669            sage: CFF
     1670            Real field with infinite precision (continued fractions)
     1671        """
     1672        return "Real field with infinite precision (continued fractions)"
     1673
     1674    def an_element(self):
     1675        r"""
     1676        Returns a continued fraction.
    1551677
    1561678        EXAMPLES::
    1571679
    158             sage: CFF._is_valid_homomorphism_(ZZ,[ZZ(1)])
    159             False
    160             sage: CFF._is_valid_homomorphism_(CFF,[CFF(1)])
    161             True
     1680            sage: CFF.an_element()
     1681            [-1; 2, 3]
    1621682        """
    163         try:
    164             return im_gens[0] == codomain._coerce_(self(1))
    165         except TypeError:
    166             return False
    167        
    168     def _repr_(self):
     1683        return self([-1,2,3])
     1684
     1685    def _element_constructor_(self, x, bits=None, nterms=None, check=True, value=None):
    1691686        """
    170         EXAMPLES::
    171        
    172             sage: CFF
    173             Field of all continued fractions
    174         """
    175         return "Field of all continued fractions"
     1687        Construct a continued fraction from the given data.
    1761688
    177     def _coerce_impl(self, x):
    178         """
    179         Anything that implicitly coerces to the rationals or a real
    180         field, implicitly coerces to the continued fraction field.
    181        
    182         EXAMPLES:
    183        
    184         The additions below call _coerce_impl implicitly::
    185        
    186             sage: a = CFF(3/5); a
    187             [0, 1, 1, 2]
    188             sage: a + 2/5
    189             [1]
    190             sage: 2/5 + a
    191             [1]
    192             sage: 1.5 + a
    193             [2, 10]
    194             sage: a  + 1.5
    195             [2, 10]
    196         """
    197         if is_RealNumber(x):
    198             return self(x)
    199         return self._coerce_try(x, [QQ, RDF])
    200 
    201     def __call__(self, x, bits=None, nterms=None):
    202         """
    2031689        INPUT:
    2041690
    2051691            - `x` -- a number
    206            
     1692
    2071693            - ``bits`` -- integer (optional) the number of bits of the
    208               input number to use when computing the continued fraction.
    209            
    210         EXAMPLES::
    211        
    212             sage: CFF(1.5)
    213             [1, 2]
    214             sage: CFF(e)
    215             [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1]
    216             sage: CFF(pi)
    217             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    218             sage: CFF([1,2,3])
    219             [1, 2, 3]
    220             sage: CFF(15/17)
    221             [0, 1, 7, 2]
    222             sage: c2 = loads(dumps(CFF))
    223             sage: c2(CFF(15/17)).parent() is c2
    224             True
     1694              input number to use when computing with continued fraction.
    2251695
    226         We illustrate varying the bits parameter::
    227        
    228             sage: CFF(pi)
    229             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    230             sage: CFF(pi, bits=20)
    231             [3, 7]
    232             sage: CFF(pi, bits=80)
    233             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1]
    234             sage: CFF(pi, bits=100)
    235             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3]
    236        
    237         And varying the nterms parameter::
    238            
    239             sage: CFF(pi, nterms=3)
    240             [3, 7, 15]
    241             sage: CFF(pi, nterms=10)
    242             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1]
    243             sage: CFF(pi, bits=10, nterms=10)
    244             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1]
    245         """
    246         return ContinuedFraction(self, x, bits, nterms)
     1696            - ``nterms`` -- integer (optional)
    2471697
    248     def __len__(self):
    249         """
    250         EXAMPLES::
    251        
    252             sage: len(CFF)
    253             Traceback (most recent call last):
    254             ...
    255             TypeError: len() of unsized object
    256         """
    257         raise TypeError, 'len() of unsized object'
    258 
    259     def gens(self):
    260         """
    261         EXAMPLES::
    262        
    263             sage: CFF.gens()
    264             ([1],)
    265         """
    266         return (self(1), )
    267 
    268     def gen(self, n=0):
    269         """
    270         EXAMPLES::
    271        
    272             sage: CFF.gen()
    273             [1]
    274             sage: CFF.0
    275             [1]       
    276         """
    277        
    278         if n == 0:
    279             return self(1)
    280         else:
    281             raise IndexError, "n must be 0"
    282 
    283     def degree(self):
    284         """
    285         EXAMPLES::
    286        
    287             sage: CFF.degree()
    288             1
    289         """
    290         return 1
    291 
    292     def ngens(self):
    293         """
    294         EXAMPLES::
    295        
    296             sage: CFF.ngens()
    297             1
    298         """
    299         return 1
    300 
    301     def is_field(self, proof = True):
    302         """
    303         Return True, since the continued fraction field is a field.
     1698            - ``check`` -- boolean (optional) -- whether or not we check the
    3041699
    3051700        EXAMPLES::
    306        
    307             sage: CFF.is_field()
     1701
     1702            sage: CFF(1.5)
     1703            [1; 2]
     1704            sage: CFF(23/17)
     1705            [1; 2, 1, 5]
     1706
     1707            sage: continued_fraction(e)
     1708            [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1, ...]
     1709            sage: continued_fraction(pi)
     1710            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2,
     1711            ...]
     1712
     1713            sage: CFF([1,2,3])
     1714            [1; 2, 3]
     1715            sage: CFF(15/17)
     1716            [0; 1, 7, 2]
     1717
     1718            sage: CFF(pi, bits=20)
     1719            [3; 7]
     1720            sage: CFF(pi, bits=80)
     1721            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 3]
     1722            sage: CFF(pi, bits=100)
     1723            [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3]
     1724            sage: CFF(pi, nterms=3)
     1725            [3; 7, 15]
     1726            sage: CFF(pi, nterms=10)
     1727            [3; 7, 15, 1, 292, 1, 1, 1, 3]
     1728        """
     1729        if isinstance(x, Element) and x.parent() is self:
     1730            return x
     1731
     1732        if isinstance(x, (list,tuple)): # a digit expansion
     1733            cls = self._element_class_periodic
     1734
     1735            if len(x) == 2 and isinstance(x[0], (list,tuple)) and isinstance(x[1], (list,tuple)):
     1736                x1 = tuple(ZZ(a) for a in x[0])
     1737                x2 = tuple(ZZ(a) for a in x[1])
     1738                x1, x2 = check_and_reduce_pair(x1, x2)
     1739            else:
     1740                x1, x2 = check_and_reduce_pair(x)
     1741            return self._element_class_periodic(self, x1, x2)
     1742
     1743        else: # a number
     1744            if x in QQ: # rational
     1745                return QQ(x).continued_fraction()
     1746            if bits is None and nterms is None and isinstance(x, Element): # an irrational
     1747                from sage.symbolic.ring import SR
     1748                if x.parent() is SR:
     1749                    # we try a conversion to the algebraic real field (in order
     1750                    # to make more accurate computations and detect purely
     1751                    # periodic expansion).
     1752                    from sage.rings.qqbar import AA
     1753                    try:
     1754                        x = AA(x)
     1755                    except StandardError:
     1756                        pass
     1757                    else:
     1758                        K,a,_ = x.as_number_field_element(minimal=True)
     1759                        # actually, the method as_number_field_element is not
     1760                        # well implemented as it is not initialized with a
     1761                        # coerce embedding...
     1762                        try:
     1763                            return a.continued_fraction()
     1764                        except AttributeError:
     1765                            pass
     1766
     1767                if x.parent() is SR or x.parent().is_exact():
     1768                    # TODO: even if ``x`` belongs to the symbolic ring, we don't
     1769                    # know if it is an exact real number (like ``pi`` or ``e``)
     1770                    # or if it is a floating point approximation...
     1771                    from sage.rings.real_mpfi import RIF
     1772                    try:
     1773                        RIF(x)
     1774                    except StandardError:
     1775                        raise ValueError("the number %s does not seem to be a real number"%x)
     1776                    return self._element_class_number(self, x)
     1777
     1778            # now work with RIF
     1779            if bits is None:
     1780                try:
     1781                    bits = x.parent().prec()
     1782                except AttributeError:
     1783                    bits = 53
     1784            x = RealIntervalField(bits)(x)
     1785            cf = []
     1786            while True:
     1787                try:
     1788                    a = x.unique_floor()
     1789                except ValueError:
     1790                    x1, x2 = check_and_reduce_pair(cf)
     1791                    cf = self._element_class_periodic(self, x1, x2)
     1792                    break
     1793                cf.append(a)
     1794                x = ~(x - a)
     1795
     1796            if nterms is None:
     1797                return cf
     1798            return self([cf.quotient(i) for i in xrange(nterms)])
     1799
     1800    def is_field(self, proof=True):
     1801        """
     1802        Return True.
     1803
     1804        EXAMPLES::
     1805
     1806            sage: CFF.is_field()
     1807            True
     1808        """
     1809        return True
     1810
     1811    def is_exact(self):
     1812        r"""
     1813        Return True.
     1814
     1815        EXAMPLES::
     1816
     1817            sage: CFF.is_exact()
    3081818            True
    3091819        """
    3101820        return True
    3111821
    3121822    def is_finite(self):
    3131823        """
    314         Return False, since the continued fraction field is not finite. 
     1824        Return False, since the continued fraction field is not finite.
    3151825
    3161826        EXAMPLES::
    317        
    318             sage: CFF.is_finite() 
     1827
     1828            sage: CFF.is_finite()
    3191829            False
    3201830        """
    3211831        return False
    class ContinuedFractionField_class(Field 
    3251835        Return 0, since the continued fraction field has characteristic 0.
    3261836
    3271837        EXAMPLES::
    328        
     1838
    3291839            sage: c = CFF.characteristic(); c
    3301840            0
    3311841            sage: parent(c)
    3321842            Integer Ring
    3331843        """
    334         return ZZ(0)
    335    
     1844        return ZZ_0
     1845
    3361846    def order(self):
    3371847        """
    3381848        EXAMPLES::
    339        
     1849
    3401850            sage: CFF.order()
    3411851            +Infinity
    3421852        """
    343         return infinity
     1853        return Infinity
    3441854
    345     def random_element(self, *args, **kwds):
     1855    def random_element(self, first_digit=0, preperiod=True, period=False,
     1856                       x=1, y=10, distribution=None):
    3461857        """
    347         EXAMPLES::
    348        
    349            sage: CFF.random_element(10,10)
    350            [0, 4]
     1858        Return a somewhat random continued fraction (the result is either
     1859        finite or ultimately periodic).
    3511860
    352         Passes extra positional or keyword arguments through::
    353        
    354            sage: [CFF.random_element(den_bound=10, num_bound=2) for x in range(4)]
    355            [[-1, 1, 3], [0, 7], [0, 3], [0, 4]]
     1861        INPUT:
    3561862
     1863        - ``first_digit`` -- an optional value for the first digit (default
     1864          ``0``)
    3571865
    358            
    359         """
    360         return self(QQ.random_element(*args, **kwds))
    361    
     1866        - ``preperiod`` -- boolean -- whether or not have a preperiod (default
     1867          is ``True``)
    3621868
    363 class ContinuedFraction(FieldElement):
    364     """
    365     A continued fraction object.
     1869        - ``period`` -- boolean -- whether or not have a period (default is
     1870          ``False``)
    3661871
    367     EXAMPLES::
    368    
    369         sage: continued_fraction(pi)
    370         [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    371         sage: CFF(pi)
    372         [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    373     """
    374     def __init__(self, parent, x, bits=None, nterms=None):
    375         """
     1872        - ``x``, ``y`` -- optional bounds for the digits (default are ``1`` and
     1873          ``10``). Note that these options are sent to ``ZZ.random_element``.
     1874
     1875        - ``distribution`` -- an optional distribution for the digits
     1876
    3761877        EXAMPLES::
    3771878
    378             sage: sage.rings.contfrac.ContinuedFraction(CFF,[1,2,3,4,1,2])
    379             [1, 2, 3, 4, 1, 2]
    380             sage: sage.rings.contfrac.ContinuedFraction(CFF,[1,2,3,4,-1,2])
    381             Traceback (most recent call last):
    382             ...
    383             ValueError: each entry except the first must be positive
     1879            sage: CFF.random_element()
     1880            [0; 4, 7]
     1881            sage: CFF.random_element()
     1882            [0; 7, 6]
     1883
     1884            sage: CFF.random_element(preperiod=False, period=True)
     1885            [0; (6)*]
     1886            sage: CFF.random_element(preperiod=False, period=True)
     1887            [0; (1, 3, 4, 5, 8)*]
    3841888        """
    385         FieldElement.__init__(self, parent)
    386         if isinstance(x, ContinuedFraction):
    387             self._x = list(x._x)
    388         elif isinstance(x, (list, tuple)):
    389             x = [ZZ(a) for a in x]
    390             for i in range(1,len(x)):
    391                 if x[i] <= 0:
    392                     raise ValueError, "each entry except the first must be positive"
    393             self._x = list(x)
    394         else:
    395             self._x = [ZZ(a) for a in continued_fraction_list(x, bits=bits, nterms=nterms)]
     1889        from sage.misc.prandom import random
     1890        x1 = [ZZ(first_digit)]
     1891        x2 = []
     1892        if preperiod:
     1893            if random() > 0.1:
     1894                x1.append(ZZ.random_element(x,y,distribution))
     1895                while random() > 0.5: # from here they may be arbitrarily many
     1896                                      # elements but the mean length is 2
     1897                    x1.append(ZZ.random_element(x,y,distribution))
     1898        if period:
     1899            x2.append(ZZ.random_element(x,y,distribution))
     1900            while random() > 0.2:
     1901                x2.append(ZZ.random_element(x,y,distribution))
    3961902
    397     def __getitem__(self, n):
    398         """
    399         Returns `n`-th term of the continued fraction.
     1903        x1,x2 = check_and_reduce_pair(x1,x2)
     1904        return self._element_class_periodic(self, x1, x2)
    4001905
    401         OUTPUT:
    402             - an integer or a a continued fraction
     1906CFF = ContinuedFractionField()
    4031907
    404         EXAMPLES::
    405        
    406             sage: a = continued_fraction(pi); a
    407             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    408             sage: a[4]
    409             292
    410             sage: a[-1]
    411             14
    412             sage: a[2:5]
    413             [15, 1, 292]
    414             sage: a[:3]
    415             [3, 7, 15]
    416             sage: a[4:]
    417             [292, 1, 1, 1, 2, 1, 3, 1, 14]
    418             sage: a[4::2]
    419             [292, 1, 2, 3, 14]
    420         """
    421         if isinstance(n, slice):
    422             start, stop, step = n.indices(len(self))
    423             return ContinuedFraction(self.parent(), self._x[start:stop:step])
    424         else:
    425             return self._x[n]
     1908def continued_fraction_list(x, type="std", partial_convergents=False, bits=None, nterms=None):
     1909    r"""
     1910    Returns the (finite) continued fraction of ``x`` as a list.
    4261911
    427     def _repr_(self):
    428         """
    429         EXAMPLES::
    430        
    431             sage: a = continued_fraction(pi); a
    432             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    433             sage: a.rename('continued fraction of pi')
    434             sage: a
    435             continued fraction of pi       
    436         """
    437         return str(self._x)
     1912    The continued fraction expansion of ``x`` are the coefficients `a_i` in
    4381913
    439     def convergents(self):
    440         """
    441         Return a list of rational numbers, which are the partial
    442         convergents of this continued fraction.
    443        
    444         OUTPUT:
     1914    .. MATH::
    4451915
    446             - list of rational numbers
     1916        x = a_0 + 1/(a_1 + 1/(...))
    4471917
    448         EXAMPLES::
    449        
    450             sage: a = CFF(pi, bits=34); a
    451             [3, 7, 15, 1, 292]
    452             sage: a.convergents()
    453             [3, 22/7, 333/106, 355/113, 103993/33102]
    454             sage: a.value()
    455             103993/33102
    456             sage: a[:-1].value()
    457             355/113       
    458         """
    459         return convergents(self._x)
     1918    with `a_0` integer and `a_1`, `...` positive integers. The Hirzebruch-Jung
     1919    continued fraction is the one for which the `+` signs are replaced with `-`
     1920    signs
    4601921
    461     def convergent(self, n):
    462         """
    463         Return the `n`-th partial convergent to self.
     1922    .. MATH::
    4641923
    465         OUTPUT:
    466        
    467             rational number
    468            
    469         EXAMPLES::
    470        
    471             sage: a = CFF(pi, bits=34); a
    472             [3, 7, 15, 1, 292]
    473             sage: a.convergents()
    474             [3, 22/7, 333/106, 355/113, 103993/33102]
    475             sage: a.convergent(0)
    476             3
    477             sage: a.convergent(1)
    478             22/7
    479             sage: a.convergent(4)
    480             103993/33102       
    481         """
    482         return convergent(self._x, n)
     1924        x = a_0 - 1/(a_1 - 1/(...))
    4831925
    484     def __len__(self):
    485         """
    486         Return the number of terms in this continued fraction.
     1926    .. SEEALSO::
    4871927
    488         EXAMPLES::
    489        
    490             sage: len(continued_fraction([1,2,3,4,5]) )
    491             5       
    492         """
    493         return len(self._x)
    494 
    495     def pn(self, n):
    496         """
    497         Return the numerator of the `n`-th partial convergent, computed
    498         using the recurrence.
    499        
    500         EXAMPLES::
    501        
    502             sage: c = continued_fraction(pi); c
    503             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    504             sage: c.pn(0), c.qn(0)
    505             (3, 1)
    506             sage: len(c)
    507             13
    508             sage: c.pn(12), c.qn(12)
    509             (80143857, 25510582)
    510         """
    511         if n < -2:
    512             raise ValueError, "n must be at least -2"
    513         if n > len(self._x):
    514             raise ValueError, "n must be at most %s"%len(self._x)
    515         try:
    516             return self.__pn[n+2]
    517         except AttributeError:
    518             self.__pn = [0, 1, self._x[0]]
    519             self.__qn = [1, 0, 1]
    520         except IndexError:
    521             pass
    522         for k in range(len(self.__pn), n+3):
    523             self.__pn.append(self._x[k-2]*self.__pn[k-1] + self.__pn[k-2])
    524             self.__qn.append(self._x[k-2]*self.__qn[k-1] + self.__qn[k-2])
    525         return self.__pn[n+2]
    526            
    527     def qn(self, n):
    528         """
    529         Return the denominator of the `n`-th partial convergent, computed
    530         using the recurrence.
    531        
    532         EXAMPLES::
    533        
    534             sage: c = continued_fraction(pi); c
    535             [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    536             sage: c.qn(0), c.pn(0)
    537             (1, 3)
    538             sage: len(c)
    539             13
    540             sage: c.pn(12), c.qn(12)
    541             (80143857, 25510582)
    542         """
    543         if n < -2:
    544             raise ValueError, "n must be at least -2"
    545         if n > len(self._x):
    546             raise ValueError, "n must be at most %s"%len(self._x)
    547         try:
    548             return self.__qn[n+2]
    549         except (AttributeError, IndexError):
    550             pass
    551         self.pn(n)
    552         return self.__qn[n+2]
    553 
    554     def _rational_(self):
    555         """
    556         EXAMPLES::
    557        
    558             sage: a = CFF(-17/389); a
    559             [-1, 1, 21, 1, 7, 2]
    560             sage: a._rational_()
    561             -17/389
    562             sage: QQ(a)
    563             -17/389
    564         """
    565         try:
    566             return self.__rational
    567         except AttributeError:
    568             r = convergents(self._x)[-1]
    569             self.__rational =r
    570             return r
    571 
    572     def value(self):
    573         """
    574         EXAMPLES::
    575        
    576             sage: a = CFF(-17/389); a
    577             [-1, 1, 21, 1, 7, 2]
    578             sage: a.value()
    579             -17/389
    580             sage: QQ(a)
    581             -17/389
    582         """
    583         return self._rational_()
    584 
    585     def numerator(self):
    586         """
    587         EXAMPLES::
    588        
    589             sage: a = CFF(-17/389); a
    590             [-1, 1, 21, 1, 7, 2]
    591             sage: a.numerator()
    592             -17
    593         """
    594         return self._rational_().numerator()
    595 
    596     def denominator(self):
    597         """
    598         EXAMPLES::
    599        
    600             sage: a = CFF(-17/389); a
    601             [-1, 1, 21, 1, 7, 2]
    602             sage: a.denominator()
    603             389
    604         """
    605         return self._rational_().denominator()
    606 
    607     def __int__(self):
    608         """
    609         EXAMPLES::
    610        
    611             sage: a = CFF(-17/389); a
    612             [-1, 1, 21, 1, 7, 2]
    613             sage: int(a)
    614             -1
    615         """
    616         return int(self._rational_())
    617 
    618     def __long__(self):
    619         """
    620         EXAMPLES::
    621        
    622             sage: a = CFF(-17/389); a
    623             [-1, 1, 21, 1, 7, 2]
    624             sage: long(a)
    625             -1L
    626         """
    627         return long(self._rational_())
    628 
    629     def __float__(self):
    630         """
    631         EXAMPLES::
    632        
    633             sage: a = CFF(-17/389); a
    634             [-1, 1, 21, 1, 7, 2]
    635             sage: float(a)
    636             -0.043701799485861184
    637         """
    638         return float(self._rational_())
    639 
    640     def _add_(self, right):
    641         """
    642         EXAMPLES::
    643        
    644             sage: a = CFF(-17/389)
    645             sage: b = CFF(1/389)
    646             sage: c = a+b; c
    647             [-1, 1, 23, 3, 5]
    648             sage: c.value()
    649             -16/389
    650         """
    651         return ContinuedFraction(self.parent(),
    652                      self._rational_() + right._rational_())
    653    
    654     def _sub_(self, right):
    655         """
    656         EXAMPLES::
    657        
    658             sage: a = CFF(-17/389)
    659             sage: b = CFF(1/389)
    660             sage: c = a - b; c
    661             [-1, 1, 20, 1, 1, 1, 1, 3]
    662             sage: c.value()
    663             -18/389
    664         """
    665         return ContinuedFraction(self.parent(),
    666                      self._rational_() - right._rational_())
    667    
    668     def _mul_(self, right):
    669         """
    670         EXAMPLES::
    671        
    672             sage: a = CFF(-17/389)
    673             sage: b = CFF(1/389)
    674             sage: c = a * b; c
    675             [-1, 1, 8900, 4, 4]
    676             sage: c.value(), (-1/389)*(17/389)
    677             (-17/151321, -17/151321)
    678         """
    679         return ContinuedFraction(self.parent(),
    680                      self._rational_() * right._rational_())
    681 
    682     def _div_(self, right):
    683         """
    684         EXAMPLES::
    685        
    686             sage: a = CFF(-17/389)
    687             sage: b = CFF(1/389)
    688             sage: c = a / b; c
    689             [-17]
    690             sage: c.value(), (17/389) / (-1/389)
    691             (-17, -17)
    692         """
    693         return ContinuedFraction(self.parent(),
    694                      self._rational_() / right._rational_())
    695 
    696     def __cmp__(self, right):
    697         """
    698         EXAMPLES::
    699        
    700             sage: a = CFF(-17/389)
    701             sage: b = CFF(1/389)
    702             sage: a < b
    703             True
    704             sage: QQ(a) < QQ(b)
    705             True
    706             sage: QQ(a)
    707             -17/389
    708             sage: QQ(b)
    709             1/389
    710         """
    711         return cmp(self._rational_(), right._rational_())
    712 
    713     def _latex_(self):
    714         """
    715         EXAMPLES::
    716        
    717             sage: a = CFF(-17/389)
    718             sage: latex(a)
    719             -1+ \frac{\displaystyle 1}{\displaystyle 1+ \frac{\displaystyle 1}{\displaystyle 21+ \frac{\displaystyle 1}{\displaystyle 1+ \frac{\displaystyle 1}{\displaystyle 7+ \frac{\displaystyle 1}{\displaystyle 2}}}}}
    720         """
    721         v = self._x
    722         if len(v) == 0:
    723             return '0'
    724         s = str(v[0])
    725         for i in range(1,len(v)):
    726             s += '+ \\frac{\\displaystyle 1}{\\displaystyle %s'%v[i]
    727         s += '}'*(len(v)-1)
    728         return s
    729 
    730     def sqrt(self, prec=53, all=False):
    731         """
    732         Return continued fraction approximation to square root of the
    733         value of this continued fraction.
    734 
    735         INPUT:
    736        
    737             - `prec` -- integer (default: 53) precision of square root
    738               that is approximated
    739 
    740             - `all` -- bool (default: False); if True, return all square
    741               roots of self, instead of just one.
    742        
    743         EXAMPLES::
    744        
    745             sage: a = CFF(4/19); a
    746             [0, 4, 1, 3]
    747             sage: b = a.sqrt(); b
    748             [0, 2, 5, 1, 1, 2, 1, 16, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1]
    749             sage: b.value()
    750             4508361/9825745
    751             sage: float(b.value()^2 - a)
    752             -5.451492525672688e-16
    753             sage: b = a.sqrt(prec=100); b
    754             [0, 2, 5, 1, 1, 2, 1, 16, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1, 16, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1, 16, 1, 2, 1, 1, 5]
    755             sage: b^2
    756             [0, 4, 1, 3, 7849253184229368265220252099, 1, 3]
    757             sage: a.sqrt(all=True)
    758             [[0, 2, 5, 1, 1, 2, 1, 16, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1],
    759              [-1, 1, 1, 5, 1, 1, 2, 1, 16, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1]]
    760             sage: a = CFF(4/25).sqrt(); a
    761             [0, 2, 2]
    762             sage: a.value()
    763             2/5
    764         """
    765         r = self._rational_()
    766         if r < 0:
    767             raise ValueError, "self must be positive"
    768         X = r.sqrt(all=all, prec=prec)
    769         if not all:
    770             return ContinuedFraction(self.parent(), X)
    771         else:
    772             return [ContinuedFraction(self.parent(), x) for x in X]
    773 
    774     def list(self):
    775         """
    776         Return copy of the underlying list of this continued fraction.
    777 
    778         EXAMPLES::
    779        
    780             sage: a = CFF(e); v = a.list(); v
    781             [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1]
    782             sage: v[0] = 5
    783             sage: a
    784             [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1]
    785         """
    786         return list(self._x)
    787 
    788     def __hash__(self):
    789         """
    790         Return hash of self, which is the same as the hash of the value
    791         of self, as a rational number.
    792 
    793         EXAMPLES::
    794        
    795             sage: a = CFF(e)
    796             sage: hash(a)
    797             19952398
    798             sage: hash(QQ(a))
    799             19952398
    800         """
    801         return hash(self._rational_())
    802 
    803     def __invert__(self):
    804         """
    805         Return the multiplicative inverse of self.
    806        
    807         EXAMPLES::
    808        
    809             sage: a = CFF(e)
    810             sage: b = ~a; b
    811             [0, 2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 2]
    812             sage: b*a
    813             [1]
    814         """
    815         return ContinuedFraction(self.parent(),
    816                      self._rational_().__invert__())
    817 
    818     def __pow__(self, n):
    819         """
    820         Return self to the power of `n`.
    821 
    822         EXAMPLES::
    823        
    824             sage: a = CFF([1,2,3]); a
    825             [1, 2, 3]
    826             sage: a^3
    827             [2, 1, 10, 1, 4, 1, 4]
    828             sage: QQ(a)^3 == QQ(a^3)
    829             True
    830             sage: a^(-3)
    831             [0, 2, 1, 10, 1, 4, 1, 4]
    832             sage: QQ(a)^(-3) == QQ(a^(-3))
    833             True
    834         """
    835         return ContinuedFraction(self.parent(),
    836                      self._rational_()**n)
    837 
    838     def __neg__(self):
    839         """
    840         Return additive inverse of self.
    841        
    842         EXAMPLES::
    843        
    844             sage: a = CFF(-17/389); a
    845             [-1, 1, 21, 1, 7, 2]
    846             sage: -a
    847             [0, 22, 1, 7, 2]
    848             sage: QQ(-a)
    849             17/389
    850         """
    851         return ContinuedFraction(self.parent(),
    852                      self._rational_().__neg__())
    853 
    854     def __abs__(self):
    855         """
    856         Return absolute value of self.
    857 
    858         EXAMPLES::
    859        
    860             sage: a = CFF(-17/389); a
    861             [-1, 1, 21, 1, 7, 2]
    862             sage: abs(a)
    863             [0, 22, 1, 7, 2]
    864             sage: QQ(abs(a))
    865             17/389       
    866         """
    867         return ContinuedFraction(self.parent(),
    868                      self._rational_().__abs__())
    869 
    870     def is_one(self):
    871         """
    872         Return True if self is one.
    873        
    874         EXAMPLES::
    875        
    876             sage: continued_fraction(1).is_one()
    877             True
    878             sage: continued_fraction(2).is_one() 
    879             False
    880         """
    881         return self._rational_().is_one()
    882 
    883     def __nonzero__(self):
    884         """
    885         Return False if self is zero.
    886        
    887         EXAMPLES::
    888        
    889             sage: continued_fraction(0).is_zero()
    890             True
    891             sage: continued_fraction(1).is_zero()
    892             False
    893         """
    894         return not self._rational_().is_zero()
    895 
    896     def _pari_(self):
    897         """
    898         Return PARI list corresponding to this continued fraction.
    899        
    900         EXAMPLES::
    901        
    902             sage: c = continued_fraction(0.12345); c
    903             [0, 8, 9, 1, 21, 1, 1]
    904             sage: pari(c)
    905             [0, 8, 9, 1, 21, 1, 1]
    906         """
    907         return pari(self._x)
    908 
    909     def _interface_init_(self, I=None):
    910         """
    911         Return list representation for other systems corresponding to
    912         this continued fraction.
    913        
    914         EXAMPLES::
    915        
    916             sage: c = continued_fraction(0.12345); c
    917             [0, 8, 9, 1, 21, 1, 1]
    918             sage: gp(c)
    919             [0, 8, 9, 1, 21, 1, 1]
    920             sage: gap(c)
    921             [ 0, 8, 9, 1, 21, 1, 1 ]
    922             sage: maxima(c)
    923             [0,8,9,1,21,1,1]
    924         """
    925         return str(self._x)
    926 
    927     def additive_order(self):
    928         """
    929         Return the additive order of this continued fraction,
    930         which we defined to be the additive order of its value.
    931        
    932         EXAMPLES::
    933 
    934             sage: CFF(-1).additive_order()
    935             +Infinity
    936             sage: CFF(0).additive_order()
    937             1
    938         """
    939         return self.value().additive_order()
    940 
    941     def multiplicative_order(self):
    942         """
    943         Return the multiplicative order of this continued fraction,
    944         which we defined to be the multiplicative order of its value.
    945        
    946         EXAMPLES::
    947 
    948             sage: CFF(-1).multiplicative_order()
    949             2
    950             sage: CFF(1).multiplicative_order()
    951             1
    952             sage: CFF(pi).multiplicative_order()
    953             +Infinity
    954         """
    955         return self.value().multiplicative_order()
    956 
    957    
    958 CFF = ContinuedFractionField_class()
    959 
    960 def ContinuedFractionField():
    961     """
    962     Return the (unique) field of all continued fractions.
    963 
    964     EXAMPLES::
    965    
    966         sage: ContinuedFractionField()
    967         Field of all continued fractions
    968     """
    969     return CFF
    970 
    971 def continued_fraction(x, bits=None, nterms=None):
    972     """
    973     Return the truncated continued fraction expansion of the real number
    974     `x`, computed with an interval floating point approximation of `x`
    975     to the given number of bits of precision. The returned continued
    976     fraction is a list-like object, with a value method and partial
    977     convergents method.
    978 
    979     If bits is not given, then use the number of valid bits of
    980     precision of `x`, if `x` is a floating point number, or 53 bits
    981     otherwise. If nterms is given, the precision is increased until
    982     the specified number of terms can be computed, if possible.
     1928        :func:`continued_fraction`
    9831929
    9841930    INPUT:
    9851931
    986         - `x` -- number
     1932    - ``x`` -- exact rational or floating-point number. The number to
     1933      compute the continued fraction of.
    9871934
    988         - ``bits`` -- None (default) or a positive integer
    989        
    990         - ``nterms`` -- None (default) or a positive integer
     1935    - ``type`` -- either "std" (default) for standard continued fractions or
     1936      "hj" for Hirzebruch-Jung ones.
     1937
     1938    - ``partial_convergents`` -- boolean. Whether to return the partial convergents.
     1939
     1940    - ``bits`` -- integer. the precision of the real interval field
     1941      that is used internally.
     1942
     1943    - ``nterms`` -- integer. The upper bound on the number of terms in
     1944      the continued fraction expansion to return.
    9911945
    9921946    OUTPUT:
    9931947
    994         - a continued fraction
     1948    A lits of integers, the coefficients in the continued fraction expansion of
     1949    ``x``. If ``partial_convergents`` is set to ``True``, then return a pair
     1950    containing the coefficient list and the partial convergents list is
     1951    returned.
    9951952
    996     EXAMPLES::
     1953     EXAMPLES::
    9971954
    998         sage: v = continued_fraction(sqrt(2)); v
    999         [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
    1000         sage: v = continued_fraction(sqrt(2), nterms=22); v
    1001         [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
    1002         sage: type(v)
    1003         <class 'sage.rings.contfrac.ContinuedFraction'>
    1004         sage: parent(v)
    1005         Field of all continued fractions
    1006         sage: v.value()
    1007         131836323/93222358
    1008         sage: RR(v.value()) == RR(sqrt(2))
    1009         True
    1010         sage: v.convergents()
    1011         [1, 3/2, 7/5, 17/12, 41/29, 99/70, 239/169,...131836323/93222358]
    1012         sage: [RR(x) for x in v.convergents()]
    1013         [1.00000000000000, 1.50000000000000, 1.40000000000000, 1.41666666666667, ...1.41421356237310]
    1014         sage: continued_fraction(sqrt(2), 10)
    1015         [1, 2, 2]
    1016         sage: v.numerator()
    1017         131836323
    1018         sage: v.denominator()
    1019         93222358
    1020         sage: [v.pn(i) for i in range(10)]
    1021         [1, 3, 7, 17, 41, 99, 239, 577, 1393, 3363]
    1022         sage: [v.qn(i) for i in range(10)]
    1023         [1, 2, 5, 12, 29, 70, 169, 408, 985, 2378]
     1955        sage: continued_fraction_list(45/19)
     1956        [2, 2, 1, 2, 2]
     1957        sage: 2 + 1/(2 + 1/(1 + 1/(2 + 1/2)))
     1958        45/19
    10241959
    1025     Here are some more examples::
    1026    
    1027         sage: continued_fraction(e, bits=20)
    1028         [2, 1, 2, 1, 1, 4, 1, 1]
    1029         sage: continued_fraction(e, bits=30)
    1030         [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8]
    1031         sage: continued_fraction(RealField(200)(e))
    1032         [2, 1, 2, 1, 1, 4, 1, 1, 6, ...36, 1, 1, 38, 1, 1]
    1033    
    1034     Initial rounding can result in incorrect trailing digits::
    1035    
     1960        sage: continued_fraction_list(45/19,type="hj")
     1961        [3, 2, 3, 2, 3]
     1962        sage: 3 - 1/(2 - 1/(3 - 1/(2 - 1/3)))
     1963        45/19
     1964
     1965    Specifying ``bits`` or ``nterms`` modify the length of the output::
     1966
     1967        sage: continued_fraction_list(e, bits=20)
     1968        [2, 1, 2, 1, 1, 4, 2]
     1969        sage: continued_fraction_list(sqrt(2)+sqrt(3), bits=30)
     1970        [3, 6, 1, 5, 7, 2]
     1971        sage: continued_fraction_list(pi, bits=53)
     1972        [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
     1973
     1974        sage: continued_fraction_list(log(3/2), nterms=15)
     1975        [0, 2, 2, 6, 1, 11, 2, 1, 2, 2, 1, 4, 3, 1, 1]
     1976        sage: continued_fraction_list(tan(sqrt(pi)), nterms=20)
     1977        [-5, 9, 4, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1, 2, 4, 3, 1, 63]
     1978
     1979    When the continued fraction is infinite (ie ``x`` is an irrational number)
     1980    and the parameters ``bits`` and ``nterms`` are not specified then a warning
     1981    is raised::
     1982
     1983        sage: continued_fraction_list(sqrt(2))
     1984        doctest:...: UserWarning: the continued fraction of sqrt(2) seems infinite, return only the first 20 terms
     1985        [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
     1986        sage: continued_fraction_list(sqrt(4/19))
     1987        doctest:...: UserWarning: the continued fraction of 2*sqrt(1/19) seems infinite, return only the first 20 terms
     1988        [0, 2, 5, 1, 1, 2, 1, 16, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1, 16]
     1989
     1990    An examples with the list of partial convergents::
     1991
     1992        sage: continued_fraction_list(RR(pi), partial_convergents=True)
     1993        ([3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3],
     1994         [(3, 1),
     1995          (22, 7),
     1996          (333, 106),
     1997          (355, 113),
     1998          (103993, 33102),
     1999          (104348, 33215),
     2000          (208341, 66317),
     2001          (312689, 99532),
     2002          (833719, 265381),
     2003          (1146408, 364913),
     2004          (4272943, 1360120),
     2005          (5419351, 1725033),
     2006          (80143857, 25510582),
     2007          (245850922, 78256779)])
     2008
     2009    TESTS::
     2010
     2011        sage: continued_fraction_list(1 + 10^-10, nterms=3)
     2012        [1, 10000000000]
     2013        sage: continued_fraction_list(1 + 10^-20 - e^-100, nterms=3)
     2014        [1, 100000000000000000000, 2688]
     2015        sage: continued_fraction_list(1 + 10^-20 - e^-100, nterms=5)
     2016        [1, 100000000000000000000, 2688, 8, 1]
     2017        sage: continued_fraction_list(1 + 10^-20 - e^-100, nterms=5)
     2018        [1, 100000000000000000000, 2688, 8, 1]
     2019    """
     2020    if type == "hj":
     2021        if bits is not None:
     2022            x = RealIntervalField(bits)(x).simplest_rational()
     2023        else:
     2024            try:
     2025                x = QQ(x)
     2026            except (TypeError,ValueError):
     2027                x = QQ(x.n())
     2028        l = x.continued_fraction_list(type="hj")
     2029        ## The C-code in sage.rings.rational is much more faster than the pure
     2030        ## Python below
     2031        ## v = []
     2032        ## while True:
     2033        ##    div, mod = divmod(x.numerator(), x.denominator())
     2034        ##    if mod == 0:
     2035        ##        v.append(div)
     2036        ##        break
     2037        ##    v.append(div+1)
     2038        ##    if nterms is not None and len(v) >= nterms:
     2039        ##        break
     2040        ##    x = 1/(div+1-x)
     2041        ## return v
     2042        if nterms is None:
     2043            return l
     2044        return l[:nterms]
     2045    if type != "std":
     2046        raise ValueError("type must be either \"std\" or \"hj\"")
     2047    if bits is not None:
     2048        x = RealIntervalField(bits)(x)
     2049    cf = CFF(x)
     2050    if nterms:
     2051        limit = min(cf.length(), nterms)
     2052    elif cf.length() != Infinity:
     2053        limit = cf.length()
     2054    else:
     2055        import warnings
     2056        warnings.warn("the continued fraction of %s seems infinite, return only the first 20 terms"%x)
     2057        limit = 20
     2058    if partial_convergents:
     2059        return [cf.quotient(i) for i in xrange(limit)], [(cf.pn(i),cf.qn(i)) for i in xrange(limit)]
     2060    return [cf.quotient(i) for i in xrange(limit)]
     2061
     2062def continued_fraction(x, bits=None, nterms=None):
     2063    r"""
     2064    Return the continued fraction of ``x``.
     2065
     2066    INPUT:
     2067
     2068        - `x` -- a number or a list of digits or two list of digits (preperiod
     2069          and period)
     2070
     2071        - ``bits`` -- None (default) or a positive integer
     2072
     2073        - ``nterms`` -- None (default) or a positive integer
     2074
     2075
     2076        - ``bits`` -- None (default) or a positive integer. If not ``None`` then
     2077          use an approximation of ``x`` with ``bits`` precision.
     2078
     2079        - ``nterms`` -- None (default) or a positive integer. If not ``None``
     2080          then return the continued fraction that consists only on the first
     2081          ``nterms`` of the continued fraction of ``x``.
     2082
     2083    EXAMPLES:
     2084
     2085    A continued fraction may be initialized by a number or by its digit
     2086    expansion::
     2087
     2088        sage: continued_fraction(12/571)
     2089        [0; 47, 1, 1, 2, 2]
     2090        sage: continued_fraction([3,2,1,4])
     2091        [3; 2, 1, 4]
     2092
     2093        sage: c = continued_fraction(golden_ratio); c
     2094        [(1)*]
     2095        sage: c.convergent(12)
     2096        377/233
     2097        sage: fibonacci(14)/fibonacci(13)
     2098        377/233
     2099
     2100        sage: c = continued_fraction(pi); c
     2101        [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     2102        sage: a = c.convergent(3); a
     2103        355/113
     2104        sage: a.n()
     2105        3.14159292035398
     2106        sage: pi.n()
     2107        3.14159265358979
     2108
     2109    When possible, it is adviced to use anything else but the symbolic ring.
     2110    Here we present an other way of dealing with the golden mean and the cube
     2111    root of 2::
     2112
     2113        sage: K.<sqrt5> = NumberField(x^2-5, embedding=2.23)
     2114        sage: my_golden_ratio = (1 + sqrt5)/2
     2115        sage: cf = continued_fraction((1+sqrt5)/2); cf
     2116        [(1)*]
     2117        sage: cf.convergent(12)
     2118        377/233
     2119        sage: cf.period()
     2120        (1,)
     2121
     2122        sage: K.<a> = NumberField(x^3-2, embedding=1.25)
     2123        sage: cf = continued_fraction(a); cf
     2124        [1; 3, 1, 5, 1, 1, 4, 1, 1, 8, 1, 14, 1, 10, 2, 1, 4, 12, 2, 3, ...]
     2125
     2126    It is possible to pass from quadratic numbers to their ultimately periodic
     2127    continued fraction expansion::
     2128
     2129        sage: c = continued_fraction(my_golden_ratio); c
     2130        [(1)*]
     2131        sage: c.preperiod()
     2132        ()
     2133        sage: c.period()
     2134        (1,)
     2135
     2136        sage: c = continued_fraction(2/3+sqrt5/5); c
     2137        [1; 8, (1, 3, 1, 1, 3, 9)*]
     2138        sage: c.preperiod()
     2139        (1, 8)
     2140        sage: c.period()
     2141        (1, 3, 1, 1, 3, 9)
     2142
     2143        sage: cf = continued_fraction([(1,1),(2,8)]); cf
     2144        [1; 1, (2, 8)*]
     2145        sage: cf.value()
     2146        2/11*sqrt5 + 14/11
     2147
     2148    Some well known mathematical constants have continued fractions with a lot
     2149    of structure::
     2150
     2151        sage: continued_fraction(e)
     2152        [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1, ...]
     2153        sage: continued_fraction(tan(1))
     2154        [1; 1, 1, 3, 1, 5, 1, 7, 1, 9, 1, 11, 1, 13, 1, 15, 1, 17, 1, 19, ...]
     2155        sage: continued_fraction(tanh(1))
     2156        [0; 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, ...]
     2157
     2158    But others do not::
     2159
     2160        sage: continued_fraction(pi)
     2161        [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     2162        sage: continued_fraction(euler_gamma)
     2163        [0; 1, 1, 2, 1, 2, 1, 4, 3, 13, 5, 1, 1, 8, 1, 2, 4, 1, 1, 40, ...]
     2164
     2165    Note that initial rounding can result in incorrect trailing digits::
     2166
    10362167        sage: continued_fraction(RealField(39)(e))
    1037         [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 2]
     2168        [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 2]
    10382169        sage: continued_fraction(RealIntervalField(39)(e))
    1039         [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10]
     2170        [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10]
    10402171    """
    10412172    return CFF(x, bits=bits, nterms=nterms)
    10422173
     2174def Hirzebruch_Jung_continued_fraction_list(x, bits=None, nterms=None):
     2175    r"""
     2176    Return the Hirzebruch-Jung continued fraction of ``x`` as a list.
    10432177
     2178    This function is deprecated since :trac:`14567`. See
     2179    :func:`continued_fraction_list` and the documentation therein.
    10442180
    1045        
     2181    INPUT:
     2182
     2183    - ``x`` -- exact rational or something that can be numerically
     2184      evaluated. The number to compute the continued fraction of.
     2185
     2186    - ``bits`` -- integer (default: the precision of ``x``). the
     2187      precision of the real interval field that is used
     2188      internally. This is only used if ``x`` is not an exact fraction.
     2189
     2190    - ``nterms`` -- integer (default: None). The upper bound on the
     2191      number of terms in the continued fraction expansion to return.
     2192    A lits of integers, the coefficients in the Hirzebruch-Jung continued
     2193    fraction expansion of ``x``.
     2194        sage: Hirzebruch_Jung_continued_fraction_list(17/11)
     2195        doctest:...: DeprecationWarning: Hirzebruch_Jung_continued_fraction_list(x) is replaced by
     2196            continued_fraction_list(x,type="hj")
     2197        or for rationals
     2198            x.continued_fraction_list(type="hj")
     2199        See http://trac.sagemath.org/14567 for details.
     2200        [2, 3, 2, 2, 2, 2]
     2201    """
     2202    from sage.misc.superseded import deprecation
     2203    deprecation(14567, 'Hirzebruch_Jung_continued_fraction_list(x) is replaced by\n\tcontinued_fraction_list(x,type="hj")\nor for rationals\n\tx.continued_fraction_list(type="hj")')
     2204    return continued_fraction_list(x, type="hj", bits=bits, nterms=nterms)
     2205
     2206# Unpickling support is needed as the file sage.rings.contfrac is renamed
     2207# sage.rings.continued_fractions in #14567
     2208from sage.structure.sage_object import register_unpickle_override
     2209register_unpickle_override('sage.rings.contfrac','ContinuedFractionField_class',ContinuedFractionField)
  • sage/rings/number_field/number_field_element_quadratic.pyx

    diff --git a/sage/rings/number_field/number_field_element_quadratic.pyx b/sage/rings/number_field/number_field_element_quadratic.pyx
    a b cdef class NumberFieldElement_quadratic( 
    896896            return 1
    897897        return -1
    898898
     899    def continued_fraction_list(self):
     900        r"""
     901        Return the preperiod and the period of the continued fraction expansion
     902        of ``self``.
     903
     904        EXAMPLES::
     905
     906            sage: K.<sqrt2> = QuadraticField(2)
     907            sage: sqrt2.continued_fraction_list()
     908            ((1,), (2,))
     909            sage: (1/2+sqrt2/3).continued_fraction_list()
     910            ((0, 1, 33), (1, 32))
     911
     912        Check that it works even for rational entries::
     913
     914            sage: K(123/567).continued_fraction_list()
     915            ((0, 4, 1, 1, 1, 1, 3, 2), (0,))
     916        """
     917        cdef NumberFieldElement_quadratic x
     918
     919        if mpz_sgn(self.b) == 0:
     920            return tuple(Rational(self).continued_fraction_list()),(0,)
     921
     922        if mpz_sgn(self.D.value) < 0:
     923            raise ValueError("the method is only available for positive discriminant")
     924
     925        x = self
     926        orbit = []
     927        quots = []
     928        while x not in orbit:
     929            quots.append(x.floor())
     930            orbit.append(x)
     931            x = ~(x - quots[-1])
     932
     933        i = orbit.index(x)
     934
     935        return tuple(quots[:i]), tuple(quots[i:])
     936
     937    def continued_fraction(self, R=None):
     938        r"""
     939        Return the (finite or ultimately periodic) continued fraction of ``self``.
     940
     941        EXAMPLES::
     942
     943            sage: K.<sqrt2> = QuadraticField(2)
     944            sage: cf = sqrt2.continued_fraction(); cf
     945            [1; (2)*]
     946            sage: cf.n()
     947            1.41421356237310
     948            sage: sqrt2.n()
     949            1.41421356237310
     950            sage: cf.value()
     951            sqrt2
     952
     953            sage: (sqrt2/3 + 1/4).continued_fraction()
     954            [0; 1, (2, 1, 1, 2, 3, 2, 1, 1, 2, 5, 1, 1, 14, 1, 1, 5)*]
     955        """
     956        if R is None:
     957            from sage.rings.continued_fractions import CFF
     958            R = CFF
     959        return R(self.continued_fraction_list())
     960
    899961#########################################################
    900962# Arithmetic
    901963#########################################################
  • sage/rings/rational.pyx

    diff --git a/sage/rings/rational.pyx b/sage/rings/rational.pyx
    a b AUTHORS: 
    2424
    2525- Travis Scrimshaw (2012-10-18): Added doctests for full coverage.
    2626
     27- Vincent Delecroix (2013): continued fraction
     28
    2729TESTS::
    2830
    2931    sage: a = -2/3
    cdef class Rational(sage.structure.eleme 
    574576        """
    575577        return [ self ]
    576578
     579    def continued_fraction_list(self, type="std"):
     580        r"""
     581        Return the list of partial quotients of this rational number.
     582
     583        INPUT:
     584
     585        - ``type`` - either "std" (the default) for the standard continued
     586          fractions or "hj" for the Hirzebruch-Jung ones.
     587
     588        EXAMPLES::
     589
     590            sage: (13/9).continued_fraction_list()
     591            [1, 2, 4]
     592            sage: 1 + 1/(2 + 1/4)
     593            13/9
     594
     595            sage: (225/157).continued_fraction_list()
     596            [1, 2, 3, 4,  5]
     597            sage: 1 + 1/(2 + 1/(3 + 1/(4 + 1/5)))
     598            225/157
     599
     600            sage: (fibonacci(20)/fibonacci(19)).continued_fraction_list()
     601            [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
     602
     603            sage: (-1/3).continued_fraction_list()
     604            [-1, 1, 2]
     605
     606        Check that the partial quotients of an integer ``n`` is simply ``[n]``::
     607
     608            sage: QQ(1).continued_fraction_list()
     609            [1]
     610            sage: QQ(0).continued_fraction_list()
     611            [0]
     612            sage: QQ(-1).continued_fraction_list()
     613            [-1]
     614
     615        Hirzebruch-Jung continued fractions::
     616
     617            sage: (11/19).continued_fraction_list("hj")
     618            [1, 3, 2, 3, 2]
     619            sage: 1 - 1/(3 - 1/(2 - 1/(3 - 1/2)))
     620            11/19
     621
     622            sage: (225/137).continued_fraction_list("hj")
     623            [2, 3, 5, 10]
     624            sage: 2 - 1/(3 - 1/(5 - 1/10))
     625            225/137
     626
     627            sage: (-23/19).continued_fraction_list("hj")
     628            [-1, 5, 4]
     629            sage: -1 - 1/(5 - 1/4)
     630            -23/19
     631        """
     632        cdef Integer z
     633        cdef mpz_t p,q,tmp
     634        cdef list res = []
     635
     636        mpz_init(tmp)
     637        mpz_init(p)
     638        mpz_init(q)
     639        mpz_set(p, mpq_numref(self.value))
     640        mpz_set(q, mpq_denref(self.value))
     641
     642        if type == "std":
     643            while mpz_sgn(q) != 0:
     644                z = PY_NEW(Integer)
     645                mpz_fdiv_qr(z.value,tmp,p,q)
     646                mpz_set(p,q)
     647                mpz_set(q,tmp)
     648                res.append(z)
     649        elif type == "hj":
     650            while mpz_sgn(q) != 0:
     651                z = PY_NEW(Integer)
     652                mpz_cdiv_qr(z.value,tmp,p,q)
     653                mpz_set(p,q)
     654                mpz_set(q,tmp)
     655                res.append(z)
     656                if mpz_sgn(q) == 0:
     657                    break
     658                z = PY_NEW(Integer)
     659                mpz_fdiv_qr(z.value,tmp,p,q)
     660                mpz_set(p,q)
     661                mpz_set(q,tmp)
     662                mpz_neg(z.value,z.value)
     663                res.append(z)
     664        else:
     665            mpz_clear(p)
     666            mpz_clear(q)
     667            mpz_clear(tmp)
     668            raise ValueError("the type must be one of 'floor', 'hj'")
     669
     670        mpz_clear(p)
     671        mpz_clear(q)
     672        mpz_clear(tmp)
     673
     674        return res
     675
     676    def continued_fraction(self, R=None):
     677        r"""
     678        Return the continued fraction of that rational.
     679
     680        EXAMPLES::
     681
     682            sage: (641/472).continued_fraction()
     683            [1; 2, 1, 3, 1, 4, 1, 5]
     684
     685            sage: a = (355/113).continued_fraction(); a
     686            [3; 7, 16]
     687            sage: a.n(digits=10)
     688            3.141592920
     689
     690        It's almost pi!
     691        """
     692        if R is None:
     693            from sage.rings.continued_fractions import CFF
     694            R = CFF
     695        return R(self.continued_fraction_list())
    577696
    578697    def __richcmp__(left, right, int op):
    579698        """
  • sage/tests/book_stein_ent.py

    diff --git a/sage/tests/book_stein_ent.py b/sage/tests/book_stein_ent.py
    a b sage: find_sqrt(5,389) # see, it 
    38938986
    390390
    391391# Several of the examples below had to be changed due to improved
    392 # behavior of the continued_fraction function #8017.
     392# behavior of the continued_fraction function #8017 and #14567.
    393393
    394 sage: continued_fraction(17/23) 
    395 [0, 1, 2, 1, 5]
     394sage: continued_fraction(17/23)
     395[0; 1, 2, 1, 5]
    396396sage: reset('e')
    397 sage: continued_fraction(e) 
    398 [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1]
     397sage: continued_fraction(e)
     398[2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1, ...]
    399399sage: continued_fraction(e, bits=21)
    400 [2, 1, 2, 1, 1, 4, 1, 1, 6]
     400[2; 1, 2, 1, 1, 4, 1, 1, 6]
    401401sage: continued_fraction(e, bits=30)
    402 [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8]
     402[2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8]
    403403sage: a = continued_fraction(17/23); a
    404 [0, 1, 2, 1, 5]
    405 sage: a.value() 
     404[0; 1, 2, 1, 5]
     405sage: a.value()
    40640617/23
    407407sage: b = continued_fraction(6/23); b
    408 [0, 3, 1, 5]
    409 sage: a + b 
     408[0; 3, 1, 5]
     409sage: a + b
    410410[1]
    411411sage: c = continued_fraction(pi,bits=35); c
    412 [3, 7, 15, 1, 292, 1]
    413 sage: c.convergents() 
    414 [3, 22/7, 333/106, 355/113, 103993/33102, 104348/33215]
    415 sage: c = continued_fraction(pi); c 
    416 [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
    417 sage: for n in range(-1, len(c)):
     412[3; 7, 15, 1, 293]
     413sage: c.convergents()
     414[3, 22/7, 333/106, 355/113, 104348/33215]
     415sage: c = continued_fraction(pi); c
     416[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
     417sage: for n in range(-1, 13):
    418418...    print c.pn(n)*c.qn(n-1) - c.qn(n)*c.pn(n-1),
    4194191 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1
    420 sage: for n in range(len(c)):
    421 ...    print c.pn(n)*c.qn(n-2) - c.qn(n)*c.pn(n-2), 
     420sage: for n in range(13):
     421...    print c.pn(n)*c.qn(n-2) - c.qn(n)*c.pn(n-2),
    4224223 -7 15 -1 292 -1 1 -1 2 -1 3 -1 14
    423 sage: c = continued_fraction([1,2,3,4,5]) 
    424 sage: c.convergents() 
     423sage: c = continued_fraction([1,2,3,4,5])
     424sage: c.convergents()
    425425[1, 3/2, 10/7, 43/30, 225/157]
    426 sage: [c.pn(n) for n in range(len(c))] 
     426sage: [c.pn(n) for n in range(len(c))]
    427427[1, 3, 10, 43, 225]
    428 sage: [c.qn(n) for n in range(len(c))] 
     428sage: [c.qn(n) for n in range(len(c))]
    429429[1, 2, 7, 30, 157]
    430430sage: c = continued_fraction([1,1,1,1,1,1,1,1])
    431431sage: v = [(i, c.pn(i)/c.qn(i)) for i in range(len(c))]
    432432sage: P = point(v, rgbcolor=(0,0,1), pointsize=40)
    433 sage: L = line(v, rgbcolor=(0.5,0.5,0.5)) 
     433sage: L = line(v, rgbcolor=(0.5,0.5,0.5))
    434434sage: L2 = line([(0,c.value()),(len(c)-1,c.value())], \
    435435...      thickness=0.5, rgbcolor=(0.7,0,0))
    436 sage: (L+L2+P).show(xmin=0,ymin=1) 
     436sage: (L+L2+P).show(xmin=0,ymin=1)
    437437sage: def cf(bits):
    438438...   x = (1 + sqrt(RealField(bits)(5))) / 2
    439439...   return continued_fraction(x)
    440440sage: cf(10)
    441 [1, 1, 1, 1, 1, 1, 1]
    442 sage: cf(30) 
    443 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    444 sage: cf(50) 
    445 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
     441[1; 1, 1, 1, 1, 1, 1, 2]
     442sage: cf(30)
     443[1; 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
     444sage: cf(50)
     445 [1; 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
    446446sage: def cf_sqrt_d(d, bits):
    447447...   x = sqrt(RealField(bits)(d))
    448448...   return continued_fraction(x)
    449449sage: cf_sqrt_d(389,50)
    450 [19, 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1, 1, 1, 2, 1, 38, 2]
     450[19; 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1, 1, 1, 2, 1, 38]
    451451sage: cf_sqrt_d(389,100)
    452 [19, 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1]
     452[19; 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 1, 1, 1, 1, 2, 1, 38, 1, 2, 2]
    453453sage: def newton_root(f, iterates=2, x0=0, prec=53):
    454454...    x = RealField(prec)(x0)
    455455...    R = PolynomialRing(ZZ,'x')
    sage: reset('x') 
    462462sage: a = newton_root(3847*x^2 - 14808904*x + 36527265); a
    4634632.46815700480740
    464464sage: cf = continued_fraction(a); cf
    465 [2, 2, 7, 2, 1, 5, 1, 1, 1, 1, 1, 1, 103, 8, 1, 2, 3]
     465[2; 2, 7, 2, 1, 5, 1, 1, 1, 1, 1, 1, 103, 8, 1, 2, 3, 2]
    466466sage: c = cf[:12]; c
    467 [2, 2, 7, 2, 1, 5, 1, 1, 1, 1, 1, 1]
     467[2; 2, 7, 2, 1, 5, 1, 1, 1, 1, 2]
    468468sage: c.value()
    4694699495/3847
    470470sage: def sum_of_two_squares_naive(n):
  • sage/tests/book_stein_modform.py

    diff --git a/sage/tests/book_stein_modform.py b/sage/tests/book_stein_modform.py
    a b sage: S = M.cuspidal_submodule() 
    7373sage: S.integral_basis()     # basis over ZZ.
    7474({-1/8, 0}, {-1/9, 0})
    7575sage: set_modsym_print_mode ('manin')    # set it back
    76 sage: convergents(4/7)
     76sage: continued_fraction(4/7).convergents()
    7777[0, 1, 1/2, 4/7]
    7878sage: M = ModularSymbols(2,2)
    7979sage: M