Commit 84d2fe34 authored by Julien Puydt's avatar Julien Puydt

New upstream version 2.15.0

parent 0167d0bb
......@@ -644,6 +644,16 @@ acb_submul_arb(acb_t z, const acb_t x, const arb_t y, slong prec)
arb_submul(acb_imagref(z), acb_imagref(x), y, prec);
}
void acb_dot_simple(acb_t res, const acb_t initial, int subtract,
acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec);
void acb_dot_precise(acb_t res, const acb_t initial, int subtract,
acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec);
void acb_dot(acb_t res, const acb_t initial, int subtract,
acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec);
void acb_approx_dot(acb_t res, const acb_t initial, int subtract,
acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec);
void acb_inv(acb_t z, const acb_t x, slong prec);
void acb_div(acb_t z, const acb_t x, const acb_t y, slong prec);
......
This diff is collapsed.
This diff is collapsed.
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_dot_precise(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep,
acb_srcptr y, slong ystep, slong len, slong prec)
{
arb_ptr tmp;
slong i;
tmp = flint_malloc(sizeof(arb_struct) * (4 * len));
for (i = 0; i < len; i++)
{
tmp[0 * len + i] = *acb_realref(x + i * xstep);
tmp[1 * len + i] = *acb_imagref(x + i * xstep);
tmp[2 * len + i] = *acb_realref(y + i * ystep);
arb_init(tmp + 3 * len + i);
arb_neg(tmp + 3 * len + i, acb_imagref(y + i * ystep));
}
arb_dot_precise(acb_realref(res), initial == NULL ? NULL : acb_realref(initial), subtract,
tmp, 1, tmp + 2 * len, 1, 2 * len, prec);
for (i = 0; i < len; i++)
arb_clear(tmp + 3 * len + i);
for (i = 0; i < len; i++)
{
tmp[0 * len + i] = *acb_realref(x + i * xstep);
tmp[1 * len + i] = *acb_imagref(x + i * xstep);
tmp[2 * len + i] = *acb_imagref(y + i * ystep);
tmp[3 * len + i] = *acb_realref(y + i * ystep);
}
arb_dot_precise(acb_imagref(res), initial == NULL ? NULL : acb_imagref(initial), subtract,
tmp, 1, tmp + 2 * len, 1, 2 * len, prec);
flint_free(tmp);
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_dot_simple(acb_t res, const acb_t initial, int subtract,
acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)
{
slong i;
if (len <= 0)
{
if (initial == NULL)
acb_zero(res);
else
acb_set_round(res, initial, prec);
return;
}
if (initial == NULL)
{
acb_mul(res, x, y, prec);
}
else
{
if (subtract)
acb_neg(res, initial);
else
acb_set(res, initial);
acb_addmul(res, x, y, prec);
}
for (i = 1; i < len; i++)
acb_addmul(res, x + i * xstep, y + i * ystep, prec);
if (subtract)
acb_neg(res, res);
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
ARB_DLL extern slong acb_dot_gauss_dot_cutoff;
int main()
{
slong iter;
flint_rand_t state;
flint_printf("approx_dot....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
{
acb_ptr x, y;
acb_t s1, s2, z;
slong i, len, prec, xbits, ybits, ebits;
int initial, subtract, revx, revy;
if (n_randint(state, 100) == 0)
len = n_randint(state, 100);
else if (n_randint(state, 10) == 0)
len = n_randint(state, 10);
else
len = n_randint(state, 3);
acb_dot_gauss_dot_cutoff = 3 + n_randint(state, 30);
if (n_randint(state, 10) != 0 || len > 10)
{
prec = 2 + n_randint(state, 500);
xbits = 2 + n_randint(state, 500);
ybits = 2 + n_randint(state, 500);
}
else
{
prec = 2 + n_randint(state, 4000);
xbits = 2 + n_randint(state, 4000);
ybits = 2 + n_randint(state, 4000);
}
if (n_randint(state, 100) == 0)
ebits = 2 + n_randint(state, 100);
else
ebits = 2 + n_randint(state, 10);
initial = n_randint(state, 2);
subtract = n_randint(state, 2);
revx = n_randint(state, 2);
revy = n_randint(state, 2);
x = _acb_vec_init(len);
y = _acb_vec_init(len);
acb_init(s1);
acb_init(s2);
acb_init(z);
switch (n_randint(state, 3))
{
case 0:
for (i = 0; i < len; i++)
{
acb_randtest(x + i, state, xbits, ebits);
acb_randtest(y + i, state, ybits, ebits);
}
break;
/* Test with cancellation */
case 1:
for (i = 0; i < len; i++)
{
if (i <= len / 2)
{
acb_randtest(x + i, state, xbits, ebits);
acb_randtest(y + i, state, ybits, ebits);
}
else
{
acb_neg(x + i, x + len - i - 1);
acb_set(y + i, y + len - i - 1);
}
}
break;
default:
for (i = 0; i < len; i++)
{
if (i <= len / 2)
{
acb_randtest(x + i, state, xbits, ebits);
acb_randtest(y + i, state, ybits, ebits);
}
else
{
acb_neg_round(x + i, x + len - i - 1, 2 + n_randint(state, 500));
acb_set_round(y + i, y + len - i - 1, 2 + n_randint(state, 500));
}
}
break;
}
acb_randtest(s1, state, 200, 100);
acb_randtest(s2, state, 200, 100);
acb_randtest(z, state, xbits, ebits);
acb_approx_dot(s1, initial ? z : NULL, subtract,
revx ? (x + len - 1) : x, revx ? -1 : 1,
revy ? (y + len - 1) : y, revy ? -1 : 1,
len, prec);
mag_zero(arb_radref(acb_realref(s1)));
mag_zero(arb_radref(acb_imagref(s1)));
/* With the fast algorithm, we expect identical results when
reversing the vectors. */
if (ebits <= 12)
{
acb_approx_dot(s2, initial ? z : NULL, subtract,
!revx ? (x + len - 1) : x, !revx ? -1 : 1,
!revy ? (y + len - 1) : y, !revy ? -1 : 1,
len, prec);
mag_zero(arb_radref(acb_realref(s2)));
mag_zero(arb_radref(acb_imagref(s2)));
if (!acb_equal(s1, s2))
{
flint_printf("FAIL (reversal)\n\n");
flint_printf("iter = %wd, len = %wd, prec = %wd, ebits = %wd\n\n", iter, len, prec, ebits);
if (initial)
{
flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
}
for (i = 0; i < len; i++)
{
flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
}
flint_printf("\n\n");
flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_abort();
}
}
/* Verify that radii are ignored */
for (i = 0; i < len; i++)
{
arb_get_mid_arb(acb_realref(x + i), acb_realref(x + i));
arb_get_mid_arb(acb_imagref(x + i), acb_imagref(x + i));
arb_get_mid_arb(acb_realref(y + i), acb_realref(y + i));
arb_get_mid_arb(acb_imagref(y + i), acb_imagref(y + i));
}
arb_get_mid_arb(acb_realref(z), acb_realref(z));
arb_get_mid_arb(acb_imagref(z), acb_imagref(z));
acb_approx_dot(s2, initial ? z : NULL, subtract,
revx ? (x + len - 1) : x, revx ? -1 : 1,
revy ? (y + len - 1) : y, revy ? -1 : 1,
len, prec);
mag_zero(arb_radref(acb_realref(s2)));
mag_zero(arb_radref(acb_imagref(s2)));
if (!acb_equal(s1, s2))
{
flint_printf("FAIL (radii)\n\n");
flint_printf("iter = %wd, len = %wd, prec = %wd, ebits = %wd\n\n", iter, len, prec, ebits);
if (initial)
{
flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
}
for (i = 0; i < len; i++)
{
flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
}
flint_printf("\n\n");
flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_abort();
}
/* Compare with acb_dot */
acb_approx_dot(s2, initial ? z : NULL, subtract,
revx ? (x + len - 1) : x, revx ? -1 : 1,
revy ? (y + len - 1) : y, revy ? -1 : 1,
len, prec);
{
mag_t err, xx, yy;
mag_init(err);
mag_init(xx);
mag_init(yy);
if (initial)
acb_get_mag(err, z);
for (i = 0; i < len; i++)
{
acb_get_mag(xx, revx ? x + len - 1 - i : x + i);
acb_get_mag(yy, revx ? y + len - 1 - i : y + i);
mag_addmul(err, xx, yy);
}
mag_mul_2exp_si(err, err, -prec + 2);
acb_add_error_mag(s2, err);
if (!acb_contains(s2, s1))
{
flint_printf("FAIL (inclusion)\n\n");
flint_printf("iter = %wd, len = %wd, prec = %wd, ebits = %wd\n\n", iter, len, prec, ebits);
if (initial)
{
flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
}
for (i = 0; i < len; i++)
{
flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
}
flint_printf("\n\n");
flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_abort();
}
mag_clear(err);
mag_clear(xx);
mag_clear(yy);
}
acb_clear(s1);
acb_clear(s2);
acb_clear(z);
_acb_vec_clear(x, len);
_acb_vec_clear(y, len);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
ARB_DLL extern slong acb_dot_gauss_dot_cutoff;
int main()
{
slong iter;
flint_rand_t state;
flint_printf("dot....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000000 * arb_test_multiplier(); iter++)
{
acb_ptr x, y;
acb_t s1, s2, z;
slong i, len, prec, xbits, ybits, ebits;
int ok, initial, subtract, revx, revy;
if (n_randint(state, 100) == 0)
len = n_randint(state, 50);
else if (n_randint(state, 10) == 0)
len = n_randint(state, 5);
else
len = n_randint(state, 3);
acb_dot_gauss_dot_cutoff = 3 + n_randint(state, 30);
if (n_randint(state, 10) != 0 || len > 10)
{
prec = 2 + n_randint(state, 500);
xbits = 2 + n_randint(state, 500);
ybits = 2 + n_randint(state, 500);
}
else
{
prec = 2 + n_randint(state, 5000);
xbits = 2 + n_randint(state, 5000);
ybits = 2 + n_randint(state, 5000);
}
if (n_randint(state, 100) == 0)
ebits = 2 + n_randint(state, 100);
else
ebits = 2 + n_randint(state, 10);
initial = n_randint(state, 2);
subtract = n_randint(state, 2);
revx = n_randint(state, 2);
revy = n_randint(state, 2);
x = _acb_vec_init(len);
y = _acb_vec_init(len);
acb_init(s1);
acb_init(s2);
acb_init(z);
switch (n_randint(state, 3))
{
case 0:
for (i = 0; i < len; i++)
{
acb_randtest(x + i, state, xbits, ebits);
acb_randtest(y + i, state, ybits, ebits);
}
break;
/* Test with cancellation */
case 1:
for (i = 0; i < len; i++)
{
if (i <= len / 2)
{
acb_randtest(x + i, state, xbits, ebits);
acb_randtest(y + i, state, ybits, ebits);
}
else
{
acb_neg(x + i, x + len - i - 1);
acb_set(y + i, y + len - i - 1);
}
}
break;
default:
for (i = 0; i < len; i++)
{
if (i <= len / 2)
{
acb_randtest(x + i, state, xbits, ebits);
acb_randtest(y + i, state, ybits, ebits);
}
else
{
acb_neg_round(x + i, x + len - i - 1, 2 + n_randint(state, 500));
acb_set_round(y + i, y + len - i - 1, 2 + n_randint(state, 500));
}
}
break;
}
acb_randtest(s1, state, 200, 100);
acb_randtest(s2, state, 200, 100);
acb_randtest(z, state, xbits, ebits);
acb_dot(s1, initial ? z : NULL, subtract,
revx ? (x + len - 1) : x, revx ? -1 : 1,
revy ? (y + len - 1) : y, revy ? -1 : 1,
len, prec);
acb_dot_precise(s2, initial ? z : NULL, subtract,
revx ? (x + len - 1) : x, revx ? -1 : 1,
revy ? (y + len - 1) : y, revy ? -1 : 1,
len, ebits <= 12 ? ARF_PREC_EXACT : 2 * prec + 100);
if (ebits <= 12)
ok = acb_contains(s1, s2);
else
ok = acb_overlaps(s1, s2);
if (!ok)
{
flint_printf("FAIL\n\n");
flint_printf("iter = %wd, len = %wd, prec = %wd, ebits = %wd, subtract = %d\n\n", iter, len, prec, ebits, subtract);
if (initial)
{
flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
}
for (i = 0; i < len; i++)
{
flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
}
flint_printf("\n\n");
flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_abort();
}
/* With the fast algorithm, we expect identical results when
reversing the vectors. */
if (ebits <= 12)
{
revx ^= 1;
revy ^= 1;
acb_dot(s2, initial ? z : NULL, subtract,
revx ? (x + len - 1) : x, revx ? -1 : 1,
revy ? (y + len - 1) : y, revy ? -1 : 1,
len, prec);
if (!acb_equal(s1, s2))
{
flint_printf("FAIL (reversal)\n\n");
flint_printf("iter = %wd, len = %wd, prec = %wd, ebits = %wd\n\n", iter, len, prec, ebits);
if (initial)
{
flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
}
for (i = 0; i < len; i++)
{
flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
}
flint_printf("\n\n");
flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
flint_abort();
}
}
acb_clear(s1);
acb_clear(s2);
acb_clear(z);
_acb_vec_clear(x, len);
_acb_vec_clear(y, len);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}
......@@ -262,9 +262,12 @@ acb_dft_bluestein_init(acb_dft_bluestein_t t, slong n, slong prec)
ACB_DFT_INLINE void
acb_dft_bluestein_clear(acb_dft_bluestein_t t)
{
_acb_vec_clear(t->z, t->n);
_acb_vec_clear(t->g, t->rad2->n);
acb_dft_rad2_clear(t->rad2);
if (t->n != 0)
{
_acb_vec_clear(t->z, t->n);
_acb_vec_clear(t->g, t->rad2->n);
acb_dft_rad2_clear(t->rad2);
}
}
void _acb_dft_crt_init(acb_dft_crt_t crt, slong dv, slong len, slong prec);
......
......@@ -88,16 +88,21 @@ _acb_dft_bluestein_init(acb_dft_bluestein_t t, slong dv, slong n, slong prec)
{
acb_ptr z, g;
slong k, n2;
int e = n_clog(2 * n - 1, 2);
int e;
t->n = n;
t->dv = dv;
if (n == 0)
return;
e = n_clog(2 * n - 1, 2);
if (DFT_VERB)
flint_printf("dft_bluestein: init z[2^%i]\n", e);
acb_dft_rad2_init(t->rad2, e, prec);
t->n = n;
t->dv = dv;
t->z = z = _acb_vec_init(n);
_acb_vec_bluestein_factors(t->z, n, prec);
......@@ -118,6 +123,9 @@ acb_dft_bluestein_precomp(acb_ptr w, acb_srcptr v, const acb_dft_bluestein_t t,
slong n = t->n, np = t->rad2->n, dv = t->dv;
acb_ptr fp;
if (n == 0)
return;
fp = _acb_vec_init(np);
_acb_vec_kronecker_mul_step(fp, t->z, v, dv, n, prec);
......
......@@ -316,6 +316,9 @@ stieltjes_mag(double n)
double pi = 3.141592653589793;
int i;
if (n <= 1.0)
return 0.0;
va = 1e-6;
vb = 0.5 * pi - 1e-6;
......@@ -754,7 +757,7 @@ acb_dirichlet_stieltjes(acb_t res, const fmpz_t n, const acb_t a, slong prec)
cutoff = FLINT_MAX(100, prec / 2);
cutoff = FLINT_MIN(cutoff, 10000);
if (acb_is_one(a) && fmpz_cmp_ui(n, cutoff) >= 0)
if (fmpz_cmp_ui(n, cutoff) >= 0)
{
acb_dirichlet_stieltjes_integral(res, n, a, prec);
}
......
......@@ -389,6 +389,7 @@ int acb_mat_solve(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec)
int acb_mat_solve_precond(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec);
void acb_mat_approx_mul(acb_mat_t C, const acb_mat_t A, const acb_mat_t B, slong prec);
void acb_mat_approx_solve_triu(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec);
void acb_mat_approx_solve_tril(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec);
int acb_mat_approx_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec);
......
......@@ -194,8 +194,7 @@ acb_mat_approx_lu_recursive(slong * P, acb_mat_t LU, const acb_mat_t A, slong pr
/* acb_mat_submul(A11, A11, A10, A01, prec); */
acb_mat_t T;
acb_mat_init(T, A10->r, A01->c);
acb_mat_mul(T, A10, A01, prec);
acb_mat_get_mid(T, T);
acb_mat_approx_mul(T, A10, A01, prec);
acb_mat_sub(A11, A11, T, prec);
acb_mat_get_mid(A11, A11);
acb_mat_clear(T);
......@@ -227,4 +226,3 @@ acb_mat_approx_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)
else
return acb_mat_approx_lu_recursive(P, LU, A, prec);
}