gcdext.c

a very popular packet of cryptography tools,it encloses the most common used algorithm and protocols

	      sub_ddmmss (Th, Tl, uh, ul, t1, t0);
	      uh = vh, ul = vl;
	      vh = Th, vl = Tl;

	      asign = ~asign;

	      add_ssaaaa (dh, dl, vh, vl, 0, C);
/*	      if (dh == 0)	should never happen
		break;	       */
	      sub_ddmmss (nh, nl, uh, ul, 0, A);
	      q1 = div2 (nh, nl, dh, dl);

	      sub_ddmmss (dh, dl, vh, vl, 0, D);
	      if (dh == 0)
		break;
	      add_ssaaaa (nh, nl, uh, ul, 0, B);
	      q2 = div2 (nh, nl, dh, dl);

	      if (q1 != q2)
		break;

	      Tac = A + q1 * C;
	      if (GMP_NAIL_BITS != 0 && Tac > GMP_NUMB_MAX)
		break;
	      Tbd = B + q1 * D;
	      if (GMP_NAIL_BITS != 0 && Tbd > GMP_NUMB_MAX)
		break;
	      A = C;
	      C = Tac;
	      B = D;
	      D = Tbd;
	      umul_ppmm (t1, t0, q1, vl);
	      t1 += q1 * vh;
	      sub_ddmmss (Th, Tl, uh, ul, t1, t0);
	      uh = vh, ul = vl;
	      vh = Th, vl = Tl;

	      asign = ~asign;
	    }
#if EXTEND
	  if (asign)
	    sign = -sign;
#endif
	}
      else /* Same, but using single-limb calculations.  */
	{
	  mp_limb_t uh, vh;
	  /* Make UH be the most significant limb of U, and make VH be
	     corresponding bits from V.  */
	  uh = up[size - 1];
	  vh = vp[size - 1];
	  count_leading_zeros (cnt, uh);
#if GMP_NAIL_BITS == 0
	  if (cnt != 0)
	    {
	      uh = (uh << cnt) | (up[size - 2] >> (GMP_LIMB_BITS - cnt));
	      vh = (vh << cnt) | (vp[size - 2] >> (GMP_LIMB_BITS - cnt));
	    }
#else
	  uh <<= cnt;
	  vh <<= cnt;
	  if (cnt < GMP_NUMB_BITS)
	    {
	      uh |= up[size - 2] >> (GMP_NUMB_BITS - cnt);
	      vh |= vp[size - 2] >> (GMP_NUMB_BITS - cnt);
	    }
	  else
	    {
	      uh |= up[size - 2] << cnt - GMP_NUMB_BITS;
	      vh |= vp[size - 2] << cnt - GMP_NUMB_BITS;
	      if (size >= 3)
		{
		  uh |= up[size - 3] >> 2 * GMP_NUMB_BITS - cnt;
		  vh |= vp[size - 3] >> 2 * GMP_NUMB_BITS - cnt;
		}
	    }
#endif

	  A = 1;
	  B = 0;
	  C = 0;
	  D = 1;

	  asign = 0;
	  for (;;)
	    {
	      mp_limb_t q, T;
	      if (vh - C == 0 || vh + D == 0)
		break;

	      q = (uh + A) / (vh - C);
	      if (q != (uh - B) / (vh + D))
		break;

	      T = A + q * C;
	      A = C;
	      C = T;
	      T = B + q * D;
	      B = D;
	      D = T;
	      T = uh - q * vh;
	      uh = vh;
	      vh = T;

	      asign = ~asign;

	      if (vh - D == 0)
		break;

	      q = (uh - A) / (vh + C);
	      if (q != (uh + B) / (vh - D))
		break;

	      T = A + q * C;
	      A = C;
	      C = T;
	      T = B + q * D;
	      B = D;
	      D = T;
	      T = uh - q * vh;
	      uh = vh;
	      vh = T;

	      asign = ~asign;
	    }
#if EXTEND
	  if (asign)
	    sign = -sign;
#endif
	}

#if RECORD
      max = MAX (A, max);  max = MAX (B, max);
      max = MAX (C, max);  max = MAX (D, max);
