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

rewrite arb_get_unique_fmpz and arb_get_interval_fmpz_2exp, reducing overhead, making them more robust with huge exponents, and documenting this more carefully

Browse source
This commit is contained in:
Fredrik Johansson 2016-02-25 16:55:21 +01:00
parent 8c5d26de65
commit 7cebb85b15
6 changed files with 288 additions and 65 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

80
arb/get_interval_fmpz_2exp.c
View file

@ -28,9 +28,79 @@
void
arb_get_interval_fmpz_2exp(fmpz_t a, fmpz_t b, fmpz_t exp, const arb_t x)
{
fmprb_t t;
fmprb_init(t);
arb_get_fmprb(t, x);
fmprb_get_interval_fmpz_2exp(a, b, exp, t);
fmprb_clear(t);
if (!arb_is_finite(x))
{
printf("arb_get_interval_fmpz_2exp: expected finite input\n");
abort();
}
else if (arb_is_exact(x))
{
arf_get_fmpz_2exp(a, exp, arb_midref(x));
fmpz_set(b, a);
}
else if (arf_is_zero(arb_midref(x)))
{
arf_t t;
arf_init_set_mag_shallow(t, arb_radref(x));
arf_get_fmpz_2exp(b, exp, t);
fmpz_neg(a, b);
}
else
{
arf_t rad;
fmpz_t tmp;
slong shift;
mp_bitcnt_t aval, bval;
fmpz_init(tmp);
arf_get_fmpz_2exp(a, exp, arb_midref(x));
arf_init_set_mag_shallow(rad, arb_radref(x));
arf_get_fmpz_2exp(b, tmp, rad);
shift = _fmpz_sub_small(exp, tmp);
if (FLINT_ABS(shift) >= WORD_MAX / 2)
{
printf("arb_get_interval_fmpz_2exp: too large shift\n");
abort();
}
if (shift >= 0)
{
fmpz_mul_2exp(a, a, shift);
fmpz_set(exp, tmp);
}
else
fmpz_mul_2exp(b, b, -shift);
fmpz_sub(tmp, a, b);
fmpz_add(b, a, b);
fmpz_swap(tmp, a);
if (fmpz_is_zero(a))
{
aval = fmpz_val2(b);
}
else if (fmpz_is_zero(b))
{
aval = fmpz_val2(a);
}
else
{
aval = fmpz_val2(a);
bval = fmpz_val2(b);
aval = FLINT_MIN(aval, bval);
}
if (aval > 0)
{
fmpz_add_ui(exp, exp, aval);
fmpz_tdiv_q_2exp(a, a, aval);
fmpz_tdiv_q_2exp(b, b, aval);
}
fmpz_clear(tmp);
}
}

91
arb/get_unique_fmpz.c
View file

@ -19,7 +19,7 @@
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2012, 2016 Fredrik Johansson
******************************************************************************/
@ -28,15 +28,90 @@
int
arb_get_unique_fmpz(fmpz_t z, const arb_t x)
{
fmprb_t t;
int ans;
if (!arb_is_finite(x))
{
return 0;
}
else if (arb_is_exact(x))
{
/* x = b*2^e, e >= 0 */
if (arf_is_int(arb_midref(x)))
{
/* arf_get_fmpz aborts on overflow */
arf_get_fmpz(z, arb_midref(x), ARF_RND_DOWN);
return 1;
}
else
{
return 0;
}
}
/* if the radius is >= 1, there are at least two integers */
else if (mag_cmp_2exp_si(arb_radref(x), 0) >= 0)
{
return 0;
}
/* there are 0 or 1 integers if the radius is < 1 */
else
{
fmpz_t a, b, exp;
int res;
fmprb_init(t);
arb_get_fmprb(t, x);
/* if the midpoint is exactly an integer, it is what we want */
if (arf_is_int(arb_midref(x)))
{
/* arf_get_fmpz aborts on overflow */
arf_get_fmpz(z, arb_midref(x), ARF_RND_DOWN);
return 1;
}
ans = fmprb_get_unique_fmpz(z, t);
fmpz_init(a);
fmpz_init(b);
fmpz_init(exp);
fmprb_clear(t);
return ans;
/* if the radius is tiny, it can't be an integer */
arf_bot(a, arb_midref(x));
if (fmpz_cmp(a, MAG_EXPREF(arb_radref(x))) > 0)
{
res = 0;
}
else
{
arb_get_interval_fmpz_2exp(a, b, exp, x);
if (COEFF_IS_MPZ(*exp))
{
flint_printf("arb_get_unique_fmpz: input too large\n");
abort();
}
if (*exp >= 0)
{
res = fmpz_equal(a, b);
if (res)
{
fmpz_mul_2exp(a, a, *exp);
fmpz_mul_2exp(b, b, *exp);
}
}
else
{
fmpz_cdiv_q_2exp(a, a, -(*exp));
fmpz_fdiv_q_2exp(b, b, -(*exp));
res = fmpz_equal(a, b);
}
if (res)
fmpz_set(z, a);
}
fmpz_clear(a);
fmpz_clear(b);
fmpz_clear(exp);
return res;
}
}

