tentative code for log

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);

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);
 }
-*/

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));
 }

mprb/add_rad_ui_2exp.c

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

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;
+}

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);

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;
+}

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);
+}

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;
 }

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;
+}

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