#endif

      if (B == 0)
	{
	  /* This is quite rare.  I.e., optimize something else!  */

	  mpn_tdiv_qr (wp, up, (mp_size_t) 0, up, size, vp, vsize);

#if EXTEND
	  MPN_COPY (tp, s0p, ssize);
	  {
	    mp_size_t qsize;
	    mp_size_t i;

	    qsize = size - vsize + 1; /* size of stored quotient from division */
	    MPN_ZERO (s1p + ssize, qsize); /* zero s1 too */

	    for (i = 0; i < qsize; i++)
	      {
		mp_limb_t cy;
		cy = mpn_addmul_1 (tp + i, s1p, ssize, wp[i]);
		tp[ssize + i] = cy;
	      }

	    ssize += qsize;
	    ssize -= tp[ssize - 1] == 0;
	  }

	  sign = -sign;
	  MP_PTR_SWAP (s0p, s1p);
	  MP_PTR_SWAP (s1p, tp);
#endif
	  size = vsize;
	  MP_PTR_SWAP (up, vp);
	}
      else
	{
#if EXTEND
	  mp_size_t tsize, wsize;
#endif
	  /* T = U*A + V*B
	     W = U*C + V*D
	     U = T
	     V = W	   */

#if STAT
	  { mp_limb_t x; x = A | B | C | D; count_leading_zeros (cnt, x);
	  arr[GMP_LIMB_BITS - cnt]++; }
#endif
	  if (A == 0)
	    {
	      /* B == 1 and C == 1 (D is arbitrary) */
	      mp_limb_t cy;
	      MPN_COPY (tp, vp, size);
	      MPN_COPY (wp, up, size);
	      mpn_submul_1 (wp, vp, size, D);
	      MP_PTR_SWAP (tp, up);
	      MP_PTR_SWAP (wp, vp);
#if EXTEND
	      MPN_COPY (tp, s1p, ssize);
	      tsize = ssize;
	      tp[ssize] = 0;	/* must zero since wp might spill below */
	      MPN_COPY (wp, s0p, ssize);
	      cy = mpn_addmul_1 (wp, s1p, ssize, D);
	      wp[ssize] = cy;
	      wsize = ssize + (cy != 0);
	      MP_PTR_SWAP (tp, s0p);
	      MP_PTR_SWAP (wp, s1p);
	      ssize = MAX (wsize, tsize);
#endif
	    }
	  else
	    {
	      mp_limb_t cy, cy1, cy2;

	      if (asign)
		{
		  mpn_mul_1 (tp, vp, size, B);
		  mpn_submul_1 (tp, up, size, A);
		  mpn_mul_1 (wp, up, size, C);
		  mpn_submul_1 (wp, vp, size, D);
		}
	      else
		{
		  mpn_mul_1 (tp, up, size, A);
		  mpn_submul_1 (tp, vp, size, B);
		  mpn_mul_1 (wp, vp, size, D);
		  mpn_submul_1 (wp, up, size, C);
		}
	      MP_PTR_SWAP (tp, up);
	      MP_PTR_SWAP (wp, vp);
#if EXTEND
	      /* Compute new s0 */
	      cy1 = mpn_mul_1 (tp, s0p, ssize, A);
	      cy2 = mpn_addmul_1 (tp, s1p, ssize, B);
	      cy = cy1 + cy2;
	      tp[ssize] = cy & GMP_NUMB_MASK;
	      tsize = ssize + (cy != 0);
#if GMP_NAIL_BITS == 0
	      if (cy < cy1)
#else
	      if (cy > GMP_NUMB_MAX)
#endif
		{
		  tp[tsize] = 1;
		  wp[tsize] = 0;
		  tsize++;
		  /* This happens just for nails, since we get more work done
		     per numb there.  */
		}

	      /* Compute new s1 */
	      cy1 = mpn_mul_1 (wp, s1p, ssize, D);
	      cy2 = mpn_addmul_1 (wp, s0p, ssize, C);
	      cy = cy1 + cy2;
	      wp[ssize] = cy & GMP_NUMB_MASK;
	      wsize = ssize + (cy != 0);
#if GMP_NAIL_BITS == 0
	      if (cy < cy1)
#else
	      if (cy > GMP_NUMB_MAX)
#endif
		{
		  wp[wsize] = 1;
		  if (wsize >= tsize)
		    tp[wsize] = 0;
		  wsize++;
		}

	      MP_PTR_SWAP (tp, s0p);
	      MP_PTR_SWAP (wp, s1p);
	      ssize = MAX (wsize, tsize);
#endif
	    }
	  size -= up[size - 1] == 0;
#if GMP_NAIL_BITS != 0
	  size -= up[size - 1] == 0;
#endif
	}

