low-level set_round/normalise code

fmpr.h

@@ -155,8 +155,7 @@ fmpr_one(fmpr_t x)
 
 /* ------------------------------------------------------------------------ */
 
-/* TODO: version for just the val2 reduction! */
-long _fmpr_normalise(fmpz_t man, fmpz_t exp, long prec, fmpr_rnd_t rnd);
+long _fmpr_normalise_naive(fmpz_t man, fmpz_t exp, long prec, fmpr_rnd_t rnd);
 
 static __inline__ void
 fmpr_set(fmpr_t y, const fmpr_t x)
@@ -168,8 +167,26 @@ fmpr_set(fmpr_t y, const fmpr_t x)
     }
 }
 
+long _fmpr_set_round(fmpz_t rman, fmpz_t rexp,
+    const fmpz_t man, const fmpz_t exp, long prec, fmpr_rnd_t rnd);
+
 static __inline__ long
-fmpr_set_round(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
+_fmpr_normalise(fmpz_t man, fmpz_t exp, long prec, fmpr_rnd_t rnd)
+{
+    if (fmpz_is_zero(man))
+    {
+        fmpz_zero(man);
+        fmpz_zero(exp);
+        return FMPR_PREC_EXACT;
+    }
+    else
+    {
+        return _fmpr_set_round(man, exp, man, exp, prec, rnd);
+    }
+}
+
+static __inline__ long
+fmpr_set_round_naive(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
 {
     fmpr_set(y, x);
     if (fmpr_is_special(y))
@@ -178,6 +195,22 @@ fmpr_set_round(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
         return _fmpr_normalise(fmpr_manref(y), fmpr_expref(y), prec, rnd);
 }
 
+static __inline__ long
+fmpr_set_round(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
+{
+    if (fmpr_is_special(x))
+    {
+        fmpr_set(y, x);
+        return FMPR_PREC_EXACT;
+    }
+    else
+    {
+        return _fmpr_set_round(fmpr_manref(y), fmpr_expref(y),
+            fmpr_manref(x), fmpr_expref(x), prec, rnd);
+    }
+}
+
+
 
 static __inline__ int
 fmpr_equal(const fmpr_t x, const fmpr_t y)

fmpr/normalise.c

@@ -25,7 +25,8 @@
 
 #include "fmpr.h"
 
-long _fmpr_normalise(fmpz_t man, fmpz_t exp, long prec, fmpr_rnd_t rnd)
+long
+_fmpr_normalise_naive(fmpz_t man, fmpz_t exp, long prec, fmpr_rnd_t rnd)
 {
     /* TODO: this should perhaps raise an exception to avoid ambiguity */
     if (fmpz_is_zero(man))

fmpr/set_round.c

@@ -0,0 +1,222 @@
+/*=============================================================================
+
+    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 "fmpr.h"
+
+static __inline__ void
+_fmpz_add_ui(fmpz_t z, const fmpz_t x, ulong y)
+{
+    fmpz f = *x;
+
+    if (!COEFF_IS_MPZ(f) && y <= COEFF_MAX)
+        fmpz_set_si(z, f + y);
+    else
+        fmpz_add_ui(z, x, y);
+}
+
+static __inline__ int
+rounds_up(fmpr_rnd_t rnd, int negative)
+{
+    if (rnd == FMPR_RND_DOWN) return 0;
+    if (rnd == FMPR_RND_UP) return 1;
+    if (rnd == FMPR_RND_FLOOR) return negative;
+    return !negative;
+}
+
+/* like mpn_scan0b, but takes an upper size */
+static __inline__ mp_bitcnt_t
+mpn_scan0b(mp_srcptr up, mp_size_t size, mp_bitcnt_t from_bit)
+{
+    mp_limb_t t;
+    long i, c;
+
+    i = from_bit / GMP_NUMB_BITS;
+    c = from_bit % FLINT_BITS;
+    t = ((~up[i]) >> c) << c;
+
+    while (t == 0)
+    {
+        i++;
+        if (i == size)
+            return size * FLINT_BITS;
+        else
+            t = ~up[i];
+    }
+
+    count_trailing_zeros(c, t);
+    return (i * FLINT_BITS) + c;
+}
+
+long
+_fmpr_set_round(fmpz_t rman, fmpz_t rexp,
+    const fmpz_t man, const fmpz_t exp, long prec, fmpr_rnd_t rnd)
+{
+    if (!COEFF_IS_MPZ(*man))
+    {
+        long lead, trail, bc, v, w, shift, ret;
+        int negative;
+
+        v = *man;
+
+        count_trailing_zeros(trail, v);
+        v >>= trail;
+        shift = trail;
+        ret = FMPR_PREC_EXACT;
+
+        /* may need to round */
+        if (prec < FLINT_BITS - 2 - trail)
+        {
+            if (v < 0)
+            {
+                negative = 1;
+                w = -v;
+            }
+            else
+            {
+                negative = 0;
+                w = v;
+            }
+
+            count_leading_zeros(lead, w);
+            bc = FLINT_BITS - lead;
+
+            /* round */
+            if (prec < bc)
+            {
+                w = (w >> (bc - prec)) + rounds_up(rnd, negative);
+                shift += bc - prec;
+                count_trailing_zeros(trail, w);
+                w >>= trail;
+                shift += trail;
+                v = negative ? -w : w;
+                ret = trail;
+            }
+        }
+
+        /* the input is small, so the output must be small too */
+        _fmpz_demote(rman);
+        *rman = v;
+        _fmpz_add_ui(rexp, exp, shift);
+        return ret;
+    }
+    else
+    {
+        long size, bc, val, val_bits, val_limbs, ret, old_val;
+        int negative, increment;
+        mp_ptr d;
+        __mpz_struct * z = COEFF_TO_PTR(*man);
+
+        size = z->_mp_size;
+        negative = size < 0;
+        size = FLINT_ABS(size);
+        d = z->_mp_d;
+
+        /* bit size */
+        count_leading_zeros(bc, d[size - 1]);
+        bc = FLINT_BITS - bc;
+        bc += (size - 1) * FLINT_BITS;
+
+        /* quick exit */
+        if (bc <= prec && (d[0] & 1))
+        {
+            if (rman != man) fmpz_set(rman, man);
+            if (rexp != exp) fmpz_set(rexp, exp);
+            return FMPR_PREC_EXACT;
+        }
+
+        /* trailing zeros */
+        val_limbs = 0;
+        while (d[val_limbs] == 0)
+            val_limbs++;
+        count_trailing_zeros(val_bits, d[val_limbs]);
+        val = val_bits + (val_limbs * FLINT_BITS);
+
+        /* no rounding necessary; just copy or shift to destination */
+        if (bc - val <= prec)
+        {
+            ret = FMPR_PREC_EXACT;
+            increment = 0;
+        }
+        else
+        {
+            /* truncation */
+            if (!rounds_up(rnd, negative))
+            {
+                old_val = val;
+                val = mpn_scan1(d, bc - prec);
+                increment = 0;
+            }
+            /* round to next higher odd mantissa */
+            else
+            {
+                old_val = val;
+                val = mpn_scan0b(d, size, bc - prec);
+
+                /* can overflow to next power of 2 */
+                if (val == bc)
+                {
+                    fmpz_set_si(rman, negative ? -1 : 1);
+                    _fmpz_add_ui(rexp, exp, bc);
+                    return prec;
+                }
+
+                /* otherwise, incrementing will not cause overflow below */
+                increment = 1;
+            }
+
+            val_limbs = val / FLINT_BITS;
+            val_bits = val % FLINT_BITS;
+            ret = prec - (bc - val);
+        }
+
+        /* the output mantissa is a small fmpz */
+        if (bc - val <= FLINT_BITS - 2)
+        {
+            mp_limb_t h;
+            if (val_limbs + 1 == size || val_bits == 0)
+                h = d[val_limbs] >> val_bits;
+            else
+                h = (d[val_limbs] >> val_bits)
+                    | (d[val_limbs + 1] << (FLINT_BITS - val_bits));
+            h += increment;
+            _fmpz_demote(rman);
+            *rman = negative ? -h : h;
+        }
+        /* the output mantissa is an mpz */
+        else
+        {
+            if (rman == man)
+                mpz_tdiv_q_2exp(z, z, val);
+            else
+                mpz_tdiv_q_2exp(_fmpz_promote(rman), z, val);
+            if (increment)
+                COEFF_TO_PTR(*rman)->_mp_d[0]++;
+        }
+
+        _fmpz_add_ui(rexp, exp, val);
+        return ret;
+    }
+}
+

fmpr/test/t-normalise.c

@@ -0,0 +1,85 @@
+/*=============================================================================
+
+    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 "fmpr.h"
+
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("normalise....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000000; iter++)
+    {
+        fmpr_t x, y;
+        long prec, ret1, ret2;
+        fmpr_rnd_t rnd;
+
+        fmpr_init(x);
+        fmpr_init(y);
+
+        fmpz_randtest(fmpr_manref(x), state, 2000);
+        fmpz_randtest(fmpr_expref(x), state, 200);
+
+        fmpz_set(fmpr_manref(y), fmpr_manref(x));
+        fmpz_set(fmpr_expref(y), fmpr_expref(x));
+
+        switch (n_randint(state, 4))
+        {
+            case 0: rnd = FMPR_RND_DOWN; break;
+            case 1: rnd = FMPR_RND_UP; break;
+            case 2: rnd = FMPR_RND_FLOOR; break;
+            default: rnd = FMPR_RND_CEIL; break;
+        }
+
+        prec = 2 + n_randint(state, 2000);
+
+        ret1 = _fmpr_normalise(fmpr_manref(x), fmpr_expref(x), prec, rnd);
+        ret2 = _fmpr_normalise_naive(fmpr_manref(y), fmpr_expref(y), prec, rnd);
+
+        if (!fmpr_equal(x, y) || ret1 != ret2)
+        {
+            printf("FAIL!\n");
+            printf("x = "); fmpr_print(x); printf("\n\n");
+            printf("y = "); fmpr_print(y); printf("\n\n");
+            printf("ret1 = %ld, ret2 = %ld\n", ret1, ret2);
+            abort();
+        }
+
+        fmpr_clear(x);
+        fmpr_clear(y);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+