diff --git a/fmpcb_poly/rsqrt_series.c b/fmpcb_poly/rsqrt_series.c index 11b068eb..d7e49963 100644 --- a/fmpcb_poly/rsqrt_series.c +++ b/fmpcb_poly/rsqrt_series.c @@ -39,6 +39,13 @@ _fmpcb_poly_rsqrt_series(fmpcb_ptr g, { _fmpcb_vec_zero(g + 1, len - 1); } + else if (len == 2) + { + fmpcb_div(g + 1, h + 1, h, prec); + fmpcb_mul(g + 1, g + 1, g, prec); + fmpcb_mul_2exp_si(g + 1, g + 1, -1); + fmpcb_neg(g + 1, g + 1); + } else { fmpcb_ptr t, u; @@ -66,10 +73,10 @@ _fmpcb_poly_rsqrt_series(fmpcb_ptr g, void fmpcb_poly_rsqrt_series(fmpcb_poly_t g, const fmpcb_poly_t h, long n, long prec) { - if (n == 0 || h->length == 0) + if (n == 0) { - printf("fmpcb_poly_rsqrt_series: require n > 0 and nonzero input\n"); - abort(); + fmpcb_poly_zero(g); + return; } if (g == h) @@ -83,7 +90,10 @@ fmpcb_poly_rsqrt_series(fmpcb_poly_t g, const fmpcb_poly_t h, long n, long prec) } fmpcb_poly_fit_length(g, n); - _fmpcb_poly_rsqrt_series(g->coeffs, h->coeffs, h->length, n, prec); + if (h->length == 0) + _fmpcb_vec_indeterminate(g->coeffs, n); + else + _fmpcb_poly_rsqrt_series(g->coeffs, h->coeffs, h->length, n, prec); _fmpcb_poly_set_length(g, n); _fmpcb_poly_normalise(g); } diff --git a/fmpcb_poly/sqrt_series.c b/fmpcb_poly/sqrt_series.c index 777e7492..75051f9d 100644 --- a/fmpcb_poly/sqrt_series.c +++ b/fmpcb_poly/sqrt_series.c @@ -29,21 +29,37 @@ void _fmpcb_poly_sqrt_series(fmpcb_ptr g, fmpcb_srcptr h, long hlen, long len, long prec) { - fmpcb_ptr t; hlen = FLINT_MIN(hlen, len); - t = _fmpcb_vec_init(len); - _fmpcb_poly_rsqrt_series(t, h, hlen, len, prec); - _fmpcb_poly_mullow(g, t, len, h, hlen, len, prec); - _fmpcb_vec_clear(t, len); + + if (hlen == 1) + { + fmpcb_sqrt(g, h, prec); + _fmpcb_vec_zero(g + 1, len - 1); + } + else if (len == 2) + { + fmpcb_sqrt(g, h, prec); + fmpcb_div(g + 1, h + 1, h, prec); + fmpcb_mul(g + 1, g + 1, g, prec); + fmpcb_mul_2exp_si(g + 1, g + 1, -1); + } + else + { + fmpcb_ptr t; + t = _fmpcb_vec_init(len); + _fmpcb_poly_rsqrt_series(t, h, hlen, len, prec); + _fmpcb_poly_mullow(g, t, len, h, hlen, len, prec); + _fmpcb_vec_clear(t, len); + } } void fmpcb_poly_sqrt_series(fmpcb_poly_t g, const fmpcb_poly_t h, long n, long prec) { - if (n == 0 || h->length == 0) + if (n == 0) { - printf("fmpcb_poly_sqrt_series: require n > 0 and nonzero input\n"); - abort(); + fmpcb_poly_zero(g); + return; } if (g == h) @@ -57,7 +73,10 @@ fmpcb_poly_sqrt_series(fmpcb_poly_t g, const fmpcb_poly_t h, long n, long prec) } fmpcb_poly_fit_length(g, n); - _fmpcb_poly_sqrt_series(g->coeffs, h->coeffs, h->length, n, prec); + if (h->length == 0) + _fmpcb_vec_indeterminate(g->coeffs, n); + else + _fmpcb_poly_sqrt_series(g->coeffs, h->coeffs, h->length, n, prec); _fmpcb_poly_set_length(g, n); _fmpcb_poly_normalise(g); } diff --git a/fmprb_poly/rsqrt_series.c b/fmprb_poly/rsqrt_series.c index 8e698da5..341da7f8 100644 --- a/fmprb_poly/rsqrt_series.c +++ b/fmprb_poly/rsqrt_series.c @@ -30,12 +30,20 @@ _fmprb_poly_rsqrt_series(fmprb_ptr g, fmprb_srcptr h, long hlen, long len, long prec) { hlen = FLINT_MIN(hlen, len); + fmprb_rsqrt(g, h, prec); if (hlen == 1) { _fmprb_vec_zero(g + 1, len - 1); } + else if (len == 2) + { + fmprb_div(g + 1, h + 1, h, prec); + fmprb_mul(g + 1, g + 1, g, prec); + fmprb_mul_2exp_si(g + 1, g + 1, -1); + fmprb_neg(g + 1, g + 1); + } else { fmprb_ptr t, u; @@ -63,10 +71,10 @@ _fmprb_poly_rsqrt_series(fmprb_ptr g, void fmprb_poly_rsqrt_series(fmprb_poly_t g, const fmprb_poly_t h, long n, long prec) { - if (n == 0 || h->length == 0) + if (n == 0) { - printf("fmprb_poly_rsqrt_series: require n > 0 and nonzero input\n"); - abort(); + fmprb_poly_zero(g); + return; } if (g == h) @@ -80,7 +88,10 @@ fmprb_poly_rsqrt_series(fmprb_poly_t g, const fmprb_poly_t h, long n, long prec) } fmprb_poly_fit_length(g, n); - _fmprb_poly_rsqrt_series(g->coeffs, h->coeffs, h->length, n, prec); + if (h->length == 0) + _fmprb_vec_indeterminate(g->coeffs, n); + else + _fmprb_poly_rsqrt_series(g->coeffs, h->coeffs, h->length, n, prec); _fmprb_poly_set_length(g, n); _fmprb_poly_normalise(g); } diff --git a/fmprb_poly/sqrt_series.c b/fmprb_poly/sqrt_series.c index b68e126e..17c2ec7e 100644 --- a/fmprb_poly/sqrt_series.c +++ b/fmprb_poly/sqrt_series.c @@ -29,21 +29,37 @@ void _fmprb_poly_sqrt_series(fmprb_ptr g, fmprb_srcptr h, long hlen, long len, long prec) { - fmprb_ptr t; hlen = FLINT_MIN(hlen, len); - t = _fmprb_vec_init(len); - _fmprb_poly_rsqrt_series(t, h, hlen, len, prec); - _fmprb_poly_mullow(g, t, len, h, hlen, len, prec); - _fmprb_vec_clear(t, len); + + if (hlen == 1) + { + fmprb_sqrt(g, h, prec); + _fmprb_vec_zero(g + 1, len - 1); + } + else if (len == 2) + { + fmprb_sqrt(g, h, prec); + fmprb_div(g + 1, h + 1, h, prec); + fmprb_mul(g + 1, g + 1, g, prec); + fmprb_mul_2exp_si(g + 1, g + 1, -1); + } + else + { + fmprb_ptr t; + t = _fmprb_vec_init(len); + _fmprb_poly_rsqrt_series(t, h, hlen, len, prec); + _fmprb_poly_mullow(g, t, len, h, hlen, len, prec); + _fmprb_vec_clear(t, len); + } } void fmprb_poly_sqrt_series(fmprb_poly_t g, const fmprb_poly_t h, long n, long prec) { - if (n == 0 || h->length == 0) + if (n == 0) { - printf("fmprb_poly_sqrt_series: require n > 0 and nonzero input\n"); - abort(); + fmprb_poly_zero(g); + return; } if (g == h) @@ -57,7 +73,10 @@ fmprb_poly_sqrt_series(fmprb_poly_t g, const fmprb_poly_t h, long n, long prec) } fmprb_poly_fit_length(g, n); - _fmprb_poly_sqrt_series(g->coeffs, h->coeffs, h->length, n, prec); + if (h->length == 0) + _fmprb_vec_indeterminate(g->coeffs, n); + else + _fmprb_poly_sqrt_series(g->coeffs, h->coeffs, h->length, n, prec); _fmprb_poly_set_length(g, n); _fmprb_poly_normalise(g); }