contrib/gmp/mpz/and.c

   1 /* mpz_and -- Logical and.
   2
   3 Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2003, 2005 Free Software
   4 Foundation, 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
  24 void
  25 mpz_and (mpz_ptr res, mpz_srcptr op1, mpz_srcptr op2)
  26 {
  27   mp_srcptr op1_ptr, op2_ptr;
  28   mp_size_t op1_size, op2_size;
  29   mp_ptr res_ptr;
  30   mp_size_t res_size;
  31   mp_size_t i;
  32   TMP_DECL;
  33
  34   TMP_MARK;
  35   op1_size = SIZ(op1);
  36   op2_size = SIZ(op2);
  37
  38   op1_ptr = PTR(op1);
  39   op2_ptr = PTR(op2);
  40   res_ptr = PTR(res);
  41
  42   if (op1_size >= 0)
  43     {
  44       if (op2_size >= 0)
  45         {
  46           res_size = MIN (op1_size, op2_size);
  47           /* First loop finds the size of the result.  */
  48           for (i = res_size - 1; i >= 0; i--)
  49             if ((op1_ptr[i] & op2_ptr[i]) != 0)
  50               break;
  51           res_size = i + 1;
  52
  53           /* Handle allocation, now then we know exactly how much space is
  54              needed for the result.  */
  55           if (UNLIKELY (ALLOC(res) < res_size))
  56             {
  57               _mpz_realloc (res, res_size);
  58               res_ptr = PTR(res);
  59               /* Don't re-read op1_ptr and op2_ptr.  Since res_size <=
  60                  MIN(op1_size, op2_size), we will not reach this code when op1
  61                  is identical to res or op2 is identical to res.  */
  62             }
  63
  64           SIZ(res) = res_size;
  65           if (LIKELY (res_size != 0))
  66             mpn_and_n (res_ptr, op1_ptr, op2_ptr, res_size);
  67           return;
  68         }
  69       else /* op2_size < 0 */
  70         {
  71           /* Fall through to the code at the end of the function.  */
  72         }
  73     }
  74   else
  75     {
  76       if (op2_size < 0)
  77         {
  78           mp_ptr opx;
  79           mp_limb_t cy;
  80           mp_size_t res_alloc;
  81
  82           /* Both operands are negative, so will be the result.
  83              -((-OP1) & (-OP2)) = -(~(OP1 - 1) & ~(OP2 - 1)) =
  84              = ~(~(OP1 - 1) & ~(OP2 - 1)) + 1 =
  85              = ((OP1 - 1) | (OP2 - 1)) + 1      */
  86
  87           /* It might seem as we could end up with an (invalid) result with
  88              a leading zero-limb here when one of the operands is of the
  89              type 1,,0,,..,,.0.  But some analysis shows that we surely
  90              would get carry into the zero-limb in this situation...  */
  91
  92           op1_size = -op1_size;
  93           op2_size = -op2_size;
  94
  95           res_alloc = 1 + MAX (op1_size, op2_size);
  96
  97           opx = TMP_ALLOC_LIMBS (op1_size);
  98           mpn_sub_1 (opx, op1_ptr, op1_size, (mp_limb_t) 1);
  99           op1_ptr = opx;
 100
 101           opx = TMP_ALLOC_LIMBS (op2_size);
 102           mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1);
 103           op2_ptr = opx;
 104
 105           if (ALLOC(res) < res_alloc)
 106             {
 107               _mpz_realloc (res, res_alloc);
 108               res_ptr = PTR(res);
 109               /* Don't re-read OP1_PTR and OP2_PTR.  They point to temporary
 110                  space--never to the space PTR(res) used to point to before
 111                  reallocation.  */
 112             }
 113
 114           if (op1_size >= op2_size)
 115             {
 116               MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size,
 117                         op1_size - op2_size);
 118               for (i = op2_size - 1; i >= 0; i--)
 119                 res_ptr[i] = op1_ptr[i] | op2_ptr[i];
 120               res_size = op1_size;
 121             }
 122           else
 123             {
 124               MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size,
 125                         op2_size - op1_size);
 126               for (i = op1_size - 1; i >= 0; i--)
 127                 res_ptr[i] = op1_ptr[i] | op2_ptr[i];
 128               res_size = op2_size;
 129             }
 130
 131           cy = mpn_add_1 (res_ptr, res_ptr, res_size, (mp_limb_t) 1);
 132           if (cy)
 133             {
 134               res_ptr[res_size] = cy;
 135               res_size++;
 136             }
 137
 138           SIZ(res) = -res_size;
 139           TMP_FREE;
 140           return;
 141         }
 142       else
 143         {
 144           /* We should compute -OP1 & OP2.  Swap OP1 and OP2 and fall
 145              through to the code that handles OP1 & -OP2.  */
 146           MPZ_SRCPTR_SWAP (op1, op2);
 147           MPN_SRCPTR_SWAP (op1_ptr,op1_size, op2_ptr,op2_size);
 148         }
 149
 150     }
 151
 152   {
 153 #if ANDNEW
 154     mp_size_t op2_lim;
 155     mp_size_t count;
 156
 157     /* OP2 must be negated as with infinite precision.
 158
 159        Scan from the low end for a non-zero limb.  The first non-zero
 160        limb is simply negated (two's complement).  Any subsequent
 161        limbs are one's complemented.  Of course, we don't need to
 162        handle more limbs than there are limbs in the other, positive
 163        operand as the result for those limbs is going to become zero
 164        anyway.  */
 165
 166     /* Scan for the least significant non-zero OP2 limb, and zero the
 167        result meanwhile for those limb positions.  (We will surely
 168        find a non-zero limb, so we can write the loop with one
 169        termination condition only.)  */
 170     for (i = 0; op2_ptr[i] == 0; i++)
 171       res_ptr[i] = 0;
 172     op2_lim = i;
 173
 174     op2_size = -op2_size;
 175
 176     if (op1_size <= op2_size)
 177       {
 178         /* The ones-extended OP2 is >= than the zero-extended OP1.
 179            RES_SIZE <= OP1_SIZE.  Find the exact size.  */
 180         for (i = op1_size - 1; i > op2_lim; i--)
 181           if ((op1_ptr[i] & ~op2_ptr[i]) != 0)
 182             break;
 183         res_size = i + 1;
 184         for (i = res_size - 1; i > op2_lim; i--)
 185           res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
 186         res_ptr[op2_lim] = op1_ptr[op2_lim] & -op2_ptr[op2_lim];
 187         /* Yes, this *can* happen!  */
 188         MPN_NORMALIZE (res_ptr, res_size);
 189       }
 190     else
 191       {
 192         /* The ones-extended OP2 is < than the zero-extended OP1.
 193            RES_SIZE == OP1_SIZE, since OP1 is normalized.  */
 194         res_size = op1_size;
 195         MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, op1_size - op2_size);
 196         for (i = op2_size - 1; i > op2_lim; i--)
 197           res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
 198         res_ptr[op2_lim] = op1_ptr[op2_lim] & -op2_ptr[op2_lim];
 199       }
 200
 201     SIZ(res) = res_size;
 202 #else
 203
 204     /* OP1 is positive and zero-extended,
 205        OP2 is negative and ones-extended.
 206        The result will be positive.
 207        OP1 & -OP2 = OP1 & ~(OP2 - 1).  */
 208
 209     mp_ptr opx;
 210
 211     op2_size = -op2_size;
 212     opx = TMP_ALLOC_LIMBS (op2_size);
 213     mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1);
 214     op2_ptr = opx;
 215
 216     if (op1_size > op2_size)
 217       {
 218         /* The result has the same size as OP1, since OP1 is normalized
 219            and longer than the ones-extended OP2.  */
 220         res_size = op1_size;
 221
 222         /* Handle allocation, now then we know exactly how much space is
 223            needed for the result.  */
 224         if (ALLOC(res) < res_size)
 225           {
 226             _mpz_realloc (res, res_size);
 227             res_ptr = PTR(res);
 228             /* Don't re-read OP1_PTR or OP2_PTR.  Since res_size = op1_size,
 229                we will not reach this code when op1 is identical to res.
 230                OP2_PTR points to temporary space.  */
 231           }
 232
 233         MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, res_size - op2_size);
 234         for (i = op2_size - 1; i >= 0; i--)
 235           res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
 236
 237         SIZ(res) = res_size;
 238       }
 239     else
 240       {
 241         /* Find out the exact result size.  Ignore the high limbs of OP2,
 242            OP1 is zero-extended and would make the result zero.  */
 243         for (i = op1_size - 1; i >= 0; i--)
 244           if ((op1_ptr[i] & ~op2_ptr[i]) != 0)
 245             break;
 246         res_size = i + 1;
 247
 248         /* Handle allocation, now then we know exactly how much space is
 249            needed for the result.  */
 250         if (ALLOC(res) < res_size)
 251           {
 252             _mpz_realloc (res, res_size);
 253             res_ptr = PTR(res);
 254             /* Don't re-read OP1_PTR.  Since res_size <= op1_size, we will
 255                not reach this code when op1 is identical to res.  */
 256             /* Don't re-read OP2_PTR.  It points to temporary space--never
 257                to the space PTR(res) used to point to before reallocation.  */
 258           }
 259
 260         for (i = res_size - 1; i >= 0; i--)
 261           res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
 262
 263         SIZ(res) = res_size;
 264       }
 265 #endif
 266   }
 267   TMP_FREE;
 268 }