continuation of last commit, adding addition

Makefile.in

@@ -102,5 +102,5 @@ build/%.lo: %.c
 build/%.o: %.c
 	$(CC) -fPIC $(CFLAGS) $(INCS) -c $< -o $@
 
-BUILD_DIRS = mpr
+BUILD_DIRS = mpr mprb
 

mpr.h

@@ -5,13 +5,16 @@
 #include "flint.h"
 #include "ulong_extras.h"
 #include "fmpz.h"
+#include "mpn_extras.h"
 #include "mpn_inlines.h"
 
 #define _MPR_BITS_TO_LIMBS(b) (((b) + FLINT_BITS - 1) / FLINT_BITS)
 
 /* Conversions / debugging ***************************************************/
-long _mpr_set_mpfr(mp_ptr y, const mpfr_t x, mp_size_t n, mpfr_rnd_t rnd);
+long _mpr_set_mpfr(mp_ptr y, const mpfr_t x, mp_size_t n);
 void _mpr_get_mpfr(mpfr_t y, mp_srcptr x, long exp, mp_size_t n, mpfr_rnd_t rnd);
+void _mpr_get_mpfr_signed(mpfr_t y, mp_srcptr x, long exp,
+    mp_size_t n, int sign, mpfr_rnd_t rnd);
 
 void _mpr_randtest(mp_ptr x, flint_rand_t state, mp_size_t n);
 

mpr/get_mpfr.c

@@ -11,3 +11,16 @@ _mpr_get_mpfr(mpfr_t y, mp_srcptr x, long exp, mp_size_t n, mpfr_rnd_t rnd)
 
     mpfr_set_z_2exp(y, &z, exp - n * FLINT_BITS, rnd);
 }
+
+void
+_mpr_get_mpfr_signed(mpfr_t y, mp_srcptr x, long exp,
+    mp_size_t n, int sign, mpfr_rnd_t rnd)
+{
+    __mpz_struct z;
+
+    z._mp_size = (sign > 0) ? n : -n;
+    z._mp_alloc = n;
+    z._mp_d = (mp_ptr) x;
+
+    mpfr_set_z_2exp(y, &z, exp - n * FLINT_BITS, rnd);
+}

mpr/randtest.c

@@ -4,6 +4,8 @@ void
 _mpr_randtest(mp_ptr x, flint_rand_t state, mp_size_t n)
 {
     __mpz_struct z;
+    mp_limb_t t;
+    long i;
 
     _flint_rand_init_gmp(state);
 
@@ -11,9 +13,15 @@ _mpr_randtest(mp_ptr x, flint_rand_t state, mp_size_t n)
     z._mp_alloc = n;
     z._mp_d = x;
 
-if (0)    mpz_urandomb(&z, state->gmp_state, n * FLINT_BITS);
+    t = n_randlimb(state);
 
-    mpn_random2(x, n);
+    if (t & 1)
+        mpz_urandomb(&z, state->gmp_state, n * FLINT_BITS);
+    else
+        mpz_rrandomb(&z, state->gmp_state, n * FLINT_BITS);
+
+    for (i = z._mp_size; i < n; i++)
+        x[i] = 0;
 
     x[n - 1] |= (1UL << (FLINT_BITS - 1));
 }

mpr/set_mpfr.c

@@ -1,19 +1,19 @@
 #include "mpr.h"
 
 long
