From 68258bd59e7c45daa7f26bd22f6d4accae22848a Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 8 Aug 2012 15:58:06 +0200 Subject: [PATCH] tentative code for log --- mpr.h | 13 +++++++ mprb.h | 48 +++++++++++++++++--------- mprb/add.c | 3 +- mprb/add_rad_ui_2exp.c | 8 ----- mprb/contains_zero.c | 42 ++++++++++++++++++++++ mprb/debug.c | 5 ++- mprb/get_lower_bound_ufloat.c | 36 +++++++++++++++++++ mprb/log_using_mpfr.c | 60 ++++++++++++++++++++++++++++++++ mprb/mul.c | 7 ++-- mprb/test/t-log.c | 65 +++++++++++++++++++++++++++++++++++ ufloat.h | 9 ++--- 11 files changed, 260 insertions(+), 36 deletions(-) delete mode 100644 mprb/add_rad_ui_2exp.c create mode 100644 mprb/contains_zero.c create mode 100644 mprb/get_lower_bound_ufloat.c create mode 100644 mprb/log_using_mpfr.c create mode 100644 mprb/test/t-log.c diff --git a/mpr.h b/mpr.h index ab8dea83..534abe23 100644 --- a/mpr.h +++ b/mpr.h @@ -10,6 +10,19 @@ #define _MPR_BITS_TO_LIMBS(b) (((b) + FLINT_BITS - 1) / FLINT_BITS) +typedef struct +{ + mp_ptr d; + long exp; + long size; + long alloc; + int sign; +} +mpr_struct; + +typedef mpr_struct mpr_t[1]; + + /* Conversions / debugging ***************************************************/ 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); diff --git a/mprb.h b/mprb.h index 883a9432..a9d1e694 100644 --- a/mprb.h +++ b/mprb.h @@ -1,16 +1,6 @@ #include "mpr.h" #include "ufloat.h" -typedef struct -{ - mp_ptr d; - long exp; - long size; - long alloc; - int sign; -} -mpr_struct; - typedef struct { mpr_struct mid; @@ -39,6 +29,7 @@ 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); + /* XXX: assumes nonzero! */ static __inline__ void _mpr_get_ufloat(ufloat_t u, mp_srcptr d, long size, long exp) @@ -55,16 +46,39 @@ _mpr_get_ufloat(ufloat_t u, mp_srcptr d, long size, long exp) } } +void mprb_get_lower_bound_ufloat(ufloat_t u, const mprb_t x); + +static __inline__ int +mprb_is_exact(const mprb_t x) +{ + return ufloat_is_zero(&x->rad); +} + +int mprb_contains_zero(const mprb_t x); + + +static __inline__ long +mprb_ulp_exp(const mprb_t x) +{ + return x->mid.exp - x->mid.size * FLINT_BITS; +} + 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); -/* + +/* Add r*2^exp to the radius of x */ static __inline__ void -_mprb_fit_to_rad(mprb_t x) +mprb_add_rad_ui_2exp(mprb_t x, mp_limb_t r, long exp) { - if (!ufloat_is_zero(&x->rad)) - { - long bot_exp; - } + ufloat_t t; + ufloat_set_ui_2exp(t, r, exp); + ufloat_add(&x->rad, &x->rad, t); +} + +/* Add 2^exp to the radius of x */ +static __inline__ void +mprb_add_rad_2exp(mprb_t x, long exp) +{ + ufloat_add_2exp(&x->rad, &x->rad, exp); } -*/ diff --git a/mprb/add.c b/mprb/add.c index 94f53271..c8ebba46 100644 --- a/mprb/add.c +++ b/mprb/add.c @@ -54,6 +54,5 @@ mprb_add(mprb_t z, const mprb_t x, const mprb_t y) z->mid.size = zsize; - /* rounding error */ - ufloat_add_2exp(&z->rad, &z->rad, z->mid.exp - z->mid.size * FLINT_BITS); + mprb_add_rad_2exp(z, mprb_ulp_exp(z)); } diff --git a/mprb/add_rad_ui_2exp.c b/mprb/add_rad_ui_2exp.c deleted file mode 100644 index 6b02d051..00000000 --- a/mprb/add_rad_ui_2exp.c +++ /dev/null @@ -1,8 +0,0 @@ -#include "mprb.h" - -void -add_rad_ui_2exp(mprb_t x, mp_limb_t rad, long exp) -{ - -} - diff --git a/mprb/contains_zero.c b/mprb/contains_zero.c new file mode 100644 index 00000000..88b06727 --- /dev/null +++ b/mprb/contains_zero.c @@ -0,0 +1,42 @@ +#include "mprb.h" + +int +mprb_contains_zero(const mprb_t x) +{ + mp_limb_t midm, radm; + long i; + + /* centered on zero */ + if (x->mid.d[x->mid.size - 1] == 0) + return 1; + + /* midpoint is larger */ + if (x->mid.exp > x->rad.exp) + return 0; + + /* radius is larger */ + if (x->mid.exp < x->rad.exp) + return 1; + + /* same exponent; compare mantissas */ + midm = x->mid.d[x->mid.size - 1]; + radm = x->rad.man << (FLINT_BITS - UFLOAT_PREC); + + /* midpoint is larger */ + if (midm > radm) + return 0; + + /* radius is larger */ + if (midm < radm) + return 1; + + /* midpoint is larger */ + for (i = x->mid.size - 2; i >= 0; i--) + { + if (x->mid.d[i] != 0) + return 0; + } + + /* equal */ + return 1; +} diff --git a/mprb/debug.c b/mprb/debug.c index d34b0e5a..3cb97a20 100644 --- a/mprb/debug.c +++ b/mprb/debug.c @@ -9,9 +9,12 @@ mprb_debug(const mprb_t x) mpz_init(t); e = mprb_get_mid_mpz_2exp(t, x); +/* printf("N: [exp=%ld] ", x->mid.exp); mpn_debug(x->mid.d, x->mid.size); - gmp_printf("Z: {%Zd * 2^%ld +- %lu * 2^%ld, size=%ld, alloc=%ld}\n", +*/ + + gmp_printf("{%Zd * 2^%ld +- %lu * 2^%ld, size=%ld, alloc=%ld}\n", t, e, x->rad.man, x->rad.exp - UFLOAT_PREC, (long) x->mid.size, (long) x->mid.alloc); mpz_clear(t); diff --git a/mprb/get_lower_bound_ufloat.c b/mprb/get_lower_bound_ufloat.c new file mode 100644 index 00000000..0e047ac3 --- /dev/null +++ b/mprb/get_lower_bound_ufloat.c @@ -0,0 +1,36 @@ +#include "mprb.h" + +/* sets u to a lower bound for mid(x)-rad(x), assuming that this is positive + +XXX: implementation is wrong if rad is zero +*/ + +void +mprb_get_lower_bound_ufloat(ufloat_t u, const mprb_t x) +{ + __mpfr_struct zs, xs, ys; + + mp_limb_t r, s; + + xs._mpfr_d = (mp_ptr) x->mid.d; + xs._mpfr_exp = x->mid.exp; + xs._mpfr_prec = x->mid.size * FLINT_BITS; + xs._mpfr_sign = 1; + + r = x->rad.man << (FLINT_BITS - UFLOAT_PREC); + + ys._mpfr_d = &r; + ys._mpfr_exp = x->rad.exp; + ys._mpfr_prec = UFLOAT_PREC; + ys._mpfr_sign = 1; + + zs._mpfr_d = &s; + zs._mpfr_exp = 0; + zs._mpfr_prec = UFLOAT_PREC; + zs._mpfr_sign = 1; + + mpfr_sub(&zs, &xs, &ys, MPFR_RNDD); + + u->man = s >> (FLINT_BITS - UFLOAT_PREC); + u->exp = zs._mpfr_exp; +} diff --git a/mprb/log_using_mpfr.c b/mprb/log_using_mpfr.c new file mode 100644 index 00000000..901e9b93 --- /dev/null +++ b/mprb/log_using_mpfr.c @@ -0,0 +1,60 @@ +#include "mprb.h" + + +/* +Let the input be [a-b, a+b] * 2^e. We require a > b >= 0 (otherwise the +interval contains zero or a negative number and the logarithm is not +defined). The error is largest at a-b, and we have + +log(a * 2^e) - log((a-b) * 2^e) = log(1 + b/(a-b)). +*/ + +void +mprb_log_error(ufloat_t err, const mprb_t x) +{ + ufloat_t t, u; + + mprb_get_lower_bound_ufloat(t, x); + ufloat_div(u, &x->rad, t); + ufloat_log1p(err, u); +} + +void +mprb_log_using_mpfr(mprb_t y, const mprb_t x) +{ + long prec; + mpfr_t t, u; + int input_approx, value_approx; + ufloat_t err; + + if (mprb_contains_zero(x) || x->mid.sign == MPRB_SIGN_MINUS) + { + printf("mprb_log: interval contains zero or negative numbers\n"); + abort(); + } + + prec = y->mid.alloc * FLINT_BITS; + + mpfr_init2(t, x->mid.size * FLINT_BITS); + mpfr_init2(u, prec); + + mprb_get_mid_mpfr(t, x, MPFR_RNDN); /* exact */ + + input_approx = !mprb_is_exact(x); + + if (input_approx) + mprb_log_error(err, x); + + /* todo: handle exact case (log(1)) */ + mpfr_log(u, t, MPFR_RNDN); + + mprb_set_mpfr(y, u); + + ufloat_set_2exp(&y->rad, mprb_ulp_exp(y)); + + if (input_approx) + ufloat_add(&y->rad, &y->rad, err); + + mpfr_clear(t); + mpfr_clear(u); +} diff --git a/mprb/mul.c b/mprb/mul.c index 4ad67fb1..af6bbc35 100644 --- a/mprb/mul.c +++ b/mprb/mul.c @@ -6,7 +6,7 @@ mprb_mul(mprb_t z, const mprb_t x, const mprb_t y) ufloat_t a, b; mp_srcptr xptr, yptr; mp_ptr zptr; - long zsize, xsize, ysize, exp, shift; + long zsize, xsize, ysize, shift; xptr = x->mid.d; yptr = y->mid.d; @@ -36,13 +36,12 @@ mprb_mul(mprb_t z, const mprb_t x, const mprb_t y) shift = _mpr_muld_using_mpfr(zptr, zsize, xptr, xsize, yptr, ysize); /* add exponents */ - exp = x->mid.exp + y->mid.exp - shift; + z->mid.exp = x->mid.exp + y->mid.exp - shift; /* rounding error */ - ufloat_add_2exp(a, a, exp - z->mid.size * FLINT_BITS); + ufloat_add_2exp(a, a, mprb_ulp_exp(z)); z->rad = *a; z->mid.size = zsize; - z->mid.exp = exp; z->mid.sign = x->mid.sign ^ y->mid.sign; } diff --git a/mprb/test/t-log.c b/mprb/test/t-log.c new file mode 100644 index 00000000..a0d98a1d --- /dev/null +++ b/mprb/test/t-log.c @@ -0,0 +1,65 @@ +#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("log...."); + fflush(stdout); + + flint_randinit(state); + _flint_rand_init_gmp(state); + + for (iter = 0; iter < 100000; iter++) + { + mprb_t x, y; + mpfr_t t, u; + + mprb_init(x, 1 + n_randint(state, 300)); + mprb_init(y, 1 + n_randint(state, 300)); + + do { mprb_randtest(x, state, -1000, 1000); } + while (x->mid.size == 0 || mprb_contains_zero(x) + || x->mid.sign == MPRB_SIGN_MINUS); + + mpfr_init2(t, 1000); + mpfr_init2(u, 1000); + + if (n_randint(state, 2)) + mprb_get_interval_mpfr(t, u, x); + else + mprb_get_interval_mpfr(u, t, x); + + mprb_log_using_mpfr(y, x); + mpfr_log(u, t, MPFR_RNDN); + + if (!mprb_contains_mpfr(y, u)) + { + printf("FAIL!\n"); + mprb_debug(x); + mprb_debug(y); + + mpfr_debug(t); + mpfr_debug(u); + + abort(); + } + + mpfr_clear(t); + mpfr_clear(u); + + mprb_clear(x); + mprb_clear(y); + } + + printf("PASS\n"); + flint_randclear(state); + return EXIT_SUCCESS; +} diff --git a/ufloat.h b/ufloat.h index 4e2f8926..be336d3d 100644 --- a/ufloat.h +++ b/ufloat.h @@ -116,9 +116,9 @@ ufloat_zero(ufloat_t u) } static __inline__ int -ufloat_is_zero(ufloat_t u) +ufloat_is_zero(const ufloat_t u) { - return u->man != 0UL; + return u->man == 0UL; } @@ -242,9 +242,10 @@ ufloat_sub(ufloat_t z, const ufloat_t x, const ufloat_t y) } static __inline__ void -ufloat_set_2exp(ufloat_t y, const ufloat_t x, long exp) +ufloat_set_2exp(ufloat_t y, long exp) { - + y->man = 1UL << (UFLOAT_PREC - 1); + y->exp = exp + 1; } static __inline__ void