yann@1250: Original patch from: mpf_eq.diff
yann@1250: 
yann@1250: -= BEGIN original header =-
yann@1250: -= END original header =-
yann@1250: 
yann@1250: diff -durN gmp-4.2.4.orig/doc/gmp.texi gmp-4.2.4/doc/gmp.texi
yann@1250: --- gmp-4.2.4.orig/doc/gmp.texi	2008-09-18 17:36:14.000000000 +0200
yann@1250: +++ gmp-4.2.4/doc/gmp.texi	2009-03-08 18:36:16.000000000 +0100
yann@1250: @@ -4849,9 +4849,12 @@
yann@1250:  equal, zero otherwise.  I.e., test if @var{op1} and @var{op2} are approximately
yann@1250:  equal.
yann@1250:  
yann@1250: -Caution: Currently only whole limbs are compared, and only in an exact
yann@1250: -fashion.  In the future values like 1000 and 0111 may be considered the same
yann@1250: -to 3 bits (on the basis that their difference is that small).
yann@1250: +Caution 1: All version of GMP up to version 4.2.4 compared just whole limbs,
yann@1250: +meaning sometimes more than @var{op3} bits, sometimes fewer.
yann@1250: +
yann@1250: +Caution 2: This function will consider XXX11...111 and XX100...000 different,
yann@1250: +even if ... is replaced by a semi-infinite number of bits.  Such numbers are
yann@1250: +really just one ulp off, and should be considered equal.
yann@1250:  @end deftypefun
yann@1250:  
yann@1250:  @deftypefun void mpf_reldiff (mpf_t @var{rop}, mpf_t @var{op1}, mpf_t @var{op2})
yann@1250: diff -durN gmp-4.2.4.orig/mpf/eq.c gmp-4.2.4/mpf/eq.c
yann@1250: --- gmp-4.2.4.orig/mpf/eq.c	2007-08-30 20:31:40.000000000 +0200
yann@1250: +++ gmp-4.2.4/mpf/eq.c	2009-03-08 18:36:16.000000000 +0100
yann@1250: @@ -1,6 +1,6 @@
yann@1250:  /* mpf_eq -- Compare two floats up to a specified bit #.
yann@1250:  
yann@1250: -Copyright 1993, 1995, 1996, 2001, 2002 Free Software Foundation, Inc.
yann@1250: +Copyright 1993, 1995, 1996, 2001, 2002, 2008 Free Software Foundation, Inc.
yann@1250:  
yann@1250:  This file is part of the GNU MP Library.
yann@1250:  
yann@1250: @@ -19,6 +19,7 @@
yann@1250:  
yann@1250:  #include "gmp.h"
yann@1250:  #include "gmp-impl.h"
yann@1250: +#include "longlong.h"
yann@1250:  
yann@1250:  int
yann@1250:  mpf_eq (mpf_srcptr u, mpf_srcptr v, unsigned long int n_bits)
yann@1250: @@ -26,6 +27,8 @@
yann@1250:    mp_srcptr up, vp;
yann@1250:    mp_size_t usize, vsize, size, i;
yann@1250:    mp_exp_t uexp, vexp;
yann@1250: +  mp_limb_t diff;
yann@1250: +  int cnt;
yann@1250:  
yann@1250:    uexp = u->_mp_exp;
yann@1250:    vexp = v->_mp_exp;
yann@1250: @@ -53,10 +56,8 @@
yann@1250:    /* U and V have the same sign and are both non-zero.  */
yann@1250:  
yann@1250:    /* 2. Are the exponents different?  */
yann@1250: -  if (uexp > vexp)
yann@1250: -    return 0;			/* ??? handle (uexp = vexp + 1)   */
yann@1250: -  if (vexp > uexp)
yann@1250: -    return 0;			/* ??? handle (vexp = uexp + 1)   */
yann@1250: +  if (uexp != vexp)
yann@1250: +    return 0;
yann@1250:  
yann@1250:    usize = ABS (usize);
yann@1250:    vsize = ABS (vsize);
yann@1250: @@ -93,17 +94,26 @@
yann@1250:        size = usize;
yann@1250:      }
yann@1250:  
yann@1250: -  if (size > (n_bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS)
yann@1250: -    size = (n_bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
yann@1250: +  up += usize;			/* point just above most significant limb */
yann@1250: +  vp += vsize;			/* point just above most significant limb */
yann@1250:  
yann@1250: -  up += usize - size;
yann@1250: -  vp += vsize - size;
yann@1250: +  count_leading_zeros (cnt, up[-1]);
yann@1250: +  if ((vp[-1] >> (GMP_LIMB_BITS - 1 - cnt)) != 1)
yann@1250: +    return 0;			/* msb positions different */
yann@1250:  
yann@1250: -  for (i = size - 1; i >= 0; i--)
yann@1250: +  n_bits += cnt - GMP_NAIL_BITS;
yann@1250: +
yann@1250: +  size = MIN (size, (n_bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS);
yann@1250: +
yann@1250: +  up -= size;			/* point at least significant relevant limb */
yann@1250: +  vp -= size;			/* point at least significant relevant limb */
yann@1250: +
yann@1250: +  for (i = size - 1; i > 0; i--)
yann@1250:      {
yann@1250:        if (up[i] != vp[i])
yann@1250:  	return 0;
yann@1250:      }
yann@1250:  
yann@1250: -  return 1;
yann@1250: +  diff = (up[0] ^ vp[0]) >> GMP_NUMB_BITS - 1 - (n_bits - 1) % GMP_NUMB_BITS;
yann@1250: +  return diff == 0;
yann@1250:  }
yann@1250: diff -durN gmp-4.2.4.orig/tests/cxx/t-prec.cc gmp-4.2.4/tests/cxx/t-prec.cc
yann@1250: --- gmp-4.2.4.orig/tests/cxx/t-prec.cc	2007-09-01 12:09:03.000000000 +0200
yann@1250: +++ gmp-4.2.4/tests/cxx/t-prec.cc	2009-03-08 18:36:16.000000000 +0100
yann@1250: @@ -1,6 +1,6 @@
yann@1250:  /* Test precision of mpf_class expressions.
yann@1250:  
yann@1250: -Copyright 2001, 2002, 2003 Free Software Foundation, Inc.
yann@1250: +Copyright 2001, 2002, 2003, 2008 Free Software Foundation, Inc.
yann@1250:  
yann@1250:  This file is part of the GNU MP Library.
yann@1250:  
yann@1250: @@ -61,7 +61,7 @@
yann@1250:      g = 1 / f;
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250:        (g, "0.11111 11111 11111 11111 11111 11111 11111 11111 11111 11111"
yann@1250: -       "     11111 11111 11111 11111 11111 11", very_large_prec);
yann@1250: +       "     11111 11111 11111 11111 11111 111", very_large_prec);
yann@1250:    }
yann@1250:    {
yann@1250:      mpf_class f(15.0, large_prec);
yann@1250: @@ -69,7 +69,7 @@
yann@1250:      g = 1 / f;
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250:        (g, "0.06666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
yann@1250: -       "     66666 66666 66666 66666 66666 67", very_large_prec);
yann@1250: +       "     66666 66666 66666 66666 66666 667", very_large_prec);
yann@1250:    }
yann@1250:  
yann@1250:    // compound expressions
yann@1250: @@ -94,14 +94,14 @@
yann@1250:      i = f / g + h;
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250:        (i, "15.33333 33333 33333 33333 33333 33333 33333 33333 33333 33333"
yann@1250: -       "      33333 33333 33333 333", very_large_prec);
yann@1250: +       "      33333 33333 33333 33333 33333 3", very_large_prec);
yann@1250:    }
yann@1250:    {
yann@1250:      mpf_class f(3.0, small_prec);
yann@1250:      mpf_class g(-(1 + f) / 3, very_large_prec);
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250:        (g, "-1.33333 33333 33333 33333 33333 33333 33333 33333 33333 33333"
yann@1250: -       "      33333 33333 33333 333", very_large_prec);
yann@1250: +       "      33333 33333 33333 33333 33333 33", very_large_prec);
yann@1250:    }
yann@1250:    {
yann@1250:      mpf_class f(9.0, medium_prec);
yann@1250: @@ -117,7 +117,7 @@
yann@1250:      g = hypot(1 + 5 / f, 1.0);
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250:        (g, "1.66666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
yann@1250: -       "     66666 66666 66666 667", very_large_prec);
yann@1250: +       "     66666 66666 66666 66666 66666 67", very_large_prec);
yann@1250:    }
yann@1250:  
yann@1250:    // compound assignments
yann@1250: @@ -142,7 +142,7 @@
yann@1250:      mpf_class g(0.0, very_large_prec);
yann@1250:      g = mpf_class(1 / f);
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250: -      (g, "0.11111 11111 11111 11111 11111 11111 11111 111", medium_prec);
yann@1250: +      (g, "0.11111 11111 11111 11111 11111 11111 11111 1111", medium_prec);
yann@1250:    }
yann@1250:    {
yann@1250:      mpf_class f(15.0, large_prec);
yann@1250: @@ -150,7 +150,7 @@
yann@1250:      g = mpf_class(1 / f);
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250:        (g, "0.06666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
yann@1250: -       "     66666 667", large_prec);
yann@1250: +       "     66666 6667", large_prec);
yann@1250:    }
yann@1250:  
yann@1250:    {
yann@1250: @@ -158,7 +158,8 @@
yann@1250:      mpf_class h(0.0, very_large_prec);
yann@1250:      h = mpf_class(f / g + 1, large_prec);
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250: -      (h, "1.33333 33333 33333 33333 33333 33333 33333 33333 33333 3333",
yann@1250: +      (h, "1.33333 33333 33333 33333 33333 33333 33333 33333 33333 33333"
yann@1250: +       "     33333 333",
yann@1250:         large_prec);
yann@1250:    }
yann@1250:  
yann@1250: @@ -170,7 +171,7 @@
yann@1250:      g = f - q;
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250:        (g, "2.66666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
yann@1250: -       "     66666 66666 66666 667", very_large_prec);
yann@1250: +       "     66666 66666 66666 66666 66666 67", very_large_prec);
yann@1250:    }
yann@1250:  
yann@1250:    {
yann@1250: @@ -179,7 +180,8 @@
yann@1250:      mpf_class g(0.0, very_large_prec);
yann@1250:      g = mpf_class(f - q, large_prec);
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250: -      (g, "2.66666 66666 66666 66666 66666 66666 66666 66666 66666 6667",
yann@1250: +      (g, "2.66666 66666 66666 66666 66666 66666 66666 66666 66666 66666"
yann@1250: +       "     66666 667",
yann@1250:         large_prec);
yann@1250:    }
yann@1250:    {
yann@1250: @@ -188,7 +190,7 @@
yann@1250:      mpf_class g(0.0, very_large_prec);
yann@1250:      g = mpf_class(f - q);
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250: -      (g, "2.66666 66666 66666 66666 66666 6667", medium_prec);
yann@1250: +      (g, "2.66666 66666 66666 66666 66666 66666 66666 667", medium_prec);
yann@1250:    }
yann@1250:    {
yann@1250:      mpf_class f(15.0, large_prec);
yann@1250: @@ -196,7 +198,8 @@
yann@1250:      mpf_class g(0.0, very_large_prec);
yann@1250:      g = mpf_class(f + q);
yann@1250:      ASSERT_ALWAYS_PREC
yann@1250: -      (g, "15.33333 33333 33333 33333 33333 33333 33333 33333 33333 3333",
yann@1250: +      (g, "15.33333 33333 33333 33333 33333 33333 33333 33333 33333 33333"
yann@1250: +       "      33333 33",
yann@1250:         large_prec);
yann@1250:    }
yann@1250:  }