Repository: sympy/sympy

Issue:
Better way to compose rational functions 
I have a degree 17/degree 17 rational function in one variable with floating point coefficients with 1000 digits of precision. I would like to compose this with the Möbius transformation `(t - c)/(t + c)`. Here `c = 51/5` but in general it will be a similar fraction. I want the final rational function in expanded form (expanded polynomial/expanded polynomial)

The problem is that if you try just about any naive and even not so naive way of doing this, you lose a ton of precision as soon as you expand the polynomials. I finally figured out the following method for doing it without losing precision. 

Write `(t - c)/(t + c)` as `1 - 2*c/(t + c)`. Then the composition corresponds to a shift by 1, a composition with `-2*c*t`, an inversion (replacing `t` with `1/t`), and a shift by c. For the inversion we can make use of the fact that `P(1/t) == t**d*P'(t)`, where `d = deg(P)` and `P'` is `P` with the terms reversed. The `t**d` will cancel in the numerator and denominator because the numerator and denominator have the same degree. In other words, we just need to reverse the terms. 

The result is something like this

```py
mpmath.mp.dps = 1000 # This is needed to make Poly use full precision. I'll open a separate issue. 

frac = map(Poly, fraction(D))
for i in range(len(frac)):
    # Shift by 1
    frac[i] = frac[i].shift(1)
    # Compose with -2*c*t
    frac[i] = frac[i].compose(Poly(2*c*t))
    # Invert
    # XXX: Is there a better way than rep.rep
    frac[i] = Poly(reversed(frac[i].rep.rep), t)
    # Shift by c
    frac[i] = frac[i].shift(c)

n, d = frac
```

The resulting rational function, from what I can tell, does not lose precision. I believe the key here is that `shift` uses the Taylor shift algorithm. 

It would be nice to have some more direct support for this somehow. 

Discussion:

Patch:
diff --git a/README.rst b/README.rst
index 303131a83a..79cd0a82bb 100644
--- a/README.rst
+++ b/README.rst
@@ -223,34 +223,18 @@ will show the top developers from the last year.
 Citation
 --------
 
