contrib/gmp/mpf/eq.c

   1 /* mpf_eq -- Compare two floats up to a specified bit #.
   2
   3 Copyright 1993, 1995, 1996, 2001, 2002, 2008, 2009 Free Software Foundation,
   4 Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MP Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  20
  21 #include "gmp.h"
  22 #include "gmp-impl.h"
  23 #include "longlong.h"
  24
  25 int
  26 mpf_eq (mpf_srcptr u, mpf_srcptr v, unsigned long int n_bits)
  27 {
  28   mp_srcptr up, vp, p;
  29   mp_size_t usize, vsize, minsize, maxsize, n_limbs, i;
  30   mp_exp_t uexp, vexp;
  31   mp_limb_t diff;
  32   int cnt;
  33
  34   uexp = u->_mp_exp;
  35   vexp = v->_mp_exp;
  36
  37   usize = u->_mp_size;
  38   vsize = v->_mp_size;
  39
  40   /* 1. Are the signs different?  */
  41   if ((usize ^ vsize) >= 0)
  42     {
  43       /* U and V are both non-negative or both negative.  */
  44       if (usize == 0)
  45         return vsize == 0;
  46       if (vsize == 0)
  47         return 0;
  48
  49       /* Fall out.  */
  50     }
  51   else
  52     {
  53       /* Either U or V is negative, but not both.  */
  54       return 0;
  55     }
  56
  57   /* U and V have the same sign and are both non-zero.  */
  58
  59   /* 2. Are the exponents different?  */
  60   if (uexp != vexp)
  61     return 0;
  62
  63   usize = ABS (usize);
  64   vsize = ABS (vsize);
  65
  66   up = u->_mp_d;
  67   vp = v->_mp_d;
  68
  69   up += usize;                  /* point just above most significant limb */
  70   vp += vsize;                  /* point just above most significant limb */
  71
  72   count_leading_zeros (cnt, up[-1]);
  73   if ((vp[-1] >> (GMP_LIMB_BITS - 1 - cnt)) != 1)
  74     return 0;                   /* msb positions different */
  75
  76   n_bits += cnt - GMP_NAIL_BITS;
  77   n_limbs = (n_bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
  78
  79   usize = MIN (usize, n_limbs);
  80   vsize = MIN (vsize, n_limbs);
  81
  82 #if 0
  83   /* Ignore zeros at the low end of U and V.  */
  84   while (up[0] == 0)
  85     up++, usize--;
  86   while (vp[0] == 0)
  87     vp++, vsize--;
  88 #endif
  89
  90   minsize = MIN (usize, vsize);
  91   maxsize = usize + vsize - minsize;
  92
  93   up -= minsize;                /* point at most significant common limb */
  94   vp -= minsize;                /* point at most significant common limb */
  95
  96   /* Compare the most significant part which has explicit limbs for U and V. */
  97   for (i = minsize - 1; i > 0; i--)
  98     {
  99       if (up[i] != vp[i])
 100         return 0;
 101     }
 102
 103   if (minsize != maxsize)
 104     {
 105       if (up[0] != vp[0])
 106         return 0;
 107     }
 108
 109   /* Now either U or V has its limbs consumed.  Check the the other operand
 110      has just zeros in the corresponding, relevant part.  */
 111
 112   if (usize > vsize)
 113     p = up + minsize - maxsize;
 114   else
 115     p = vp + minsize - maxsize;
 116
 117   for (i = maxsize - minsize - 1; i > 0; i--)
 118     {
 119       if (p[i] != 0)
 120         return 0;
 121     }
 122
 123   n_bits -= (maxsize - 1) * GMP_NUMB_BITS;
 124
 125   if (minsize != maxsize)
 126     diff = p[0];
 127   else
 128     diff = up[0] ^ vp[0];
 129
 130   if (n_bits < GMP_NUMB_BITS)
 131     diff >>= GMP_NUMB_BITS - n_bits;
 132
 133   return diff == 0;
 134 }