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

add mag_log, mag_log_lower, mag_neg_log, mag_neg_log_lower

Browse source
This commit is contained in:
Fredrik Johansson 2018-01-04 22:47:13 +01:00
parent 948cd473e2
commit da07b2ea20
8 changed files with 730 additions and 29 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

18
doc/source/mag.rst
View file

@ -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)`.

5
mag.h
View file

@ -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

220
mag/log.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}
}
}

42
mag/log1p.c
View file

@ -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);
}
}
}

116
mag/test/t-log.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

116
mag/test/t-log_lower.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

121
mag/test/t-neg_log.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

121
mag/test/t-neg_log_lower.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}