Merge pull request #2672 from fredrik-johansson/series

Overhaul tan-like functions of power series

Author: fredrik-johansson

doc/source/gr_poly.rst

@@ -1178,10 +1178,26 @@ Power series special functions
               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`,
+    Compute `s = \sin(h)`, `c = \cos(h)` as power series truncated to length `n`,
     or `s = \sin(\pi h)`, `c = \cos(\pi h)` for the ``pi`` variants.
     The underscore methods allow aliasing.
 
+.. function:: int _gr_poly_tan_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_tan_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_tanh_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_tanh_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_cot_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_cot_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_coth_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_coth_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_tan_pi_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_tan_pi_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+              int _gr_poly_cot_pi_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+              int gr_poly_cot_pi_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+
+    Compute the respective trigonometric and hyperbolic functions of power series
+    truncated to length `n`. 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)
@@ -1207,15 +1223,42 @@ Power series special functions
 
     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)
-              int _gr_poly_tan_series_newton(gr_ptr f, gr_srcptr h, slong hlen, slong n, slong cutoff, gr_ctx_t ctx)
-              int gr_poly_tan_series_newton(gr_poly_t f, const gr_poly_t h, slong n, slong cutoff, gr_ctx_t ctx)
-              int _gr_poly_tan_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
-              int gr_poly_tan_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+    The *cutoff* parameter specifies the initial length to compute with the
+    basecase algorithm.
+
+.. function:: int _gr_poly_tan_series_basecase(gr_ptr f, gr_srcptr h, slong hlen, slong n, int func, gr_ctx_t ctx)
+              int gr_poly_tan_series_basecase(gr_poly_t f, const gr_poly_t h, slong n, int func, gr_ctx_t ctx)
+              int _gr_poly_tan_series_newton(gr_ptr f, gr_srcptr h, slong hlen, slong n, slong cutoff, int func, gr_ctx_t ctx)
+              int gr_poly_tan_series_newton(gr_poly_t f, const gr_poly_t h, slong n, slong cutoff, int func, gr_ctx_t ctx)
+              int _gr_poly_tan_series_sine_cosine(gr_ptr f, gr_srcptr h, slong hlen, slong n, int func, gr_ctx_t ctx)
+              int gr_poly_tan_series_sine_cosine(gr_poly_t res, const gr_poly_t h, slong len, int func, gr_ctx_t ctx)
+              int _gr_poly_tan_series_exponential(gr_ptr f, gr_srcptr h, slong hlen, slong n, int func, gr_ctx_t ctx)
+              int gr_poly_tan_series_exponential(gr_poly_t res, const gr_poly_t h, slong len, int func, gr_ctx_t ctx)
+
+    Various algorithms to compute tangent-like functions of power series.
+    The *func* parameter specifies the function as follows:
+
+    * 0 - tan
+    * 1 - tanh
+    * 2 - cot
+    * 3 - coth
+    * 4 - tan_pi
+    * 5 - tanh_pi  (not currently implemented)
+    * 6 - cot_pi
+    * 7 - coth_pi  (not currently implemented)
+
+    The *basecase* algorithm uses an `O(n^2)` recurrence for the coefficients.
+    The *newton* algorithm uses Newton iteration to invert the integral
+    defining the corresponding inverse function; the *cutoff* parameter
+    specifies the initial length to compute with the basecase algorithm.
+    The *sine_cosine* algorithm computes the function as a quotient of
+    a sine and cosine; for the hyperbolic functions this algorithm requires that
+    the ring contains the imaginary unit.
+    The *exponential* algorithm computes the function as a quotient of
+    `\pm 1` plus an exponential, where the sign of the exponent is chosen
+    (if possible) so that the exponential is small. For the trigonometric
+    functions this algorithm requires that the ring contains the imaginary unit.
+
 
 Modular arithmetic and composition
 --------------------------------------------------------------------------------

doc/source/gr_series.rst

@@ -128,6 +128,11 @@ directly.
               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_tanh(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_cot(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_coth(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_tan_pi(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
+              int gr_series_mod_cot_pi(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)
               int gr_series_mod_atan(gr_poly_t res, const gr_poly_t x, gr_ctx_t ctx)
@@ -260,6 +265,11 @@ directly.
               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_tanh(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_cot(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_coth(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_tan_pi(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)
+              int gr_series_cot_pi(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)
               int gr_series_atan(gr_series_t res, const gr_series_t x, gr_ctx_t ctx)

doc/source/nmod_poly.rst

@@ -2451,9 +2451,7 @@ of polynomial multiplication.
 
 .. function:: void _nmod_poly_tan_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)
 
-    Set `g = \operatorname{tan}(h) + O(x^n)`. Assumes `n > 0` and that `h`
-    is zero-padded as necessary to length `n`. Aliasing of `g` and `h` is
-    not allowed. Uses Newton iteration to invert the atan function.
+    Set `g = \operatorname{tan}(h) + O(x^n)`. Assumes `n > 0` and `hlen > 0`.
 
 .. function:: void nmod_poly_tan_series(nmod_poly_t g, const nmod_poly_t h, slong n)
 
@@ -2480,11 +2478,9 @@ of polynomial multiplication.
 
     Set `g = \operatorname{cosh}(h) + O(x^n)`.
 
-.. function:: void _nmod_poly_tanh_series(nn_ptr g, nn_srcptr h, slong n, nmod_t mod)
+.. function:: void _nmod_poly_tanh_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)
 
-    Set `g = \operatorname{tanh}(h) + O(x^n)`. Assumes `n > 0` and that `h`
-    is zero-padded as necessary to length `n`. Uses the identity
-    `\tanh(x) = (e^{2x}-1)/(e^{2x}+1)`.
+    Set `g = \operatorname{tanh}(h) + O(x^n)`. Assumes `n > 0` and `hlen > 0`.
 
 .. function:: void nmod_poly_tanh_series(nmod_poly_t g, const nmod_poly_t h, slong n)
 

src/acb_poly/cot_pi_series.c

@@ -9,6 +9,7 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "gr_poly.h"
 #include "acb_poly.h"
 
 void
@@ -23,51 +24,10 @@ _acb_poly_cot_pi_series(acb_ptr g, acb_srcptr h, slong hlen, slong len, slong pr
     }
     else
     {
-        acb_ptr t, u;
-
-        t = _acb_vec_init(len);
-        u = _acb_vec_init(len);
-
-        if (arf_cmpabs_2exp_si(arb_midref(acb_imagref(h)), 0) < 0)
-        {
-            _acb_poly_sin_cos_pi_series(t, u, h, hlen, len, prec);
-            _acb_poly_div_series(g, u, len, t, len, len, prec);
-        }
-        else
-        {
-            _acb_vec_scalar_mul_2exp_si(t, h, hlen, 1);
-
-            if (arf_sgn(arb_midref(acb_imagref(h))) > 0)
-            {
-                acb_const_pi(u, prec);
-                acb_mul_onei(u, u);
-                _acb_vec_scalar_mul(t, t, hlen, u, prec);
-                _acb_poly_exp_series(t, t, hlen, len, prec);
-                acb_sub_ui(u, t, 1, prec);
-                _acb_vec_set(u + 1, t + 1, len - 1);
-                _acb_poly_div_series(g, t, len, u, len, len, prec);
-                _acb_vec_scalar_mul_2exp_si(g, g, len, 1);
-                acb_sub_ui(g, g, 1, prec);
-                _acb_vec_scalar_mul_onei(g, g, len);
-            }
-            else
-            {
-                acb_const_pi(u, prec);
-                acb_div_onei(u, u);
-                _acb_vec_scalar_mul(t, t, hlen, u, prec);
-                _acb_poly_exp_series(t, t, hlen, len, prec);
-                acb_sub_ui(u, t, 1, prec);
-                _acb_vec_set(u + 1, t + 1, len - 1);
-                _acb_poly_div_series(g, t, len, u, len, len, prec);
-                _acb_vec_scalar_mul_2exp_si(g, g, len, 1);
-                acb_sub_ui(g, g, 1, prec);
-                _acb_vec_scalar_mul_onei(g, g, len);
-                _acb_vec_neg(g, g, len);
-            }
-        }
-
-        _acb_vec_clear(t, len);
-        _acb_vec_clear(u, len);
+        gr_ctx_t ctx;
+        gr_ctx_init_complex_acb(ctx, prec);
+        if (_gr_poly_cot_pi_series(g, h, hlen, len, ctx) != GR_SUCCESS)
+            _acb_vec_indeterminate(g, len);
     }
 }
 

src/acb_poly/tan_series.c

@@ -12,17 +12,12 @@
 #include "acb_poly.h"
 #include "gr_poly.h"
 
-#define TAN_NEWTON_CUTOFF 20
-
 void
 _acb_poly_tan_series(acb_ptr res, acb_srcptr h, slong hlen, slong len, slong prec)
 {
     gr_ctx_t ctx;
     gr_ctx_init_complex_acb(ctx, prec);
-
-    hlen = FLINT_MIN(hlen, len);
-
-    if (_gr_poly_tan_series_newton(res, h, hlen, len, TAN_NEWTON_CUTOFF, ctx) != GR_SUCCESS)
+    if (_gr_poly_tan_series(res, h, hlen, len, ctx) != GR_SUCCESS)
         _acb_vec_indeterminate(res, len);
 }
 

src/arb_poly/cot_pi_series.c

@@ -9,6 +9,7 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "gr_poly.h"
 #include "arb_poly.h"
 
 void
@@ -23,16 +24,10 @@ _arb_poly_cot_pi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong pr
     }
     else
     {
-        arb_ptr t, u;
-
-        t = _arb_vec_init(len);
-        u = _arb_vec_init(len);
-
-        _arb_poly_sin_cos_pi_series(t, u, h, hlen, len, prec);
-        _arb_poly_div_series(g, u, len, t, len, len, prec);
-
-        _arb_vec_clear(t, len);
-        _arb_vec_clear(u, len);
+        gr_ctx_t ctx;
+        gr_ctx_init_real_arb(ctx, prec);
+        if (_gr_poly_cot_pi_series(g, h, hlen, len, ctx) != GR_SUCCESS)
+            _arb_vec_indeterminate(g, len);
     }
 }
 

src/arb_poly/tan_series.c

@@ -12,17 +12,12 @@
 #include "arb_poly.h"
 #include "gr_poly.h"
 
-#define TAN_NEWTON_CUTOFF 20
-
 void
 _arb_poly_tan_series(arb_ptr res, arb_srcptr h, slong hlen, slong len, slong prec)
 {
     gr_ctx_t ctx;
     gr_ctx_init_real_arb(ctx, prec);
-
-    hlen = FLINT_MIN(hlen, len);
-
-    if (_gr_poly_tan_series_newton(res, h, hlen, len, TAN_NEWTON_CUTOFF, ctx) != GR_SUCCESS)
+    if (_gr_poly_tan_series(res, h, hlen, len, ctx) != GR_SUCCESS)
         _arb_vec_indeterminate(res, len);
 }
 

src/gr/acb.c

@@ -1188,6 +1188,8 @@ DEF_FUNC(cos)
 DEF_FUNC(cos_pi)
 DEF_FUNC_SING(tan)
 DEF_FUNC_SING(cot)
+DEF_FUNC_SING(tan_pi)
+DEF_FUNC_SING(cot_pi)
 
 DEF_FUNC(sinc)
 DEF_FUNC(sinc_pi)
@@ -2343,7 +2345,9 @@ gr_method_tab_input _acb_methods_input[] =
     {GR_METHOD_COS_PI,          (gr_funcptr) _gr_acb_cos_pi},
     {GR_METHOD_SIN_COS_PI,      (gr_funcptr) _gr_acb_sin_cos_pi},
     {GR_METHOD_TAN,             (gr_funcptr) _gr_acb_tan},
+    {GR_METHOD_TAN_PI,          (gr_funcptr) _gr_acb_tan_pi},
     {GR_METHOD_COT,             (gr_funcptr) _gr_acb_cot},
+    {GR_METHOD_COT_PI,          (gr_funcptr) _gr_acb_cot_pi},
     {GR_METHOD_SINC,            (gr_funcptr) _gr_acb_sinc},
     {GR_METHOD_SINC_PI,         (gr_funcptr) _gr_acb_sinc_pi},
     {GR_METHOD_SINH,            (gr_funcptr) _gr_acb_sinh},

src/gr_generic/generic.c

@@ -2874,6 +2874,7 @@ const gr_method_tab_input _gr_generic_methods[] =
     {GR_METHOD_COS,                     (gr_funcptr) gr_generic_cos},
     {GR_METHOD_SIN_COS,                 (gr_funcptr) gr_generic_sin_cos},
     {GR_METHOD_TAN,                     (gr_funcptr) gr_generic_tan},
+    {GR_METHOD_TANH,                    (gr_funcptr) gr_generic_tanh},
     {GR_METHOD_ASIN,                    (gr_funcptr) gr_generic_asin},
     {GR_METHOD_ATAN,                    (gr_funcptr) gr_generic_atan},
     {GR_METHOD_ASINH,                   (gr_funcptr) gr_generic_asinh},

src/gr_poly.h

@@ -588,12 +588,28 @@ WARN_UNUSED_RESULT int gr_poly_sin_pi_series(gr_poly_t s, const gr_poly_t h, slo
 WARN_UNUSED_RESULT int gr_poly_cos_series(gr_poly_t c, const gr_poly_t h, slong n, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_poly_cos_pi_series(gr_poly_t c, const gr_poly_t h, slong n, gr_ctx_t ctx);
 
-WARN_UNUSED_RESULT int _gr_poly_tan_series_basecase(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx);
-WARN_UNUSED_RESULT int gr_poly_tan_series_basecase(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx);
-WARN_UNUSED_RESULT int _gr_poly_tan_series_newton(gr_ptr f, gr_srcptr h, slong hlen, slong n, slong cutoff, gr_ctx_t ctx);
-WARN_UNUSED_RESULT int gr_poly_tan_series_newton(gr_poly_t f, const gr_poly_t h, slong n, slong cutoff, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_tan_series_basecase(gr_ptr f, gr_srcptr h, slong hlen, slong n, int func, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_tan_series_basecase(gr_poly_t f, const gr_poly_t h, slong n, int func, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_tan_series_newton(gr_ptr f, gr_srcptr h, slong hlen, slong n, slong cutoff, int func, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_tan_series_newton(gr_poly_t f, const gr_poly_t h, slong n, slong cutoff, int func, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_tan_series_sine_cosine(gr_ptr f, gr_srcptr h, slong hlen, slong n, int func, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_tan_series_sine_cosine(gr_poly_t res, const gr_poly_t h, slong len, int func, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_tan_series_exponential(gr_ptr f, gr_srcptr h, slong hlen, slong n, int func, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_tan_series_exponential(gr_poly_t res, const gr_poly_t h, slong len, int func, gr_ctx_t ctx);
+
 WARN_UNUSED_RESULT int _gr_poly_tan_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_tanh_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_cot_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_coth_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_tan_pi_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_cot_pi_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_poly_tan_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_tanh_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_cot_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_coth_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_tan_pi_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_cot_pi_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx);
+
 
 /* Modular arithmetic and composition */