-_mpr_set_mpfr(mp_ptr y, const mpfr_t x, mp_size_t n, mpfr_rnd_t rnd)
+_mpr_set_mpfr(mp_ptr y, const mpfr_t x, mp_size_t n)
 {
     long xlimbs, exp, m;
 
-    if (!mpfr_regular_p(x) || mpfr_sgn(x) != 1)
+    if (!mpfr_regular_p(x))
     {
-        printf("_mpr_set_mpfr: input must be a positive regular number\n");
+        printf("_mpr_set_mpfr: input must be a regular number\n");
         abort();
     }
 
     xlimbs = _MPR_BITS_TO_LIMBS(mpfr_get_prec(x));
 
-    if (n >= xlimbs || rnd == MPFR_RNDD || rnd == MPFR_RNDZ)
+    if (n >= xlimbs)
     {
         m = FLINT_MIN(xlimbs, n);
         mpn_copyi(y + n - m, x->_mpfr_d + m - xlimbs, m);
@@ -24,7 +24,7 @@ _mpr_set_mpfr(mp_ptr y, const mpfr_t x, mp_size_t n, mpfr_rnd_t rnd)
     {
         mpfr_t t;
         mpfr_init2(t, n * FLINT_BITS);
-        mpfr_set(t, x, rnd);
+        mpfr_set(t, x, MPFR_RNDZ);
         mpn_copyi(y, t->_mpfr_d, n);
         exp = mpfr_get_exp(t);
         mpfr_clear(t);

mpr/test/t-addd_n.c

@@ -0,0 +1,76 @@
+#include "mpr.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("addd_n....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000000; iter++)
+    {
+        mp_ptr x, y, z;
+        mpfr_t t, u, v, w;
+        long n, xexp, yexp, zexp;
+
+        n = 1 + n_randint(state, 40);
+
+        x = malloc(n * sizeof(mp_limb_t));
+        y = malloc(n * sizeof(mp_limb_t));
+        z = malloc(n * sizeof(mp_limb_t));
+
+        mpfr_init2(t, n * FLINT_BITS);
+        mpfr_init2(u, n * FLINT_BITS);
+        mpfr_init2(v, n * FLINT_BITS);
+        mpfr_init2(w, n * FLINT_BITS);
+
+        _mpr_randtest(x, state, n);
+        _mpr_randtest(y, state, n);
+        _mpr_randtest(z, state, n);
+
+        xexp = n_randint(state, 10) - 20;
+        yexp = n_randint(state, 10) - 20;
+
+        _mpr_get_mpfr(t, x, xexp, n, MPFR_RNDN);
+        _mpr_get_mpfr(u, y, yexp, n, MPFR_RNDN);
+
+        if (xexp >= yexp)
+            zexp = xexp + _mpr_addd_n(z, x, y, xexp - yexp, n);
+        else
+            zexp = yexp + _mpr_addd_n(z, y, x, yexp - xexp, n);
+
+        _mpr_get_mpfr(w, z, zexp, n, MPFR_RNDN);
+
+        mpfr_add(v, t, u, MPFR_RNDD);
+
+        if (mpfr_cmp(v, w) != 0)
+        {
+            printf("FAIL!\n");
+            _mpr_printr(x, n);
+            _mpr_printr(y, n);
+            _mpr_printr(z, n);
+
+            mpfr_printf("\n%.80Rf\n", t);
+            mpfr_printf("\n%.80Rf\n", u);
+            mpfr_printf("\n%.80Rf\n", v);
+            mpfr_printf("\n%.80Rf\n", w);
+            abort();
+        }
+
+        mpfr_clear(t);
+        mpfr_clear(u);
+        mpfr_clear(v);
+        mpfr_clear(w);
+        free(x);
+        free(y);
+        free(z);
+    }
+
+    printf("PASS\n");
+
+    flint_randclear(state);
+    return EXIT_SUCCESS;
+}

mpr/test/t-muld_n.c

@@ -0,0 +1,71 @@
+#include "mpr.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("muld_n....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        mp_ptr x, y, z;
+        mpfr_t t, u, v, w;
+        long n, xexp, yexp, zexp;
+
+        n = 1 + n_randint(state, 30);
+
+        x = malloc(n * sizeof(mp_limb_t));
+        y = malloc(n * sizeof(mp_limb_t));
+        z = malloc(n * sizeof(mp_limb_t));
+
+        mpfr_init2(t, n * FLINT_BITS);
+        mpfr_init2(u, n * FLINT_BITS);
+        mpfr_init2(v, n * FLINT_BITS);
+        mpfr_init2(w, n * FLINT_BITS);
+
+        _mpr_randtest(x, state, n);
+        _mpr_randtest(y, state, n);
+        _mpr_randtest(z, state, n);
+
+        xexp = n_randint(state, 10) - 20;
+        yexp = n_randint(state, 10) - 20;
+
+        _mpr_get_mpfr(t, x, xexp, n, MPFR_RNDN);
+        _mpr_get_mpfr(u, y, yexp, n, MPFR_RNDN);
+
+        zexp = xexp + yexp - _mpr_muld_n(z, x, y, n);
+        _mpr_get_mpfr(w, z, zexp, n, MPFR_RNDN);
+
+        mpfr_mul(v, t, u, MPFR_RNDD);
+
+        if (mpfr_cmp(v, w) != 0)
+        {
+            printf("FAIL!\n");
+            _mpr_printr(x, n);
+            _mpr_printr(y, n);
+            _mpr_printr(z, n);
+
+            mpfr_printf("\n%.80Rf\n", t);
+            mpfr_printf("\n%.80Rf\n", u);
+            mpfr_printf("\n%.80Rf\n", v);
+            mpfr_printf("\n%.80Rf\n", w);
+            abort();
+        }
+
+        mpfr_clear(t);
+        mpfr_clear(u);
+        mpfr_clear(v);
+        mpfr_clear(w);
+        free(x);
+        free(y);
+        free(z);
+    }
+
+    printf("PASS\n");
+    flint_randclear(state);
+    return EXIT_SUCCESS;
+}

