From 13bbbd92c68637fb2e47471f353ffcba799af8ca Mon Sep 17 00:00:00 2001
From: Fredrik Johansson <fredrik.johansson@gmail.com>
Date: Fri, 29 Apr 2016 13:50:05 +0200
Subject: [PATCH] add arf_frexp

---
 arf.h              |  2 ++
 arf/frexp.c        | 23 ++++++++++++
 arf/test/t-frexp.c | 87 ++++++++++++++++++++++++++++++++++++++++++++++
 doc/source/arf.rst |  6 ++++
 4 files changed, 118 insertions(+)
 create mode 100644 arf/frexp.c
 create mode 100644 arf/test/t-frexp.c

diff --git a/arf.h b/arf.h
index e16d8343..ecd4c4be 100644
--- a/arf.h
+++ b/arf.h
@@ -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);
diff --git a/arf/frexp.c b/arf/frexp.c
new file mode 100644
index 00000000..32af9e72
--- /dev/null
+++ b/arf/frexp.c
@@ -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));
+}
+
diff --git a/arf/test/t-frexp.c b/arf/test/t-frexp.c
new file mode 100644
index 00000000..d079c164
--- /dev/null
+++ b/arf/test/t-frexp.c
@@ -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;
+}
+
diff --git a/doc/source/arf.rst b/doc/source/arf.rst
index fd70d0ea..5e727b47 100644
--- a/doc/source/arf.rst
+++ b/doc/source/arf.rst
@@ -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*.