From 8eddc357e2018b75976abe16879a779730a7ab5b Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sun, 27 May 2012 19:38:26 +0200 Subject: [PATCH] continuation of last commit, adding addition --- Makefile.in | 2 +- mpr.h | 5 ++- mpr/adds.c | 0 mpr/get_mpfr.c | 13 +++++++ mpr/randtest.c | 12 +++++- mpr/set_mpfr.c | 10 ++--- mpr/test/t-addd_n.c | 76 ++++++++++++++++++++++++++++++++++++++ mpr/test/t-muld_n.c | 71 ++++++++++++++++++++++++++++++++++++ mpr/test/t-set_mpfr.c | 54 +++++++++++++++++++++++++++ mprb.h | 36 ++++++++++++++++++ mprb/add.c | 85 +++++++++++++++++++++++++++++++++++++++++++ mprb/test/t-add.c | 79 ++++++++++++++++++++++++++++++++++++++++ 12 files changed, 434 insertions(+), 9 deletions(-) create mode 100644 mpr/adds.c create mode 100644 mpr/test/t-addd_n.c create mode 100644 mpr/test/t-muld_n.c create mode 100644 mpr/test/t-set_mpfr.c create mode 100644 mprb.h create mode 100644 mprb/add.c create mode 100644 mprb/test/t-add.c diff --git a/Makefile.in b/Makefile.in index 27a39da3..3488683f 100644 --- a/Makefile.in +++ b/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 diff --git a/mpr.h b/mpr.h index a1ff0f59..ed8bd998 100644 --- a/mpr.h +++ b/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); diff --git a/mpr/adds.c b/mpr/adds.c new file mode 100644 index 00000000..e69de29b diff --git a/mpr/get_mpfr.c b/mpr/get_mpfr.c index 59ab88c2..2d3d4831 100644 --- a/mpr/get_mpfr.c +++ b/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); +} diff --git a/mpr/randtest.c b/mpr/randtest.c index 09c66f20..84bc9988 100644 --- a/mpr/randtest.c +++ b/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)); } diff --git a/mpr/set_mpfr.c b/mpr/set_mpfr.c index 54327c9e..588b0486 100644 --- a/mpr/set_mpfr.c +++ b/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); diff --git a/mpr/test/t-addd_n.c b/mpr/test/t-addd_n.c new file mode 100644 index 00000000..ef85ae0d --- /dev/null +++ b/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; +} diff --git a/mpr/test/t-muld_n.c b/mpr/test/t-muld_n.c new file mode 100644 index 00000000..993839c7 --- /dev/null +++ b/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; +} diff --git a/mpr/test/t-set_mpfr.c b/mpr/test/t-set_mpfr.c new file mode 100644 index 00000000..c84cace0 --- /dev/null +++ b/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; +} \ No newline at end of file diff --git a/mprb.h b/mprb.h new file mode 100644 index 00000000..172bc639 --- /dev/null +++ b/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); diff --git a/mprb/add.c b/mprb/add.c new file mode 100644 index 00000000..8056b8d3 --- /dev/null +++ b/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; + } + } +} diff --git a/mprb/test/t-add.c b/mprb/test/t-add.c new file mode 100644 index 00000000..4c7965c3 --- /dev/null +++ b/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; +} \ No newline at end of file