xref: /dflybsd-src/contrib/gmp/mpz/and.c (revision 86d7f5d305c6adaa56ff4582ece9859d73106103)
186d7f5d3SJohn Marino /* mpz_and -- Logical and.
286d7f5d3SJohn Marino 
386d7f5d3SJohn Marino Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2003, 2005 Free Software
486d7f5d3SJohn Marino Foundation, Inc.
586d7f5d3SJohn Marino 
686d7f5d3SJohn Marino This file is part of the GNU MP Library.
786d7f5d3SJohn Marino 
886d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
986d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
1086d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
1186d7f5d3SJohn Marino option) any later version.
1286d7f5d3SJohn Marino 
1386d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
1486d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1586d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
1686d7f5d3SJohn Marino License for more details.
1786d7f5d3SJohn Marino 
1886d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
1986d7f5d3SJohn Marino along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
2086d7f5d3SJohn Marino 
2186d7f5d3SJohn Marino #include "gmp.h"
2286d7f5d3SJohn Marino #include "gmp-impl.h"
2386d7f5d3SJohn Marino 
2486d7f5d3SJohn Marino void
mpz_and(mpz_ptr res,mpz_srcptr op1,mpz_srcptr op2)2586d7f5d3SJohn Marino mpz_and (mpz_ptr res, mpz_srcptr op1, mpz_srcptr op2)
2686d7f5d3SJohn Marino {
2786d7f5d3SJohn Marino   mp_srcptr op1_ptr, op2_ptr;
2886d7f5d3SJohn Marino   mp_size_t op1_size, op2_size;
2986d7f5d3SJohn Marino   mp_ptr res_ptr;
3086d7f5d3SJohn Marino   mp_size_t res_size;
3186d7f5d3SJohn Marino   mp_size_t i;
3286d7f5d3SJohn Marino   TMP_DECL;
3386d7f5d3SJohn Marino 
3486d7f5d3SJohn Marino   TMP_MARK;
3586d7f5d3SJohn Marino   op1_size = SIZ(op1);
3686d7f5d3SJohn Marino   op2_size = SIZ(op2);
3786d7f5d3SJohn Marino 
3886d7f5d3SJohn Marino   op1_ptr = PTR(op1);
3986d7f5d3SJohn Marino   op2_ptr = PTR(op2);
4086d7f5d3SJohn Marino   res_ptr = PTR(res);
4186d7f5d3SJohn Marino 
4286d7f5d3SJohn Marino   if (op1_size >= 0)
4386d7f5d3SJohn Marino     {
4486d7f5d3SJohn Marino       if (op2_size >= 0)
4586d7f5d3SJohn Marino 	{
4686d7f5d3SJohn Marino 	  res_size = MIN (op1_size, op2_size);
4786d7f5d3SJohn Marino 	  /* First loop finds the size of the result.  */
4886d7f5d3SJohn Marino 	  for (i = res_size - 1; i >= 0; i--)
4986d7f5d3SJohn Marino 	    if ((op1_ptr[i] & op2_ptr[i]) != 0)
5086d7f5d3SJohn Marino 	      break;
5186d7f5d3SJohn Marino 	  res_size = i + 1;
5286d7f5d3SJohn Marino 
5386d7f5d3SJohn Marino 	  /* Handle allocation, now then we know exactly how much space is
5486d7f5d3SJohn Marino 	     needed for the result.  */
5586d7f5d3SJohn Marino 	  if (UNLIKELY (ALLOC(res) < res_size))
5686d7f5d3SJohn Marino 	    {
5786d7f5d3SJohn Marino 	      _mpz_realloc (res, res_size);
5886d7f5d3SJohn Marino 	      res_ptr = PTR(res);
5986d7f5d3SJohn Marino 	      /* Don't re-read op1_ptr and op2_ptr.  Since res_size <=
6086d7f5d3SJohn Marino 		 MIN(op1_size, op2_size), we will not reach this code when op1
6186d7f5d3SJohn Marino 		 is identical to res or op2 is identical to res.  */
6286d7f5d3SJohn Marino 	    }
6386d7f5d3SJohn Marino 
6486d7f5d3SJohn Marino 	  SIZ(res) = res_size;
6586d7f5d3SJohn Marino           if (LIKELY (res_size != 0))
6686d7f5d3SJohn Marino             mpn_and_n (res_ptr, op1_ptr, op2_ptr, res_size);
6786d7f5d3SJohn Marino 	  return;
6886d7f5d3SJohn Marino 	}
6986d7f5d3SJohn Marino       else /* op2_size < 0 */
7086d7f5d3SJohn Marino 	{
7186d7f5d3SJohn Marino 	  /* Fall through to the code at the end of the function.  */
7286d7f5d3SJohn Marino 	}
7386d7f5d3SJohn Marino     }
7486d7f5d3SJohn Marino   else
7586d7f5d3SJohn Marino     {
7686d7f5d3SJohn Marino       if (op2_size < 0)
7786d7f5d3SJohn Marino 	{
7886d7f5d3SJohn Marino 	  mp_ptr opx;
7986d7f5d3SJohn Marino 	  mp_limb_t cy;
8086d7f5d3SJohn Marino 	  mp_size_t res_alloc;
8186d7f5d3SJohn Marino 
8286d7f5d3SJohn Marino 	  /* Both operands are negative, so will be the result.
8386d7f5d3SJohn Marino 	     -((-OP1) & (-OP2)) = -(~(OP1 - 1) & ~(OP2 - 1)) =
8486d7f5d3SJohn Marino 	     = ~(~(OP1 - 1) & ~(OP2 - 1)) + 1 =
8586d7f5d3SJohn Marino 	     = ((OP1 - 1) | (OP2 - 1)) + 1      */
8686d7f5d3SJohn Marino 
8786d7f5d3SJohn Marino 	  /* It might seem as we could end up with an (invalid) result with
8886d7f5d3SJohn Marino 	     a leading zero-limb here when one of the operands is of the
8986d7f5d3SJohn Marino 	     type 1,,0,,..,,.0.  But some analysis shows that we surely
9086d7f5d3SJohn Marino 	     would get carry into the zero-limb in this situation...  */
9186d7f5d3SJohn Marino 
9286d7f5d3SJohn Marino 	  op1_size = -op1_size;
9386d7f5d3SJohn Marino 	  op2_size = -op2_size;
9486d7f5d3SJohn Marino 
9586d7f5d3SJohn Marino 	  res_alloc = 1 + MAX (op1_size, op2_size);
9686d7f5d3SJohn Marino 
9786d7f5d3SJohn Marino 	  opx = TMP_ALLOC_LIMBS (op1_size);
9886d7f5d3SJohn Marino 	  mpn_sub_1 (opx, op1_ptr, op1_size, (mp_limb_t) 1);
9986d7f5d3SJohn Marino 	  op1_ptr = opx;
10086d7f5d3SJohn Marino 
10186d7f5d3SJohn Marino 	  opx = TMP_ALLOC_LIMBS (op2_size);
10286d7f5d3SJohn Marino 	  mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1);
10386d7f5d3SJohn Marino 	  op2_ptr = opx;
10486d7f5d3SJohn Marino 
10586d7f5d3SJohn Marino 	  if (ALLOC(res) < res_alloc)
10686d7f5d3SJohn Marino 	    {
10786d7f5d3SJohn Marino 	      _mpz_realloc (res, res_alloc);
10886d7f5d3SJohn Marino 	      res_ptr = PTR(res);
10986d7f5d3SJohn Marino 	      /* Don't re-read OP1_PTR and OP2_PTR.  They point to temporary
11086d7f5d3SJohn Marino 		 space--never to the space PTR(res) used to point to before
11186d7f5d3SJohn Marino 		 reallocation.  */
11286d7f5d3SJohn Marino 	    }
11386d7f5d3SJohn Marino 
11486d7f5d3SJohn Marino 	  if (op1_size >= op2_size)
11586d7f5d3SJohn Marino 	    {
11686d7f5d3SJohn Marino 	      MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size,
11786d7f5d3SJohn Marino 			op1_size - op2_size);
11886d7f5d3SJohn Marino 	      for (i = op2_size - 1; i >= 0; i--)
11986d7f5d3SJohn Marino 		res_ptr[i] = op1_ptr[i] | op2_ptr[i];
12086d7f5d3SJohn Marino 	      res_size = op1_size;
12186d7f5d3SJohn Marino 	    }
12286d7f5d3SJohn Marino 	  else
12386d7f5d3SJohn Marino 	    {
12486d7f5d3SJohn Marino 	      MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size,
12586d7f5d3SJohn Marino 			op2_size - op1_size);
12686d7f5d3SJohn Marino 	      for (i = op1_size - 1; i >= 0; i--)
12786d7f5d3SJohn Marino 		res_ptr[i] = op1_ptr[i] | op2_ptr[i];
12886d7f5d3SJohn Marino 	      res_size = op2_size;
12986d7f5d3SJohn Marino 	    }
13086d7f5d3SJohn Marino 
13186d7f5d3SJohn Marino 	  cy = mpn_add_1 (res_ptr, res_ptr, res_size, (mp_limb_t) 1);
13286d7f5d3SJohn Marino 	  if (cy)
13386d7f5d3SJohn Marino 	    {
13486d7f5d3SJohn Marino 	      res_ptr[res_size] = cy;
13586d7f5d3SJohn Marino 	      res_size++;
13686d7f5d3SJohn Marino 	    }
13786d7f5d3SJohn Marino 
13886d7f5d3SJohn Marino 	  SIZ(res) = -res_size;
13986d7f5d3SJohn Marino 	  TMP_FREE;
14086d7f5d3SJohn Marino 	  return;
14186d7f5d3SJohn Marino 	}
14286d7f5d3SJohn Marino       else
14386d7f5d3SJohn Marino 	{
14486d7f5d3SJohn Marino 	  /* We should compute -OP1 & OP2.  Swap OP1 and OP2 and fall
14586d7f5d3SJohn Marino 	     through to the code that handles OP1 & -OP2.  */
14686d7f5d3SJohn Marino           MPZ_SRCPTR_SWAP (op1, op2);
14786d7f5d3SJohn Marino           MPN_SRCPTR_SWAP (op1_ptr,op1_size, op2_ptr,op2_size);
14886d7f5d3SJohn Marino 	}
14986d7f5d3SJohn Marino 
15086d7f5d3SJohn Marino     }
15186d7f5d3SJohn Marino 
15286d7f5d3SJohn Marino   {
15386d7f5d3SJohn Marino #if ANDNEW
15486d7f5d3SJohn Marino     mp_size_t op2_lim;
15586d7f5d3SJohn Marino     mp_size_t count;
15686d7f5d3SJohn Marino 
15786d7f5d3SJohn Marino     /* OP2 must be negated as with infinite precision.
15886d7f5d3SJohn Marino 
15986d7f5d3SJohn Marino        Scan from the low end for a non-zero limb.  The first non-zero
16086d7f5d3SJohn Marino        limb is simply negated (two's complement).  Any subsequent
16186d7f5d3SJohn Marino        limbs are one's complemented.  Of course, we don't need to
16286d7f5d3SJohn Marino        handle more limbs than there are limbs in the other, positive
16386d7f5d3SJohn Marino        operand as the result for those limbs is going to become zero
16486d7f5d3SJohn Marino        anyway.  */
16586d7f5d3SJohn Marino 
16686d7f5d3SJohn Marino     /* Scan for the least significant non-zero OP2 limb, and zero the
16786d7f5d3SJohn Marino        result meanwhile for those limb positions.  (We will surely
16886d7f5d3SJohn Marino        find a non-zero limb, so we can write the loop with one
16986d7f5d3SJohn Marino        termination condition only.)  */
17086d7f5d3SJohn Marino     for (i = 0; op2_ptr[i] == 0; i++)
17186d7f5d3SJohn Marino       res_ptr[i] = 0;
17286d7f5d3SJohn Marino     op2_lim = i;
17386d7f5d3SJohn Marino 
17486d7f5d3SJohn Marino     op2_size = -op2_size;
17586d7f5d3SJohn Marino 
17686d7f5d3SJohn Marino     if (op1_size <= op2_size)
17786d7f5d3SJohn Marino       {
17886d7f5d3SJohn Marino 	/* The ones-extended OP2 is >= than the zero-extended OP1.
17986d7f5d3SJohn Marino 	   RES_SIZE <= OP1_SIZE.  Find the exact size.  */
18086d7f5d3SJohn Marino 	for (i = op1_size - 1; i > op2_lim; i--)
18186d7f5d3SJohn Marino 	  if ((op1_ptr[i] & ~op2_ptr[i]) != 0)
18286d7f5d3SJohn Marino 	    break;
18386d7f5d3SJohn Marino 	res_size = i + 1;
18486d7f5d3SJohn Marino 	for (i = res_size - 1; i > op2_lim; i--)
18586d7f5d3SJohn Marino 	  res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
18686d7f5d3SJohn Marino 	res_ptr[op2_lim] = op1_ptr[op2_lim] & -op2_ptr[op2_lim];
18786d7f5d3SJohn Marino 	/* Yes, this *can* happen!  */
18886d7f5d3SJohn Marino 	MPN_NORMALIZE (res_ptr, res_size);
18986d7f5d3SJohn Marino       }
19086d7f5d3SJohn Marino     else
19186d7f5d3SJohn Marino       {
19286d7f5d3SJohn Marino 	/* The ones-extended OP2 is < than the zero-extended OP1.
19386d7f5d3SJohn Marino 	   RES_SIZE == OP1_SIZE, since OP1 is normalized.  */
19486d7f5d3SJohn Marino 	res_size = op1_size;
19586d7f5d3SJohn Marino 	MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, op1_size - op2_size);
19686d7f5d3SJohn Marino 	for (i = op2_size - 1; i > op2_lim; i--)
19786d7f5d3SJohn Marino 	  res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
19886d7f5d3SJohn Marino 	res_ptr[op2_lim] = op1_ptr[op2_lim] & -op2_ptr[op2_lim];
19986d7f5d3SJohn Marino       }
20086d7f5d3SJohn Marino 
20186d7f5d3SJohn Marino     SIZ(res) = res_size;
20286d7f5d3SJohn Marino #else
20386d7f5d3SJohn Marino 
20486d7f5d3SJohn Marino     /* OP1 is positive and zero-extended,
20586d7f5d3SJohn Marino        OP2 is negative and ones-extended.
20686d7f5d3SJohn Marino        The result will be positive.
20786d7f5d3SJohn Marino        OP1 & -OP2 = OP1 & ~(OP2 - 1).  */
20886d7f5d3SJohn Marino 
20986d7f5d3SJohn Marino     mp_ptr opx;
21086d7f5d3SJohn Marino 
21186d7f5d3SJohn Marino     op2_size = -op2_size;
21286d7f5d3SJohn Marino     opx = TMP_ALLOC_LIMBS (op2_size);
21386d7f5d3SJohn Marino     mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1);
21486d7f5d3SJohn Marino     op2_ptr = opx;
21586d7f5d3SJohn Marino 
21686d7f5d3SJohn Marino     if (op1_size > op2_size)
21786d7f5d3SJohn Marino       {
21886d7f5d3SJohn Marino 	/* The result has the same size as OP1, since OP1 is normalized
21986d7f5d3SJohn Marino 	   and longer than the ones-extended OP2.  */
22086d7f5d3SJohn Marino 	res_size = op1_size;
22186d7f5d3SJohn Marino 
22286d7f5d3SJohn Marino 	/* Handle allocation, now then we know exactly how much space is
22386d7f5d3SJohn Marino 	   needed for the result.  */
22486d7f5d3SJohn Marino 	if (ALLOC(res) < res_size)
22586d7f5d3SJohn Marino 	  {
22686d7f5d3SJohn Marino 	    _mpz_realloc (res, res_size);
22786d7f5d3SJohn Marino 	    res_ptr = PTR(res);
22886d7f5d3SJohn Marino 	    /* Don't re-read OP1_PTR or OP2_PTR.  Since res_size = op1_size,
22986d7f5d3SJohn Marino 	       we will not reach this code when op1 is identical to res.
23086d7f5d3SJohn Marino 	       OP2_PTR points to temporary space.  */
23186d7f5d3SJohn Marino 	  }
23286d7f5d3SJohn Marino 
23386d7f5d3SJohn Marino 	MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, res_size - op2_size);
23486d7f5d3SJohn Marino 	for (i = op2_size - 1; i >= 0; i--)
23586d7f5d3SJohn Marino 	  res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
23686d7f5d3SJohn Marino 
23786d7f5d3SJohn Marino 	SIZ(res) = res_size;
23886d7f5d3SJohn Marino       }
23986d7f5d3SJohn Marino     else
24086d7f5d3SJohn Marino       {
24186d7f5d3SJohn Marino 	/* Find out the exact result size.  Ignore the high limbs of OP2,
24286d7f5d3SJohn Marino 	   OP1 is zero-extended and would make the result zero.  */
24386d7f5d3SJohn Marino 	for (i = op1_size - 1; i >= 0; i--)
24486d7f5d3SJohn Marino 	  if ((op1_ptr[i] & ~op2_ptr[i]) != 0)
24586d7f5d3SJohn Marino 	    break;
24686d7f5d3SJohn Marino 	res_size = i + 1;
24786d7f5d3SJohn Marino 
24886d7f5d3SJohn Marino 	/* Handle allocation, now then we know exactly how much space is
24986d7f5d3SJohn Marino 	   needed for the result.  */
25086d7f5d3SJohn Marino 	if (ALLOC(res) < res_size)
25186d7f5d3SJohn Marino 	  {
25286d7f5d3SJohn Marino 	    _mpz_realloc (res, res_size);
25386d7f5d3SJohn Marino 	    res_ptr = PTR(res);
25486d7f5d3SJohn Marino 	    /* Don't re-read OP1_PTR.  Since res_size <= op1_size, we will
25586d7f5d3SJohn Marino 	       not reach this code when op1 is identical to res.  */
25686d7f5d3SJohn Marino 	    /* Don't re-read OP2_PTR.  It points to temporary space--never
25786d7f5d3SJohn Marino 	       to the space PTR(res) used to point to before reallocation.  */
25886d7f5d3SJohn Marino 	  }
25986d7f5d3SJohn Marino 
26086d7f5d3SJohn Marino 	for (i = res_size - 1; i >= 0; i--)
26186d7f5d3SJohn Marino 	  res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
26286d7f5d3SJohn Marino 
26386d7f5d3SJohn Marino 	SIZ(res) = res_size;
26486d7f5d3SJohn Marino       }
26586d7f5d3SJohn Marino #endif
26686d7f5d3SJohn Marino   }
26786d7f5d3SJohn Marino   TMP_FREE;
26886d7f5d3SJohn Marino }
269