cfinite_sequence.py

# -*- coding: utf-8 -*-
r"""
C-Finite Sequences

C-finite infinite sequences satisfy homogenous linear recurrences with constant coefficients:

.. MATH::

    a_{n+d} = c_0a_n + c_1a_{n+1} + \cdots + c_{d-1}a_{n+d-1}, \quad d>0.

CFiniteSequences are completely defined by their ordinary generating function (o.g.f., which
is always a :mod:`fraction <sage.rings.fraction_field_element>` of
:mod:`polynomials <sage.rings.polynomial.polynomial_element>` over `\mathbb{Z}` or `\mathbb{Q}` ).

EXAMPLES::

    sage: fibo = CFiniteSequence(x/(1-x-x^2))        # the Fibonacci sequence
    sage: fibo
    C-finite sequence, generated by x/(-x^2 - x + 1)
    sage: fibo.parent()
    The ring of C-Finite sequences in x over Rational Field
    sage: fibo.parent().category()
    Category of commutative rings
    sage: C.<x> = CFiniteSequences(QQ);
    sage: fibo.parent() == C
    True
    sage: C
    The ring of C-Finite sequences in x over Rational Field
    sage: C(x/(1-x-x^2))
    C-finite sequence, generated by x/(-x^2 - x + 1)
    sage: C(x/(1-x-x^2)) == fibo
    True
    sage: var('y')
    y
    sage: CFiniteSequence(y/(1-y-y^2))
    C-finite sequence, generated by y/(-y^2 - y + 1)
    sage: CFiniteSequence(y/(1-y-y^2)) == fibo
    False

Finite subsets of the sequence are accessible via python slices::

    sage: fibo[137]                 #the 137th term of the Fibonacci sequence
    19134702400093278081449423917
    sage: fibo[137] == fibonacci(137)
    True
    sage: fibo[0:12]
    [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]
    sage: fibo[14:4:-2]
    [377, 144, 55, 21, 8]

They can be created also from the coefficients and start values of a recurrence::

    sage: r = C.from_recurrence([1,1],[0,1])
    sage: r == fibo
    True

Given enough values, the o.g.f. of a C-finite sequence
can be guessed::

    sage: r = C.guess([0,1,1,2,3,5,8])
    sage: r == fibo
    True

.. SEEALSO::

    :func:`fibonacci`, :class:`BinaryRecurrenceSequence`

AUTHORS:

- Ralf Stephan (2014): initial version

REFERENCES:

.. [GK82] Greene, Daniel H.; Knuth, Donald E. (1982), "2.1.1 Constant
   coefficients - A) Homogeneous equations", Mathematics for the Analysis
   of Algorithms (2nd ed.), Birkhauser, p. 17.
.. [SZ94] Bruno Salvy and Paul Zimmermann. - Gfun: a Maple package for
   the manipulation of generating and holonomic functions in one variable.
   - Acm transactions on mathematical software, 20.2:163-177, 1994.
.. [Z11] Zeilberger, Doron. "The C-finite ansatz." The Ramanujan Journal
   (2011): 1-10.
"""

#*****************************************************************************
#       Copyright (C) 2014 Ralf Stephan <gtrwst9@gmail.com>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#                  http://www.gnu.org/licenses/
#*****************************************************************************


from sage.categories.fields import Fields
from sage.misc.inherit_comparison import InheritComparisonClasscallMetaclass
from sage.rings.ring import CommutativeRing
from sage.rings.integer import Integer
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.arith.all import gcd
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.polynomial.polynomial_ring import PolynomialRing_general
from sage.rings.laurent_series_ring import LaurentSeriesRing
from sage.rings.power_series_ring import PowerSeriesRing
from sage.rings.fraction_field import FractionField
from sage.structure.element import FieldElement
from sage.structure.unique_representation import UniqueRepresentation

from sage.interfaces.gp import Gp
from sage.misc.all import sage_eval

_gp = None

def CFiniteSequences(base_ring, names = None, category = None):
    r"""
    Return the ring of C-Finite sequences.

    The ring is defined over a base ring (`\mathbb{Z}` or `\mathbb{Q}` )
    and each element is represented by its ordinary generating function (ogf)
    which is a rational function over the base ring.

    INPUT:

    - ``base_ring`` -- the base ring to construct the fraction field
      representing the C-Finite sequences
    - ``names`` -- (optional) the list of variables.

    EXAMPLES::

        sage: C.<x> = CFiniteSequences(QQ)
        sage: C
        The ring of C-Finite sequences in x over Rational Field
        sage: C.an_element()
        C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
        sage: C.category()
        Category of commutative rings
        sage: C.one()
        Finite sequence [1], offset = 0
        sage: C.zero()
        Constant infinite sequence 0.
        sage: C(x)
        Finite sequence [1], offset = 1
        sage: C(1/x)
        Finite sequence [1], offset = -1
        sage: C((-x + 2)/(-x^2 - x + 1))
        C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)

    TESTS::

        sage: TestSuite(C).run()
    """
    if isinstance(base_ring, PolynomialRing_general):
        polynomial_ring = base_ring
        base_ring = polynomial_ring.base_ring()
    if names is None:
        names = ['x']
    elif len(names)>1:
        raise NotImplementedError("Multidimensional o.g.f. not implemented.")
    if category is None:
        category = Fields()
    if not(base_ring in (QQ, ZZ)):
        raise ValueError("O.g.f. base not rational.")
    polynomial_ring = PolynomialRing(base_ring, names)
    return CFiniteSequences_generic(polynomial_ring, category)


