more functions and test functions

arf.h

@@ -117,6 +117,9 @@ arf_rnd_to_mpfr(arf_rnd_t rnd)
 /* Encoding for special values. */
 #define ARF_IS_SPECIAL(x) (ARF_XSIZE(x) == 0)
 
+/* Value is +/- a power of two */
+#define ARF_IS_POW2(x) (ARF_SIZE(x) == 1) && (ARF_NOPTR_D(x)[0] == LIMB_TOP)
+
 /* To set a special value, first call this and then set the exponent. */
 #define ARF_MAKE_SPECIAL(x)         \
     do {                            \
@@ -494,6 +497,8 @@ int arf_neg_round(arf_t y, const arf_t x, long prec, arf_rnd_t rnd);
 
 void arf_get_fmpr(fmpr_t y, const arf_t x);
 
+void arf_set_fmpr(arf_t y, const fmpr_t x);
+
 int arf_get_mpfr(mpfr_t x, const arf_t y, mpfr_rnd_t rnd);
 
 void arf_set_mpfr(arf_t x, const mpfr_t y);
@@ -639,6 +644,38 @@ arf_set_round_fmpz_2exp(arf_t y, const fmpz_t x, const fmpz_t exp, long prec, ar
     }
 }
 
+static __inline__ void
+arf_abs_bound_lt_2exp_fmpz(fmpz_t b, const arf_t x)
+{
+    if (arf_is_special(x))
+        fmpz_zero(b);
+    else
+        fmpz_set(b, ARF_EXPREF(x));
+}
+
+static __inline__ void
+arf_abs_bound_le_2exp_fmpz(fmpz_t b, const arf_t x)
+{
+    if (arf_is_special(x))
+        fmpz_zero(b);
+    else if (ARF_IS_POW2(x))
+        fmpz_sub_ui(b, ARF_EXPREF(x), 1);
+    else
+        fmpz_set(b, ARF_EXPREF(x));
+}
+
+long arf_abs_bound_lt_2exp_si(const arf_t x);
+
+void arf_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const arf_t x);
+
+static __inline__ void
+arf_set_fmpz_2exp(arf_t x, const fmpz_t man, const fmpz_t exp)
+{
+    arf_set_fmpz(x, man);
+    if (!arf_is_zero(x))
+        fmpz_add_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp);
+}
+
 void arf_debug(const arf_t x);
 
 #define arf_print arf_debug

arf/abs_bound_lt_2exp_si.c

@@ -0,0 +1,56 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+long
+arf_abs_bound_lt_2exp_si(const arf_t x)
+{
+    long res;
+
+    if (arf_is_special(x))
+    {
+        if (arf_is_zero(x))
+            return -ARF_PREC_EXACT;
+        else
+            return ARF_PREC_EXACT;
+    }
+
+    if (!COEFF_IS_MPZ(ARF_EXP(x)))
+        return ARF_EXP(x);
+
+    if (fmpz_fits_si(ARF_EXPREF(x)))
+        res = fmpz_get_si(ARF_EXPREF(x));
+    else
+        res = fmpz_sgn(ARF_EXPREF(x)) < 0 ? -ARF_PREC_EXACT : ARF_PREC_EXACT;
+
+    if (res < -ARF_PREC_EXACT)
+        res = -ARF_PREC_EXACT;
+    if (res > ARF_PREC_EXACT)
+        res = ARF_PREC_EXACT;
+
+    return res;
+}
+

arf/cmp_2exp_si.c

