summaryrefslogtreecommitdiff
path: root/patches/mpfr/3.1.3/150-divhigh-basecase.patch
blob: ef8d8912446faa53dd8316d35f22a503a260a3cf (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
diff -Naurd mpfr-3.1.3-a/PATCHES mpfr-3.1.3-b/PATCHES
--- mpfr-3.1.3-a/PATCHES	2015-10-29 13:47:46.735901185 +0000
+++ mpfr-3.1.3-b/PATCHES	2015-10-29 13:47:46.763900609 +0000
@@ -0,0 +1 @@
+divhigh-basecase
diff -Naurd mpfr-3.1.3-a/VERSION mpfr-3.1.3-b/VERSION
--- mpfr-3.1.3-a/VERSION	2015-07-17 08:58:21.118986898 +0000
+++ mpfr-3.1.3-b/VERSION	2015-10-29 13:47:46.763900609 +0000
@@ -1 +1 @@
-3.1.3-p4
+3.1.3-p5
diff -Naurd mpfr-3.1.3-a/src/mpfr.h mpfr-3.1.3-b/src/mpfr.h
--- mpfr-3.1.3-a/src/mpfr.h	2015-07-17 08:58:21.114986979 +0000
+++ mpfr-3.1.3-b/src/mpfr.h	2015-10-29 13:47:46.759900692 +0000
@@ -27,7 +27,7 @@
 #define MPFR_VERSION_MAJOR 3
 #define MPFR_VERSION_MINOR 1
 #define MPFR_VERSION_PATCHLEVEL 3
-#define MPFR_VERSION_STRING "3.1.3-p4"
+#define MPFR_VERSION_STRING "3.1.3-p5"
 
 /* Macros dealing with MPFR VERSION */
 #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c))
diff -Naurd mpfr-3.1.3-a/src/mulders.c mpfr-3.1.3-b/src/mulders.c
--- mpfr-3.1.3-a/src/mulders.c	2015-06-19 19:55:10.000000000 +0000
+++ mpfr-3.1.3-b/src/mulders.c	2015-10-29 13:47:46.751900855 +0000
@@ -236,9 +236,10 @@
          that in addition to the limb np[n-1] to reduce, we have at least 2
          extra limbs, thus accessing np[n-3] is valid. */
 
-      /* warning: we can have np[n-1]=d1 and np[n-2]=d0, but since {np,n} < D,
-         the largest possible partial quotient is B-1 */
-      if (MPFR_UNLIKELY(np[n - 1] == d1 && np[n - 2] == d0))
+      /* Warning: we can have np[n-1]>d1 or (np[n-1]=d1 and np[n-2]>=d0) here,
+         since we truncate the divisor at each step, but since {np,n} < D
+         originally, the largest possible partial quotient is B-1. */
+      if (MPFR_UNLIKELY(np[n-1] > d1 || (np[n-1] == d1 && np[n-2] >= d0)))
         q2 = ~ (mp_limb_t) 0;
       else
         udiv_qr_3by2 (q2, q1, q0, np[n - 1], np[n - 2], np[n - 3],
diff -Naurd mpfr-3.1.3-a/src/version.c mpfr-3.1.3-b/src/version.c
--- mpfr-3.1.3-a/src/version.c	2015-07-17 08:58:21.118986898 +0000
+++ mpfr-3.1.3-b/src/version.c	2015-10-29 13:47:46.763900609 +0000
@@ -25,5 +25,5 @@
 const char *
 mpfr_get_version (void)
 {
-  return "3.1.3-p4";
+  return "3.1.3-p5";
 }
diff -Naurd mpfr-3.1.3-a/tests/tdiv.c mpfr-3.1.3-b/tests/tdiv.c
--- mpfr-3.1.3-a/tests/tdiv.c	2015-06-19 19:55:10.000000000 +0000
+++ mpfr-3.1.3-b/tests/tdiv.c	2015-10-29 13:47:46.751900855 +0000
@@ -1099,6 +1099,69 @@
   mpfr_set_emax (old_emax);
 }
 
+/* Bug in mpfr_divhigh_n_basecase when all limbs of q (except the most
+   significant one) are B-1 where B=2^GMP_NUMB_BITS. Since we truncate
+   the divisor at each step, it might happen at some point that
+   (np[n-1],np[n-2]) > (d1,d0), and not only the equality.
+   Reported by Ricky Farr
+   <https://sympa.inria.fr/sympa/arc/mpfr/2015-10/msg00023.html>
+   To get a failure, a MPFR_DIVHIGH_TAB entry below the MPFR_DIV_THRESHOLD
+   limit must have a value 0. With most mparam.h files, this cannot occur. */
+static void
+test_20151023 (void)
+{
+  mpfr_prec_t p;
+  mpfr_t n, d, q, q0;
+  int inex, i;
+
+  for (p = GMP_NUMB_BITS; p <= 2000; p++)
+    {
+      mpfr_init2 (n, 2*p);
+      mpfr_init2 (d, p);
+      mpfr_init2 (q, p);
+      mpfr_init2 (q0, GMP_NUMB_BITS);
+
+      /* generate a random divisor of p bits */
+      mpfr_urandomb (d, RANDS);
+      /* generate a random quotient of GMP_NUMB_BITS bits */
+      mpfr_urandomb (q0, RANDS);
+      /* zero-pad the quotient to p bits */
+      inex = mpfr_prec_round (q0, p, MPFR_RNDN);
+      MPFR_ASSERTN(inex == 0);
+
+      for (i = 0; i < 3; i++)
+        {
+          /* i=0: try with the original quotient xxx000...000
+             i=1: try with the original quotient minus one ulp
+             i=2: try with the original quotient plus one ulp */
+          if (i == 1)
+            mpfr_nextbelow (q0);
+          else if (i == 2)
+            {
+              mpfr_nextabove (q0);
+              mpfr_nextabove (q0);
+            }
+
+          inex = mpfr_mul (n, d, q0, MPFR_RNDN);
+          MPFR_ASSERTN(inex == 0);
+          mpfr_nextabove (n);
+          mpfr_div (q, n, d, MPFR_RNDN);
+          MPFR_ASSERTN(mpfr_cmp (q, q0) == 0);
+
+          inex = mpfr_mul (n, d, q0, MPFR_RNDN);
+          MPFR_ASSERTN(inex == 0);
+          mpfr_nextbelow (n);
+          mpfr_div (q, n, d, MPFR_RNDN);
+          MPFR_ASSERTN(mpfr_cmp (q, q0) == 0);
+        }
+
+      mpfr_clear (n);
+      mpfr_clear (d);
+      mpfr_clear (q);
+      mpfr_clear (q0);
+    }
+}
+
 #define TEST_FUNCTION test_div
 #define TWO_ARGS
 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
@@ -1219,6 +1282,7 @@
   consistency ();
   test_20070603 ();
   test_20070628 ();
+  test_20151023 ();
   test_generic (2, 800, 50);
   test_extreme ();