preliminary version of a ball type, with multiplication implemented

mprb/Makefile

@@ -0,0 +1,35 @@
+SOURCES = $(wildcard *.c)
+
+OBJS = $(patsubst %.c, $(BUILD_DIR)/%.o, $(SOURCES))
+
+LIB_OBJS = $(patsubst %.c, $(BUILD_DIR)/%.lo, $(SOURCES))
+
+TEST_SOURCES = $(wildcard test/*.c)
+
+PROF_SOURCES = $(wildcard profile/*.c)
+
+TESTS = $(patsubst %.c, %, $(TEST_SOURCES))
+
+PROFS = $(patsubst %.c, %, $(PROF_SOURCES))
+
+all: $(OBJS)
+
+library: $(LIB_OBJS)
+
+profile:
+	$(foreach prog, $(PROFS), $(CC) -O2 -std=c99 $(INCS) $(prog).c ../profiler.o -o $(BUILD_DIR)/$(prog) $(LIBS) -lflint -larb;)
+        
+$(BUILD_DIR)/%.o: %.c
+	$(CC) $(CFLAGS) -c $(INCS) $< -o $@
+
+$(BUILD_DIR)/%.lo: %.c
+	$(CC) -fPIC $(CFLAGS) $(INCS) -c $< -o $@
+
+clean:
+	rm -rf $(BUILD_DIR)	
+
+check: library
+	$(foreach prog, $(TESTS), $(CC) $(CFLAGS) $(INCS) $(prog).c -o $(BUILD_DIR)/$(prog) $(LIBS) -lflint -larb;)
+	$(foreach prog, $(TESTS), $(BUILD_DIR)/$(prog);)
+
+.PHONY: profile clean check all

mprb/add_rad_ui_2exp.c

@@ -0,0 +1,8 @@
+#include "mprb.h"
+
+void
+add_rad_ui_2exp(mprb_t x, mp_limb_t rad, long exp)
+{
+
+}
+

mprb/clear.c

@@ -0,0 +1,7 @@
+#include "mprb.h"
+
+void
+mprb_clear(mprb_t x)
+{
+    free(x->d);
+}

mprb/contains_mpfr.c

@@ -0,0 +1,21 @@
+#include "mprb.h"
+
+int
+mprb_contains_mpfr(const mprb_t x, const mpfr_t f)
+{
+    mpfr_t a, b;
+    int result;
+
+    /* 1 extra bit allows adding radius exactly */
+    mpfr_init2(a, FLINT_BITS * x->size + 1);
+    mpfr_init2(b, FLINT_BITS * x->size + 1);
+
+    mprb_get_interval_mpfr(a, b, x);
+
+    result = (mpfr_cmp(a, f) <= 0) && (mpfr_cmp(f, b) <= 0);
+
+    mpfr_clear(a);
+    mpfr_clear(b);
+
+    return result;
+}

mprb/debug.c

@@ -0,0 +1,18 @@
+#include "mprb.h"
+
+void
+mprb_debug(const mprb_t x)
+{
+    long e;
+    mpz_t t;
+
+    mpz_init(t);
+    e = mprb_get_mid_mpz_2exp(t, x);
+
+    printf("N: [exp=%ld] ", x->exp);
+    mpn_debug(x->d, x->size);
+    gmp_printf("Z: {mid=%Zx, exp=%ld, rad=%lu, size=%ld, alloc=%ld}\n",
+        t, e, x->rad, (long) x->size, (long) x->alloc);
+
+    mpz_clear(t);
+}

mprb/get_interval_mpfr.c

@@ -0,0 +1,30 @@
+#include "mprb.h"
+
+void
+mprb_get_interval_mpfr(mpfr_t a, mpfr_t b, const mprb_t x)
+{
+    mpfr_t r;
+    mpfr_init2(r, FLINT_BITS);
+
+    mprb_get_mid_mpfr(a, x, MPFR_RNDD);
+    mprb_get_mid_mpfr(b, x, MPFR_RNDU);
+
+    mprb_get_rad_mpfr(r, x);
+
+#if 0
+    printf("GETTING INTRVLA\n");
+    mpfr_printf("a = %Rg [xsign = %ld]\n", a, x->sign);
+    mpfr_printf("b = %Rg\n", b);
+    mpfr_printf("r = %Rg\n", r);
+#endif
+
+    mpfr_sub(a, a, r, MPFR_RNDD);
+    mpfr_add(b, b, r, MPFR_RNDU);
+
+#if 0
+    mpfr_printf("a2 = %Rg\n", a);
+    mpfr_printf("b2 = %Rg\n", b);
+#endif
+
+    mpfr_clear(r);
+}

