Merge pull request #2365 from d0sboots/set_str

Make arb_set_str round correctly almost always

Author: fredrik-johansson

src/acb_hypgeom/test/t-beta_lower.c

@@ -188,7 +188,7 @@ TEST_FUNCTION_START(acb_hypgeom_beta_lower, state)
         acb_hypgeom_beta_lower(res, a, b, z, 1, prec);
         arb_set_str(acb_realref(res2), "0.9999999999999999999999999999999999999999999999999999999999999999998474866919 +/- 2.07e-77", prec);
 
-        if (!acb_overlaps(res, res2) || acb_rel_accuracy_bits(res) < acb_rel_accuracy_bits(res2) - 20)
+        if (!acb_overlaps(res, res2) || acb_rel_accuracy_bits(res) < prec - 20)
         {
             flint_printf("FAIL: test case (4)\n\n");
             flint_printf("res = "); acb_printd(res, 100); flint_printf("\n\n");

src/arb/set_str.c

@@ -127,26 +127,42 @@ arb_set_float_str(arb_t res, const char * inp, slong prec)
     }
     else if (fmpz_is_zero(exp))
     {
-        arb_set_round_fmpz(res, man, prec);
+        int inexact;
+        inexact = arf_set_round_fmpz(arb_midref(res), man, prec, ARF_RND_NEAR);
+
+        /* nearest rounding guarantees half-ulp accuracy */
+        if (inexact)
+            arf_mag_set_ulp(arb_radref(res), arb_midref(res), prec+1);
+        else
+            mag_zero(arb_radref(res));
     }
     else
     {
+        slong wp;
         arb_t t;
         arb_init(t);
         arb_set_ui(t, 10);
         arb_set_fmpz(res, man);
+        wp = arb_bits(res);
+        wp = (wp > prec ? wp : prec) + 10;
 
         if (fmpz_sgn(exp) > 0)
         {
-            arb_pow_fmpz_binexp(t, t, exp, prec + 4);
-            arb_mul(res, res, t, prec);
+            arb_pow_fmpz_binexp(t, t, exp, wp);
+            arb_mul(res, res, t, prec + 10);
         }
         else
         {
             fmpz_neg(exp, exp);
-            arb_pow_fmpz_binexp(t, t, exp, prec + 4);
-            arb_div(res, res, t, prec);
+            arb_pow_fmpz_binexp(t, t, exp, wp);
+            arb_div(res, res, t, prec + 10);
         }
+        /* There doesn't seem to be a faster way to do this in one shot, since
+         * all arb methods round down. */
+        arf_set_round(arb_midref(t), arb_midref(res), prec, ARF_RND_NEAR);
+        arf_swap(arb_midref(t), arb_midref(res));
+        arf_sub(arb_midref(t), arb_midref(t), arb_midref(res), 30, ARF_RND_UP);
+        arb_add_error_arf(res, arb_midref(t));
 
         arb_clear(t);
     }

src/arb/test/t-set_str.c