mpr/test/t-set_mpfr.c

@@ -0,0 +1,54 @@
+#include "mpr.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("set_mpfr....");
+    fflush(stdout);
+
+    flint_randinit(state);
+    _flint_rand_init_gmp(state);
+
+    /* test exact roundtrip (todo: test rounding) */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        mp_ptr x;
+        long exp;
+        mpfr_t t, u;
+        long n;
+        mpfr_rnd_t rnd;
+
+        n = 1 + n_randint(state, 10);
+        x = malloc(n * sizeof(mp_limb_t));
+        mpfr_init2(t, n * FLINT_BITS);
+        mpfr_init2(u, n * FLINT_BITS);
+
+        mpfr_urandom(t, state->gmp_state, MPFR_RNDN);
+        mpfr_mul_2exp(t, t, n_randint(state, 20), MPFR_RNDN);
+        mpfr_div_2exp(t, t, n_randint(state, 20), MPFR_RNDN);
+
+        rnd = MPFR_RNDD;
+        exp = _mpr_set_mpfr(x, t, n);
+        _mpr_get_mpfr(u, x, exp, n, rnd);
+
+        if (mpfr_cmp(t, u) != 0)
+        {
+            printf("FAIL!\n");
+            mpfr_printf("\n%.200Rf\n", t);
+            mpfr_printf("\n%.200Rf\n", u);
+            abort();
+        }
+
+        mpfr_clear(t);
+        mpfr_clear(u);
+        free(x);
+    }
+
+    printf("PASS\n");
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    return EXIT_SUCCESS;
+}

mprb.h

@@ -0,0 +1,36 @@
+#include "mpr.h"
+
+typedef struct
+{
+    mp_ptr d;
+    mp_limb_t rad;
+    long size;
+    long alloc;
+    long sign;
+    long exp;
+}
+mprb_struct;
+
+typedef mprb_struct mprb_t[1];
+typedef mprb_struct * mprb_ptr;
+typedef const mprb_struct * mprb_srcptr;
+
+#define MPRB_SIGN_PLUS 0
+#define MPRB_SIGN_MINUS 1
+
+void mprb_init(mprb_t x, long bits);
+void mprb_clear(mprb_t x);
+
+void mprb_debug(const mprb_t x);
+void mprb_randtest(mprb_t x, flint_rand_t state, long emin, long emax);
+
+long mprb_get_mid_mpz_2exp(mpz_t a, const mprb_t x);
+
+void mprb_set_mpfr(mprb_t x, const mpfr_t v);
+void mprb_get_rad_mpfr(mpfr_t r, const mprb_t x);
+void mprb_get_mid_mpfr(mpfr_t x, const mprb_t v, mpfr_rnd_t rnd);
+void mprb_get_interval_mpfr(mpfr_t a, mpfr_t b, const mprb_t x);
+int mprb_contains_mpfr(const mprb_t x, const mpfr_t f);
+
+void mprb_add(mprb_t z, const mprb_t x, const mprb_t y);
+void mprb_mul(mprb_t z, const mprb_t x, const mprb_t y);

mprb/add.c

