| [d3318e9] | 1 | Submitted By: Jim Gifford (jim at cross-lfs dot org)
|
|---|
| 2 | Date: 2009-01-03
|
|---|
| 3 | Initial Package Version: 4.2.4
|
|---|
| 4 | Origin: GMP Website
|
|---|
| 5 | Upstream Status: Fixed
|
|---|
| 6 | Description: See http://gmplib.org Website Under Curent Status
|
|---|
| 7 |
|
|---|
| 8 | diff -Naur gmp-4.2.4.orig/doc/gmp.texi gmp-4.2.4/doc/gmp.texi
|
|---|
| 9 | --- gmp-4.2.4.orig/doc/gmp.texi 2008-09-18 11:36:14.000000000 -0400
|
|---|
| 10 | +++ gmp-4.2.4/doc/gmp.texi 2009-01-03 13:48:39.498471376 -0500
|
|---|
| 11 | @@ -4849,9 +4849,12 @@
|
|---|
| 12 | equal, zero otherwise. I.e., test if @var{op1} and @var{op2} are approximately
|
|---|
| 13 | equal.
|
|---|
| 14 |
|
|---|
| 15 | -Caution: Currently only whole limbs are compared, and only in an exact
|
|---|
| 16 | -fashion. In the future values like 1000 and 0111 may be considered the same
|
|---|
| 17 | -to 3 bits (on the basis that their difference is that small).
|
|---|
| 18 | +Caution 1: All version of GMP up to version 4.2.4 compared just whole limbs,
|
|---|
| 19 | +meaning sometimes more than @var{op3} bits, sometimes fewer.
|
|---|
| 20 | +
|
|---|
| 21 | +Caution 2: This function will consider XXX11...111 and XX100...000 different,
|
|---|
| 22 | +even if ... is replaced by a semi-infinite number of bits. Such numbers are
|
|---|
| 23 | +really just one ulp off, and should be considered equal.
|
|---|
| 24 | @end deftypefun
|
|---|
| 25 |
|
|---|
| 26 | @deftypefun void mpf_reldiff (mpf_t @var{rop}, mpf_t @var{op1}, mpf_t @var{op2})
|
|---|
| 27 | diff -Naur gmp-4.2.4.orig/mpf/eq.c gmp-4.2.4/mpf/eq.c
|
|---|
| 28 | --- gmp-4.2.4.orig/mpf/eq.c 2007-08-30 14:31:40.000000000 -0400
|
|---|
| 29 | +++ gmp-4.2.4/mpf/eq.c 2009-01-03 13:48:39.498471376 -0500
|
|---|
| 30 | @@ -1,6 +1,6 @@
|
|---|
| 31 | /* mpf_eq -- Compare two floats up to a specified bit #.
|
|---|
| 32 |
|
|---|
| 33 | -Copyright 1993, 1995, 1996, 2001, 2002 Free Software Foundation, Inc.
|
|---|
| 34 | +Copyright 1993, 1995, 1996, 2001, 2002, 2008 Free Software Foundation, Inc.
|
|---|
| 35 |
|
|---|
| 36 | This file is part of the GNU MP Library.
|
|---|
| 37 |
|
|---|
| 38 | @@ -19,6 +19,7 @@
|
|---|
| 39 |
|
|---|
| 40 | #include "gmp.h"
|
|---|
| 41 | #include "gmp-impl.h"
|
|---|
| 42 | +#include "longlong.h"
|
|---|
| 43 |
|
|---|
| 44 | int
|
|---|
| 45 | mpf_eq (mpf_srcptr u, mpf_srcptr v, unsigned long int n_bits)
|
|---|
| 46 | @@ -26,6 +27,8 @@
|
|---|
| 47 | mp_srcptr up, vp;
|
|---|
| 48 | mp_size_t usize, vsize, size, i;
|
|---|
| 49 | mp_exp_t uexp, vexp;
|
|---|
| 50 | + mp_limb_t diff;
|
|---|
| 51 | + int cnt;
|
|---|
| 52 |
|
|---|
| 53 | uexp = u->_mp_exp;
|
|---|
| 54 | vexp = v->_mp_exp;
|
|---|
| 55 | @@ -53,10 +56,8 @@
|
|---|
| 56 | /* U and V have the same sign and are both non-zero. */
|
|---|
| 57 |
|
|---|
| 58 | /* 2. Are the exponents different? */
|
|---|
| 59 | - if (uexp > vexp)
|
|---|
| 60 | - return 0; /* ??? handle (uexp = vexp + 1) */
|
|---|
| 61 | - if (vexp > uexp)
|
|---|
| 62 | - return 0; /* ??? handle (vexp = uexp + 1) */
|
|---|
| 63 | + if (uexp != vexp)
|
|---|
| 64 | + return 0;
|
|---|
| 65 |
|
|---|
| 66 | usize = ABS (usize);
|
|---|
| 67 | vsize = ABS (vsize);
|
|---|
| 68 | @@ -93,17 +94,26 @@
|
|---|
| 69 | size = usize;
|
|---|
| 70 | }
|
|---|
| 71 |
|
|---|
| 72 | - if (size > (n_bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS)
|
|---|
| 73 | - size = (n_bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
|
|---|
| 74 | + up += usize; /* point just above most significant limb */
|
|---|
| 75 | + vp += vsize; /* point just above most significant limb */
|
|---|
| 76 |
|
|---|
| 77 | - up += usize - size;
|
|---|
| 78 | - vp += vsize - size;
|
|---|
| 79 | + count_leading_zeros (cnt, up[-1]);
|
|---|
| 80 | + if ((vp[-1] >> (GMP_LIMB_BITS - 1 - cnt)) != 1)
|
|---|
| 81 | + return 0; /* msb positions different */
|
|---|
| 82 |
|
|---|
| 83 | - for (i = size - 1; i >= 0; i--)
|
|---|
| 84 | + n_bits += cnt - GMP_NAIL_BITS;
|
|---|
| 85 | +
|
|---|
| 86 | + size = MIN (size, (n_bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS);
|
|---|
| 87 | +
|
|---|
| 88 | + up -= size; /* point at least significant relevant limb */
|
|---|
| 89 | + vp -= size; /* point at least significant relevant limb */
|
|---|
| 90 | +
|
|---|
| 91 | + for (i = size - 1; i > 0; i--)
|
|---|
| 92 | {
|
|---|
| 93 | if (up[i] != vp[i])
|
|---|
| 94 | return 0;
|
|---|
| 95 | }
|
|---|
| 96 |
|
|---|
| 97 | - return 1;
|
|---|
| 98 | + diff = (up[0] ^ vp[0]) >> GMP_NUMB_BITS - 1 - (n_bits - 1) % GMP_NUMB_BITS;
|
|---|
| 99 | + return diff == 0;
|
|---|
| 100 | }
|
|---|
| 101 | diff -Naur gmp-4.2.4.orig/mpf/set_str.c gmp-4.2.4/mpf/set_str.c
|
|---|
| 102 | --- gmp-4.2.4.orig/mpf/set_str.c 2008-08-25 10:11:37.000000000 -0400
|
|---|
| 103 | +++ gmp-4.2.4/mpf/set_str.c 2009-01-03 13:48:18.493274358 -0500
|
|---|
| 104 | @@ -137,7 +137,12 @@
|
|---|
| 105 | c = (unsigned char) *++str;
|
|---|
| 106 | }
|
|---|
| 107 |
|
|---|
| 108 | + /* Default base to decimal. */
|
|---|
| 109 | + if (base == 0)
|
|---|
| 110 | + base = 10;
|
|---|
| 111 | +
|
|---|
| 112 | exp_base = base;
|
|---|
| 113 | +
|
|---|
| 114 | if (base < 0)
|
|---|
| 115 | {
|
|---|
| 116 | exp_base = 10;
|
|---|
| 117 | @@ -165,10 +170,6 @@
|
|---|
| 118 | return -1;
|
|---|
| 119 | }
|
|---|
| 120 |
|
|---|
| 121 | - /* Default base to decimal. */
|
|---|
| 122 | - if (base == 0)
|
|---|
| 123 | - base = 10;
|
|---|
| 124 | -
|
|---|
| 125 | /* Locate exponent part of the input. Look from the right of the string,
|
|---|
| 126 | since the exponent is usually a lot shorter than the mantissa. */
|
|---|
| 127 | expptr = NULL;
|
|---|
| 128 | diff -Naur gmp-4.2.4.orig/mpz/perfpow.c gmp-4.2.4/mpz/perfpow.c
|
|---|
| 129 | --- gmp-4.2.4.orig/mpz/perfpow.c 2007-08-30 14:31:41.000000000 -0400
|
|---|
| 130 | +++ gmp-4.2.4/mpz/perfpow.c 2009-01-03 13:47:51.611742467 -0500
|
|---|
| 131 | @@ -1,7 +1,7 @@
|
|---|
| 132 | /* mpz_perfect_power_p(arg) -- Return non-zero if ARG is a perfect power,
|
|---|
| 133 | zero otherwise.
|
|---|
| 134 |
|
|---|
| 135 | -Copyright 1998, 1999, 2000, 2001, 2005 Free Software Foundation, Inc.
|
|---|
| 136 | +Copyright 1998, 1999, 2000, 2001, 2005, 2008 Free Software Foundation, Inc.
|
|---|
| 137 |
|
|---|
| 138 | This file is part of the GNU MP Library.
|
|---|
| 139 |
|
|---|
| 140 | @@ -59,6 +59,8 @@
|
|---|
| 141 | #define SMALLEST_OMITTED_PRIME 1009
|
|---|
| 142 |
|
|---|
| 143 |
|
|---|
| 144 | +#define POW2P(a) (((a) & ((a) - 1)) == 0)
|
|---|
| 145 | +
|
|---|
| 146 | int
|
|---|
| 147 | mpz_perfect_power_p (mpz_srcptr u)
|
|---|
| 148 | {
|
|---|
| 149 | @@ -72,16 +74,13 @@
|
|---|
| 150 | mp_size_t usize = SIZ (u);
|
|---|
| 151 | TMP_DECL;
|
|---|
| 152 |
|
|---|
| 153 | - if (usize == 0)
|
|---|
| 154 | - return 1; /* consider 0 a perfect power */
|
|---|
| 155 | + if (mpz_cmpabs_ui (u, 1) <= 0)
|
|---|
| 156 | + return 1; /* -1, 0, and +1 are perfect powers */
|
|---|
| 157 |
|
|---|
| 158 | n2 = mpz_scan1 (u, 0);
|
|---|
| 159 | if (n2 == 1)
|
|---|
| 160 | return 0; /* 2 divides exactly once. */
|
|---|
| 161 |
|
|---|
| 162 | - if (n2 != 0 && (n2 & 1) == 0 && usize < 0)
|
|---|
| 163 | - return 0; /* 2 has even multiplicity with negative U */
|
|---|
| 164 | -
|
|---|
| 165 | TMP_MARK;
|
|---|
| 166 |
|
|---|
| 167 | uns = ABS (usize) - n2 / BITS_PER_MP_LIMB;
|
|---|
| 168 | @@ -89,6 +88,14 @@
|
|---|
| 169 | MPZ_TMP_INIT (u2, uns);
|
|---|
| 170 |
|
|---|
| 171 | mpz_tdiv_q_2exp (u2, u, n2);
|
|---|
| 172 | + mpz_abs (u2, u2);
|
|---|
| 173 | +
|
|---|
| 174 | + if (mpz_cmp_ui (u2, 1) == 0)
|
|---|
| 175 | + {
|
|---|
| 176 | + TMP_FREE;
|
|---|
| 177 | + /* factoring completed; consistent power */
|
|---|
| 178 | + return ! (usize < 0 && POW2P(n2));
|
|---|
| 179 | + }
|
|---|
| 180 |
|
|---|
| 181 | if (isprime (n2))
|
|---|
| 182 | goto n2prime;
|
|---|
| 183 | @@ -97,6 +104,9 @@
|
|---|
| 184 | {
|
|---|
| 185 | prime = primes[i];
|
|---|
| 186 |
|
|---|
| 187 | + if (mpz_cmp_ui (u2, prime) < 0)
|
|---|
| 188 | + break;
|
|---|
| 189 | +
|
|---|
| 190 | if (mpz_divisible_ui_p (u2, prime)) /* divisible by this prime? */
|
|---|
| 191 | {
|
|---|
| 192 | rem = mpz_tdiv_q_ui (q, u2, prime * prime);
|
|---|
| 193 | @@ -115,12 +125,6 @@
|
|---|
| 194 | n++;
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 | - if ((n & 1) == 0 && usize < 0)
|
|---|
| 198 | - {
|
|---|
| 199 | - TMP_FREE;
|
|---|
| 200 | - return 0; /* even multiplicity with negative U, reject */
|
|---|
| 201 | - }
|
|---|
| 202 | -
|
|---|
| 203 | n2 = gcd (n2, n);
|
|---|
| 204 | if (n2 == 1)
|
|---|
| 205 | {
|
|---|
| 206 | @@ -128,10 +132,11 @@
|
|---|
| 207 | return 0; /* we have multiplicity 1 of some factor */
|
|---|
| 208 | }
|
|---|
| 209 |
|
|---|
| 210 | - if (mpz_cmpabs_ui (u2, 1) == 0)
|
|---|
| 211 | + if (mpz_cmp_ui (u2, 1) == 0)
|
|---|
| 212 | {
|
|---|
| 213 | TMP_FREE;
|
|---|
| 214 | - return 1; /* factoring completed; consistent power */
|
|---|
| 215 | + /* factoring completed; consistent power */
|
|---|
| 216 | + return ! (usize < 0 && POW2P(n2));
|
|---|
| 217 | }
|
|---|
| 218 |
|
|---|
| 219 | /* As soon as n2 becomes a prime number, stop factoring.
|
|---|
| 220 | @@ -169,6 +174,10 @@
|
|---|
| 221 | else
|
|---|
| 222 | {
|
|---|
| 223 | unsigned long int nth;
|
|---|
| 224 | +
|
|---|
| 225 | + if (usize < 0 && POW2P(n2))
|
|---|
| 226 | + return 0;
|
|---|
| 227 | +
|
|---|
| 228 | /* We found some factors above. We just need to consider values of n
|
|---|
| 229 | that divides n2. */
|
|---|
| 230 | for (nth = 2; nth <= n2; nth++)
|
|---|
| 231 | @@ -184,8 +193,11 @@
|
|---|
| 232 | exact = mpz_root (q, u2, nth);
|
|---|
| 233 | if (exact)
|
|---|
| 234 | {
|
|---|
| 235 | - TMP_FREE;
|
|---|
| 236 | - return 1;
|
|---|
| 237 | + if (! (usize < 0 && POW2P(nth)))
|
|---|
| 238 | + {
|
|---|
| 239 | + TMP_FREE;
|
|---|
| 240 | + return 1;
|
|---|
| 241 | + }
|
|---|
| 242 | }
|
|---|
| 243 | if (mpz_cmp_ui (q, SMALLEST_OMITTED_PRIME) < 0)
|
|---|
| 244 | {
|
|---|
| 245 | @@ -199,6 +211,9 @@
|
|---|
| 246 | }
|
|---|
| 247 |
|
|---|
| 248 | n2prime:
|
|---|
| 249 | + if (usize < 0 && POW2P(n2))
|
|---|
| 250 | + return 0;
|
|---|
| 251 | +
|
|---|
| 252 | exact = mpz_root (NULL, u2, n2);
|
|---|
| 253 | TMP_FREE;
|
|---|
| 254 | return exact;
|
|---|
| 255 | diff -Naur gmp-4.2.4.orig/tests/cxx/t-prec.cc gmp-4.2.4/tests/cxx/t-prec.cc
|
|---|
| 256 | --- gmp-4.2.4.orig/tests/cxx/t-prec.cc 2007-09-01 06:09:03.000000000 -0400
|
|---|
| 257 | +++ gmp-4.2.4/tests/cxx/t-prec.cc 2009-01-03 13:48:39.498471376 -0500
|
|---|
| 258 | @@ -1,6 +1,6 @@
|
|---|
| 259 | /* Test precision of mpf_class expressions.
|
|---|
| 260 |
|
|---|
| 261 | -Copyright 2001, 2002, 2003 Free Software Foundation, Inc.
|
|---|
| 262 | +Copyright 2001, 2002, 2003, 2008 Free Software Foundation, Inc.
|
|---|
| 263 |
|
|---|
| 264 | This file is part of the GNU MP Library.
|
|---|
| 265 |
|
|---|
| 266 | @@ -61,7 +61,7 @@
|
|---|
| 267 | g = 1 / f;
|
|---|
| 268 | ASSERT_ALWAYS_PREC
|
|---|
| 269 | (g, "0.11111 11111 11111 11111 11111 11111 11111 11111 11111 11111"
|
|---|
| 270 | - " 11111 11111 11111 11111 11111 11", very_large_prec);
|
|---|
| 271 | + " 11111 11111 11111 11111 11111 111", very_large_prec);
|
|---|
| 272 | }
|
|---|
| 273 | {
|
|---|
| 274 | mpf_class f(15.0, large_prec);
|
|---|
| 275 | @@ -69,7 +69,7 @@
|
|---|
| 276 | g = 1 / f;
|
|---|
| 277 | ASSERT_ALWAYS_PREC
|
|---|
| 278 | (g, "0.06666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
|
|---|
| 279 | - " 66666 66666 66666 66666 66666 67", very_large_prec);
|
|---|
| 280 | + " 66666 66666 66666 66666 66666 667", very_large_prec);
|
|---|
| 281 | }
|
|---|
| 282 |
|
|---|
| 283 | // compound expressions
|
|---|
| 284 | @@ -94,14 +94,14 @@
|
|---|
| 285 | i = f / g + h;
|
|---|
| 286 | ASSERT_ALWAYS_PREC
|
|---|
| 287 | (i, "15.33333 33333 33333 33333 33333 33333 33333 33333 33333 33333"
|
|---|
| 288 | - " 33333 33333 33333 333", very_large_prec);
|
|---|
| 289 | + " 33333 33333 33333 33333 33333 3", very_large_prec);
|
|---|
| 290 | }
|
|---|
| 291 | {
|
|---|
| 292 | mpf_class f(3.0, small_prec);
|
|---|
| 293 | mpf_class g(-(1 + f) / 3, very_large_prec);
|
|---|
| 294 | ASSERT_ALWAYS_PREC
|
|---|
| 295 | (g, "-1.33333 33333 33333 33333 33333 33333 33333 33333 33333 33333"
|
|---|
| 296 | - " 33333 33333 33333 333", very_large_prec);
|
|---|
| 297 | + " 33333 33333 33333 33333 33333 33", very_large_prec);
|
|---|
| 298 | }
|
|---|
| 299 | {
|
|---|
| 300 | mpf_class f(9.0, medium_prec);
|
|---|
| 301 | @@ -117,7 +117,7 @@
|
|---|
| 302 | g = hypot(1 + 5 / f, 1.0);
|
|---|
| 303 | ASSERT_ALWAYS_PREC
|
|---|
| 304 | (g, "1.66666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
|
|---|
| 305 | - " 66666 66666 66666 667", very_large_prec);
|
|---|
| 306 | + " 66666 66666 66666 66666 66666 67", very_large_prec);
|
|---|
| 307 | }
|
|---|
| 308 |
|
|---|
| 309 | // compound assignments
|
|---|
| 310 | @@ -142,7 +142,7 @@
|
|---|
| 311 | mpf_class g(0.0, very_large_prec);
|
|---|
| 312 | g = mpf_class(1 / f);
|
|---|
| 313 | ASSERT_ALWAYS_PREC
|
|---|
| 314 | - (g, "0.11111 11111 11111 11111 11111 11111 11111 111", medium_prec);
|
|---|
| 315 | + (g, "0.11111 11111 11111 11111 11111 11111 11111 1111", medium_prec);
|
|---|
| 316 | }
|
|---|
| 317 | {
|
|---|
| 318 | mpf_class f(15.0, large_prec);
|
|---|
| 319 | @@ -150,7 +150,7 @@
|
|---|
| 320 | g = mpf_class(1 / f);
|
|---|
| 321 | ASSERT_ALWAYS_PREC
|
|---|
| 322 | (g, "0.06666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
|
|---|
| 323 | - " 66666 667", large_prec);
|
|---|
| 324 | + " 66666 6667", large_prec);
|
|---|
| 325 | }
|
|---|
| 326 |
|
|---|
| 327 | {
|
|---|
| 328 | @@ -158,7 +158,8 @@
|
|---|
| 329 | mpf_class h(0.0, very_large_prec);
|
|---|
| 330 | h = mpf_class(f / g + 1, large_prec);
|
|---|
| 331 | ASSERT_ALWAYS_PREC
|
|---|
| 332 | - (h, "1.33333 33333 33333 33333 33333 33333 33333 33333 33333 3333",
|
|---|
| 333 | + (h, "1.33333 33333 33333 33333 33333 33333 33333 33333 33333 33333"
|
|---|
| 334 | + " 33333 333",
|
|---|
| 335 | large_prec);
|
|---|
| 336 | }
|
|---|
| 337 |
|
|---|
| 338 | @@ -170,7 +171,7 @@
|
|---|
| 339 | g = f - q;
|
|---|
| 340 | ASSERT_ALWAYS_PREC
|
|---|
| 341 | (g, "2.66666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
|
|---|
| 342 | - " 66666 66666 66666 667", very_large_prec);
|
|---|
| 343 | + " 66666 66666 66666 66666 66666 67", very_large_prec);
|
|---|
| 344 | }
|
|---|
| 345 |
|
|---|
| 346 | {
|
|---|
| 347 | @@ -179,7 +180,8 @@
|
|---|
| 348 | mpf_class g(0.0, very_large_prec);
|
|---|
| 349 | g = mpf_class(f - q, large_prec);
|
|---|
| 350 | ASSERT_ALWAYS_PREC
|
|---|
| 351 | - (g, "2.66666 66666 66666 66666 66666 66666 66666 66666 66666 6667",
|
|---|
| 352 | + (g, "2.66666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
|
|---|
| 353 | + " 66666 667",
|
|---|
| 354 | large_prec);
|
|---|
| 355 | }
|
|---|
| 356 | {
|
|---|
| 357 | @@ -188,7 +190,7 @@
|
|---|
| 358 | mpf_class g(0.0, very_large_prec);
|
|---|
| 359 | g = mpf_class(f - q);
|
|---|
| 360 | ASSERT_ALWAYS_PREC
|
|---|
| 361 | - (g, "2.66666 66666 66666 66666 66666 6667", medium_prec);
|
|---|
| 362 | + (g, "2.66666 66666 66666 66666 66666 66666 66666 667", medium_prec);
|
|---|
| 363 | }
|
|---|
| 364 | {
|
|---|
| 365 | mpf_class f(15.0, large_prec);
|
|---|
| 366 | @@ -196,7 +198,8 @@
|
|---|
| 367 | mpf_class g(0.0, very_large_prec);
|
|---|
| 368 | g = mpf_class(f + q);
|
|---|
| 369 | ASSERT_ALWAYS_PREC
|
|---|
| 370 | - (g, "15.33333 33333 33333 33333 33333 33333 33333 33333 33333 3333",
|
|---|
| 371 | + (g, "15.33333 33333 33333 33333 33333 33333 33333 33333 33333 33333"
|
|---|
| 372 | + " 33333 33",
|
|---|
| 373 | large_prec);
|
|---|
| 374 | }
|
|---|
| 375 | }
|
|---|