@@ -100,6 +100,56 @@ const char * testdata_floats[] = {
     "-12345.125",
     "+12345.125",
 
+    "1.0625",
+    "1.03125",
+    "1.015625",
+    "1.0078125",
+    "1.00390625",
+    "1.001953125",
+    "1.0009765625",
+    "1.00048828125",
+    "1.000244140625",
+    "1.0001220703125",
+    "1.00006103515625",
+    "1.000030517578125",
+    "1.0000152587890625",
+    "1.00000762939453125",
+    "1.000003814697265625",
+    "1.0000019073486328125",
+    "1.00000095367431640625",
+    "1.000000476837158203125",
+    "1.0000002384185791015625",
+    "1.00000011920928955078125",
+    "1.000000059604644775390625",
+    "1.0000000298023223876953125",
+    "1.00000001490116119384765625",
+    "1.000000007450580596923828125",
+    "1.0000000037252902984619140625",
+    "1.00000000186264514923095703125",
+    "1.000000000931322574615478515625",
+    "1.0000000004656612873077392578125",
+    "1.00000000023283064365386962890625",
+    "1.000000000116415321826934814453125",
+    "1.0000000000582076609134674072265625",
+    "1.00000000002910383045673370361328125",
+    "1.000000000014551915228366851806640625",
+    "1.0000000000072759576141834259033203125",
+    "1.00000000000363797880709171295166015625",
+    "1.000000000001818989403545856475830078125",
+    "1.0000000000009094947017729282379150390625",
+    "1.00000000000045474735088646411895751953125",
+    "1.000000000000227373675443232059478759765625",
+    "1.0000000000001136868377216160297393798828125",
+    "1.00000000000005684341886080801486968994140625",
+    "1.000000000000028421709430404007434844970703125",
+    "1.0000000000000142108547152020037174224853515625",
+    "1.00000000000000710542735760100185871124267578125",
+    "1.000000000000003552713678800500929355621337890625",
+    "1.0000000000000017763568394002504646778106689453125",
+    "1.00000000000000088817841970012523233890533447265625",
+    "1.000000000000000444089209850062616169452667236328125",
+    "1.0000000000000002220446049250313080847263336181640625",
+
     "inf",
     "-inf",
     "+inf",
@@ -149,6 +199,31 @@ const char * testdata_invalid[] = {
     NULL,
 };
 
+const char * testdata_rounding[] = {
+    "1.1",
+    "-1.1",
+    "1.01",
+    "1.001",
+    "1.0001",
+    "1.00001",
+    "1.000001",
+    "1.0000001",
+    "1.00000001",
+    "1.000000001",
+    "1.0000000001",
+    "1.00000000001",
+    "1.000000000001",
+    "1.0000000000001",
+    "1.00000000000001",
+    "1.000000000000001",
+    "1.0000000000000001",
+    "1.00000000000000001",
+    "100000000000000001",
+    "1000000000000000001",
+    "10000000000000000001",
+    NULL,
+};
+
 TEST_FUNCTION_START(arb_set_str, state)
 {
     arb_t t, u, v;
@@ -291,6 +366,47 @@ TEST_FUNCTION_START(arb_set_str, state)
         }
     }
 
+    for (i = 0; testdata_rounding[i] != NULL; i++)
+    {
+        arb_const_pi(t, 53);
+
+        error = arb_set_str(t, testdata_rounding[i], 53);
+
+        /* It is not *guaranteed* that the rounding will be to half-ulp
+         * accuracy, just like it is not guaranteed that all inputs will be
+         * correctly rounded. But it is very likely for a random string, and
+         * it *is* guaranteed if the input is the result of get_str performed
+         * with a sufficient number of digits.
+         * Thus, we are justified in testing for these things here. */
+        x = strtod(testdata_rounding[i], NULL);
+        arf_set_d(arb_midref(u), x);
+        arf_mag_set_ulp(arb_radref(u), arb_midref(u), 54);
+
+        if (error != 0 || !arf_equal(arb_midref(t), arb_midref(u)) ||
+            mag_cmp(arb_radref(t), arb_radref(u)) > 0)
+        {
+            flint_printf("FAIL (valid input): %s\n", testdata_rounding[i]);
+            arb_printd(t, 19); flint_printf("\n");
+            arb_printd(u, 19); flint_printf("\n");
+            flint_abort();
+        }
+
+        error = arb_set_str(t, testdata_rounding[i], 24);
+
+        x = strtof(testdata_rounding[i], NULL);
+        arf_set_d(arb_midref(u), x);
+        arf_mag_set_ulp(arb_radref(u), arb_midref(u), 25);
+
+        if (error != 0 || !arf_equal(arb_midref(t), arb_midref(u)) ||
+            mag_cmp(arb_radref(t), arb_radref(u)) > 0)
+        {
+            flint_printf("FAIL (valid input): %s\n", testdata_rounding[i]);
+            arb_printd(t, 19); flint_printf("\n");
+            arb_printd(u, 19); flint_printf("\n");
+            flint_abort();
+        }
+    }
+
     arb_clear(t);
     arb_clear(u);
     arb_clear(v);

src/python/flint_ctypes.py

@@ -4407,9 +4407,9 @@ def expm1(self):
         Exponential function minus 1.
 
             >>> RR("1e-10").expm1()
-            [1.000000000050000e-10 +/- 3.86e-26]
+            [1.000000000050000e-10 +/- 1.69e-26]
             >>> CC(RR("1e-10")).expm1()
-            [1.000000000050000e-10 +/- 3.86e-26]
+            [1.000000000050000e-10 +/- 1.69e-26]
             >>> RF("1e-10").expm1()
             1.000000000050000e-10
             >>> CF(RF("1e-10")).expm1()