@@ -0,0 +1,85 @@
+#include "mprb.h"
+
+void
+mprb_add(mprb_t z, const mprb_t x, const mprb_t y)
+{
+    mp_limb_t xrad, yrad, t, u;
+    mp_srcptr xptr, yptr;
+    long size, xexp, yexp, shift, carry, exp;
+
+    if (x->sign != y->sign)
+    {
+        printf("mprb_add: different signs not supported\n");
+        abort();
+    }
+
+    size = FLINT_MIN(x->size, y->size);
+    size = FLINT_MIN(size, z->alloc);
+    xrad = (size == x->size) ? x->rad : 2UL;
+    yrad = (size == y->size) ? y->rad : 2UL;
+    xptr = x->d + x->size - size;
+    yptr = y->d + y->size - size;
+
+    xexp = x->exp;
+    yexp = y->exp;
+
+    if (xexp >= yexp)
+    {
+        shift = xexp - yexp;
+        carry = _mpr_addd_n(z->d, xptr, yptr, shift, size);
+        exp = xexp;
+    }
+    else
+    {
+        shift = yexp - xexp;
+        carry = _mpr_addd_n(z->d, yptr, xptr, shift, size);
+        t = xrad;
+        xrad = yrad;
+        yrad = t;
+        exp = yexp;
+    }
+
+    xrad = (xrad >> carry) + carry;
+    yrad = (yrad >> carry) + carry;
+
+    if (shift >= FLINT_BITS)
+    {
+        add_ssaaaa(u, t, 0, xrad, 0, 3);
+    }
+    else if (shift == 0)
+    {
+        add_ssaaaa(u, t, 0, xrad, 0, yrad);
+        add_ssaaaa(u, t, u, t, 0, carry);
+    }
+    else
+    {
+        add_ssaaaa(u, t, 0, xrad, 0, (yrad >> shift) + 3);
+    }
+
+    if (u == 0)
+    {
+        z->rad = t;
+        z->size = size;
+        z->sign = x->sign;
+        z->exp = exp + carry;
+    }
+    else
+    {
+        if (size == 1)
+        {
+            z->rad = 1UL << (FLINT_BITS - 1);
+            z->d[0] = 0;
+            z->size = 1;
+            z->sign = x->sign;
+            z->exp = exp + carry + 3;
+        }
+        else
+        {
+            z->rad = 3UL;
+            mpn_copyi(z->d, z->d + 1, size - 1);
+            z->size = size - 1;
+            z->sign = x->sign;
+            z->exp = exp + carry;
+        }
+    }
+}

mprb/test/t-add.c

@@ -0,0 +1,79 @@
+#include "mprb.h"
+
+void mpfr_debug(mpfr_t t)
+{
+    mpfr_printf("\n%Rf\n", t);
+    mpn_debug(t->_mpfr_d, (t->_mpfr_prec + FLINT_BITS - 1) / FLINT_BITS);
+}
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("add....");
+    fflush(stdout);
+
+    flint_randinit(state);
+    _flint_rand_init_gmp(state);
+
+    for (iter = 0; iter < 1000000; iter++)
+    {
+        mprb_t x, y, z;
+        mpfr_t t, u, v;
+
+        mprb_init(x, 1 + n_randint(state, 300));
+        mprb_init(y, 1 + n_randint(state, 300));
+        mprb_init(z, 1 + n_randint(state, 300));
+
+        do { mprb_randtest(x, state, -10, 10); } while (x->size == 0);
+        do { mprb_randtest(y, state, -10, 10); } while (y->size == 0);
+        do { mprb_randtest(z, state, -10, 10); } while (z->size == 0);
+
+        y->sign = x->sign;
+
+        mpfr_init2(t, 1000);
+        mpfr_init2(u, 1000);
+        mpfr_init2(v, 1000);
+
+        if (n_randint(state, 2))
+            mprb_get_interval_mpfr(t, v, x);
+        else
+            mprb_get_interval_mpfr(v, t, x);
+
+        if (n_randint(state, 2))
+            mprb_get_interval_mpfr(u, v, y);
+        else
+            mprb_get_interval_mpfr(v, u, y);
+
+        mprb_add(z, x, y);
+        mpfr_add(v, t, u, MPFR_RNDN);
+
+        if (!mprb_contains_mpfr(z, v))
+        {
+            printf("FAIL!\n");
+            mprb_debug(x);
+            mprb_debug(y);
+            mprb_debug(z);
+
+            mpfr_debug(t);
+            mpfr_debug(u);
+            mpfr_debug(v);
+
+            abort();
+        }
+
+        mpfr_clear(t);
+        mpfr_clear(u);
+        mpfr_clear(v);
+
+        mprb_clear(x);
+        mprb_clear(y);
+        mprb_clear(z);
+    }
+
+    printf("PASS\n");
+    flint_randclear(state);
+    _fmpz_cleanup();
+    return EXIT_SUCCESS;
+}