-To cite SymPy in publications use
-
-    Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A,
-    Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP,
-    Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR,
-    Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy:
-    symbolic computing in Python. *PeerJ Computer Science* 3:e103
-    https://doi.org/10.7717/peerj-cs.103
-
-A BibTeX entry for LaTeX users is
-
-.. code-block:: none
-
-    @article{10.7717/peerj-cs.103,
-     title = {SymPy: symbolic computing in Python},
-     author = {Meurer, Aaron and Smith, Christopher P. and Paprocki, Mateusz and Čertík, Ondřej and Kirpichev, Sergey B. and Rocklin, Matthew and Kumar, AMiT and Ivanov, Sergiu and Moore, Jason K. and Singh, Sartaj and Rathnayake, Thilina and Vig, Sean and Granger, Brian E. and Muller, Richard P. and Bonazzi, Francesco and Gupta, Harsh and Vats, Shivam and Johansson, Fredrik and Pedregosa, Fabian and Curry, Matthew J. and Terrel, Andy R. and Roučka, Štěpán and Saboo, Ashutosh and Fernando, Isuru and Kulal, Sumith and Cimrman, Robert and Scopatz, Anthony},
-     year = 2017,
-     month = jan,
-     keywords = {Python, Computer algebra system, Symbolics},
-     abstract = {
-                SymPy is an open source computer algebra system written in pure Python. It is built with a focus on extensibility and ease of use, through both interactive and programmatic applications. These characteristics have led SymPy to become a popular symbolic library for the scientific Python ecosystem. This paper presents the architecture of SymPy, a description of its features, and a discussion of select submodules. The supplementary material provide additional examples and further outline details of the architecture and features of SymPy.
-             },
-     volume = 3,
-     pages = {e103},
-     journal = {PeerJ Computer Science},
-     issn = {2376-5992},
-     url = {https://doi.org/10.7717/peerj-cs.103},
-     doi = {10.7717/peerj-cs.103}
+To cite SymPy in publications use::
+
+    SymPy Development Team (2016). SymPy: Python library for symbolic mathematics
+    URL http://www.sympy.org.
+
+A BibTeX entry for LaTeX users is::
+
+    @Manual{,
+    title = {SymPy: Python library for symbolic mathematics},
+    author = {{SymPy Development Team}},
+    year = {2016},
+    url = {http://www.sympy.org},
     }
 
 SymPy is BSD licensed, so you are free to use it whatever you like, be it
diff --git a/doc/src/citing.rst b/doc/src/citing.rst
deleted file mode 100644
index ec6a62031d..0000000000
--- a/doc/src/citing.rst
+++ /dev/null
@@ -1,40 +0,0 @@
-Citing SymPy
-============
-
-To cite SymPy in publications use
-
-    Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A,
-    Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP,
-    Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR,
-    Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy:
-    symbolic computing in Python. *PeerJ Computer Science* 3:e103
-    https://doi.org/10.7717/peerj-cs.103
-
-A BibTeX entry for LaTeX users is
-
-.. code-block:: none
-
-    @article{10.7717/peerj-cs.103,
-     title = {SymPy: symbolic computing in Python},
-     author = {Meurer, Aaron and Smith, Christopher P. and Paprocki, Mateusz and Čertík, Ondřej and Kirpichev, Sergey B. and Rocklin, Matthew and Kumar, AMiT and Ivanov, Sergiu and Moore, Jason K. and Singh, Sartaj and Rathnayake, Thilina and Vig, Sean and Granger, Brian E. and Muller, Richard P. and Bonazzi, Francesco and Gupta, Harsh and Vats, Shivam and Johansson, Fredrik and Pedregosa, Fabian and Curry, Matthew J. and Terrel, Andy R. and Roučka, Štěpán and Saboo, Ashutosh and Fernando, Isuru and Kulal, Sumith and Cimrman, Robert and Scopatz, Anthony},
-     year = 2017,
-     month = jan,
-     keywords = {Python, Computer algebra system, Symbolics},
-     abstract = {
-                SymPy is an open source computer algebra system written in pure Python. It is built with a focus on extensibility and ease of use, through both interactive and programmatic applications. These characteristics have led SymPy to become a popular symbolic library for the scientific Python ecosystem. This paper presents the architecture of SymPy, a description of its features, and a discussion of select submodules. The supplementary material provide additional examples and further outline details of the architecture and features of SymPy.
-             },
-     volume = 3,
-     pages = {e103},
-     journal = {PeerJ Computer Science},
-     issn = {2376-5992},
-     url = {https://doi.org/10.7717/peerj-cs.103},
-     doi = {10.7717/peerj-cs.103}
-    }
-
-SymPy is BSD licensed, so you are free to use it whatever you like, be it
-academic, commercial, creating forks or derivatives, as long as you copy the
-BSD statement if you redistribute it (see the LICENSE file for details).  That
-said, although not required by the SymPy license, if it is convenient for you,
-please cite SymPy when using it in your work and also consider contributing
-all your changes back, so that we can incorporate it and all of us will
-benefit in the end.
diff --git a/doc/src/index.rst b/doc/src/index.rst
index f695b3f608..fcbdd8690b 100644
--- a/doc/src/index.rst
+++ b/doc/src/index.rst
@@ -23,9 +23,8 @@ Contents:
    gotchas.rst
    guide.rst
    modules/index.rst
-   special_topics/index.rst
    python-comparisons.rst
    wiki.rst
    outreach.rst
    aboutus.rst
-   citing.rst
+   special_topics/index.rst
diff --git a/sympy/geometry/ellipse.py b/sympy/geometry/ellipse.py
index e8f0b4fd12..5a1b215f05 100644
--- a/sympy/geometry/ellipse.py
+++ b/sympy/geometry/ellipse.py
@@ -581,7 +581,7 @@ def foci(self):
 
     @property
     def focus_distance(self):
-        """The focal distance of the ellipse.
+        """The focale distance of the ellipse.
 
         The distance between the center and one focus.
 
@@ -1324,7 +1324,7 @@ class Circle(Ellipse):
     >>> c1.hradius, c1.vradius, c1.radius
     (5, 5, 5)
 
-    >>> # a circle constructed from three points
+    >>> # a circle costructed from three points
     >>> c2 = Circle(Point(0, 0), Point(1, 1), Point(1, 0))
     >>> c2.hradius, c2.vradius, c2.radius, c2.center
     (sqrt(2)/2, sqrt(2)/2, sqrt(2)/2, Point2D(1/2, 1/2))
diff --git a/sympy/polys/polyclasses.py b/sympy/polys/polyclasses.py
index 96417efcf0..b853060a3e 100644
--- a/sympy/polys/polyclasses.py
+++ b/sympy/polys/polyclasses.py
@@ -99,6 +99,7 @@ def _perify_factors(per, result, include):
     dmp_compose,
     dup_decompose,
     dup_shift,
+    dup_transform,
     dmp_lift)
 
 from sympy.polys.euclidtools import (
@@ -721,6 +722,20 @@ def shift(f, a):
         else:
             raise ValueError('univariate polynomial expected')
 
+    def transform(f, p, q):
+        """Evaluate functional transformation ``q**n * f(p/q)``."""
+        if f.lev:
+            raise ValueError('univariate polynomial expected')
+
+        lev, dom, per, P, Q = p.unify(q)
+        lev, dom, per, F, P = f.unify(per(P, dom, lev))
+        lev, dom, per, F, Q = per(F, dom, lev).unify(per(Q, dom, lev))
+
+        if not lev:
+            return per(dup_transform(F, P, Q, dom))
+        else:
+            raise ValueError('univariate polynomial expected')
+
     def sturm(f):
         """Computes the Sturm sequence of ``f``. """
         if not f.lev:
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index 5f7f175948..dbab1ecb72 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -2944,6 +2944,32 @@ def shift(f, a):
 
         return f.per(result)
 
+    def transform(f, p, q):
+        """
+        Efficiently evaluate the functional transformation ``q**n * f(p/q)``.
+
+
+        Examples
+        ========
+
+        >>> from sympy import Poly
+        >>> from sympy.abc import x
+
+        >>> Poly(x**2 - 2*x + 1, x).transform(Poly(x + 1, x), Poly(x - 1, x))
+        Poly(4, x, domain='ZZ')
+
+        """
+        P, Q = p.unify(q)
+        F, P = f.unify(P)
+        F, Q = F.unify(Q)
+
+        if hasattr(F.rep, 'transform'):
+            result = F.rep.transform(P.rep, Q.rep)
+        else:  # pragma: no cover
+            raise OperationNotSupported(F, 'transform')
+
+        return F.per(result)
+
     def sturm(self, auto=True):
         """
         Computes the Sturm sequence of ``f``.
diff --git a/sympy/series/fourier.py b/sympy/series/fourier.py
index 4a2a70a556..2196cea554 100644
--- a/sympy/series/fourier.py
+++ b/sympy/series/fourier.py
@@ -10,7 +10,7 @@
 from sympy.core.symbol import Dummy, Symbol
 from sympy.core.compatibility import is_sequence
 from sympy.core.containers import Tuple
-from sympy.functions.elementary.trigonometric import sin, cos, sinc
+from sympy.functions.elementary.trigonometric import sin, cos
 from sympy.sets.sets import Interval
 from sympy.series.series_class import SeriesBase
 from sympy.series.sequences import SeqFormula
@@ -154,34 +154,9 @@ def _eval_subs(self, old, new):
             return self
 
     def truncate(self, n=3):
-        """
-        Return the first n nonzero terms of the series.
-
-        If n is None return an iterator.
-
-        Parameters
-        ==========
-        n : int or None
-            Amount of non-zero terms in approximation or None.
-
-        Returns
-        =======
-        Expr or iterator
-            Approximation of function expanded into Fourier series.
+        """Returns the first n (non-zero)terms of the series.
 
-        Examples
-        ========
-
-        >>> from sympy import fourier_series, pi
-        >>> from sympy.abc import x
-        >>> s = fourier_series(x, (x, -pi, pi))
-        >>> s.truncate(4)
-        2*sin(x) - sin(2*x) + 2*sin(3*x)/3 - sin(4*x)/2
-
-        See Also
-        ========
-
-        sympy.series.fourier.FourierSeries.sigma_approximation
+        If n is ``None`` returns an iterator.
         """
         if n is None:
             return iter(self)
@@ -195,71 +170,6 @@ def truncate(self, n=3):
 
         return Add(*terms)
 
-    def sigma_approximation(self, n=3):
-        """
-        Return :math:`\sigma`-approximation of Fourier series with respect
-        to order n.
-
-        Sigma approximation adjusts a Fourier summation to eliminate the Gibbs
-        phenomenon which would otherwise occur at discontinuities.
-        A sigma-approximated summation for a Fourier series of a T-periodical
-        function can be written as
-
-        .. math::
-            s(\\theta) = \\frac{1}{2} a_0 + \\sum _{k=1}^{m-1}
-            \\operatorname{sinc} \\Bigl( \\frac{k}{m} \\Bigr) \\cdot
-            \\left[ a_k \\cos \\Bigl( \\frac{2\\pi k}{T} \\theta \\Bigr)
-            + b_k \\sin \\Bigl( \\frac{2\\pi k}{T} \\theta \\Bigr) \\right],
-
-        where :math:`a_0, a_k, b_k, k=1,\ldots,{m-1}` are standard Fourier
-        series coefficients and
-        :math:`\\operatorname{sinc} \\Bigl( \\frac{k}{m} \\Bigr)` is a Lanczos
-        :math:`\sigma` factor (expressed in terms of normalized
-        :math:`\operatorname{sinc}` function).
-
-        Parameters
-        ==========
-        n : int
-            Highest order of the terms taken into account in approximation.
-
-        Returns
-        =======
-        Expr
-            Sigma approximation of function expanded into Fourier series.
-
-        Examples
-        ========
-
-        >>> from sympy import fourier_series, pi
-        >>> from sympy.abc import x
-        >>> s = fourier_series(x, (x, -pi, pi))
-        >>> s.sigma_approximation(4)
-        2*sin(x)*sinc(pi/4) - 2*sin(2*x)/pi + 2*sin(3*x)*sinc(3*pi/4)/3
-
-        See Also
-        ========
-
-        sympy.series.fourier.FourierSeries.truncate
-
-        Notes
-        =====
-
-        The behaviour of
-        :meth:`~sympy.series.fourier.FourierSeries.sigma_approximation`
-        is different from :meth:`~sympy.series.fourier.FourierSeries.truncate`
-        - it takes all nonzero terms of degree smaller than n, rather than
-        first n nonzero ones.
-
-        References
-        ==========
-
-        .. [1] https://en.wikipedia.org/wiki/Gibbs_phenomenon
-        .. [2] https://en.wikipedia.org/wiki/Sigma_approximation
-        """
-        terms = [sinc(pi * i / n) * t for i, t in enumerate(self[:n])
-                 if t is not S.Zero]
-        return Add(*terms)
-
     def shift(self, s):
         """Shift the function by a term independent of x.
 


Test patch:
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index 6db0499a24..727ad8ddb8 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -2119,6 +2119,35 @@ def test_compose():
 def test_shift():
     assert Poly(x**2 - 2*x + 1, x).shift(2) == Poly(x**2 + 2*x + 1, x)
 
+def test_transform():
+    # Also test that 3-way unification is done correctly
+    assert Poly(x**2 - 2*x + 1, x).transform(Poly(x + 1), Poly(x - 1)) == \
+        Poly(4, x) == \
+        cancel((x - 1)**2*(x**2 - 2*x + 1).subs(x, (x + 1)/(x - 1)))
+
+    assert Poly(x**2 - x/2 + 1, x).transform(Poly(x + 1), Poly(x - 1)) == \
+        Poly(3*x**2/2 + S(5)/2, x) == \
+        cancel((x - 1)**2*(x**2 - x/2 + 1).subs(x, (x + 1)/(x - 1)))
+
+    assert Poly(x**2 - 2*x + 1, x).transform(Poly(x + S(1)/2), Poly(x - 1)) == \
+        Poly(S(9)/4, x) == \
+        cancel((x - 1)**2*(x**2 - 2*x + 1).subs(x, (x + S(1)/2)/(x - 1)))
+
+    assert Poly(x**2 - 2*x + 1, x).transform(Poly(x + 1), Poly(x - S(1)/2)) == \
+        Poly(S(9)/4, x) == \
+        cancel((x - S(1)/2)**2*(x**2 - 2*x + 1).subs(x, (x + 1)/(x - S(1)/2)))
+
+    # Unify ZZ, QQ, and RR
+    assert Poly(x**2 - 2*x + 1, x).transform(Poly(x + 1.0), Poly(x - S(1)/2)) == \
+        Poly(9/4, x) == \
+        cancel((x - S(1)/2)**2*(x**2 - 2*x + 1).subs(x, (x + 1.0)/(x - S(1)/2)))
+
+    raises(ValueError, lambda: Poly(x*y).transform(Poly(x + 1), Poly(x - 1)))
+    raises(ValueError, lambda: Poly(x).transform(Poly(y + 1), Poly(x - 1)))
+    raises(ValueError, lambda: Poly(x).transform(Poly(x + 1), Poly(y - 1)))
+    raises(ValueError, lambda: Poly(x).transform(Poly(x*y + 1), Poly(x - 1)))
+    raises(ValueError, lambda: Poly(x).transform(Poly(x + 1), Poly(x*y - 1)))
+
 
 def test_sturm():
     f, F = x, Poly(x, domain='QQ')
diff --git a/sympy/series/tests/test_fourier.py b/sympy/series/tests/test_fourier.py
index 6ff4722b00..c1b3acb87e 100644
--- a/sympy/series/tests/test_fourier.py
+++ b/sympy/series/tests/test_fourier.py
@@ -1,4 +1,4 @@
-from sympy import (symbols, pi, Piecewise, sin, cos, sinc, Rational,
+from sympy import (symbols, pi, Piecewise, sin, cos, Rational,
                    oo, fourier_series, Add)
 from sympy.series.fourier import FourierSeries
 from sympy.utilities.pytest import raises
@@ -65,17 +65,6 @@ def test_FourierSeries_2():
                             4*cos(pi*x / 2) / pi**2 + Rational(1, 2))
 
 
-def test_fourier_series_square_wave():
-    """Test if fourier_series approximates discontinuous function correctly."""
-    square_wave = Piecewise((1, x < pi), (-1, True))
-    s = fourier_series(square_wave, (x, 0, 2*pi))
-
-    assert s.truncate(3) == 4 / pi * sin(x) + 4 / (3 * pi) * sin(3 * x) + \
-        4 / (5 * pi) * sin(5 * x)
-    assert s.sigma_approximation(4) == 4 / pi * sin(x) * sinc(pi / 4) + \
-        4 / (3 * pi) * sin(3 * x) * sinc(3 * pi / 4)
-
-
 def test_FourierSeries__operations():
     fes = fe.scale(-1).shift(pi**2)
     assert fes.truncate() == 4*cos(x) - cos(2*x) + 2*pi**2 / 3