make exp_series work with an arbitrary constant term; some improvements

arb.h

@@ -32,6 +32,8 @@
 #include "flint.h"
 #include "fmpz.h"
 #include "fmpq.h"
+#include "ufloat.h"
+
 
 typedef struct
 {
@@ -138,7 +140,10 @@ void arb_sqrt_ui(arb_t x, ulong n);
 void arb_log_ui(arb_t x, ulong n);
 void arb_log(arb_t y, const arb_t x);
 
+void arb_exp(arb_t y, const arb_t x);
+
 void arb_add_error_2exp(arb_t x, long c);
+void _arb_rad_add_ufloat(arb_t y, const ufloat_t err);
 
 void arb_const_pi_chudnovsky(arb_t x);
 void arb_const_euler_brent_mcmillan(arb_t x);
@@ -165,5 +170,4 @@ void _fmpz_addmul_abs(fmpz_t, const fmpz_t, const fmpz_t);
 
 void _fmpz_addmul_abs_ui(fmpz_t, const fmpz_t, ulong x);
 
-
 #endif

arb/exp.c

@@ -0,0 +1,82 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb.h"
+#include "ufloat.h"
+
+/* TODO: check for overflow (both mid and rad) */
+
+
+/*
+exp((a+b)*2^r) - exp(a*2^r) = exp(a*2^r) * (exp(b*2^r)-1)
+*/
+void
+arb_exp_error(ufloat_t err, const fmpz_t mid, const fmpz_t rad, const fmpz_t exp)
+{
+    mpfr_t a, b;
+
+    mpfr_init2(a, 32);
+    mpfr_init2(b, 32);
+
+    _arb_get_mpfr(a, mid, exp, MPFR_RNDU);
+    _arb_get_mpfr(b, rad, exp, MPFR_RNDU);
+
+    mpfr_exp(a, a, MPFR_RNDU);
+    mpfr_expm1(b, b, MPFR_RNDU);
+    mpfr_mul(a, a, b, MPFR_RNDU);
+
+    ufloat_set_mpfr(err, a);
+
+    mpfr_clear(a);
+    mpfr_clear(b);
+}
+
+void
+arb_exp(arb_t y, const arb_t x)
+{
+    long prec;
+    mpfr_t t, u;
+    int input_approx;
+    ufloat_t err;
+
+    prec = FLINT_MAX(2, arb_prec(y));
+
+    mpfr_init2(t, 2 + fmpz_bits(arb_midref(x)));
+    mpfr_init2(u, prec);
+    arb_get_mpfr(t, x, MPFR_RNDN);  /* exact */
+
+    input_approx = !fmpz_is_zero(arb_radref(x));
+    if (input_approx)
+        arb_exp_error(err, arb_midref(x), arb_radref(x), arb_expref(x));
+
+    mpfr_exp(u, t, MPFR_RNDN);
+    arb_set_mpfr(y, u, 1);
+
+    if (input_approx)
+        _arb_rad_add_ufloat(y, err);
+
+    mpfr_clear(t);
+    mpfr_clear(u);
+}

arb/log.c

@@ -51,17 +51,6 @@ arb_log_error(ufloat_t err, const fmpz_t mid, const fmpz_t rad)
     fmpz_clear(d);
 }
 
-static void
-_arb_rad_add_ufloat(arb_t y, ufloat_t err)
-{
-    ufloat_t w;
-
-    err->exp -= *arb_expref(y);
-    ufloat_set_fmpz(w, arb_radref(y));
-    ufloat_add(w, w, err);
-    ufloat_get_fmpz(arb_radref(y), w);
-}
-
 void
 arb_log(arb_t y, const arb_t x)
 {
@@ -84,6 +73,7 @@ arb_log(arb_t y, const arb_t x)
     arb_get_mpfr(t, x, MPFR_RNDN);  /* exact */
 
     input_approx = !fmpz_is_zero(arb_radref(x));
+
     if (input_approx)
         arb_log_error(err, arb_midref(x), arb_radref(x));
 

arb/rad_add_ufloat.c

@@ -0,0 +1,62 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb.h"
+
+void
+_arb_rad_add_ufloat(arb_t y, const ufloat_t err)
+{
+    ufloat_t t, w;
+
+    if (fmpz_is_zero(arb_radref(y)))
+    {
+        long yexp, eexp;
+
+        yexp = *arb_expref(y);
+        eexp = err->exp;
+
+        if (yexp >= eexp)
+        {
+            fmpz_mul_2exp(arb_midref(y), arb_midref(y), yexp - eexp);
+            fmpz_set_ui(arb_radref(y), err->man);
+        }
+        else
+        {
+            fmpz_tdiv_q_2exp(arb_midref(y), arb_midref(y), eexp - yexp);
+            fmpz_set_ui(arb_radref(y), err->man + 1);
+        }
+
+        fmpz_set_si(arb_expref(y), eexp);
+    }
+    else
+    {
+        t->man = err->man;
+        t->exp = err->exp - *arb_expref(y);
+
+        ufloat_set_fmpz(w, arb_radref(y));
+        ufloat_add(w, w, t);
+        ufloat_get_fmpz(arb_radref(y), w);
+    }
+}

arb/test/t-exp.c

@@ -0,0 +1,88 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("exp....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        arb_t r, s;
+        fmpq_t q;
+        mpfr_t x, y;
+        long wp;
+
+        arb_init(r, 1 + n_randint(state, 600));
+        arb_init(s, 1 + n_randint(state, 600));
+
+        wp = FLINT_MAX(arb_prec(r), arb_prec(s)) + 100;
+
+        mpfr_init2(x, wp);
+        mpfr_init2(y, wp);
+
+        fmpq_init(q);
+
+        do {
+            arb_randtest(r, state, 5);
+        } while (FLINT_MAX(fmpz_bits(arb_midref(r)),
+            fmpz_bits(arb_radref(r))) + *arb_expref(r) > 13);
+
+        arb_get_rand_fmpq(q, state, r);
+        fmpq_get_mpfr(x, q, MPFR_RNDN);
+
+        arb_exp(s, r);
+        mpfr_exp(y, x, MPFR_RNDN);
+
+        if (!arb_contains_mpfr(s, y))
+        {
+            printf("FAIL: containment\n\n");
+            printf("r = "); arb_debug(r); printf("\n\n");
+            printf("s = "); arb_debug(s); printf("\n\n");
+            abort();
+        }
+
+        arb_clear(r);
+        arb_clear(s);
+
+        fmpq_clear(q);
+
+        mpfr_clear(x);
+        mpfr_clear(y);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    mpfr_free_cache();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}