add arf_frexp

This commit is contained in:
Fredrik Johansson 2016-04-29 13:50:05 +02:00
parent 4c05069ba9
commit 13bbbd92c6
4 changed files with 118 additions and 0 deletions

2
arf.h
View file

@ -740,6 +740,8 @@ arf_abs_bound_le_2exp_fmpz(fmpz_t b, const arf_t x)
slong arf_abs_bound_lt_2exp_si(const arf_t x);
void arf_frexp(arf_t man, fmpz_t exp, const arf_t x);
void arf_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const arf_t x);
int _arf_get_integer_mpn(mp_ptr y, mp_srcptr x, mp_size_t xn, slong exp);

23
arf/frexp.c Normal file
View file

@ -0,0 +1,23 @@
/*
Copyright (C) 2016 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 "arf.h"
void
arf_frexp(arf_t man, fmpz_t exp, const arf_t x)
{
arf_set(man, x);
fmpz_zero(exp);
if (!arf_is_special(man))
fmpz_swap(exp, ARF_EXPREF(man));
}

87
arf/test/t-frexp.c Normal file
View file

@ -0,0 +1,87 @@
/*
Copyright (C) 2016 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 "arf.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("frexp....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
{
arf_t x, y, z;
fmpz_t e, f;
int ok;
arf_init(x);
arf_init(y);
arf_init(z);
fmpz_init(e);
fmpz_init(f);
arf_randtest_special(x, state, 200, 100);
arf_randtest_special(y, state, 200, 100);
fmpz_randtest(e, state, 100);
fmpz_randtest(f, state, 100);
arf_frexp(y, e, x);
arf_mul_2exp_fmpz(z, y, e);
ok = 1;
if (!arf_equal(z, x)) ok = 0;
if (arf_is_special(x) && !fmpz_is_zero(e)) ok = 0;
if (!arf_is_special(x) && !(arf_cmpabs_2exp_si(y, 0) < 0
&& arf_cmpabs_2exp_si(y, -1) >= 0)) ok = 0;
if (!ok)
{
printf("FAIL\n");
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); arf_print(y); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("e = "); fmpz_print(e); printf("\n\n");
printf("f = "); fmpz_print(f); printf("\n\n");
abort();
}
arf_frexp(x, f, x);
if (!arf_equal(x, y) || !fmpz_equal(e, f))
{
printf("FAIL (aliasing)\n");
printf("x = "); arf_print(x); printf("\n\n");
printf("y = "); arf_print(y); printf("\n\n");
printf("z = "); arf_print(z); printf("\n\n");
printf("e = "); fmpz_print(e); printf("\n\n");
printf("f = "); fmpz_print(f); printf("\n\n");
abort();
}
arf_clear(x);
arf_clear(y);
arf_clear(z);
fmpz_clear(e);
fmpz_clear(f);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -246,6 +246,12 @@ Assignment, rounding and conversions
If *x* is zero, both *m* and *e* are set to zero. If *x* is
infinite or NaN, the result is undefined.
.. function:: void arf_frexp(arf_t m, fmpz_t e, const arf_t x)
Writes *x* as `m \times 2^e`, where `1/2 \le |m| < 1` if *x* is a normal
value. If *x* is a special value, copies this to *m* and sets *e* to zero.
Note: for the inverse operation (*ldexp*), use :func:`arf_mul_2exp_fmpz`.
.. function:: double arf_get_d(const arf_t x, arf_rnd_t rnd)
Returns *x* rounded to a double in the direction specified by *rnd*.