hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

use arb_dot in arb_poly series sin_cos and sinh_cosh + improve cutoffs

Browse source
This commit is contained in:
fredrik 2018-08-22 02:09:56 +02:00
parent 20579719ef
commit 1465e566dd
6 changed files with 55 additions and 37 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
arb_poly/exp_series.c
View file

@ -112,10 +112,6 @@ _arb_poly_exp_series(arb_ptr f, arb_srcptr h, slong hlen, slong n, slong prec)
_arb_vec_zero(f + j - d + 1, n - (j - d + 1)); _arb_vec_zero(f + j - d + 1, n - (j - d + 1));
arb_clear(t); arb_clear(t);
} }
else if (hlen <= arb_poly_newton_exp_cutoff)
{
_arb_poly_exp_series_basecase(f, h, hlen, n, prec);
}
else else
{ {
slong cutoff; slong cutoff;

19
arb_poly/sin_cos_pi_series.c
View file

@ -11,8 +11,6 @@
#include "arb_poly.h" #include "arb_poly.h"
#define TANGENT_CUTOFF 240
void void
_arb_poly_sin_cos_pi_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec) _arb_poly_sin_cos_pi_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)
{ {
@ -36,10 +34,25 @@ _arb_poly_sin_cos_pi_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slon
arb_mul(c + 1, s, t, prec); arb_mul(c + 1, s, t, prec);
arb_clear(t); arb_clear(t);
} }
else if (hlen < TANGENT_CUTOFF) else
{
slong cutoff;
if (prec <= 128)
{
cutoff = 1400;
}
else
{
cutoff = 100000 / pow(log(prec), 3);
cutoff = FLINT_MIN(cutoff, 700);
}
if (hlen < cutoff)
_arb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec, 1); _arb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec, 1);
else else
_arb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec, 1); _arb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec, 1);
}
} }
void void

19
arb_poly/sin_cos_series.c
View file

@ -11,8 +11,6 @@
#include "arb_poly.h" #include "arb_poly.h"
#define TANGENT_CUTOFF 240
void void
_arb_poly_sin_cos_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec) _arb_poly_sin_cos_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)
{ {
@ -35,10 +33,25 @@ _arb_poly_sin_cos_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n
arb_mul(c + 1, s, t, prec); arb_mul(c + 1, s, t, prec);
arb_clear(t); arb_clear(t);
} }
else if (hlen < TANGENT_CUTOFF) else
{
slong cutoff;
if (prec <= 128)
{
cutoff = 1400;
}
else
{
cutoff = 100000 / pow(log(prec), 3);
cutoff = FLINT_MIN(cutoff, 700);
}
if (hlen < cutoff)
_arb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec, 0); _arb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec, 0);
else else
_arb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec, 0); _arb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec, 0);
}
} }
void void

13
arb_poly/sin_cos_series_basecase.c
View file

@ -15,7 +15,7 @@ void
_arb_poly_sin_cos_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, _arb_poly_sin_cos_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen,
slong n, slong prec, int times_pi) slong n, slong prec, int times_pi)
{ {
slong j, k, alen = FLINT_MIN(n, hlen); slong k, alen = FLINT_MIN(n, hlen);
arb_ptr a; arb_ptr a;
arb_t t, u; arb_t t, u;
@ -46,15 +46,8 @@ _arb_poly_sin_cos_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen
for (k = 1; k < n; k++) for (k = 1; k < n; k++)
{ {
arb_zero(t); arb_dot(t, NULL, 1, a + 1, 1, s + k - 1, -1, FLINT_MIN(k, hlen - 1), prec);
arb_zero(u); arb_dot(u, NULL, 0, a + 1, 1, c + k - 1, -1, FLINT_MIN(k, hlen - 1), prec);
for (j = 1; j < FLINT_MIN(k + 1, hlen); j++)
{
arb_submul(t, a + j, s + k - j, prec);
arb_addmul(u, a + j, c + k - j, prec);
}
arb_div_ui(c + k, t, k, prec); arb_div_ui(c + k, t, k, prec);
arb_div_ui(s + k, u, k, prec); arb_div_ui(s + k, u, k, prec);
} }

12
arb_poly/sinh_cosh_series.c
View file

@ -32,10 +32,20 @@ _arb_poly_sinh_cosh_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong
arb_mul(c + 1, s, t, prec); arb_mul(c + 1, s, t, prec);
arb_clear(t); arb_clear(t);
} }
else if (hlen < 60 || n < 120) else
{
slong cutoff;
if (prec <= 128)
cutoff = 400;
else
cutoff = 30000 / pow(log(prec), 3);
if (hlen < cutoff)
_arb_poly_sinh_cosh_series_basecase(s, c, h, hlen, n, prec); _arb_poly_sinh_cosh_series_basecase(s, c, h, hlen, n, prec);
else else
_arb_poly_sinh_cosh_series_exponential(s, c, h, hlen, n, prec); _arb_poly_sinh_cosh_series_exponential(s, c, h, hlen, n, prec);
}
} }
void void

13
arb_poly/sinh_cosh_series_basecase.c
View file

@ -15,7 +15,7 @@ void
_arb_poly_sinh_cosh_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, _arb_poly_sinh_cosh_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen,
slong n, slong prec) slong n, slong prec)
{ {
slong j, k, alen = FLINT_MIN(n, hlen); slong k, alen = FLINT_MIN(n, hlen);
arb_ptr a; arb_ptr a;
arb_t t, u; arb_t t, u;
@ -37,15 +37,8 @@ _arb_poly_sinh_cosh_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, slong hl
for (k = 1; k < n; k++) for (k = 1; k < n; k++)
{ {
arb_zero(t); arb_dot(t, NULL, 0, a + 1, 1, s + k - 1, -1, FLINT_MIN(k, hlen - 1), prec);
arb_zero(u); arb_dot(u, NULL, 0, a + 1, 1, c + k - 1, -1, FLINT_MIN(k, hlen - 1), prec);
for (j = 1; j < FLINT_MIN(k + 1, hlen); j++)
{
arb_addmul(t, a + j, s + k - j, prec);
arb_addmul(u, a + j, c + k - j, prec);
}
arb_div_ui(c + k, t, k, prec); arb_div_ui(c + k, t, k, prec);
arb_div_ui(s + k, u, k, prec); arb_div_ui(s + k, u, k, prec);
} }