15
arb/test/t-get_interval_fmpz_2exp.c
View file

@ -53,7 +53,7 @@ int main()
arf_set_fmpz_2exp(y, a, exp);
if (!arb_contains_arf(x, y))
{
flint_printf("FAIL:\n\n");
flint_printf("FAIL (a):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("a = "); fmpz_print(a);
flint_printf(" exp = "); fmpz_print(exp); flint_printf("\n\n");
@ -63,13 +63,24 @@ int main()
arf_set_fmpz_2exp(y, b, exp);
if (!arb_contains_arf(x, y))
{
flint_printf("FAIL:\n\n");
flint_printf("FAIL (b):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("b = "); fmpz_print(b);
flint_printf(" exp = "); fmpz_print(exp); flint_printf("\n\n");
abort();
}
if (fmpz_is_even(a) && fmpz_is_even(b) &&
!(fmpz_is_zero(a) && fmpz_is_zero(b)))
{
flint_printf("FAIL (odd):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("a = "); fmpz_print(a); flint_printf("\n\n");
flint_printf("b = "); fmpz_print(b); flint_printf("\n\n");
flint_printf(" exp = "); fmpz_print(exp); flint_printf("\n\n");
abort();
}
arb_clear(x);
arf_clear(y);
fmpz_clear(a);

130
arb/test/t-get_unique_fmpz.c
View file

@ -19,7 +19,7 @@
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2012, 2016 Fredrik Johansson
******************************************************************************/
@ -49,63 +49,103 @@ int main()
arb_randtest(x, state, 2000, 10);
unique = arb_get_unique_fmpz(y, x);
arb_get_interval_fmpz_2exp(a, b, exp, x);
if (fmpz_sgn(exp) >= 0)
/* generate tiny and huge radii */
if (iter % 2 == 0)
{
fmpz_mul_2exp(a, a, fmpz_get_si(exp));
fmpz_mul_2exp(b, b, fmpz_get_si(exp));
mag_randtest_special(arb_radref(x), state, 100);
unique = arb_get_unique_fmpz(y, x);
arf_get_fmpz(a, arb_midref(x), ARF_RND_FLOOR);
fmpz_add_ui(b, a, 1);
if (unique)
{
if (arb_contains_fmpz(x, a) == arb_contains_fmpz(x, b))
{
flint_printf("FAIL (1):\n\n");
flint_printf("x = "); arb_printd(x, 100); flint_printf("\n\n");
flint_printf("unique = %d\n\n", unique);
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
flint_printf("a = "); fmpz_print(a); flint_printf("\n\n");
flint_printf("b = "); fmpz_print(b); flint_printf("\n\n");
abort();
}
}
else
{
if (arb_contains_fmpz(x, a) != arb_contains_fmpz(x, b))
{
flint_printf("FAIL (2):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("unique = %d\n\n", unique);
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
flint_printf("a = "); fmpz_print(a); flint_printf("\n\n");
flint_printf("b = "); fmpz_print(b); flint_printf("\n\n");
abort();
}
}
}
else
{
fmpz_cdiv_q_2exp(a, a, -fmpz_get_si(exp));
fmpz_fdiv_q_2exp(b, b, -fmpz_get_si(exp));
}
unique2 = fmpz_equal(a, b);
unique = arb_get_unique_fmpz(y, x);
if ((unique != unique2) || (unique && !fmpz_equal(y, a)))
{
flint_printf("FAIL:\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("unique = %d, unique2 = %d\n\n", unique, unique2);
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
flint_printf("a = "); fmpz_print(a); flint_printf("\n\n");
flint_printf("b = "); fmpz_print(b); flint_printf("\n\n");
flint_printf(" exp = "); fmpz_print(exp); flint_printf("\n\n");
abort();
}
if (unique)
{
arb_set_fmpz(z, y);
arb_set_round(z, z, 2 + n_randint(state, 1000));
if (!arb_overlaps(x, z))
arb_get_interval_fmpz_2exp(a, b, exp, x);
if (fmpz_sgn(exp) >= 0)
{
flint_printf("FAIL (overlap):\n\n");
fmpz_mul_2exp(a, a, fmpz_get_si(exp));
fmpz_mul_2exp(b, b, fmpz_get_si(exp));
}
else
{
fmpz_cdiv_q_2exp(a, a, -fmpz_get_si(exp));
fmpz_fdiv_q_2exp(b, b, -fmpz_get_si(exp));
}
unique2 = fmpz_equal(a, b);
if ((unique != unique2) || (unique && !fmpz_equal(y, a)))
{
flint_printf("FAIL:\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("unique = %d, unique2 = %d\n\n", unique, unique2);
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
flint_printf("z = "); arb_print(z); flint_printf("\n\n");
flint_printf("a = "); fmpz_print(a); flint_printf("\n\n");
flint_printf("b = "); fmpz_print(b); flint_printf("\n\n");
flint_printf(" exp = "); fmpz_print(exp); flint_printf("\n\n");
abort();
}
fmpz_add_ui(b, y, 1);
if (arb_contains_fmpz(x, b))
if (unique)
{
flint_printf("FAIL (contains a + 1):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
abort();
}
arb_set_fmpz(z, y);
arb_set_round(z, z, 2 + n_randint(state, 1000));
fmpz_sub_ui(b, y, 1);
if (arb_contains_fmpz(x, b))
{
flint_printf("FAIL (contains a - 1):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
abort();
if (!arb_overlaps(x, z))
{
flint_printf("FAIL (overlap):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
flint_printf("z = "); arb_print(z); flint_printf("\n\n");
abort();
}
fmpz_add_ui(b, y, 1);
if (arb_contains_fmpz(x, b))
{
flint_printf("FAIL (contains a + 1):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
abort();
}
fmpz_sub_ui(b, y, 1);
if (arb_contains_fmpz(x, b))
{
flint_printf("FAIL (contains a - 1):\n\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
flint_printf("y = "); fmpz_print(y); flint_printf("\n\n");
abort();
}
}
}

25
doc/source/arb.rst
View file

@ -368,10 +368,19 @@ Radius and interval operations
Computes the exact interval represented by *x*, in the form of an integer
interval multiplied by a power of two, i.e. `x = [a, b] \times 2^{\text{exp}}`.
The result is normalized by removing common trailing zeros
from *a* and *b*.
The outcome is undefined if the midpoint or radius of *x* is non-finite,
or if the difference in magnitude between the midpoint and radius
is so large that representing the endpoints exactly would cause overflows.
This method aborts if *x* is infinite or NaN, or if the difference between
the exponents of the midpoint and the radius is so large that allocating
memory for the result fails.
Warning: this method will allocate a huge amount of memory to store
the result if the exponent difference is huge. Memory allocation could
succeed even if the required space is far larger than the physical
memory available on the machine, resulting in swapping. It is recommended
to check that the midpoint and radius of *x* both are within a
reasonable range before calling this method.
.. function:: void arb_set_interval_arf(arb_t x, const arf_t a, const arf_t b, slong prec)
@ -420,6 +429,16 @@ Radius and interval operations
nonzero. Otherwise (if *x* represents no integers or more than one integer),
returns zero.
This method aborts if there is a unique integer but that integer
is so large that allocating memory for the result fails.
Warning: this method will allocate a huge amount of memory to store
the result if there is a unique integer and that integer is huge.
Memory allocation could succeed even if the required space is far
larger than the physical memory available on the machine, resulting
in swapping. It is recommended to check that the midpoint of *x* is
within a reasonable range before calling this method.
.. function:: void arb_floor(arb_t y, const arb_t x, slong prec)
.. function:: void arb_ceil(arb_t y, const arb_t x, slong prec)

12
doc/source/arf.rst
View file

@ -258,8 +258,16 @@ Assignment, rounding and conversions
Sets *z* to *x* rounded to the nearest integer in the direction
specified by *rnd*. If rnd is *ARF_RND_NEAR*, rounds to the nearest
even integer in case of a tie. Aborts if *x* is infinite, NaN or if the
exponent is unreasonably large.
even integer in case of a tie.
This method aborts if *x* is infinite or NaN, or if the exponent of *x*
is so large that allocating memory for the result fails.
Warning: this method will allocate a huge amount of memory to store
the result if the exponent of *x* is huge. Memory allocation could
succeed even if the required space is far larger than the physical
memory available on the machine, resulting in swapping. It is recommended
to check that *x* is within a reasonable range before calling this method.
.. function:: slong arf_get_si(const arf_t x, arf_rnd_t rnd)