diff --git a/src/sage/data_structures/stream.py b/src/sage/data_structures/stream.py index 8838a4f1754..c67c2b3de04 100644 --- a/src/sage/data_structures/stream.py +++ b/src/sage/data_structures/stream.py @@ -596,7 +596,6 @@ class Stream_exact(Stream): The convention for ``order`` is different to the one in :class:`sage.rings.lazy_series_ring.LazySeriesRing`, where the input is shifted to have the prescribed order. - """ def __init__(self, initial_coefficients, constant=None, degree=None, order=None): """ @@ -1038,6 +1037,175 @@ def __eq__(self, other): return isinstance(other, type(self)) and self.get_coefficient == other.get_coefficient +class Stream_taylor(Stream_inexact): + r""" + Class that returns a stream for the Taylor series of a function. + + INPUT: + + - ``function`` -- a function that has a ``derivative`` method and takes + a single input + - ``is_sparse`` -- boolean; specifies whether the stream is sparse + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_taylor + sage: g(x) = sin(x) + sage: f = Stream_taylor(g, False) + sage: f[3] + -1/6 + sage: [f[i] for i in range(10)] + [0, 1, 0, -1/6, 0, 1/120, 0, -1/5040, 0, 1/362880] + + sage: g(y) = cos(y) + sage: f = Stream_taylor(g, False) + sage: n = f.iterate_coefficients() + sage: [next(n) for _ in range(10)] + [1, 0, -1/2, 0, 1/24, 0, -1/720, 0, 1/40320, 0] + + sage: g(z) = 1 / (1 - 2*z) + sage: f = Stream_taylor(g, True) + sage: [f[i] for i in range(4)] + [1, 2, 4, 8] + """ + def __init__(self, function, is_sparse): + """ + Initialize. + + TESTS:: + + sage: from sage.data_structures.stream import Stream_taylor + sage: f = Stream_taylor(polygen(QQ, 'x')^3, False) + sage: TestSuite(f).run(skip="_test_pickling") + """ + from sage.symbolic.ring import SR + from sage.structure.element import parent + if parent(function) is SR: + self._is_symbolic = True + if function.number_of_arguments() != 1: + raise NotImplementedError("the function can only take a single input") + self._arg = function.arguments()[0] + else: + self._is_symbolic = False + self._func = function + super().__init__(is_sparse, False) + self._approximate_order = 0 + + def __hash__(self): + """ + Return the hash of ``self``. + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_taylor + sage: x = polygen(QQ, 'x') + sage: f = Stream_taylor(x^3 + x - 2, True) + sage: g = Stream_taylor(x^3 + x - 2, False) + sage: hash(f) == hash(g) + True + """ + # We don't hash the function as it might not be hashable. + return hash(type(self)) + + def __eq__(self, other): + r""" + Return whether ``self`` and ``other`` are known to be equal. + + INPUT: + + - ``other`` -- a stream + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_taylor + sage: x = polygen(QQ, 'x') + sage: fun = x^2 + sage: f = Stream_taylor(x^2, True) + sage: g = Stream_taylor(x^2, False) + sage: h = Stream_taylor((x^3 + x^2) / (x + 1), False) + sage: f == g + True + sage: f == h + True + + sage: F(y) = cos(y)^2 + sin(y)^2 + y + sage: G(y) = y + 1 + sage: f = Stream_taylor(F, True) + sage: g = Stream_taylor(G, False) + sage: f == g + True + """ + # The bool call is needed when the functions are in SR + return isinstance(other, type(self)) and bool(self._func == other._func) + + def get_coefficient(self, n): + r""" + Return the ``n``-th coefficient of ``self``. + + INPUT: + + - ``n`` -- integer; the degree for the coefficient + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_taylor + sage: g(x) = exp(x) + sage: f = Stream_taylor(g, True) + sage: f.get_coefficient(5) + 1/120 + + sage: from sage.data_structures.stream import Stream_taylor + sage: y = SR.var('y') + sage: f = Stream_taylor(sin(y), True) + sage: f.get_coefficient(0) + 0 + sage: f.get_coefficient(5) + 1/120 + """ + if n == 0: + if self._is_symbolic: + return self._func.subs({self._arg: ZZ.zero()}) + return self._func(ZZ.zero()) + + from sage.functions.other import factorial + if self._is_symbolic: + num = self._func.derivative(n).subs({self._arg: ZZ.zero()}) + else: + num = self._func.derivative(n)(ZZ.zero()) + return num / factorial(n) + + def iterate_coefficients(self): + r""" + A generator for the coefficients of ``self``. + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_taylor + sage: x = polygen(QQ, 'x') + sage: f = Stream_taylor(x^3, False) + sage: it = f.iterate_coefficients() + sage: [next(it) for _ in range(10)] + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0] + + sage: y = SR.var('y') + sage: f = Stream_taylor(y^3, False) + sage: it = f.iterate_coefficients() + sage: [next(it) for _ in range(10)] + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0] + """ + cur = self._func + n = 0 + denom = 1 + while True: + if self._is_symbolic: + yield cur({self._arg: ZZ.zero()}) / denom + else: + yield cur(ZZ.zero()) / denom + cur = cur.derivative() + n += 1 + denom *= n + + class Stream_uninitialized(Stream_inexact): r""" Coefficient stream for an uninitialized series. @@ -1051,7 +1219,7 @@ class Stream_uninitialized(Stream_inexact): .. TODO:: - shouldn't instances of this class share the cache with its + Should instances of this class share the cache with its ``_target``? EXAMPLES:: @@ -1247,7 +1415,6 @@ class Stream_binary(Stream_inexact): sage: [h[i] for i in range(10)] [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] """ - def __init__(self, left, right, is_sparse): """ Initialize ``self``. @@ -3371,7 +3538,7 @@ def __getitem__(self, n): sage: [f2[i] for i in range(-1, 4)] [0, 2, 6, 12, 20] """ - return (prod(n + k for k in range(1, self._shift + 1)) + return (ZZ.prod(range(n + 1, n + self._shift + 1)) * self._series[n + self._shift]) def __hash__(self): @@ -3433,6 +3600,141 @@ def is_nonzero(self): return self._series.is_nonzero() +class Stream_integral(Stream_unary): + """ + Operator for taking integrals of a non-exact stream. + + INPUT: + + - ``series`` -- a :class:`Stream` + - ``integration_constants`` -- a list of integration constants + - ``is_sparse`` -- boolean + """ + def __init__(self, series, integration_constants, is_sparse): + """ + Initialize ``self``. + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_exact, Stream_integral + sage: f = Stream_exact([1, 2, 3]) + sage: f2 = Stream_integral(f, [-2], True) + sage: TestSuite(f2).run() + """ + self._shift = len(integration_constants) + self._integration_constants = tuple(integration_constants) + super().__init__(series, is_sparse, False) + + @lazy_attribute + def _approximate_order(self): + """ + Compute and return the approximate order of ``self``. + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_function, Stream_integral + sage: f = Stream_function(lambda n: (n+1)*(n+2), True, 2) + sage: h = Stream_integral(f, [0, 0], True) + sage: h._approximate_order + 0 + sage: [h[i] for i in range(10)] + [0, 0, 0, 0, 1, 1, 1, 1, 1, 1] + sage: h._approximate_order + 4 + """ + # this is generally not the true order + return min(self._series._approximate_order + self._shift, 0) + + def get_coefficient(self, n): + r""" + Return the ``n``-th coefficient of ``self``. + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_function, Stream_integral + sage: f = Stream_function(lambda n: n + 1, True, -3) + sage: [f[i] for i in range(-3, 4)] + [-2, -1, 0, 1, 2, 3, 4] + sage: f2 = Stream_integral(f, [0], True) + sage: [f2.get_coefficient(i) for i in range(-3, 5)] + [0, 1, 1, 0, 1, 1, 1, 1] + + sage: f = Stream_function(lambda n: (n + 1)*(n+2), True, 2) + sage: [f[i] for i in range(-1, 4)] + [0, 0, 0, 12, 20] + sage: f2 = Stream_integral(f, [-1, -1, -1], True) + sage: [f2.get_coefficient(i) for i in range(-1, 7)] + [0, -1, -1, -1/2, 0, 0, 1/5, 1/6] + """ + if 0 <= n < self._shift: + return (self._integration_constants[n] / ZZ.prod(range(2, n + 1))) + return (self._series[n - self._shift] / + ZZ.prod(range(n - self._shift + 1, n + 1))) + + def __hash__(self): + """ + Return the hash of ``self``. + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_function, Stream_integral + sage: a = Stream_function(lambda n: 2*n, False, 1) + sage: f = Stream_integral(a, [0, 1], True) + sage: g = Stream_integral(a, [0, 2], True) + sage: hash(f) == hash(f) + True + sage: hash(f) == hash(g) + False + """ + return hash((type(self), self._series, self._integration_constants)) + + def __eq__(self, other): + """ + Return whether ``self`` and ``other`` are known to be equal. + + INPUT: + + - ``other`` -- a stream + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_function, Stream_integral + sage: a = Stream_function(lambda n: 2*n, False, 1) + sage: f = Stream_integral(a, [1], True) + sage: g = Stream_integral(a, [2], True) + sage: f == g + False + sage: f == Stream_integral(a, [1], True) + True + """ + return (isinstance(other, type(self)) + and self._integration_constants == other._integration_constants + and self._series == other._series) + + def is_nonzero(self): + r""" + Return ``True`` if and only if this stream is known + to be non-zero. + + EXAMPLES:: + + sage: from sage.data_structures.stream import Stream_function, Stream_integral + sage: f = Stream_function(lambda n: 2*n, True, 1) + sage: f[1] + 2 + sage: f.is_nonzero() + True + sage: Stream_integral(f, [0], True).is_nonzero() + True + sage: f = Stream_function(lambda n: 0, False, 1) + sage: Stream_integral(f, [0, 0, 0], False).is_nonzero() + False + sage: Stream_integral(f, [0, 2], False).is_nonzero() + True + """ + return self._series.is_nonzero() or any(self._integration_constants) + + class Stream_infinite_operator(Stream): r""" Stream defined by applying an operator an infinite number of times. diff --git a/src/sage/rings/lazy_series.py b/src/sage/rings/lazy_series.py index 5fdfb64cea4..16ffe2a2f10 100644 --- a/src/sage/rings/lazy_series.py +++ b/src/sage/rings/lazy_series.py @@ -257,6 +257,7 @@ Stream_truncated, Stream_function, Stream_derivative, + Stream_integral, Stream_dirichlet_convolve, Stream_dirichlet_invert, Stream_plethysm @@ -4533,6 +4534,172 @@ def derivative(self, *args): P.is_sparse()) return P.element_class(P, coeff_stream) + def integral(self, variable=None, *, constants=None): + r""" + Return the integral of ``self`` with respect to ``variable``. + + INPUT: + + - ``variable`` -- (optional) the variable to integrate + - ``constants`` -- (optional; keyword-only) list of integration + constants for the integrals of ``self`` (the last constant + corresponds to the first integral) + + If the first argument is a list, then this method iterprets it as + integration constants. If it is a positive integer, the method + interprets it as the number of times to integrate the function. + If ``variable`` is not the variable of the Laurent series, then + the coefficients are integrated with respect to ``variable``. + + If the integration constants are not specified, they are considered + to be `0`. + + EXAMPLES:: + + sage: L. = LazyLaurentSeriesRing(QQ) + sage: f = t^-3 + 2 + 3*t + t^5 + sage: f.integral() + -1/2*t^-2 + 2*t + 3/2*t^2 + 1/6*t^6 + sage: f.integral([-2, -2]) + 1/2*t^-1 - 2 - 2*t + t^2 + 1/2*t^3 + 1/42*t^7 + sage: f.integral(t) + -1/2*t^-2 + 2*t + 3/2*t^2 + 1/6*t^6 + sage: f.integral(2) + 1/2*t^-1 + t^2 + 1/2*t^3 + 1/42*t^7 + sage: L.zero().integral() + 0 + sage: L.zero().integral([0, 1, 2, 3]) + t + t^2 + 1/2*t^3 + + We solve the ODE `f' = a f` by integrating both sides and + the recursive definition:: + + sage: R. = QQ[] + sage: L. = LazyLaurentSeriesRing(R) + sage: f = L.undefined(0) + sage: f.define((a*f).integral(constants=[C])) + sage: f + C + a*C*x + 1/2*a^2*C*x^2 + 1/6*a^3*C*x^3 + 1/24*a^4*C*x^4 + + 1/120*a^5*C*x^5 + 1/720*a^6*C*x^6 + O(x^7) + sage: C * exp(a*x) + C + a*C*x + 1/2*a^2*C*x^2 + 1/6*a^3*C*x^3 + 1/24*a^4*C*x^4 + + 1/120*a^5*C*x^5 + 1/720*a^6*C*x^6 + O(x^7) + + We can integrate both the series and coefficients:: + + sage: R. = QQ[] + sage: L. = LazyLaurentSeriesRing(R) + sage: f = (x*t^2 + y*t^-2 + z)^2; f + y^2*t^-4 + 2*y*z*t^-2 + (2*x*y + z^2) + 2*x*z*t^2 + x^2*t^4 + sage: f.integral(x) + x*y^2*t^-4 + 2*x*y*z*t^-2 + (x^2*y + x*z^2) + x^2*z*t^2 + 1/3*x^3*t^4 + sage: f.integral(t) + -1/3*y^2*t^-3 - 2*y*z*t^-1 + (2*x*y + z^2)*t + 2/3*x*z*t^3 + 1/5*x^2*t^5 + sage: f.integral(y, constants=[x*y*z]) + -1/9*y^3*t^-3 - y^2*z*t^-1 + x*y*z + (x*y^2 + y*z^2)*t + 2/3*x*y*z*t^3 + 1/5*x^2*y*t^5 + + TESTS:: + + sage: L. = LazyLaurentSeriesRing(QQ) + sage: f = t^-2 + sage: f.integral(t, constants=[0, 0, 0]) + Traceback (most recent call last): + ... + ValueError: cannot integrate 3 times the series t^-2 + sage: f = t^-5 + t^-2 + sage: f.integral(3) + Traceback (most recent call last): + ... + ValueError: cannot integrate 3 times the series t^-5 + t^-2 + sage: f.integral([0, 1], constants=[0, 1]) + Traceback (most recent call last): + ... + ValueError: integration constants given twice + sage: f.integral(4, constants=[0, 1]) + Traceback (most recent call last): + ... + ValueError: the number of integrations does not match the number of integration constants + """ + P = self.parent() + zero = P.base_ring().zero() + if variable is None: + if constants is None: + constants = [zero] + elif variable != P.gen(): + if isinstance(variable, (list, tuple)): + if constants is not None: + raise ValueError("integration constants given twice") + constants = tuple(variable) + variable = None + elif variable in ZZ and ZZ(variable) >= 0: + if constants is None: + constants = [zero] * ZZ(variable) + elif ZZ(variable) != len(constants): + raise ValueError("the number of integrations does not match" + " the number of integration constants") + variable = None + if constants is None: + constants = [] + else: + if constants is None: + constants = [zero] + variable = None + + nints = len(constants) + + coeff_stream = self._coeff_stream + if isinstance(coeff_stream, Stream_zero): + if any(constants): + coeff_stream = Stream_exact([c / ZZ.prod(k for k in range(1, i+1)) + for i, c in enumerate(constants)], + order=0, + constant=zero) + return P.element_class(P, coeff_stream) + return self + + if (isinstance(coeff_stream, Stream_exact) and not coeff_stream._constant): + coeffs = [c / ZZ.prod(k for k in range(1, i+1)) + for i, c in enumerate(constants)] + if coeff_stream._approximate_order < 0: + ic = coeff_stream._initial_coefficients + ao = coeff_stream._approximate_order + if nints > -ao or any(ic[-ao-nints:-ao]): + raise ValueError(f"cannot integrate {nints} times the series {self}") + if variable is not None: + coeffs = [c.integral(variable) / ZZ.prod(i+k for k in range(1, nints+1)) + for i, c in enumerate(ic[:-ao-nints], ao)] + coeffs + else: + coeffs = [c / ZZ.prod(i+k for k in range(1, nints+1)) + for i, c in enumerate(ic[:-ao-nints], ao)] + coeffs + + ic = ic[-ao:] + val = ao + nints + ao = 0 + else: + coeffs += [zero] * coeff_stream._approximate_order + ic = coeff_stream._initial_coefficients + val = 0 + ao = coeff_stream._approximate_order + if variable: + coeffs += [c.integral(variable) / ZZ.prod(i+k for k in range(1, nints+1)) + for i, c in enumerate(ic, ao)] + else: + coeffs += [c / ZZ.prod(i+k for k in range(1, nints+1)) + for i, c in enumerate(ic, ao)] + if not any(coeffs): + return P.zero() + coeff_stream = Stream_exact(coeffs, order=val, constant=zero) + return P.element_class(P, coeff_stream) + + if nints: + coeff_stream = Stream_integral(coeff_stream, constants, P.is_sparse()) + + if variable is not None: + coeff_stream = Stream_map_coefficients(coeff_stream, + lambda c: c.integral(variable), + P.is_sparse()) + return P.element_class(P, coeff_stream) + def approximate_series(self, prec, name=None): r""" Return the Laurent series with absolute precision ``prec`` approximated @@ -5301,14 +5468,18 @@ def derivative(self, *args): sage: T. = LazyPowerSeriesRing(ZZ) sage: z.derivative() 1 - sage: (1+z+z^2).derivative(3) + sage: (1 + z + z^2).derivative(3) 0 - sage: (1/(1-z)).derivative() + sage: (z^2 + z^4 + z^10).derivative(3) + 24*z + 720*z^7 + sage: (1 / (1-z)).derivative() 1 + 2*z + 3*z^2 + 4*z^3 + 5*z^4 + 6*z^5 + 7*z^6 + O(z^7) + sage: T([1, 1, 1], constant=4).derivative() + 1 + 2*z + 12*z^2 + 16*z^3 + 20*z^4 + 24*z^5 + 28*z^6 + O(z^7) sage: R. = QQ[] sage: L. = LazyPowerSeriesRing(R) - sage: f = 1/(1-q*x+y); f + sage: f = 1 / (1-q*x+y); f 1 + (q*x-y) + (q^2*x^2+(-2*q)*x*y+y^2) + (q^3*x^3+(-3*q^2)*x^2*y+3*q*x*y^2-y^3) + (q^4*x^4+(-4*q^3)*x^3*y+6*q^2*x^2*y^2+(-4*q)*x*y^3+y^4) @@ -5322,6 +5493,53 @@ def derivative(self, *args): + (6*q^5*x^6+(-30*q^4)*x^5*y+60*q^3*x^4*y^2+(-60*q^2)*x^3*y^3+30*q*x^2*y^4+(-6)*x*y^5) + O(x,y)^7 + Multivariate:: + + sage: L. = LazyPowerSeriesRing(QQ) + sage: f = (x + y^2 + z)^3; f + (x^3+3*x^2*z+3*x*z^2+z^3) + (3*x^2*y^2+6*x*y^2*z+3*y^2*z^2) + (3*x*y^4+3*y^4*z) + y^6 + sage: f.derivative(x) + (3*x^2+6*x*z+3*z^2) + (6*x*y^2+6*y^2*z) + 3*y^4 + sage: f.derivative(y, 5) + 720*y + sage: f.derivative(z, 5) + 0 + sage: f.derivative(x, y, z) + 12*y + + sage: f = (1 + x + y^2 + z)^-1 + sage: f.derivative(x) + -1 + (2*x+2*z) + (-3*x^2+2*y^2-6*x*z-3*z^2) + ... + O(x,y,z)^6 + sage: f.derivative(y, 2) + -2 + (4*x+4*z) + (-6*x^2+12*y^2-12*x*z-6*z^2) + ... + O(x,y,z)^5 + sage: f.derivative(x, y) + 4*y + (-12*x*y-12*y*z) + (24*x^2*y-12*y^3+48*x*y*z+24*y*z^2) + + (-40*x^3*y+48*x*y^3-120*x^2*y*z+48*y^3*z-120*x*y*z^2-40*y*z^3) + O(x,y,z)^5 + sage: f.derivative(x, y, z) + (-12*y) + (48*x*y+48*y*z) + (-120*x^2*y+48*y^3-240*x*y*z-120*y*z^2) + O(x,y,z)^4 + + sage: R. = QQ[] + sage: L. = LazyPowerSeriesRing(R) + sage: f = ((t^2-3)*x + t*y^2 - t*z)^2 + sage: f.derivative(t,x,t,y) + 24*t*y + sage: f.derivative(t, 2) + ((12*t^2-12)*x^2+(-12*t)*x*z+2*z^2) + (12*t*x*y^2+(-4)*y^2*z) + 2*y^4 + sage: f.derivative(z, t) + ((-6*t^2+6)*x+4*t*z) + ((-4*t)*y^2) + sage: f.derivative(t, 10) + 0 + + sage: f = (1 + t*(x + y + z))^-1 + sage: f.derivative(x, t, y) + 4*t + ((-18*t^2)*x+(-18*t^2)*y+(-18*t^2)*z) + + (48*t^3*x^2+96*t^3*x*y+48*t^3*y^2+96*t^3*x*z+96*t^3*y*z+48*t^3*z^2) + + ... + O(x,y,z)^5 + sage: f.derivative(t, 2) + (2*x^2+4*x*y+2*y^2+4*x*z+4*y*z+2*z^2) + ... + O(x,y,z)^7 + sage: f.derivative(x, y, z, t) + (-18*t^2) + (96*t^3*x+96*t^3*y+96*t^3*z) + ... + O(x,y,z)^4 + TESTS: Check that :issue:`36154` is fixed:: @@ -5341,7 +5559,7 @@ def derivative(self, *args): if x is None: order += 1 elif x in V: - gen_vars.append(x) + gen_vars.append(x._coeff_stream[1]) else: vars.append(x) @@ -5357,6 +5575,16 @@ def derivative(self, *args): if P._arity > 1: v = gen_vars + vars d = -len(gen_vars) + + if isinstance(coeff_stream, Stream_exact): # the constant should be 0 + ao = coeff_stream._approximate_order + val = max(ao + d, 0) + coeffs = [R(c).derivative(v) for c in coeff_stream._initial_coefficients[val-(ao+d):]] + if any(coeffs): + coeff_stream = Stream_exact(coeffs, order=val, constant=coeff_stream._constant) + return P.element_class(P, coeff_stream) + return P.zero() + coeff_stream = Stream_map_coefficients(coeff_stream, lambda c: R(c).derivative(v), P.is_sparse()) @@ -5390,6 +5618,248 @@ def derivative(self, *args): P.is_sparse()) return P.element_class(P, coeff_stream) + def integral(self, variable=None, *, constants=None): + r""" + Return the integral of ``self`` with respect to ``variable``. + + INPUT: + + - ``variable`` -- (optional) the variable to integrate + - ``constants`` -- (optional; keyword-only) list of integration + constants for the integrals of ``self`` (the last constant + corresponds to the first integral) + + For multivariable series, only ``variable`` should be + specified; the integration constant is taken to be `0`. + + Now we assume the series is univariate. If the first argument is a + list, then this method iterprets it as integration constants. If it + is a positive integer, the method interprets it as the number of times + to integrate the function. If ``variable`` is not the variable of + the power series, then the coefficients are integrated with respect + to ``variable``. If the integration constants are not specified, + they are considered to be `0`. + + EXAMPLES:: + + sage: L. = LazyPowerSeriesRing(QQ) + sage: f = 2 + 3*t + t^5 + sage: f.integral() + 2*t + 3/2*t^2 + 1/6*t^6 + sage: f.integral([-2, -2]) + -2 - 2*t + t^2 + 1/2*t^3 + 1/42*t^7 + sage: f.integral(t) + 2*t + 3/2*t^2 + 1/6*t^6 + sage: f.integral(2) + t^2 + 1/2*t^3 + 1/42*t^7 + sage: (t^3 + t^5).integral() + 1/4*t^4 + 1/6*t^6 + sage: L.zero().integral() + 0 + sage: L.zero().integral([0, 1, 2, 3]) + t + t^2 + 1/2*t^3 + sage: L([1, 2 ,3], constant=4).integral() + t + t^2 + t^3 + t^4 + 4/5*t^5 + 2/3*t^6 + O(t^7) + + We solve the ODE `f'' - f' - 2 f = 0` by solving for `f''`, then + integrating and applying a recursive definition:: + + sage: R. = QQ[] + sage: L. = LazyPowerSeriesRing(R) + sage: f = L.undefined() + sage: f.define((f.derivative() + 2*f).integral(constants=[C, D])) + sage: f + C + D*x + ((C+1/2*D)*x^2) + ((1/3*C+1/2*D)*x^3) + + ((1/4*C+5/24*D)*x^4) + ((1/12*C+11/120*D)*x^5) + + ((11/360*C+7/240*D)*x^6) + O(x^7) + sage: f.derivative(2) - f.derivative() - 2*f + O(x^7) + + We compare this with the answer we get from the + characteristic polynomial:: + + sage: g = C * exp(-x) + D * exp(2*x); g + (C+D) + ((-C+2*D)*x) + ((1/2*C+2*D)*x^2) + ((-1/6*C+4/3*D)*x^3) + + ((1/24*C+2/3*D)*x^4) + ((-1/120*C+4/15*D)*x^5) + + ((1/720*C+4/45*D)*x^6) + O(x^7) + sage: g.derivative(2) - g.derivative() - 2*g + O(x^7) + + Note that ``C`` and ``D`` are playing different roles, so we need + to perform a substitution to the coefficients of ``f`` to recover + the solution ``g``:: + + sage: fp = f.map_coefficients(lambda c: c(C=C+D, D=2*D-C)); fp + (C+D) + ((-C+2*D)*x) + ((1/2*C+2*D)*x^2) + ((-1/6*C+4/3*D)*x^3) + + ((1/24*C+2/3*D)*x^4) + ((-1/120*C+4/15*D)*x^5) + + ((1/720*C+4/45*D)*x^6) + O(x^7) + sage: fp - g + O(x^7) + + We can integrate both the series and coefficients:: + + sage: R. = QQ[] + sage: L. = LazyPowerSeriesRing(R) + sage: f = (x*t^2 + y*t + z)^2; f + z^2 + 2*y*z*t + ((y^2+2*x*z)*t^2) + 2*x*y*t^3 + x^2*t^4 + sage: f.integral(x) + x*z^2 + 2*x*y*z*t + ((x*y^2+x^2*z)*t^2) + x^2*y*t^3 + 1/3*x^3*t^4 + sage: f.integral(t) + z^2*t + y*z*t^2 + ((1/3*y^2+2/3*x*z)*t^3) + 1/2*x*y*t^4 + 1/5*x^2*t^5 + sage: f.integral(y, constants=[x*y*z]) + x*y*z + y*z^2*t + 1/2*y^2*z*t^2 + ((1/9*y^3+2/3*x*y*z)*t^3) + 1/4*x*y^2*t^4 + 1/5*x^2*y*t^5 + + We can integrate multivariate power series:: + + sage: R. = QQ[] + sage: L. = LazyPowerSeriesRing(R) + sage: f = ((t^2 + t) - t * y^2 + t^2 * (y + z))^2; f + (t^4+2*t^3+t^2) + ((2*t^4+2*t^3)*y+(2*t^4+2*t^3)*z) + + ((t^4-2*t^3-2*t^2)*y^2+2*t^4*y*z+t^4*z^2) + + ((-2*t^3)*y^3+(-2*t^3)*y^2*z) + t^2*y^4 + sage: g = f.integral(x); g + ((t^4+2*t^3+t^2)*x) + ((2*t^4+2*t^3)*x*y+(2*t^4+2*t^3)*x*z) + + ((t^4-2*t^3-2*t^2)*x*y^2+2*t^4*x*y*z+t^4*x*z^2) + + ((-2*t^3)*x*y^3+(-2*t^3)*x*y^2*z) + t^2*x*y^4 + sage: g[0] + 0 + sage: g[1] + (t^4 + 2*t^3 + t^2)*x + sage: g[2] + (2*t^4 + 2*t^3)*x*y + (2*t^4 + 2*t^3)*x*z + sage: f.integral(z) + ((t^4+2*t^3+t^2)*z) + ((2*t^4+2*t^3)*y*z+(t^4+t^3)*z^2) + + ((t^4-2*t^3-2*t^2)*y^2*z+t^4*y*z^2+1/3*t^4*z^3) + + ((-2*t^3)*y^3*z+(-t^3)*y^2*z^2) + t^2*y^4*z + sage: f.integral(t) + (1/5*t^5+1/2*t^4+1/3*t^3) + ((2/5*t^5+1/2*t^4)*y+(2/5*t^5+1/2*t^4)*z) + + ((1/5*t^5-1/2*t^4-2/3*t^3)*y^2+2/5*t^5*y*z+1/5*t^5*z^2) + + ((-1/2*t^4)*y^3+(-1/2*t^4)*y^2*z) + 1/3*t^3*y^4 + + sage: L. = LazyPowerSeriesRing(QQ) + sage: (x + y - z^2).integral(z) + (x*z+y*z) + (-1/3*z^3) + + TESTS:: + + sage: L. = LazyPowerSeriesRing(QQ) + sage: f = t^2 + sage: f.integral([0, 1], constants=[0, 1]) + Traceback (most recent call last): + ... + ValueError: integration constants given twice + sage: f.integral(4, constants=[0, 1]) + Traceback (most recent call last): + ... + ValueError: the number of integrations does not match the number of integration constants + + sage: L. = LazyPowerSeriesRing(QQ) + sage: x.integral(y, constants=[2]) + Traceback (most recent call last): + ... + ValueError: integration constants must not be given for multivariate series + sage: x.integral() + Traceback (most recent call last): + ... + ValueError: the integration variable must be specified + """ + P = self.parent() + coeff_stream = self._coeff_stream + R = P._laurent_poly_ring + + if P._arity > 1: + if constants is not None: + raise ValueError("integration constants must not be given for multivariate series") + if variable is None: + raise ValueError("the integration variable must be specified") + + if isinstance(coeff_stream, Stream_zero): + return self + + if variable in P.gens(): + variable = variable._coeff_stream[1] + shift = 1 + else: + shift = 0 + + if isinstance(coeff_stream, Stream_exact): # constant is 0 because arity is at least 2 + ao = coeff_stream._approximate_order + coeffs = [R(c).integral(variable) for c in coeff_stream._initial_coefficients] + coeff_stream = Stream_exact(coeffs, order=ao+shift, constant=coeff_stream._constant) + return P.element_class(P, coeff_stream) + + coeff_stream = Stream_map_coefficients(coeff_stream, + lambda c: c.integral(variable), + P.is_sparse()) + if shift: + coeff_stream = Stream_shift(coeff_stream, 1) + return P.element_class(P, coeff_stream) + + # the univariate case + + zero = P.base_ring().zero() + # This is copied from the LazyLaurentSeries.integral + if variable is None: + if constants is None: + constants = [zero] + elif variable != P.gen(): + if isinstance(variable, (list, tuple)): + if constants is not None: + raise ValueError("integration constants given twice") + constants = tuple(variable) + variable = None + elif variable in ZZ and ZZ(variable) >= 0: + if constants is None: + constants = [zero] * ZZ(variable) + elif ZZ(variable) != len(constants): + raise ValueError("the number of integrations does not match" + " the number of integration constants") + variable = None + if constants is None: + constants = [] + else: + if constants is None: + constants = [zero] + variable = None + + nints = len(constants) + + if isinstance(coeff_stream, Stream_zero): + if any(constants): + coeff_stream = Stream_exact([c / ZZ.prod(k for k in range(1, i+1)) + for i, c in enumerate(constants)], + order=0, + constant=zero) + return P.element_class(P, coeff_stream) + + return self + + if (isinstance(coeff_stream, Stream_exact) and not coeff_stream._constant): + coeffs = [c / ZZ.prod(k for k in range(1, i+1)) + for i, c in enumerate(constants)] + coeffs += [zero] * coeff_stream._approximate_order + ic = coeff_stream._initial_coefficients + ao = coeff_stream._approximate_order + if variable: + coeffs += [c.integral(variable) / ZZ.prod(i+k for k in range(1, nints+1)) + for i, c in enumerate(ic, ao)] + else: + coeffs += [c / ZZ.prod(i+k for k in range(1, nints+1)) + for i, c in enumerate(ic, ao)] + if not any(coeffs): + return P.zero() + coeff_stream = Stream_exact(coeffs, order=0, constant=zero) + return P.element_class(P, coeff_stream) + + if nints: + coeff_stream = Stream_integral(coeff_stream, constants, P.is_sparse()) + + if variable is not None: + coeff_stream = Stream_map_coefficients(coeff_stream, + lambda c: c.integral(variable), + P.is_sparse()) + return P.element_class(P, coeff_stream) + def _format_series(self, formatter, format_strings=False): """ Return nonzero ``self`` formatted by ``formatter``. diff --git a/src/sage/rings/lazy_series_ring.py b/src/sage/rings/lazy_series_ring.py index 07e08938696..ba444712df8 100644 --- a/src/sage/rings/lazy_series_ring.py +++ b/src/sage/rings/lazy_series_ring.py @@ -87,7 +87,8 @@ Stream_function, Stream_iterator, Stream_exact, - Stream_uninitialized + Stream_uninitialized, + Stream_taylor ) from types import GeneratorType @@ -1792,6 +1793,43 @@ def residue_field(self): raise TypeError("the base ring is not a field") return R + def taylor(self, f): + r""" + Return the Taylor expansion around `0` of the function ``f``. + + INPUT: + + - ``f`` -- a function such that one of the following works: + + * the substitution `f(z)`, where `z` is a generator of ``self`` + * `f` is a function of a single variable with no poles at `0` + and has a ``derivative`` method + + EXAMPLES:: + + sage: L. = LazyLaurentSeriesRing(QQ) + sage: x = SR.var('x') + sage: f(x) = (1 + x) / (1 - x^2) + sage: L.taylor(f) + 1 + z + z^2 + z^3 + z^4 + z^5 + z^6 + O(z^7) + + For inputs as symbolic functions/expressions, the function must + not have any poles at `0`:: + + sage: f(x) = (1 + x^2) / sin(x^2) + sage: L.taylor(f) + + sage: def g(a): return (1 + a^2) / sin(a^2) + sage: L.taylor(g) + z^-2 + 1 + 1/6*z^2 + 1/6*z^4 + O(z^5) + """ + try: + return f(self.gen()) + except (ValueError, TypeError): + pass + stream = Stream_taylor(f, self.is_sparse()) + return self.element_class(self, stream) + # === special functions === def q_pochhammer(self, q=None): @@ -2503,6 +2541,87 @@ def some_elements(self): elts.extend([(z-3)*(2+z)**2, (1 - 2*z**3)/(1 - z + 3*z**2), self(lambda n: sum_gens**n)]) return elts + def taylor(self, f): + r""" + Return the Taylor expansion around `0` of the function ``f``. + + INPUT: + + - ``f`` -- a function such that one of the following works: + + * the substitution `f(z_1, \ldots, z_n)`, where `(z_1, \ldots, z_n)` + are the generators of ``self`` + * `f` is a function with no poles at `0` and has a ``derivative`` + method + + .. WARNING:: + + For inputs as symbolic functions/expressions, this does not check + that the function does not have poles at `0`. + + EXAMPLES:: + + sage: L. = LazyPowerSeriesRing(QQ) + sage: x = SR.var('x') + sage: f(x) = (1 + x) / (1 - x^3) + sage: L.taylor(f) + 1 + z + z^3 + z^4 + z^6 + O(z^7) + sage: (1 + z) / (1 - z^3) + 1 + z + z^3 + z^4 + z^6 + O(z^7) + sage: f(x) = cos(x + pi/2) + sage: L.taylor(f) + -z + 1/6*z^3 - 1/120*z^5 + O(z^7) + + For inputs as symbolic functions/expressions, the function must + not have any poles at `0`:: + + sage: L. = LazyPowerSeriesRing(QQ, sparse=True) + sage: f = 1 / sin(x) + sage: L.taylor(f) + + + Different multivariate inputs:: + + sage: L. = LazyPowerSeriesRing(QQ) + sage: def f(x, y): return (1 + x) / (1 + y) + sage: L.taylor(f) + 1 + (a-b) + (-a*b+b^2) + (a*b^2-b^3) + (-a*b^3+b^4) + (a*b^4-b^5) + (-a*b^5+b^6) + O(a,b)^7 + sage: g(w, z) = (1 + w) / (1 + z) + sage: L.taylor(g) + 1 + (a-b) + (-a*b+b^2) + (a*b^2-b^3) + (-a*b^3+b^4) + (a*b^4-b^5) + (-a*b^5+b^6) + O(a,b)^7 + sage: y = SR.var('y') + sage: h = (1 + x) / (1 + y) + sage: L.taylor(h) + 1 + (a-b) + (-a*b+b^2) + (a*b^2-b^3) + (-a*b^3+b^4) + (a*b^4-b^5) + (-a*b^5+b^6) + O(a,b)^7 + """ + try: + return f(*self.gens()) + except (ValueError, TypeError): + pass + + if self._arity != 1: + R = self._laurent_poly_ring + BR = R.base_ring() + args = f.arguments() + subs = {str(va): ZZ.zero() for va in args} + gens = R.gens() + ell = len(subs) + from sage.combinat.integer_vector import integer_vectors_nk_fast_iter + from sage.arith.misc import factorial + + def taylor_expand(deg): + if deg == 0: + return BR(f(**subs)) + return R.sum(BR(f.diff(*sum(([g] * e for g, e in zip(args, al)), []))(**subs) + / ZZ.prod(factorial(a) for a in al)) + * R.monomial(*al) for al in integer_vectors_nk_fast_iter(deg, ell)) + + coeff_stream = Stream_function(taylor_expand, self._sparse, self._minimal_valuation) + else: + coeff_stream = Stream_taylor(f, self._sparse) + return self.element_class(self, coeff_stream) + + ######################################################################