mprb/get_mid_mpfr.c

@@ -0,0 +1,12 @@
+#include "mprb.h"
+
+/* XXX: signs! */
+void
+mprb_get_mid_mpfr(mpfr_t x, const mprb_t v, mpfr_rnd_t rnd)
+{
+    if ((v->size == 1) && (v->d[0] == 0))
+        mpfr_set_ui(x, 0, MPFR_RNDD);
+    else
+        _mpr_get_mpfr_signed(x, v->d, v->exp, v->size,
+            (v->sign == MPRB_SIGN_PLUS) ? 1 : -1, rnd);
+}

mprb/get_mid_mpz_2exp.c

@@ -0,0 +1,19 @@
+#include "mprb.h"
+
+/* todo: handle 0 */
+long
+mprb_get_mid_mpz_2exp(mpz_t a, const mprb_t x)
+{
+    if ((x->size == 1) && (x->d[0] == 0))
+    {
+        mpz_set_ui(a, 0);
+        return 0;
+    }
+    else
+    {
+        mpz_realloc2(a, x->size * FLINT_BITS);
+        mpn_copyi(a->_mp_d, x->d, x->size);
+        a->_mp_size = (x->sign == MPRB_SIGN_PLUS) ? x->size : -(x->size);
+        return x->exp - x->size * FLINT_BITS;
+    }
+}

mprb/get_rad_mpfr.c

@@ -0,0 +1,7 @@
+#include "mprb.h"
+
+void
+mprb_get_rad_mpfr(mpfr_t r, const mprb_t x)
+{
+    mpfr_set_ui_2exp(r, x->rad, x->exp - FLINT_BITS * x->size, MPFR_RNDU);
+}

mprb/init.c

@@ -0,0 +1,14 @@
+#include "mprb.h"
+
+void
+mprb_init(mprb_t x, long bits)
+{
+    long limbs = _MPR_BITS_TO_LIMBS(bits);
+
+    x->d = calloc(limbs, sizeof(mp_limb_t));
+    x->exp = 0;
+    x->sign = MPRB_SIGN_PLUS;
+    x->alloc = limbs;
+    x->size = limbs;
+    x->rad = 0UL;
+}

mprb/mul.c

@@ -0,0 +1,114 @@
+#include "mprb.h"
+
+/* throw away to avoid overflow */
+#define RAD_MUL_SHIFT 8
+
+void
+mprb_mul(mprb_t z, const mprb_t x, const mprb_t y)
+{
+    mp_limb_t xrad, yrad;
+    mp_srcptr xptr, yptr;
+    long size;
+
+    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;
+
+    if (size == 1)
+    {
+        mp_limb_t u1, u0, a2, a1, a0, b1, b0;
+
+        umul_ppmm(a1, a0, xptr[0], yrad);
+        umul_ppmm(b1, b0, yptr[0], xrad);
+        add_sssaaaaaa(a2, a1, a0, 0, a1, a0, 0, b1, b0);
+        umul_ppmm(b1, b0, xrad, yrad);
+        add_sssaaaaaa(a2, a1, a0, a2, a1, a0, 0, b1, b0);
+
+        umul_ppmm(u1, u0, xptr[0], yptr[0]);
+
+        if ((u1 >> (FLINT_BITS - 1)) == 0)
+        {
+            u1 = (u1 << 1) | (u0 >> (FLINT_BITS - 1));
+            u0 = (u0 << 1);
+            a2 = (a2 << 1) | (a1 >> (FLINT_BITS - 1));
+            a1 = (a1 << 1) | (a0 >> (FLINT_BITS - 1));
+            a0 = (a0 << 1);
+            z->exp = x->exp + y->exp - 1;
+        }
+        else
+        {
+            z->exp = x->exp + y->exp;
+        }
+
+        /* round up */
+        add_ssaaaa(a2, a1, a2, a1, 0, (u0 != 0) + (a0 != 0));
+
+        if (a2 == 0)
+        {
+            z->rad = a1;
+            z->d[0] = u1;
+        }
+        else
+        {
+            z->rad = -1UL;
+            z->d[0] = (1UL << (FLINT_BITS - 1));
+            z->exp += 2;
+        }
+
+        z->size = 1;
+        z->sign = x->sign ^ y->sign;
+    }
+    else
+    {
+        mp_limb_t t0, t1, u0, u1, rad;
+        mp_limb_t xtop, ytop;
+        int shift, need_limb_shift;
+
+        /* upper bound for x*yrad + y*xrad */
+        xtop = (xptr[size-1] >> RAD_MUL_SHIFT) + 1UL;
+        ytop = (yptr[size-1] >> RAD_MUL_SHIFT) + 1UL;
+        umul_ppmm(t1, t0, xtop, y->rad);
+        umul_ppmm(u1, u0, ytop, x->rad);
+        add_ssaaaa(t1, t0, t1, t0, u1, u0);
+
+        /* scale back radius to the true magnitude */
+        t0 = (t0 >> (FLINT_BITS - RAD_MUL_SHIFT)) | (t1 << RAD_MUL_SHIFT);
+        t1 = (t1 >> (FLINT_BITS - RAD_MUL_SHIFT));
+
+        /* xrad * yrad < 1 ulp
+            + 1 ulp error in the main multiplication
+            + 1 ulp error from the shift above */
+        add_ssaaaa(t1, t0, t1, t0, 0, 3UL);
+
+        shift = _mpr_muld_n(z->d, xptr, yptr, size);
+
+        /* adjust radius */
+        t1 = (t1 << shift) | ((t0 >> (FLINT_BITS - 1)) & shift);
+        t0 = t0 << shift;
+
+        /* remove one limb */
+        if (t1 != 0)
+        {
+            need_limb_shift = 1;
+
+            /* +1 ulp error for the rounding of t1 */
+            /* +1 ulp error in the main product, which doubles
+                if we perform a shift */
+            rad = t1 + 2 + shift;
+            mpn_copyi(z->d, z->d + 1, size - 1);
+        }
+        else
+        {
+            need_limb_shift = 0;
+            rad = t0;
+        }
+
+        z->rad = rad;
+        z->exp = x->exp + y->exp - shift;
+        z->size = size - need_limb_shift;
+        z->sign = x->sign ^ y->sign;
+    }
+}