class CFiniteSequence(FieldElement):
    r"""
    Create a C-finite sequence given its ordinary generating function.

    INPUT:

    - ``ogf`` -- a rational function, the ordinary generating function
      (can be a an element from the symbolic ring, fraction field or polynomial
      ring)

    OUTPUT:

    - A CFiniteSequence object

    EXAMPLES::

        sage: CFiniteSequence((2-x)/(1-x-x^2))     # the Lucas sequence
        C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
        sage: CFiniteSequence(x/(1-x)^3)           # triangular numbers
        C-finite sequence, generated by x/(-x^3 + 3*x^2 - 3*x + 1)

    Polynomials are interpreted as finite sequences, or recurrences of degree 0::

        sage: CFiniteSequence(x^2-4*x^5)
        Finite sequence [1, 0, 0, -4], offset = 2
        sage: CFiniteSequence(1)
        Finite sequence [1], offset = 0

    This implementation allows any polynomial fraction as o.g.f. by interpreting
    any power of `x` dividing the o.g.f. numerator or denominator as a right or left shift
    of the sequence offset::

        sage: CFiniteSequence(x^2+3/x)
        Finite sequence [3, 0, 0, 1], offset = -1
        sage: CFiniteSequence(1/x+4/x^3)
        Finite sequence [4, 0, 1], offset = -3
        sage: P = LaurentPolynomialRing(QQ.fraction_field(), 'X')
        sage: X=P.gen()
        sage: CFiniteSequence(1/(1-X))
        C-finite sequence, generated by 1/(-X + 1)

    The o.g.f. is always normalized to get a denominator constant coefficient of `+1`::

        sage: CFiniteSequence(1/(x-2))
        C-finite sequence, generated by -1/2/(-1/2*x + 1)

    The given ``ogf`` is used to create an appropriate parent: it can
    be a symbolic expression, a polynomial , or a fraction field element
    as long as it can be coerced into a proper fraction field over the
    rationals::

        sage: var('x')
        x
        sage: f1 = CFiniteSequence((2-x)/(1-x-x^2))
        sage: P.<x> = QQ[]
        sage: f2 = CFiniteSequence((2-x)/(1-x-x^2))
        sage: f1 == f2
        True
        sage: f1.parent()
        The ring of C-Finite sequences in x over Rational Field
        sage: f1.ogf().parent()
        Fraction Field of Univariate Polynomial Ring in x over Rational Field
        sage: CFiniteSequence(log(x))
        Traceback (most recent call last):
        ...
        TypeError: unable to convert log(x) to a rational

    TESTS::

        sage: P.<x> = QQ[]
        sage: CFiniteSequence(0.1/(1-x))
        C-finite sequence, generated by 1/10/(-x + 1)
        sage: CFiniteSequence(pi/(1-x))
        Traceback (most recent call last):
        ...
        TypeError: unable to convert -pi to a rational
        sage: P.<x,y> = QQ[]
        sage: CFiniteSequence(x*y)
        Traceback (most recent call last):
        ...
        NotImplementedError: Multidimensional o.g.f. not implemented.
    """

    __metaclass__ = InheritComparisonClasscallMetaclass

    @staticmethod
    def __classcall_private__(cls, ogf):
        r"""
        Ensures that elements created by :class:`CFiniteSequence` have the same
        parent than the ones created by the parent itself and follow the category
        framework (they should be instance of :class:`CFiniteSequences` automatic
        element class).

        This method is called before the ``__init__`` method, it checks the
        o.g.f to create the appropriate parent.

        INPUT:

        - ``ogf`` - a rational function

        TESTS::

            sage: f1 = CFiniteSequence((2-x)/(1-x-x^2))
            sage: f1
            C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
            sage: C.<x> = CFiniteSequences(QQ);
            sage: f2 = CFiniteSequence((2-x)/(1-x-x^2))
            sage: f2
            C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
            sage: f3 = C((2-x)/(1-x-x^2))
            sage: f3
            C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
            sage: f1 == f2 and f2 == f3
            True
            sage: f1.parent() == f2.parent() and f2.parent() == f3.parent()
            True
            sage: type(f1)
            <class 'sage.rings.cfinite_sequence.CFiniteSequences_generic_with_category.element_class'>
            sage: type(f1) == type(f2) and type(f2) == type(f3)
            True
            sage: CFiniteSequence(log(x))
            Traceback (most recent call last):
            ...
            TypeError: unable to convert log(x) to a rational
            sage: CFiniteSequence(pi)
            Traceback (most recent call last):
            ...
            TypeError: Unable to coerce pi (<class 'sage.symbolic.constants.Pi'>) to Rational
            sage: var('y')
            y
            sage: f4 = CFiniteSequence((2-y)/(1-y-y^2))
            sage: f4
            C-finite sequence, generated by (-y + 2)/(-y^2 - y + 1)
            sage: f4 == f1
            False
            sage: f4.parent() == f1.parent()
            False
            sage: f4.parent()
            The ring of C-Finite sequences in y over Rational Field
        """

        br = ogf.base_ring()
        if not(br in (QQ, ZZ)):
            br = QQ  # if the base ring of the o.g.f is not QQ, we force it to QQ and see if the o.g.f converts nicely

        # trying to figure out the ogf variables
        variables = []
        if not ogf in br:
            if hasattr(ogf, 'variables'):
                variables =  ogf.variables()
            elif hasattr(ogf.parent(), 'gens'):
                variables = ogf.parent().gens()
                # for some reason, fraction field elements don't have the variables
                # method, but symbolic elements don't have the gens method so we check both

        if len(variables)==0:
            parent = CFiniteSequences(QQ) # if we cannot find variables, we create the default parent (with x)
        else:
            parent = CFiniteSequences(QQ, variables)
        return parent(ogf) # if ogf cannot be converted to a fraction field, this will break and raise the proper error

    def __init__(self, parent, ogf):
        r"""
        Initialize the C-Finite sequence.

        The ``__init__`` method can only be called by the :class:`CFiniteSequences`
        class. By Default, a class call reaches the ``__classcall_private__``
        which first creates a proper parent and then call the ``__init__``.

        INPUT:

        - ``ogf`` -- the ordinary generating function, a fraction of polynomials over the rationals
        - ``parent`` --  the parent of the C-Finite sequence, an occurence of :class:`CFiniteSequences`

        OUTPUT:

        - A CFiniteSequence object

        TESTS::

            sage: C.<x> = CFiniteSequences(QQ);
            sage: C((2-x)/(1-x-x^2))            # indirect doctest
            C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
        """

        br = parent.base_ring()
        ogf = parent.fraction_field()(ogf)
        P = parent.polynomial_ring()
        num = ogf.numerator()
        den = ogf.denominator()

        FieldElement.__init__(self, parent)

        if den == 1:
            self._c = []
            self._off = num.valuation()
            self._deg = 0
            if ogf == 0:
                self._a = [0]
            else:
                self._a = P((num / (P.gen()) ** self._off)).list()
        else:

            # Transform the ogf numerator and denominator to canonical form
            # to get the correct offset, degree, and recurrence coeffs and
            # start values.
            self._off = 0
            self._deg = 0
            x = P.gen()
            if num.constant_coefficient() == 0:
                self._off = num.valuation()
                num = P(num / x ** self._off)
            elif den.constant_coefficient() == 0:
                self._off = -den.valuation()
                den = P(den * x ** self._off)
            f = den.constant_coefficient()
            num = P(num / f)
            den = P(den / f)
            f = gcd(num, den)
            num = P(num / f)
            den = P(den / f)
            self._deg = den.degree()
            self._c = [-den.list()[i] for i in range(1, self._deg + 1)]
            if self._off >= 0:
                num = x ** self._off * num
            else:
                den = x ** (-self._off) * den

            # determine start values (may be different from _get_item_ values)
            alen = max(self._deg, num.degree() + 1)
            R = LaurentSeriesRing(br, parent.variable_name(), default_prec=alen)
            rem = num % den
            if den != 1:
                self._a = R(num / den).list()
                self._aa = R(rem / den).list()[:self._deg]  # needed for _get_item_
            else:
                self._a = num.list()
            if len(self._a) < alen:
                self._a.extend([0] * (alen - len(self._a)))

            ogf = num / den

        self._ogf = ogf


    def _repr_(self):
        """
        Return textual definition of sequence.

        TESTS::

            sage: CFiniteSequence(1/x^5)
            Finite sequence [1], offset = -5
            sage: CFiniteSequence(x^3)
            Finite sequence [1], offset = 3
        """
        if self._deg == 0:
            if self.ogf() == 0:
                return 'Constant infinite sequence 0.'
            else:
                return 'Finite sequence ' + str(self._a) + ', offset = ' + str(self._off)
        else:
            return 'C-finite sequence, generated by ' + str(self.ogf())

    def __hash__(self):
        r"""
        Hash value for C finite sequence.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: hash(C((2-x)/(1-x-x^2)))  # random
            42
        """
        return hash(self.parent()) ^ hash(self._ogf)

    def _add_(self, other):
        """
        Addition of C-finite sequences.

        TESTS::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: r = C(1/(1-2*x))
            sage: r[0:5]                  # a(n) = 2^n
            [1, 2, 4, 8, 16]
            sage: s = C.from_recurrence([1],[1])
            sage: (r + s)[0:5]            # a(n) = 2^n + 1
            [2, 3, 5, 9, 17]
            sage: r + 0 == r
            True
            sage: (r + x^2)[0:5]
            [1, 2, 5, 8, 16]
            sage: (r + 3/x)[-1]
            3
            sage: r = CFiniteSequence(x)
            sage: r + 0 == r
            True
            sage: CFiniteSequence(0) + CFiniteSequence(0)
            Constant infinite sequence 0.
        """
        return CFiniteSequence(self.ogf() + other.numerator() / other.denominator())

    def _sub_(self, other):
        """
        Subtraction of C-finite sequences.

        TESTS::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: r = C(1/(1-2*x))
            sage: r[0:5]                  # a(n) = 2^n
            [1, 2, 4, 8, 16]
            sage: s = C.from_recurrence([1],[1])
            sage: (r - s)[0:5]            # a(n) = 2^n + 1
            [0, 1, 3, 7, 15]
        """
        return CFiniteSequence(self.ogf() - other.numerator() / other.denominator())

    def _mul_(self, other):
        """
        Multiplication of C-finite sequences.

        TESTS::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: r = C.guess([1,2,3,4,5,6])
            sage: (r*r)[0:6]                  # self-convolution
            [1, 4, 10, 20, 35, 56]
            sage: r = C(x)
            sage: r*1 == r
            True
            sage: r*-1
            Finite sequence [-1], offset = 1
            sage: C(0) * C(1)
            Constant infinite sequence 0.
        """
        return CFiniteSequence(self.ogf() * other.numerator() / other.denominator())

    def _div_(self, other):
        """
        Division of C-finite sequences.

        TESTS::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: r = C.guess([1,2,3,4,5,6])
            sage: (r/2)[0:6]
            [1/2, 1, 3/2, 2, 5/2, 3]
            sage: s = C(x)
            sage: s/(s*-1 + 1)
            C-finite sequence, generated by x/(-x + 1)
        """
        return CFiniteSequence(self.ogf() / (other.numerator() / other.denominator()))

    def coefficients(self):
        """
        Return the coefficients of the recurrence representation of the
        C-finite sequence.

        OUTPUT:

        - A list of values

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: lucas = C((2-x)/(1-x-x^2))   # the Lucas sequence
            sage: lucas.coefficients()
            [1, 1]
        """
        return self._c

    def __eq__(self, other):
        """
        Compare two CFiniteSequences.

        EXAMPLES::

            sage: f = CFiniteSequence((2-x)/(1-x-x^2))
            sage: f2 = CFiniteSequence((2-x)/(1-x-x^2))
            sage: f == f2
            True
            sage: f == (2-x)/(1-x-x^2)
            False
            sage: (2-x)/(1-x-x^2) == f
            False
            sage: C.<x> = CFiniteSequences(QQ)
            sage: r = C.from_recurrence([1,1],[2,1])
            sage: s = C.from_recurrence([-1],[1])
            sage: r == s
            False
            sage: r = C.from_recurrence([-1],[1])
            sage: s = C(1/(1+x))
            sage: r == s
            True
        """
        if not isinstance(other, CFiniteSequence):
            return False
        return self.ogf() == other.ogf()

    def __getitem__(self, key):
        r"""
        Return a slice of the sequence.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: r = C.from_recurrence([3,3],[2,1])
            sage: r[2]
            9
            sage: r[101]
            16158686318788579168659644539538474790082623100896663971001
            sage: r = C(1/(1-x))
            sage: r[5]
            1
            sage: r = C(x)
            sage: r[0]
            0
            sage: r[1]
            1
            sage: r = C(0)
            sage: r[66]
            0
            sage: lucas = C.from_recurrence([1,1],[2,1])
            sage: lucas[5:10]
            [11, 18, 29, 47, 76]
            sage: r = C((2-x)/x/(1-x-x*x))
            sage: r[0:4]
            [1, 3, 4, 7]
            sage: r = C(1-2*x^2)
            sage: r[0:4]
            [1, 0, -2, 0]
            sage: r[-1:4]             # not tested, python will not allow this!
            [0, 1, 0 -2, 0]
            sage: r = C((-2*x^3 + x^2 + 1)/(-2*x + 1))
            sage: r[0:5]              # handle ogf > 1
            [1, 2, 5, 8, 16]
            sage: r[-2]
            0
            sage: r = C((-2*x^3 + x^2 - x + 1)/(2*x^2 - 3*x + 1))
            sage: r[0:5]
            [1, 2, 5, 9, 17]
            sage: s=C((1-x)/(-x^2 - x + 1))
            sage: s[0:5]
            [1, 0, 1, 1, 2]
            sage: s=C((1+x^20+x^40)/(1-x^12)/(1-x^30))
            sage: s[0:20]
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]
            sage: s=C(1/((1-x^2)*(1-x^6)*(1-x^8)*(1-x^12)))
            sage: s[999998]
            289362268629630
        """
        if isinstance(key, slice):
            m = max(key.start, key.stop)
            return [self[ii] for ii in xrange(*key.indices(m + 1))]
        elif isinstance(key, (int, Integer)):
            from sage.matrix.constructor import Matrix
            d = self._deg
            if (self._off <= key and key < self._off + len(self._a)):
                return self._a[key - self._off]
            elif d == 0:
                return 0
            (quo, rem) = self.numerator().quo_rem(self.denominator())
            wp = quo[key - self._off]
            if key < self._off:
                return wp
            A = Matrix(QQ, 1, d, self._c)
            B = Matrix.identity(QQ, d - 1)
            C = Matrix(QQ, d - 1, 1, 0)
            if quo == 0:
                V = Matrix(QQ, d, 1, self._a[:d][::-1])
            else:
                V = Matrix(QQ, d, 1, self._aa[:d][::-1])
            M = Matrix.block([[A], [B, C]], subdivide=False)

            return wp + list(M ** (key - self._off) * V)[d - 1][0]
        else:
            raise TypeError("invalid argument type")

    def ogf(self):
        """
        Return the ordinary generating function associated with
        the CFiniteSequence.

        This is always a fraction of polynomials in the base ring.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: r = C.from_recurrence([2],[1])
            sage: r.ogf()
            1/(-2*x + 1)
            sage: C(0).ogf()
            0
        """
        return self._ogf

    def numerator(self):
        r"""
        Return the numerator of the o.g.f of ``self``.

        EXAMPLES::

            sage: f = CFiniteSequence((2-x)/(1-x-x^2)); f
            C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
            sage: f.numerator()
            -x + 2
        """
        return self.ogf().numerator()

    def denominator(self):
        r"""
        Return the numerator of the o.g.f of ``self``.

        EXAMPLES::

            sage: f = CFiniteSequence((2-x)/(1-x-x^2)); f
            C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
            sage: f.denominator()
            -x^2 - x + 1
        """
        return self.ogf().denominator()

    def recurrence_repr(self):
        """
        Return a string with the recurrence representation of
        the C-finite sequence.

        OUTPUT:

        - A string

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: C((2-x)/(1-x-x^2)).recurrence_repr()
            'Homogenous linear recurrence with constant coefficients of degree 2: a(n+2) = a(n+1) + a(n), starting a(0...) = [2, 1]'
            sage: C(x/(1-x)^3).recurrence_repr()
            'Homogenous linear recurrence with constant coefficients of degree 3: a(n+3) = 3*a(n+2) - 3*a(n+1) + a(n), starting a(1...) = [1, 3, 6]'
            sage: C(1).recurrence_repr()
            'Finite sequence [1], offset 0'
            sage: r = C((-2*x^3 + x^2 - x + 1)/(2*x^2 - 3*x + 1))
            sage: r.recurrence_repr()
            'Homogenous linear recurrence with constant coefficients of degree 2: a(n+2) = 3*a(n+1) - 2*a(n), starting a(0...) = [1, 2, 5, 9]'
            sage: r = CFiniteSequence(x^3/(1-x-x^2))
            sage: r.recurrence_repr()
            'Homogenous linear recurrence with constant coefficients of degree 2: a(n+2) = a(n+1) + a(n), starting a(3...) = [1, 1, 2, 3]'
        """
        if self._deg == 0:
            return 'Finite sequence %s, offset %d' % (str(self._a), self._off)
        else:
            if self._c[0] == 1:
                cstr = 'a(n+%d) = a(n+%d)' % (self._deg, self._deg - 1)
            elif self._c[0] == -1:
                cstr = 'a(n+%d) = -a(n+%d)' % (self._deg, self._deg - 1)
            else:
                cstr = 'a(n+%d) = %s*a(n+%d)' % (self._deg, str(self._c[0]), self._deg - 1)
            for i in range(1, self._deg):
                j = self._deg - i - 1
                if self._c[i] < 0:
                    if self._c[i] == -1:
                        cstr = cstr + ' - a(n+%d)' % (j,)
                    else:
                        cstr = cstr + ' - %d*a(n+%d)' % (-(self._c[i]), j)
                elif self._c[i] > 0:
                    if self._c[i] == 1:
                        cstr = cstr + ' + a(n+%d)' % (j,)
                    else:
                        cstr = cstr + ' + %d*a(n+%d)' % (self._c[i], j)
            cstr = cstr.replace('+0', '')
        astr = ', starting a(%s...) = [' % str(self._off)
        maxwexp = self.numerator().quo_rem(self.denominator())[0].degree() + 1
        for i in range(maxwexp + self._deg):
            astr = astr + str(self[self._off + i]) + ', '
        astr = astr[:-2] + ']'
        return 'Homogenous linear recurrence with constant coefficients of degree ' + str(self._deg) + ': ' + cstr + astr

    def series(self, n):
        """
        Return the Laurent power series associated with the
        CFiniteSequence, with precision `n`.

        INPUT:

        - `n` -- a nonnegative integer

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: r = C.from_recurrence([-1,2],[0,1])
            sage: s = r.series(4); s
            x + 2*x^2 + 3*x^3 + 4*x^4 + O(x^5)
            sage: type(s)
            <type 'sage.rings.laurent_series_ring_element.LaurentSeries'>
        """
        R = LaurentSeriesRing(QQ, 'x', default_prec=n)
        return R(self.ogf())