#if WANT_GCDEXT_ONE_STEP
      TMP_FREE (mark);
      return 0;
#endif
    }

#if RECORD
  printf ("max: %lx\n", max);
#endif

#if STAT
 {int i; for (i = 0; i <= GMP_LIMB_BITS; i++) printf ("%d:%d\n", i, arr[i]);}
#endif

  if (vsize == 0)
    {
      if (gp != up && gp != 0)
	MPN_COPY (gp, up, size);
#if EXTEND
      MPN_NORMALIZE (s0p, ssize);
      if (orig_s0p != s0p)
	MPN_COPY (orig_s0p, s0p, ssize);
      *s0size = sign >= 0 ? ssize : -ssize;
#endif
      TMP_FREE (mark);
      return size;
    }
  else
    {
      mp_limb_t vl, ul, t;
#if EXTEND
      mp_size_t qsize, i;
#endif
      vl = vp[0];
#if EXTEND
      t = mpn_divmod_1 (wp, up, size, vl);

      MPN_COPY (tp, s0p, ssize);

      qsize = size - (wp[size - 1] == 0); /* size of quotient from division */
      if (ssize < qsize)
	{
	  MPN_ZERO (tp + ssize, qsize - ssize);
	  MPN_ZERO (s1p + ssize, qsize); /* zero s1 too */
	  for (i = 0; i < ssize; i++)
	    {
	      mp_limb_t cy;
	      cy = mpn_addmul_1 (tp + i, wp, qsize, s1p[i]);
	      tp[qsize + i] = cy;
	    }
	}
      else
	{
	  MPN_ZERO (s1p + ssize, qsize); /* zero s1 too */
	  for (i = 0; i < qsize; i++)
	    {
	      mp_limb_t cy;
	      cy = mpn_addmul_1 (tp + i, s1p, ssize, wp[i]);
	      tp[ssize + i] = cy;
	    }
	}
      ssize += qsize;
      ssize -= tp[ssize - 1] == 0;

      sign = -sign;
      MP_PTR_SWAP (s0p, s1p);
      MP_PTR_SWAP (s1p, tp);
#else
      t = mpn_mod_1 (up, size, vl);
#endif
      ul = vl;
      vl = t;
      while (vl != 0)
	{
	  mp_limb_t t;
#if EXTEND
	  mp_limb_t q;
	  q = ul / vl;
	  t = ul - q * vl;

	  MPN_COPY (tp, s0p, ssize);

	  MPN_ZERO (s1p + ssize, 1); /* zero s1 too */

	  {
	    mp_limb_t cy;
	    cy = mpn_addmul_1 (tp, s1p, ssize, q);
	    tp[ssize] = cy;
	  }

	  ssize += 1;
	  ssize -= tp[ssize - 1] == 0;

	  sign = -sign;
	  MP_PTR_SWAP (s0p, s1p);
	  MP_PTR_SWAP (s1p, tp);
#else
	  t = ul % vl;
#endif
	  ul = vl;
	  vl = t;
	}
      if (gp != 0)
	gp[0] = ul;
#if EXTEND
      MPN_NORMALIZE (s0p, ssize);
      if (orig_s0p != s0p)
	MPN_COPY (orig_s0p, s0p, ssize);
      *s0size = sign >= 0 ? ssize : -ssize;
#endif
      TMP_FREE (mark);
      return 1;
    }
}