mprb/randtest.c

@@ -0,0 +1,16 @@
+#include "mprb.h"
+
+void
+mprb_randtest(mprb_t x, flint_rand_t state, long emin, long emax)
+{
+    long n;
+
+    n = n_randint(state, x->alloc) + 1;
+
+    _mpr_randtest(x->d, state, n);
+
+    x->rad = n_randtest(state);
+    x->size = n;
+    x->sign = n_randint(state, 2) ? MPRB_SIGN_PLUS : MPRB_SIGN_MINUS;
+    x->exp = emin + n_randint(state, emax - emin + 1);
+}

mprb/set_mpfr.c

@@ -0,0 +1,16 @@
+#include "mprb.h"
+
+/* todo: handle zero */
+void
+mprb_set_mpfr(mprb_t x, const mpfr_t v)
+{
+    long xprec, vprec, prec;
+
+    xprec = x->alloc;
+    vprec = _MPR_BITS_TO_LIMBS(v->_mpfr_prec);
+
+    x->exp = _mpr_set_mpfr(x->d, v, xprec);
+    x->rad = (xprec >= vprec) ? 0UL : 1UL;
+    x->sign = (v->_mpfr_sign == 1) ? MPRB_SIGN_PLUS : MPRB_SIGN_MINUS;
+    x->size = x->alloc;
+}

mprb/test/t-mul.c

@@ -0,0 +1,77 @@
+#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("mul....");
+    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);
+
+        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_mul(z, x, y);
+        mpfr_mul(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;
+}

mprb/test/t-set_mpfr.c

@@ -0,0 +1,50 @@
+#include "mprb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("set_mpfr....");
+    fflush(stdout);
+
+    flint_randinit(state);
+    _flint_rand_init_gmp(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        mprb_t x;
+        mpfr_t t;
+
+        mprb_init(x, 1 + n_randint(state, 400));
+        /* mprb_randtest(x, state); */
+
+        mpfr_init2(t, 2 + n_randint(state, 400));
+
+        mpfr_urandom(t, state->gmp_state, MPFR_RNDN);
+        if (n_randint(state, 2))
+            mpfr_neg(t, t, MPFR_RNDN);
+
+        mpfr_mul_2exp(t, t, n_randint(state, 20), MPFR_RNDN);
+        mpfr_div_2exp(t, t, n_randint(state, 20), MPFR_RNDN);
+
+        mprb_set_mpfr(x, t);
+
+        if (!mprb_contains_mpfr(x, t))
+        {
+            printf("FAIL!\n");
+            mprb_debug(x);
+            mpfr_printf("\n%Rf\n", t);
+            mpn_debug(t->_mpfr_d, (t->_mpfr_prec + FLINT_BITS - 1) / FLINT_BITS);
+            abort();
+        }
+
+        mpfr_clear(t);
+        mprb_clear(x);
+    }
+
+    printf("PASS\n");
+    flint_randclear(state);
+    _fmpz_cleanup();
+    return EXIT_SUCCESS;
+}