class CFiniteSequences_generic(CommutativeRing, UniqueRepresentation):
    r"""
    The class representing the ring of C-Finite Sequences

    TESTS::

        sage: C.<x> = CFiniteSequences(QQ)
        sage: from sage.rings.cfinite_sequence import CFiniteSequences_generic
        sage: isinstance(C,CFiniteSequences_generic)
        True
        sage: type(C)
        <class 'sage.rings.cfinite_sequence.CFiniteSequences_generic_with_category'>
        sage: C
        The ring of C-Finite sequences in x over Rational Field
    """

    Element = CFiniteSequence
    def __init__(self, polynomial_ring, category):
        r"""
        Create the ring of CFiniteSequences over ``base_ring``

        INPUT:

        - ``base_ring`` -- the base ring for the o.g.f (either ``QQ`` or ``ZZ``)
        - ``names`` -- an iterable of variables (shuould contain only one variable)
        - ``category`` -- the category of the ring (default: ``Fields()``)

        TESTS::

            sage: C.<y> = CFiniteSequences(QQ); C
            The ring of C-Finite sequences in y over Rational Field
            sage: C.<x> = CFiniteSequences(QQ); C
            The ring of C-Finite sequences in x over Rational Field
            sage: C.<x> = CFiniteSequences(ZZ); C
            The ring of C-Finite sequences in x over Integer Ring
            sage: C.<x,y> = CFiniteSequences(ZZ)
            Traceback (most recent call last):
            ...
            NotImplementedError: Multidimensional o.g.f. not implemented.
            sage: C.<x> = CFiniteSequences(CC)
            Traceback (most recent call last):
            ...
            ValueError: O.g.f. base not rational.
        """
        base_ring = polynomial_ring.base_ring()
        self._polynomial_ring = polynomial_ring
        self._fraction_field = FractionField(self._polynomial_ring)
        CommutativeRing.__init__(self,base_ring, self._polynomial_ring.gens(), category)

    def _repr_(self):
        r"""
        Return the string representation of ``self``

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: C
            The ring of C-Finite sequences in x over Rational Field
        """
        return "The ring of C-Finite sequences in {} over {}".format(self.gen(), self.base_ring())

    def _element_constructor_(self, ogf):
        r"""
        Construct a C-Finite Sequence

        INPUT:

         - ``ogf`` -- the ordinary generating function, a fraction of polynomials over the rationals

        TESTS::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: C((2-x)/(1-x-x^2))
            C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
            sage: C(x/(1-x)^3)
            C-finite sequence, generated by x/(-x^3 + 3*x^2 - 3*x + 1)
            sage: C(x^2-4*x^5)
            Finite sequence [1, 0, 0, -4], offset = 2
            sage: C(x^2+3/x)
            Finite sequence [3, 0, 0, 1], offset = -1
            sage: C(1/x + 4/x^3)
            Finite sequence [4, 0, 1], offset = -3
            sage: P = LaurentPolynomialRing(QQ.fraction_field(), 'X')
            sage: X = P.gen()
            sage: C(1/(1-X))
            C-finite sequence, generated by 1/(-x + 1)
            sage: C = CFiniteSequences(QQ)
            sage: C(x)
            Finite sequence [1], offset = 1
        """
        ogf = self.fraction_field()(ogf)
        return self.element_class(self, ogf)

    def ngens(self):
        r"""
        Return the number of generators of ``self``

        EXAMPLES::

            sage: from sage.rings.cfinite_sequence import CFiniteSequences
            sage: C.<x> = CFiniteSequences(QQ);
            sage: C.ngens()
            1
        """
        return 1

    def gen(self,i=0):
        r"""
        Return the i-th generator of ``self``.

        INPUT:

        - ``i`` -- an integer (default:0)

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ);
            sage: C.gen()
            x
            sage: x == C.gen()
            True

        TESTS::

            sage: C.gen(2)
            Traceback (most recent call last):
            ...
            ValueError: The ring of C-Finite sequences in x over Rational Field has only one generator (i=0)
        """
        if i!= 0:
            raise ValueError("{} has only one generator (i=0)".format(self))
        return self.polynomial_ring().gen()

    def an_element(self):
        r"""
        Return an element of C-Finite Sequences.

        OUTPUT:

        The Lucas sequence.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ);
            sage: C.an_element()
            C-finite sequence, generated by (-x + 2)/(-x^2 - x + 1)
        """
        x = self.gen()
        return self((2-x)/(1-x-x**2))


    def __contains__(self, x):
        """
        Return True if x is an element of ``CFinteSequences`` or
        canonically coerces to this ring.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ);
            sage: x in C
            True
            sage: 1/x in C
            True
            sage: 5 in C
            True
            sage: pi in C
            False
            sage: Cy.<y> = CFiniteSequences(QQ);
            sage: y in C
            False
            sage: y in Cy
            True
        """
        if x.parent() == self:
            return True
        try:
            self._coerce_(x)
        except TypeError:
            return False
        return True


    def fraction_field(self):
        r"""
        Return the faction field used to represent the elements of ``self``.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ);
            sage: C.fraction_field()
            Fraction Field of Univariate Polynomial Ring in x over Rational Field
        """
        return self._fraction_field

    def polynomial_ring(self):
        r"""
        Return the polynomial ring used to represent the elements of ``self``.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ);
            sage: C.polynomial_ring()
            Univariate Polynomial Ring in x over Rational Field
        """
        return self._polynomial_ring

    def _coerce_map_from_(self, S):
        """
        A coercion from `S` exists, if `S` coerces into ``self``'s fraction
        field.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ);
            sage: C.has_coerce_map_from(C.fraction_field())
            True
            sage: C.has_coerce_map_from(QQ)
            True
            sage: C.has_coerce_map_from(QQ[x])
            True
            sage: C.has_coerce_map_from(ZZ)
            True
        """
        if self.fraction_field().has_coerce_map_from(S):
            return True

    def from_recurrence(self, coefficients, values):
        """
        Create a C-finite sequence given the coefficients $c$ and
        starting values $a$ of a homogenous linear recurrence.

        .. MATH::

            a_{n+d} = c_0a_n + c_1a_{n+1} + \cdots + c_{d-1}a_{n+d-1}, \quad d\ge0.

        INPUT:

        - ``coefficients`` -- a list of rationals
        - ``values`` -- start values, a list of rationals

        OUTPUT:

        - A CFiniteSequence object

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: C.from_recurrence([1,1],[0,1])   # Fibonacci numbers
            C-finite sequence, generated by x/(-x^2 - x + 1)
            sage: C.from_recurrence([-1,2],[0,1])    # natural numbers
            C-finite sequence, generated by x/(x^2 - 2*x + 1)
            sage: r = C.from_recurrence([-1],[1])
            sage: s = C.from_recurrence([-1],[1,-1])
            sage: r == s
            True
            sage: r = C(x^3/(1-x-x^2))
            sage: s = C.from_recurrence([1,1],[0,0,0,1,1])
            sage: r == s
            True
            sage: C.from_recurrence(1,1)
            Traceback (most recent call last):
            ...
            ValueError: Wrong type for recurrence coefficient list.
        """
        if not isinstance(coefficients, list):
            raise ValueError("Wrong type for recurrence coefficient list.")
        if not isinstance(values, list):
            raise ValueError("Wrong type for recurrence start value list.")
        deg = len(coefficients)

        co = coefficients[::-1]
        co.extend([0] * (len(values) - deg))
        R = self.polynomial_ring()
        x = R.gen()
        den = -1 + sum([x ** (n + 1) * co[n] for n in range(deg)])
        num = -values[0] + sum([x ** n * (-values[n]
                                          + sum([values[k] * co[n - 1 - k]
                                                 for k in range(n)]))
                                for n in range(1, len(values))])
        return self(num / den)

    def guess(self, sequence, algorithm='sage'):
        """
        Return the minimal CFiniteSequence that generates the sequence.

        Assume the first value has index 0.

        INPUT:

        - ``sequence`` -- list of integers
        - ``algorithm`` -- string
            - 'sage' - the default is to use Sage's matrix kernel function
            - 'pari' - use Pari's implementation of LLL
            - 'bm' - use Sage's Berlekamp-Massey algorithm

        OUTPUT:

        - a CFiniteSequence, or 0 if none could be found

        With the default kernel method, trailing zeroes are chopped
        off before a guessing attempt. This may reduce the data
        below the accepted length of six values.

        EXAMPLES::

            sage: C.<x> = CFiniteSequences(QQ)
            sage: C.guess([1,2,4,8,16,32])
            C-finite sequence, generated by 1/(-2*x + 1)
            sage: r = C.guess([1,2,3,4,5])
            Traceback (most recent call last):
            ...
            ValueError: Sequence too short for guessing.

        With Berlekamp-Massey, if an odd number of values is given, the last one is dropped.
        So with an odd number of values the result may not generate the last value::

            sage: r = C.guess([1,2,4,8,9], algorithm='bm'); r
            C-finite sequence, generated by 1/(-2*x + 1)
            sage: r[0:5]
            [1, 2, 4, 8, 16]
        """
        S = self.polynomial_ring()
        if algorithm == 'bm':
            from sage.matrix.berlekamp_massey import berlekamp_massey
            if len(sequence) < 2:
                raise ValueError('Sequence too short for guessing.')
            R = PowerSeriesRing(QQ, 'x')
            if len(sequence) % 2 == 1:
                sequence = sequence[:-1]
            l = len(sequence) - 1
            denominator = S(berlekamp_massey(sequence).list()[::-1])
            numerator = R(S(sequence) * denominator, prec=l).truncate()

            return CFiniteSequence(numerator / denominator)
        elif algorithm == 'pari':
            global _gp
            if len(sequence) < 6:
                raise ValueError('Sequence too short for guessing.')
            if _gp is None:
                _gp = Gp()
                _gp("ggf(v)=local(l,m,p,q,B);l=length(v);B=floor(l/2);\
                if(B<3,return(0));m=matrix(B,B,x,y,v[x-y+B+1]);\
                q=qflll(m,4)[1];if(length(q)==0,return(0));\
                p=sum(k=1,B,x^(k-1)*q[k,1]);\
                q=Pol(Pol(vector(l,n,v[l-n+1]))*p+O(x^(B+1)));\
                if(polcoeff(p,0)<0,q=-q;p=-p);q=q/p;p=Ser(q+O(x^(l+1)));\
                for(m=1,l,if(polcoeff(p,m-1)!=v[m],return(0)));q")
            _gp.set('gf', sequence)
            _gp("gf=ggf(gf)")
            num = S(sage_eval(_gp.eval("Vec(numerator(gf))"))[::-1])
            den = S(sage_eval(_gp.eval("Vec(denominator(gf))"))[::-1])
            if num == 0:
                return 0
            else:
                return CFiniteSequence(num / den)
        else:
            from sage.matrix.constructor import matrix
            from sage.functions.other import floor, ceil
            from numpy import trim_zeros
            l = len(sequence)
            while l > 0 and sequence[l-1] == 0:
                l -= 1
            sequence = sequence[:l]
            if l == 0:
                return 0
            if l < 6:
                raise ValueError('Sequence too short for guessing.')

            hl = ceil(ZZ(l)/2)
            A = matrix([sequence[k:k+hl] for k in range(hl)])
            K = A.kernel()
            if K.dimension() == 0:
                return 0
            R = PolynomialRing(QQ, 'x')
            den = R(trim_zeros(K.basis()[-1].list()[::-1]))
            if den == 1:
                return 0
            offset = next((i for i, x in enumerate(sequence) if x!=0), None)
            S = PowerSeriesRing(QQ, 'x', default_prec=l-offset)
            num = S(R(sequence)*den).add_bigoh(floor(ZZ(l)/2+1)).truncate()
            if num == 0 or sequence != S(num/den).list():
                return 0
            else:
                return CFiniteSequence(num / den)
"""
.. TODO::

    sage: CFiniteSequence(x+x^2+x^3+x^4+x^5+O(x^6)) # not implemented
    sage: latex(r)        # not implemented
    \big\{a_{n\ge0}\big|a_{n+2}=\sum_{i=0}^{1}c_ia_{n+i}, c=\{1,1\}, a_{n<2}=\{0,0,0,1\}\big\}
    sage: r.egf()      # not implemented
    exp(2*x)
    sage: r = CFiniteSequence(1/(1-y-x*y), x) # not implemented
"""