@@ -28,8 +28,6 @@
 int
 arf_cmp_2exp_si(const arf_t x, long e)
 {
-    int pow2;
-
     if (arf_is_special(x))
     {
         if (arf_is_zero(x)) return -1;
@@ -41,18 +39,16 @@ arf_cmp_2exp_si(const arf_t x, long e)
     if (ARF_SGNBIT(x))
         return -1;
 
-    pow2 = (ARF_SIZE(x) == 1) && (ARF_NOPTR_D(x)[0] == LIMB_TOP);
-
     /* Fast path. */
     if (!COEFF_IS_MPZ(ARF_EXP(x)))
     {
-        if (pow2 && (ARF_EXP(x) - 1 == e))
+        if (ARF_IS_POW2(x) && (ARF_EXP(x) - 1 == e))
             return 0;
         else
             return (ARF_EXP(x) <= e) ? -1 : 1;
     }
 
-    if (pow2)
+    if (ARF_IS_POW2(x))
     {
         fmpz_t t;
         fmpz_init(t);

arf/cmpabs_2exp_si.c

@@ -28,8 +28,6 @@
 int
 arf_cmpabs_2exp_si(const arf_t x, long e)
 {
-    int pow2;
-
     if (arf_is_special(x))
     {
         if (arf_is_zero(x)) return -1;
@@ -37,18 +35,16 @@ arf_cmpabs_2exp_si(const arf_t x, long e)
         return 0;
     }
 
-    pow2 = (ARF_SIZE(x) == 1) && (ARF_NOPTR_D(x)[0] == LIMB_TOP);
-
     /* Fast path. */
     if (!COEFF_IS_MPZ(ARF_EXP(x)))
     {
-        if (pow2 && (ARF_EXP(x) - 1 == e))
+        if (ARF_IS_POW2(x) && (ARF_EXP(x) - 1 == e))
             return 0;
         else
             return (ARF_EXP(x) <= e) ? -1 : 1;
     }
 
-    if (pow2)
+    if (ARF_IS_POW2(x))
     {
         fmpz_t t;
         fmpz_init(t);

arf/get_fmpr.c

@@ -41,16 +41,7 @@ arf_get_fmpr(fmpr_t y, const arf_t x)
     }
     else
     {
-        mp_srcptr xptr;
-        mp_size_t xn;
-        long shift;
-
-        ARF_GET_MPN_READONLY(xptr, xn, x);
-
-        _fmpr_set_round_mpn(&shift, fmpr_manref(y),
-            xptr, xn, ARF_SGNBIT(x), FMPR_PREC_EXACT, FMPR_RND_DOWN);
-
-        fmpz_add_si(fmpr_expref(y), ARF_EXPREF(x), shift - xn * FLINT_BITS);
+        arf_get_fmpz_2exp(fmpr_manref(y), fmpr_expref(y), x);
     }
 }
 

arf/get_fmpz_2exp.c

@@ -0,0 +1,72 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+void
+arf_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const arf_t x)
+{
+    if (arf_is_special(x))
+    {
+        fmpz_zero(man);
+        fmpz_zero(exp);
+    }
+    else
+    {
+        mp_srcptr xptr;
+        mp_size_t xn;
+        int shift;
+
+        ARF_GET_MPN_READONLY(xptr, xn, x);
+
+        count_trailing_zeros(shift, xptr[0]);
+
+        fmpz_sub_ui(exp, ARF_EXPREF(x), xn * FLINT_BITS - shift);
+
+        if (xn == 1)
+        {
+            if (ARF_SGNBIT(x))
+                fmpz_neg_ui(man, xptr[0] >> shift);
+            else
+                fmpz_set_ui(man, xptr[0] >> shift);
+        }
+        else
+        {
+            __mpz_struct * z = _fmpz_promote(man);
+
+            if (z->_mp_alloc < xn)
+                mpz_realloc(z, xn);
+
+            if (shift == 0)
+                flint_mpn_copyi(z->_mp_d, xptr, xn);
+            else
+                mpn_rshift(z->_mp_d, xptr, xn, shift);
+
+            /* top limb cannot be zero */
+            z->_mp_size = ARF_SGNBIT(x) ? -xn : xn;
+        }
+    }
+}
+

arf/test/t-abs_bound_le_2exp_fmpz.c

@@ -0,0 +1,82 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("abs_bound_le_2exp_fmpz....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arf_t x, y;
+        fmpz_t b;
+        int cmp1, cmp2;
+
+        arf_init(x);
+        arf_init(y);
+        fmpz_init(b);
+
+        arf_randtest_not_zero(x, state, 2 + n_randint(state, 1000), 100);
+        arf_abs_bound_le_2exp_fmpz(b, x);
+
+        arf_one(y);
+        arf_mul_2exp_fmpz(y, y, b);
+
+        cmp1 = (arf_cmpabs(x, y) <= 0);
+
+        arf_mul_2exp_si(y, y, -1);
+
+        cmp2 = (arf_cmpabs(y, x) < 0);
+
+        arf_mul_2exp_si(y, y, 1);
+
+        if (!cmp1 || !cmp2)
+        {
+            printf("FAIL\n\n");
+            printf("x = "); arf_print(x); printf("\n\n");
+            printf("y = "); arf_print(y); printf("\n\n");
+            printf("b = "); fmpz_print(b); printf("\n\n");
+            printf("cmp1 = %d, cmp2 = %d\n\n", cmp1, cmp2);
+            abort();
+        }
+
+        arf_clear(x);
+        arf_clear(y);
+        fmpz_clear(b);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arf/test/t-abs_bound_lt_2exp_fmpz.c

@@ -0,0 +1,82 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("abs_bound_lt_2exp_fmpz....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arf_t x, y;
+        fmpz_t b;
+        int cmp1, cmp2;
+
+        arf_init(x);
+        arf_init(y);
+        fmpz_init(b);
+
+        arf_randtest_not_zero(x, state, 2 + n_randint(state, 1000), 100);
+        arf_abs_bound_lt_2exp_fmpz(b, x);
+
+        arf_one(y);
+        arf_mul_2exp_fmpz(y, y, b);
+
+        cmp1 = (arf_cmpabs(x, y) < 0);
+
+        arf_mul_2exp_si(y, y, -1);
+
+        cmp2 = (arf_cmpabs(y, x) <= 0);
+
+        arf_mul_2exp_si(y, y, 1);
+
+        if (!cmp1 || !cmp2)
+        {
+            printf("FAIL\n\n");
+            printf("x = "); arf_print(x); printf("\n\n");
+            printf("y = "); arf_print(y); printf("\n\n");
+            printf("b = "); fmpz_print(b); printf("\n\n");
+            printf("cmp1 = %d, cmp2 = %d\n\n", cmp1, cmp2);
+            abort();
+        }
+
+        arf_clear(x);
+        arf_clear(y);
+        fmpz_clear(b);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arf/test/t-abs_bound_lt_2exp_si.c

@@ -0,0 +1,103 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("abs_bound_lt_2exp_si....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arf_t x;
+        fmpz_t b;
+        long c;
+
+        arf_init(x);
+        fmpz_init(b);
+
+        if (iter < 5)
+        {
+            arf_set_ui_2exp_si(x, 1, LONG_MIN + 2);
+            arf_mul_2exp_si(x, x, -iter);
+        }
+        else if (iter < 10)
+        {
+            arf_set_ui_2exp_si(x, 1, LONG_MAX - 2);
+            arf_mul_2exp_si(x, x, iter - 5);
+        }
+        else
+        {
+            arf_randtest_special(x, state, 2 + n_randint(state, 1000), 100);
+        }
+
+        arf_abs_bound_lt_2exp_fmpz(b, x);
+        c = arf_abs_bound_lt_2exp_si(x);
+
+        if (c == -ARF_PREC_EXACT)
+        {
+            if (!(arf_is_zero(x) || fmpz_cmp_si(b, -ARF_PREC_EXACT) <= 0))
+            {
+                printf("FAIL (small/zero)\n\n");
+                abort();
+            }
+        }
+        else if (c == ARF_PREC_EXACT)
+        {
+            if (!(arf_is_inf(x) || arf_is_nan(x) ||
+                fmpz_cmp_si(b, ARF_PREC_EXACT) >= 0))
+            {
+                printf("FAIL (large/inf/nan)\n\n");
+                abort();
+            }
+        }
+        else
+        {
+            if (fmpz_cmp_si(b, c) != 0)
+            {
+                printf("FAIL (normal)\n\n");
+                printf("x = "); arf_print(x); printf("\n\n");
+                printf("b = "); fmpz_print(b); printf("\n\n");
+                printf("c = %ld\n\n", c);
+                abort();
+            }
+        }
+
+        arf_clear(x);
+        fmpz_clear(b);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arf/test/t-set_fmpr.c

@@ -0,0 +1,110 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("set_fmpr....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* test exact roundtrip R -> Q -> R */
+    for (iter = 0; iter < 100000; iter++)
+    {
+        long bits;
+        arf_t x, z;
+        fmpr_t y;
+
+        bits = 2 + n_randint(state, 1000);
+
+        arf_init(x);
+        arf_init(z);
+        fmpr_init(y);
+
+        arf_randtest_special(x, state, bits, 1 + n_randint(state, 100));
+        arf_randtest_special(z, state, bits, 1 + n_randint(state, 100));
+
+        arf_get_fmpr(y, x);
+        arf_set_fmpr(z, y);
+
+        if (!arf_equal(x, z))
+        {
+            printf("FAIL\n\n");
+            printf("bits: %ld\n", bits);
+            printf("x = "); arf_print(x); printf("\n\n");
+            printf("y = "); fmpr_print(y); printf("\n\n");
+            printf("z = "); arf_print(z); printf("\n\n");
+            abort();
+        }
+
+        arf_clear(x);
+        arf_clear(z);
+        fmpr_clear(y);
+    }
+
+    /* test exact roundtrip Q -> R -> Q */
+    for (iter = 0; iter < 100000; iter++)
+    {
+        long bits;
+        fmpr_t x, z;
+        arf_t y;
+
+        bits = 2 + n_randint(state, 1000);
+
+        fmpr_init(x);
+        fmpr_init(z);
+        arf_init(y);
+
+        fmpr_randtest_special(x, state, bits, 1 + n_randint(state, 100));
+        fmpr_randtest_special(z, state, bits, 1 + n_randint(state, 100));
+
+        arf_set_fmpr(y, x);
+        arf_get_fmpr(z, y);
+
+        if (!fmpr_equal(x, z))
+        {
+            printf("FAIL (2)\n\n");
+            printf("bits: %ld\n", bits);
+            printf("x = "); fmpr_print(x); printf("\n\n");
+            printf("y = "); arf_print(y); printf("\n\n");
+            printf("z = "); fmpr_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmpr_clear(x);
+        fmpr_clear(z);
+        arf_clear(y);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arf/test/t-set_fmpz_2exp.c

@@ -0,0 +1,79 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("set_fmpz_2exp....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* test exact roundtrip R -> Q -> R */
+    for (iter = 0; iter < 100000; iter++)
+    {
+        long bits;
+        arf_t x, z;
+        fmpz_t y, e;
+
+        bits = 2 + n_randint(state, 200);
+
+        arf_init(x);
+        arf_init(z);
+        fmpz_init(y);
+        fmpz_init(e);
+
+        arf_randtest(x, state, bits, 1 + n_randint(state, 100));
+        arf_randtest(z, state, bits, 1 + n_randint(state, 100));
+
+        arf_get_fmpz_2exp(y, e, x);
+        arf_set_fmpz_2exp(z, y, e);
+
+        if (!arf_equal(x, z))
+        {
+            printf("FAIL\n\n");
+            printf("bits: %ld\n", bits);
+            printf("x = "); arf_print(x); printf("\n\n");
+            printf("y = "); fmpz_print(y); printf("\n\n");
+            printf("e = "); fmpz_print(e); printf("\n\n");
+            printf("z = "); arf_print(z); printf("\n\n");
+            abort();
+        }
+
+        arf_clear(x);
+        arf_clear(z);
+        fmpz_clear(y);
+        fmpz_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}