Merge pull request #2663 from fredrik-johansson/series

Improve sine and cosine of power series

Author: fredrik-johansson

doc/source/gr_poly.rst

@@ -1165,10 +1165,34 @@ Power series special functions
               int _gr_poly_exp_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
               int gr_poly_exp_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
 
+.. function:: int _gr_poly_sin_cos_series(gr_ptr s, gr_ptr c, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int _gr_poly_sin_cos_pi_series(gr_ptr s, gr_ptr c, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_sin_cos_series(gr_poly_t s, gr_poly_t c, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int gr_poly_sin_cos_pi_series(gr_poly_t s, gr_poly_t c, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_sin_series(gr_ptr s, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_sin_series(gr_poly_t s, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_sin_pi_series(gr_ptr s, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_sin_pi_series(gr_poly_t s, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_cos_series(gr_ptr c, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_cos_series(gr_poly_t c, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_cos_pi_series(gr_ptr c, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_cos_pi_series(gr_poly_t c, const gr_poly_t h, slong n, gr_ctx_t ctx);
+
+    Compute `s = \sin(h)`, `c = \cos(h)` as power series truncated to length `m`,
+    or `s = \sin(\pi h)`, `c = \cos(\pi h)` for the ``pi`` variants.
+    The underscore methods allow aliasing.
+
 .. function:: int _gr_poly_sin_cos_series_basecase(gr_ptr s, gr_ptr c, gr_srcptr h, slong hlen, slong n, int times_pi, gr_ctx_t ctx)
               int gr_poly_sin_cos_series_basecase(gr_poly_t s, gr_poly_t c, const gr_poly_t h, slong n, int times_pi, gr_ctx_t ctx)
               int _gr_poly_sin_cos_series_tangent(gr_ptr s, gr_ptr c, gr_srcptr h, slong hlen, slong n, int times_pi, gr_ctx_t ctx)
               int gr_poly_sin_cos_series_tangent(gr_poly_t s, gr_poly_t c, const gr_poly_t h, slong n, int times_pi, gr_ctx_t ctx)
+              int _gr_poly_sin_cos_series_newton(gr_ptr s, gr_ptr c, gr_srcptr h, slong hlen, slong n, slong cutoff, int times_pi, gr_ctx_t ctx)
+              int gr_poly_sin_cos_series_newton(gr_poly_t s, gr_poly_t c, const gr_poly_t h, slong n, slong cutoff, int times_pi, gr_ctx_t ctx)
+
+    Various algorithms to compute sine and cosine of power series.
+    The default functions (:func:`_gr_poly_sin_cos_series` et al.)
+    choose an algorithm automatically.
+    The *times_pi* flag specifies that the input is to be multiplied by `\pi`.
 
     The *basecase* version uses a simple recurrence for the coefficients,
     requiring `O(nm)` operations where `m` is the length of `h`.
@@ -1181,8 +1205,10 @@ Power series special functions
     `\sin(h_0 + h_1) = \cos(h_0) \sin(h_1) + \sin(h_0) \cos(h_1)`,
     `\cos(h_0 + h_1) = \cos(h_0) \cos(h_1) - \sin(h_0) \sin(h_1)`.
 
-    The *basecase* and *tangent* versions take a flag *times_pi*
-    specifying that the input is to be multiplied by `\pi`.
+    The *newton* version uses Newton iteration for `\exp(ih)`. The complex
+    parts are represented formally; complex arithmetic is not required.
+    The *cutoff* parameter specifies the cutoff for using the basecase
+    algorithm.
 
 .. function:: int _gr_poly_tan_series_basecase(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
               int gr_poly_tan_series_basecase(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)

doc/source/gr_series.rst

@@ -117,10 +117,16 @@ directly.
               int gr_series_mod_mul(gr_poly_t res, const gr_poly_t x, const gr_poly_t y, gr_ctx_t ctx)
               int gr_series_mod_inv(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
               int gr_series_mod_div(gr_poly_t res, const gr_poly_t x, const gr_poly_t y, gr_ctx_t ctx)
-              int gr_series_mod_exp(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
-              int gr_series_mod_log(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
               int gr_series_mod_sqrt(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
               int gr_series_mod_rsqrt(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_exp(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_log(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_sin(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_cos(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_sin_pi(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_cos_pi(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_sin_cos(gr_poly_t res1, gr_poly_t res2, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_sin_cos_pi(gr_poly_t res1, gr_poly_t res2, const gr_poly_t x, gr_ctx_t ctx)
               int gr_series_mod_tan(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
               int gr_series_mod_asin(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
               int gr_series_mod_acos(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
@@ -247,6 +253,12 @@ directly.
               int gr_series_rsqrt(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
               int gr_series_exp(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
               int gr_series_log(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_sin(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_cos(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_sin_pi(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_cos_pi(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_sin_cos(gr_series_t res1, gr_series_t res2, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_sin_cos_pi(gr_series_t res1, gr_series_t res2, const gr_series_t x, gr_ctx_t ctx)
               int gr_series_tan(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
               int gr_series_asin(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
               int gr_series_acos(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)

doc/source/nmod_poly.rst

@@ -2431,23 +2431,19 @@ of polynomial multiplication.
 
     Set `g = \operatorname{asinh}(h) + O(x^n)`.
 
-.. function:: void _nmod_poly_sin_series(nn_ptr g, nn_srcptr h, slong n, nmod_t mod)
+.. function:: void _nmod_poly_sin_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)
 
-    Set `g = \operatorname{sin}(h) + O(x^n)`. Assumes `n > 0` and that `h`
-    is zero-padded as necessary to length `n`. Aliasing of `g` and `h` is
-    allowed. The value is computed using the identity
-    `\sin(x) = 2 \tan(x/2)) / (1 + \tan^2(x/2)).`
+    Set `g = \operatorname{sin}(h) + O(x^n)`. Assumes `n > 0` and `hlen > 0`.
+    Aliasing of `g` and `h` is allowed.
 
 .. function:: void nmod_poly_sin_series(nmod_poly_t g, const nmod_poly_t h, slong n)
 
     Set `g = \operatorname{sin}(h) + O(x^n)`.
 
-.. function:: void _nmod_poly_cos_series(nn_ptr g, nn_srcptr h, slong n, nmod_t mod)
+.. function:: void _nmod_poly_cos_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)
 
-    Set `g = \operatorname{cos}(h) + O(x^n)`. Assumes `n > 0` and that `h`
-    is zero-padded as necessary to length `n`. Aliasing of `g` and `h` is
-    allowed. The value is computed using the identity
-    `\cos(x) = (1-\tan^2(x/2)) / (1 + \tan^2(x/2)).`
+    Set `g = \operatorname{cos}(h) + O(x^n)`. Assumes `n > 0` and `hlen > 0`.
+    Aliasing of `g` and `h` is allowed.`
 
 .. function:: void nmod_poly_cos_series(nmod_poly_t g, const nmod_poly_t h, slong n)