diff --git a/doc/source/mag.rst b/doc/source/mag.rst
index ae50fd05..8656a062 100644
--- a/doc/source/mag.rst
+++ b/doc/source/mag.rst
@@ -354,6 +354,24 @@ Powers and logarithms
Sets *z* to an upper bound for `\log(n)`.
+.. function:: void mag_log(mag_t z, const mag_t x)
+
+ Sets *z* to an upper bound for `\log(\max(1,x))`.
+
+.. function:: void mag_log_lower(mag_t z, const mag_t x)
+
+ Sets *z* to a lower bound for `\log(\max(1,x))`.
+
+.. function:: void mag_neg_log(mag_t z, const mag_t x)
+
+ Sets *z* to an upper bound for `-\log(\min(1,x))`, i.e. an upper
+ bound for `|\log(x)|` for `x \le 1`.
+
+.. function:: void mag_neg_log_lower(mag_t z, const mag_t x)
+
+ Sets *z* to a lower bound for `-\log(\min(1,x))`, i.e. a lower
+ bound for `|\log(x)|` for `x \le 1`.
+
.. function:: void mag_exp(mag_t z, const mag_t x)
Sets *z* to an upper bound for `\exp(x)`.
diff --git a/mag.h b/mag.h
index 4c700ee3..1c582e2b 100644
--- a/mag.h
+++ b/mag.h
@@ -607,6 +607,11 @@ void mag_log1p(mag_t z, const mag_t x);
void mag_log_ui(mag_t t, ulong n);
+void mag_log(mag_t z, const mag_t x);
+void mag_log_lower(mag_t z, const mag_t x);
+void mag_neg_log(mag_t z, const mag_t x);
+void mag_neg_log_lower(mag_t z, const mag_t x);
+
void mag_exp_maglim(mag_t y, const mag_t x, slong maglim);
MAG_INLINE void
diff --git a/mag/log.c b/mag/log.c
new file mode 100644
index 00000000..83b83148
--- /dev/null
+++ b/mag/log.c
@@ -0,0 +1,220 @@
+/*
+ Copyright (C) 2017 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 .
+*/
+
+#include "mag.h"
+
+void
+mag_neg_log(mag_t z, const mag_t x)
+{
+ if (mag_is_special(x))
+ {
+ if (mag_is_zero(x))
+ mag_inf(z);
+ else
+ mag_zero(z);
+ }
+ else
+ {
+ double t;
+ fmpz exp = MAG_EXP(x);
+
+ if (!COEFF_IS_MPZ(exp))
+ {
+ if (exp >= 1)
+ {
+ mag_zero(z);
+ }
+ else if (exp > -(1000 - MAG_BITS))
+ {
+ t = ldexp(MAG_MAN(x), exp - MAG_BITS);
+ t = -mag_d_log_lower_bound(t);
+ mag_set_d(z, t);
+ }
+ else
+ {
+ /* -log(2^(exp-1) * (2*v)) = -exp*log(2) - log(2*v), 1 <= 2v < 2 */
+ t = MAG_MAN(x) * ldexp(1.0, 1 - MAG_BITS);
+ t = -mag_d_log_lower_bound(t);
+ t -= (exp - 1) * 0.69314718055994530942;
+ t *= (1 + 1e-13);
+ mag_set_d(z, t);
+ }
+ }
+ else if (fmpz_sgn(MAG_EXPREF(x)) > 0)
+ {
+ mag_zero(z);
+ }
+ else
+ {
+ mag_inv(z, x);
+ mag_log(z, z);
+ }
+ }
+}
+
+void
+mag_neg_log_lower(mag_t z, const mag_t x)
+{
+ if (mag_is_special(x))
+ {
+ if (mag_is_zero(x))
+ mag_inf(z);
+ else
+ mag_zero(z);
+ }
+ else
+ {
+ double t;
+ fmpz exp = MAG_EXP(x);
+
+ if (!COEFF_IS_MPZ(exp))
+ {
+ if (exp >= 1)
+ {
+ mag_zero(z);
+ }
+ else if (exp > -(1000 - MAG_BITS))
+ {
+ t = ldexp(MAG_MAN(x), exp - MAG_BITS);
+ t = -mag_d_log_upper_bound(t);
+ mag_set_d_lower(z, t);
+ }
+ else
+ {
+ /* -log(2^(exp-1) * (2*v)) = -exp*log(2) - log(2*v), 1 <= 2v < 2 */
+ t = MAG_MAN(x) * ldexp(1.0, 1 - MAG_BITS);
+ t = -mag_d_log_upper_bound(t);
+ t -= (exp - 1) * 0.69314718055994530942;
+ t *= (1 - 1e-13);
+ mag_set_d_lower(z, t);
+ }
+ }
+ else if (fmpz_sgn(MAG_EXPREF(x)) > 0)
+ {
+ mag_zero(z);
+ }
+ else
+ {
+ mag_inv_lower(z, x);
+ mag_log_lower(z, z);
+ }
+ }
+}
+
+void
+mag_log(mag_t z, const mag_t x)
+{
+ if (mag_is_special(x))
+ {
+ if (mag_is_zero(x))
+ mag_zero(z);
+ else
+ mag_inf(z);
+ }
+ else
+ {
+ double t;
+ fmpz exp = MAG_EXP(x);
+
+ if (!COEFF_IS_MPZ(exp))
+ {
+ if (exp <= 0 || (exp == 1 && MAG_MAN(x) == MAG_ONE_HALF))
+ {
+ mag_zero(z);
+ }
+ else if (exp < 1000) /* safely within the double range */
+ {
+ t = ldexp(MAG_MAN(x), exp - MAG_BITS);
+ t = mag_d_log_upper_bound(t);
+ mag_set_d(z, t);
+ }
+ else
+ {
+ /* log(2^(exp-1)*(2v)) = (exp-1)*log(2) + log(2v), 1 <= 2v < 2 */
+ t = MAG_MAN(x) * ldexp(1.0, 1 - MAG_BITS); /* exact */
+ t = mag_d_log_upper_bound(t);
+ t += (exp - 1.0) * 0.69314718055994530942;
+ t *= (1 + 1e-13); /* last step could have rounded down */
+ mag_set_d(z, t);
+ }
+ }
+ else if (fmpz_sgn(MAG_EXPREF(x)) < 0)
+ {
+ mag_zero(z);
+ }
+ else
+ {
+ /* log(2^exp) = exp*log(2), log(2) < 744261118/2^30 */
+ fmpz_t t;
+ fmpz_init(t);
+ fmpz_mul_ui(t, MAG_EXPREF(x), 744261118);
+ mag_set_fmpz(z, t);
+ mag_mul_2exp_si(z, z, -30);
+ fmpz_clear(t);
+ }
+ }
+}
+
+void
+mag_log_lower(mag_t z, const mag_t x)
+{
+ if (mag_is_special(x))
+ {
+ if (mag_is_zero(x))
+ mag_zero(z);
+ else
+ mag_inf(z);
+ }
+ else
+ {
+ double t;
+ fmpz exp = MAG_EXP(x);
+
+ if (!COEFF_IS_MPZ(exp))
+ {
+ if (exp <= 0 || (exp == 1 && MAG_MAN(x) == MAG_ONE_HALF))
+ {
+ mag_zero(z);
+ }
+ else if (exp < 1000)
+ {
+ t = ldexp(MAG_MAN(x), exp - MAG_BITS);
+ t = mag_d_log_lower_bound(t);
+ mag_set_d_lower(z, t);
+ }
+ else
+ {
+ /* log(2^(exp-1) * (2*v)) = exp*log(2) + log(2*v), 1 <= 2v < 2 */
+ t = MAG_MAN(x) * ldexp(1.0, 1 - MAG_BITS);
+ t = mag_d_log_lower_bound(t);
+ t += (exp - 1) * 0.69314718055994530942;
+ t *= (1 - 1e-13);
+ mag_set_d_lower(z, t);
+ }
+ }
+ else if (fmpz_sgn(MAG_EXPREF(x)) < 0)
+ {
+ mag_zero(z);
+ }
+ else
+ {
+ /* log(2^exp) = exp*log(2), log(2) > 744261117/2^30 */
+ fmpz_t t;
+ fmpz_init(t);
+ fmpz_sub_ui(t, MAG_EXPREF(x), 1); /* lower bound for x */
+ fmpz_mul_ui(t, t, 744261117);
+ mag_set_fmpz_lower(z, t);
+ mag_mul_2exp_si(z, z, -30);
+ fmpz_clear(t);
+ }
+ }
+}
+
diff --git a/mag/log1p.c b/mag/log1p.c
index 4298fd9a..bbc9565b 100644
--- a/mag/log1p.c
+++ b/mag/log1p.c
@@ -24,6 +24,7 @@ mag_log1p(mag_t z, const mag_t x)
}
else
{
+ double t;
fmpz exp = MAG_EXP(x);
if (!COEFF_IS_MPZ(exp))
@@ -32,50 +33,33 @@ mag_log1p(mag_t z, const mag_t x)
if (exp < -10)
{
mag_set(z, x);
- return;
}
else if (exp < 1000)
{
- double t;
t = ldexp(MAG_MAN(x), exp - MAG_BITS);
t = (1.0 + t) * (1 + 1e-14);
t = mag_d_log_upper_bound(t);
mag_set_d(z, t);
- return;
+ }
+ else
+ {
+ /* log(2^(exp-1) * (2*v)) = exp*log(2) + log(2*v), 1 <= 2v < 2 */
+ t = (MAG_MAN(x) + 1) * ldexp(1.0, 1 - MAG_BITS);
+ t = mag_d_log_upper_bound(t);
+ t += (exp - 1) * 0.69314718055994530942;
+ t *= (1 + 1e-13);
+ mag_set_d(z, t);
}
}
else if (fmpz_sgn(MAG_EXPREF(x)) < 0)
{
/* Quick bound by x */
mag_set(z, x);
- return;
}
-
- /* Now we must have x >= 2^1000 */
- /* Use log(2^(exp-1) * (2*v)) = exp*log(2) + log(2*v) */
+ else
{
- double t;
- fmpz_t b;
- mag_t u;
-
- mag_init(u);
- fmpz_init(b);
-
- /* incrementing the mantissa gives an upper bound for x+1 */
- t = ldexp(MAG_MAN(x) + 1, 1 - MAG_BITS);
- t = mag_d_log_upper_bound(t);
- mag_set_d(u, t);
-
- /* log(2) < 744261118/2^30 */
- _fmpz_add_fast(b, MAG_EXPREF(x), -1);
- fmpz_mul_ui(b, b, 744261118);
- mag_set_fmpz(z, b);
- _fmpz_add_fast(MAG_EXPREF(z), MAG_EXPREF(z), -30);
-
- mag_add(z, z, u);
-
- mag_clear(u);
- fmpz_clear(b);
+ mag_add_ui(z, x, 1);
+ mag_log(z, z);
}
}
}
diff --git a/mag/test/t-log.c b/mag/test/t-log.c
new file mode 100644
index 00000000..7f7327b2
--- /dev/null
+++ b/mag/test/t-log.c
@@ -0,0 +1,116 @@
+/*
+ Copyright (C) 2014 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 .
+*/
+
+#include "mag.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("log....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
+ {
+ fmpr_t x, y, z, z2;
+ fmpz_t n;
+ mag_t xb, yb;
+
+ fmpr_init(x);
+ fmpr_init(y);
+ fmpr_init(z);
+ fmpr_init(z2);
+
+ mag_init(xb);
+ mag_init(yb);
+
+ fmpz_init(n);
+
+ mag_randtest_special(xb, state, 25);
+ mag_randtest_special(yb, state, 25);
+
+ if (n_randint(state, 2))
+ {
+ mag_log(yb, xb);
+
+ mag_get_fmpr(x, xb);
+ if (mag_cmp_2exp_si(xb, 0) <= 0)
+ fmpr_zero(z);
+ else
+ fmpr_log(z, x, MAG_BITS, FMPR_RND_UP);
+ }
+ else
+ {
+ mag_get_fmpr(x, xb);
+
+ fmpz_randtest(n, state, 100);
+ mag_mul_2exp_fmpz(xb, xb, n);
+ mag_log(yb, xb);
+
+ if (mag_cmp_2exp_si(xb, 0) <= 0)
+ {
+ fmpr_zero(z);
+ }
+ else
+ {
+ fmpr_log(z, x, 2 * MAG_BITS, FMPR_RND_UP);
+ fmpr_set_ui(z2, 2);
+ fmpr_log(z2, z2, 2 * MAG_BITS, FMPR_RND_UP);
+ fmpr_addmul_fmpz(z, z2, n, 2 * MAG_BITS, FMPR_RND_UP);
+ }
+ }
+
+ mag_get_fmpr(y, yb);
+
+ fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
+ fmpr_mul_2exp_si(z2, z2, -10);
+
+ MAG_CHECK_BITS(xb)
+ MAG_CHECK_BITS(yb)
+
+ if (!(fmpr_cmpabs(z, y) <= 0 && fmpr_cmpabs(y, z2) <= 0))
+ {
+ flint_printf("FAIL\n\n");
+ flint_printf("x = "); fmpr_print(x); flint_printf("\n\n");
+ flint_printf("y = "); fmpr_print(y); flint_printf("\n\n");
+ flint_printf("z = "); fmpr_print(z); flint_printf("\n\n");
+ flint_printf("z2 = "); fmpr_print(z2); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ mag_log(xb, xb);
+
+ if (!mag_equal(xb, yb))
+ {
+ flint_printf("FAIL (aliasing)\n\n");
+ flint_abort();
+ }
+
+ fmpr_clear(x);
+ fmpr_clear(y);
+ fmpr_clear(z);
+ fmpr_clear(z2);
+
+ mag_clear(xb);
+ mag_clear(yb);
+
+ fmpz_clear(n);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
+
diff --git a/mag/test/t-log_lower.c b/mag/test/t-log_lower.c
new file mode 100644
index 00000000..cae60d2a
--- /dev/null
+++ b/mag/test/t-log_lower.c
@@ -0,0 +1,116 @@
+/*
+ Copyright (C) 2014 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 .
+*/
+
+#include "mag.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("log_lower....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
+ {
+ fmpr_t x, y, z, z2;
+ fmpz_t n;
+ mag_t xb, yb;
+
+ fmpr_init(x);
+ fmpr_init(y);
+ fmpr_init(z);
+ fmpr_init(z2);
+
+ mag_init(xb);
+ mag_init(yb);
+
+ fmpz_init(n);
+
+ mag_randtest_special(xb, state, 25);
+ mag_randtest_special(yb, state, 25);
+
+ if (n_randint(state, 2))
+ {
+ mag_log_lower(yb, xb);
+
+ mag_get_fmpr(x, xb);
+ if (mag_cmp_2exp_si(xb, 0) <= 0)
+ fmpr_zero(z);
+ else
+ fmpr_log(z, x, MAG_BITS, FMPR_RND_DOWN);
+ }
+ else
+ {
+ mag_get_fmpr(x, xb);
+
+ fmpz_randtest(n, state, 100);
+ mag_mul_2exp_fmpz(xb, xb, n);
+ mag_log_lower(yb, xb);
+
+ if (mag_cmp_2exp_si(xb, 0) <= 0)
+ {
+ fmpr_zero(z);
+ }
+ else
+ {
+ fmpr_log(z, x, 2 * MAG_BITS, FMPR_RND_DOWN);
+ fmpr_set_ui(z2, 2);
+ fmpr_log(z2, z2, 2 * MAG_BITS, FMPR_RND_DOWN);
+ fmpr_addmul_fmpz(z, z2, n, 2 * MAG_BITS, FMPR_RND_DOWN);
+ }
+ }
+
+ mag_get_fmpr(y, yb);
+
+ fmpr_mul_ui(z2, z, 1023, MAG_BITS, FMPR_RND_DOWN);
+ fmpr_mul_2exp_si(z2, z2, -10);
+
+ MAG_CHECK_BITS(xb)
+ MAG_CHECK_BITS(yb)
+
+ if (!(fmpr_cmpabs(z2, y) <= 0 && fmpr_cmpabs(y, z) <= 0))
+ {
+ flint_printf("FAIL\n\n");
+ flint_printf("x = "); fmpr_print(x); flint_printf("\n\n");
+ flint_printf("y = "); fmpr_print(y); flint_printf("\n\n");
+ flint_printf("z = "); fmpr_print(z); flint_printf("\n\n");
+ flint_printf("z2 = "); fmpr_print(z2); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ mag_log_lower(xb, xb);
+
+ if (!mag_equal(xb, yb))
+ {
+ flint_printf("FAIL (aliasing)\n\n");
+ flint_abort();
+ }
+
+ fmpr_clear(x);
+ fmpr_clear(y);
+ fmpr_clear(z);
+ fmpr_clear(z2);
+
+ mag_clear(xb);
+ mag_clear(yb);
+
+ fmpz_clear(n);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
+
diff --git a/mag/test/t-neg_log.c b/mag/test/t-neg_log.c
new file mode 100644
index 00000000..20c7cbe8
--- /dev/null
+++ b/mag/test/t-neg_log.c
@@ -0,0 +1,121 @@
+/*
+ Copyright (C) 2014 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 .
+*/
+
+#include "mag.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("neg_log....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
+ {
+ fmpr_t x, y, z, z2;
+ fmpz_t n;
+ mag_t xb, yb;
+
+ fmpr_init(x);
+ fmpr_init(y);
+ fmpr_init(z);
+ fmpr_init(z2);
+
+ mag_init(xb);
+ mag_init(yb);
+
+ fmpz_init(n);
+
+ mag_randtest_special(xb, state, 25);
+ mag_randtest_special(yb, state, 25);
+
+ if (n_randint(state, 2))
+ {
+ mag_neg_log(yb, xb);
+
+ mag_get_fmpr(x, xb);
+ if (mag_cmp_2exp_si(xb, 0) >= 0)
+ fmpr_zero(z);
+ else
+ {
+ fmpr_log(z, x, MAG_BITS, FMPR_RND_UP);
+ fmpr_neg(z, z);
+ }
+ }
+ else
+ {
+ mag_get_fmpr(x, xb);
+
+ fmpz_randtest(n, state, 100);
+ mag_mul_2exp_fmpz(xb, xb, n);
+ mag_neg_log(yb, xb);
+
+ if (mag_cmp_2exp_si(xb, 0) >= 0)
+ {
+ fmpr_zero(z);
+ }
+ else
+ {
+ fmpr_log(z, x, 2 * MAG_BITS, FMPR_RND_UP);
+ fmpr_neg(z, z);
+ fmpr_set_ui(z2, 2);
+ fmpr_log(z2, z2, 2 * MAG_BITS, FMPR_RND_UP);
+ fmpr_submul_fmpz(z, z2, n, 2 * MAG_BITS, FMPR_RND_UP);
+ }
+ }
+
+ mag_get_fmpr(y, yb);
+
+ fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
+ fmpr_mul_2exp_si(z2, z2, -10);
+
+ MAG_CHECK_BITS(xb)
+ MAG_CHECK_BITS(yb)
+
+ if (!(fmpr_cmpabs(z, y) <= 0 && fmpr_cmpabs(y, z2) <= 0))
+ {
+ flint_printf("FAIL\n\n");
+ flint_printf("n = "); fmpz_print(n); flint_printf("\n\n");
+ flint_printf("x = "); fmpr_print(x); flint_printf("\n\n");
+ flint_printf("y = "); fmpr_print(y); flint_printf("\n\n");
+ flint_printf("z = "); fmpr_print(z); flint_printf("\n\n");
+ flint_printf("z2 = "); fmpr_print(z2); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ mag_neg_log(xb, xb);
+
+ if (!mag_equal(xb, yb))
+ {
+ flint_printf("FAIL (aliasing)\n\n");
+ flint_abort();
+ }
+
+ fmpr_clear(x);
+ fmpr_clear(y);
+ fmpr_clear(z);
+ fmpr_clear(z2);
+
+ mag_clear(xb);
+ mag_clear(yb);
+
+ fmpz_clear(n);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
+
diff --git a/mag/test/t-neg_log_lower.c b/mag/test/t-neg_log_lower.c
new file mode 100644
index 00000000..1a2bff12
--- /dev/null
+++ b/mag/test/t-neg_log_lower.c
@@ -0,0 +1,121 @@
+/*
+ Copyright (C) 2014 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 .
+*/
+
+#include "mag.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("neg_log_lower....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
+ {
+ fmpr_t x, y, z, z2;
+ fmpz_t n;
+ mag_t xb, yb;
+
+ fmpr_init(x);
+ fmpr_init(y);
+ fmpr_init(z);
+ fmpr_init(z2);
+
+ mag_init(xb);
+ mag_init(yb);
+
+ fmpz_init(n);
+
+ mag_randtest_special(xb, state, 25);
+ mag_randtest_special(yb, state, 25);
+
+ if (n_randint(state, 2))
+ {
+ mag_neg_log_lower(yb, xb);
+
+ mag_get_fmpr(x, xb);
+ if (mag_cmp_2exp_si(xb, 0) >= 0)
+ fmpr_zero(z);
+ else
+ {
+ fmpr_log(z, x, MAG_BITS, FMPR_RND_DOWN);
+ fmpr_neg(z, z);
+ }
+ }
+ else
+ {
+ mag_get_fmpr(x, xb);
+
+ fmpz_randtest(n, state, 100);
+ mag_mul_2exp_fmpz(xb, xb, n);
+ mag_neg_log_lower(yb, xb);
+
+ if (mag_cmp_2exp_si(xb, 0) >= 0)
+ {
+ fmpr_zero(z);
+ }
+ else
+ {
+ fmpr_log(z, x, 2 * MAG_BITS, FMPR_RND_DOWN);
+ fmpr_neg(z, z);
+ fmpr_set_ui(z2, 2);
+ fmpr_log(z2, z2, 2 * MAG_BITS, FMPR_RND_DOWN);
+ fmpr_submul_fmpz(z, z2, n, 2 * MAG_BITS, FMPR_RND_DOWN);
+ }
+ }
+
+ mag_get_fmpr(y, yb);
+
+ fmpr_mul_ui(z2, z, 1023, MAG_BITS, FMPR_RND_DOWN);
+ fmpr_mul_2exp_si(z2, z2, -10);
+
+ MAG_CHECK_BITS(xb)
+ MAG_CHECK_BITS(yb)
+
+ if (!(fmpr_cmpabs(z2, y) <= 0 && fmpr_cmpabs(y, z) <= 0))
+ {
+ flint_printf("FAIL\n\n");
+ flint_printf("n = "); fmpz_print(n); flint_printf("\n\n");
+ flint_printf("x = "); fmpr_print(x); flint_printf("\n\n");
+ flint_printf("y = "); fmpr_print(y); flint_printf("\n\n");
+ flint_printf("z = "); fmpr_print(z); flint_printf("\n\n");
+ flint_printf("z2 = "); fmpr_print(z2); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ mag_neg_log_lower(xb, xb);
+
+ if (!mag_equal(xb, yb))
+ {
+ flint_printf("FAIL (aliasing)\n\n");
+ flint_abort();
+ }
+
+ fmpr_clear(x);
+ fmpr_clear(y);
+ fmpr_clear(z);
+ fmpr_clear(z2);
+
+ mag_clear(xb);
+ mag_clear(yb);
+
+ fmpz_clear(n);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
+