divide (mpn_t a, mpn_t b, mpn_t *q)
{
  /* Algorithm:
     First normalise a and b: a=[a[m-1],...,a[0]], b=[b[n-1],...,b[0]]
     with m>=0 and n>0 (in base beta = 2^GMP_LIMB_BITS).
     If m<n, then q:=0 and r:=a.
     If m>=n=1, perform a single-precision division:
       r:=0, j:=m,
       while j>0 do
         {Here (q[m-1]*beta^(m-1)+...+q[j]*beta^j) * b[0] + r*beta^j =
               = a[m-1]*beta^(m-1)+...+a[j]*beta^j und 0<=r<b[0]<beta}
         j:=j-1, r:=r*beta+a[j], q[j]:=floor(r/b[0]), r:=r-b[0]*q[j].
       Normalise [q[m-1],...,q[0]], yields q.
     If m>=n>1, perform a multiple-precision division:
       We have a/b < beta^(m-n+1).
       s:=intDsize-1-(highest bit in b[n-1]), 0<=s<intDsize.
       Shift a and b left by s bits, copying them. r:=a.
       r=[r[m],...,r[0]], b=[b[n-1],...,b[0]] with b[n-1]>=beta/2.
       For j=m-n,...,0: {Here 0 <= r < b*beta^(j+1).}
         Compute q* :
           q* := floor((r[j+n]*beta+r[j+n-1])/b[n-1]).
           In case of overflow (q* >= beta) set q* := beta-1.
           Compute c2 := ((r[j+n]*beta+r[j+n-1]) - q* * b[n-1])*beta + r[j+n-2]
           and c3 := b[n-2] * q*.
           {We have 0 <= c2 < 2*beta^2, even 0 <= c2 < beta^2 if no overflow
            occurred.  Furthermore 0 <= c3 < beta^2.
            If there was overflow and
            r[j+n]*beta+r[j+n-1] - q* * b[n-1] >= beta, i.e. c2 >= beta^2,
            the next test can be skipped.}
           While c3 > c2, {Here 0 <= c2 < c3 < beta^2}
             Put q* := q* - 1, c2 := c2 + b[n-1]*beta, c3 := c3 - b[n-2].
           If q* > 0:
             Put r := r - b * q* * beta^j. In detail:
               [r[n+j],...,r[j]] := [r[n+j],...,r[j]] - q* * [b[n-1],...,b[0]].
               hence: u:=0, for i:=0 to n-1 do
                              u := u + q* * b[i],
                              r[j+i]:=r[j+i]-(u mod beta) (+ beta, if carry),
                              u:=u div beta (+ 1, if carry in subtraction)
                      r[n+j]:=r[n+j]-u.
               {Since always u = (q* * [b[i-1],...,b[0]] div beta^i) + 1
                               < q* + 1 <= beta,
                the carry u does not overflow.}
             If a negative carry occurs, put q* := q* - 1
               and [r[n+j],...,r[j]] := [r[n+j],...,r[j]] + [0,b[n-1],...,b[0]].
         Set q[j] := q*.
       Normalise [q[m-n],..,q[0]]; this yields the quotient q.
       Shift [r[n-1],...,r[0]] right by s bits and normalise; this yields the
       rest r.
       The room for q[j] can be allocated at the memory location of r[n+j].
     Finally, round-to-even:
       Shift r left by 1 bit.
       If r > b or if r = b and q[0] is odd, q := q+1.
   */
  const mp_limb_t *a_ptr = a.limbs;
  size_t a_len = a.nlimbs;
  const mp_limb_t *b_ptr = b.limbs;
  size_t b_len = b.nlimbs;
  mp_limb_t *roomptr;
  mp_limb_t *tmp_roomptr = NULL;
  mp_limb_t *q_ptr;
  size_t q_len;
  mp_limb_t *r_ptr;
  size_t r_len;

  /* Allocate room for a_len+2 digits.
     (Need a_len+1 digits for the real division and 1 more digit for the
     final rounding of q.)  */
  roomptr = (mp_limb_t *) malloc ((a_len + 2) * sizeof (mp_limb_t));
  if (roomptr == NULL)
    return NULL;

  /* Normalise a.  */
  while (a_len > 0 && a_ptr[a_len - 1] == 0)
    a_len--;

  /* Normalise b.  */
  for (;;)
    {
      if (b_len == 0)
        /* Division by zero.  */
        abort ();
      if (b_ptr[b_len - 1] == 0)
        b_len--;
      else
        break;
    }

  /* Here m = a_len >= 0 and n = b_len > 0.  */

  if (a_len < b_len)
    {
      /* m<n: trivial case.  q=0, r := copy of a.  */
      r_ptr = roomptr;
      r_len = a_len;
      memcpy (r_ptr, a_ptr, a_len * sizeof (mp_limb_t));
      q_ptr = roomptr + a_len;
      q_len = 0;
    }
  else if (b_len == 1)
    {
      /* n=1: single precision division.
         beta^(m-1) <= a < beta^m  ==>  beta^(m-2) <= a/b < beta^m  */
      r_ptr = roomptr;
      q_ptr = roomptr + 1;
      {
        mp_limb_t den = b_ptr[0];
        mp_limb_t remainder = 0;
        const mp_limb_t *sourceptr = a_ptr + a_len;
        mp_limb_t *destptr = q_ptr + a_len;
        size_t count;
        for (count = a_len; count > 0; count--)
          {
            mp_twolimb_t num =
              ((mp_twolimb_t) remainder << GMP_LIMB_BITS) | *--sourceptr;
            *--destptr = num / den;
            remainder = num % den;
          }
        /* Normalise and store r.  */
        if (remainder > 0)
          {
            r_ptr[0] = remainder;
            r_len = 1;
          }
        else
          r_len = 0;
        /* Normalise q.  */
        q_len = a_len;
        if (q_ptr[q_len - 1] == 0)
          q_len--;
      }
    }
  else
    {
      /* n>1: multiple precision division.
         beta^(m-1) <= a < beta^m, beta^(n-1) <= b < beta^n  ==>
         beta^(m-n-1) <= a/b < beta^(m-n+1).  */
      /* Determine s.  */
      size_t s;
      {
        mp_limb_t msd = b_ptr[b_len - 1]; /* = b[n-1], > 0 */
        /* Determine s = GMP_LIMB_BITS - integer_length (msd).
           Code copied from gnulib's integer_length.c.  */
# if __GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4)
        s = __builtin_clz (msd);
# else
#  if defined DBL_EXPBIT0_WORD && defined DBL_EXPBIT0_BIT
        if (GMP_LIMB_BITS <= DBL_MANT_BIT)
          {
            /* Use 'double' operations.
               Assumes an IEEE 754 'double' implementation.  */
#   define DBL_EXP_MASK ((DBL_MAX_EXP - DBL_MIN_EXP) | 7)
#   define DBL_EXP_BIAS (DBL_EXP_MASK / 2 - 1)
#   define NWORDS \
     ((sizeof (double) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
            union { double value; unsigned int word[NWORDS]; } m;

            /* Use a single integer to floating-point conversion.  */
            m.value = msd;

            s = GMP_LIMB_BITS
                - (((m.word[DBL_EXPBIT0_WORD] >> DBL_EXPBIT0_BIT) & DBL_EXP_MASK)
                   - DBL_EXP_BIAS);
          }
        else
#   undef NWORDS
#  endif
          {
            s = 31;
            if (msd >= 0x10000)
              {
                msd = msd >> 16;
                s -= 16;
              }
            if (msd >= 0x100)
              {
                msd = msd >> 8;
                s -= 8;
              }
            if (msd >= 0x10)
              {
                msd = msd >> 4;
                s -= 4;
              }
            if (msd >= 0x4)
              {
                msd = msd >> 2;
                s -= 2;
              }
            if (msd >= 0x2)
              {
                msd = msd >> 1;
                s -= 1;
              }
          }
# endif
      }
      /* 0 <= s < GMP_LIMB_BITS.
         Copy b, shifting it left by s bits.  */
      if (s > 0)
        {
          tmp_roomptr = (mp_limb_t *) malloc (b_len * sizeof (mp_limb_t));
          if (tmp_roomptr == NULL)
            {
              free (roomptr);
              return NULL;
            }
          {
            const mp_limb_t *sourceptr = b_ptr;
            mp_limb_t *destptr = tmp_roomptr;
            mp_twolimb_t accu = 0;
            size_t count;
            for (count = b_len; count > 0; count--)
              {
                accu += (mp_twolimb_t) *sourceptr++ << s;
                *destptr++ = (mp_limb_t) accu;
                accu = accu >> GMP_LIMB_BITS;
              }
            /* accu must be zero, since that was how s was determined.  */
            if (accu != 0)
              abort ();
          }
          b_ptr = tmp_roomptr;
        }
      /* Copy a, shifting it left by s bits, yields r.
         Memory layout:
         At the beginning: r = roomptr[0..a_len],
         at the end: r = roomptr[0..b_len-1], q = roomptr[b_len..a_len]  */
      r_ptr = roomptr;
      if (s == 0)
        {
          memcpy (r_ptr, a_ptr, a_len * sizeof (mp_limb_t));
          r_ptr[a_len] = 0;
        }
      else
        {
          const mp_limb_t *sourceptr = a_ptr;
          mp_limb_t *destptr = r_ptr;
          mp_twolimb_t accu = 0;
          size_t count;
          for (count = a_len; count > 0; count--)
            {
              accu += (mp_twolimb_t) *sourceptr++ << s;
              *destptr++ = (mp_limb_t) accu;
              accu = accu >> GMP_LIMB_BITS;
            }
          *destptr++ = (mp_limb_t) accu;
        }
      q_ptr = roomptr + b_len;
      q_len = a_len - b_len + 1; /* q will have m-n+1 limbs */
      {
        size_t j = a_len - b_len; /* m-n */
        mp_limb_t b_msd = b_ptr[b_len - 1]; /* b[n-1] */
        mp_limb_t b_2msd = b_ptr[b_len - 2]; /* b[n-2] */
        mp_twolimb_t b_msdd = /* b[n-1]*beta+b[n-2] */
          ((mp_twolimb_t) b_msd << GMP_LIMB_BITS) | b_2msd;
        /* Division loop, traversed m-n+1 times.
           j counts down, b is unchanged, beta/2 <= b[n-1] < beta.  */
        for (;;)
          {
            mp_limb_t q_star;
            mp_limb_t c1;
            if (r_ptr[j + b_len] < b_msd) /* r[j+n] < b[n-1] ? */
              {
                /* Divide r[j+n]*beta+r[j+n-1] by b[n-1], no overflow.  */
                mp_twolimb_t num =
                  ((mp_twolimb_t) r_ptr[j + b_len] << GMP_LIMB_BITS)
                  | r_ptr[j + b_len - 1];
                q_star = num / b_msd;
                c1 = num % b_msd;
              }
            else
              {
                /* Overflow, hence r[j+n]*beta+r[j+n-1] >= beta*b[n-1].  */
                q_star = (mp_limb_t)~(mp_limb_t)0; /* q* = beta-1 */
                /* Test whether r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] >= beta
                   <==> r[j+n]*beta+r[j+n-1] + b[n-1] >= beta*b[n-1]+beta
                   <==> b[n-1] < floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta)
                        {<= beta !}.
                   If yes, jump directly to the subtraction loop.
                   (Otherwise, r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] < beta
                    <==> floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta) = b[n-1] ) */
                if (r_ptr[j + b_len] > b_msd
                    || (c1 = r_ptr[j + b_len - 1] + b_msd) < b_msd)
                  /* r[j+n] >= b[n-1]+1 or
                     r[j+n] = b[n-1] and the addition r[j+n-1]+b[n-1] gives a
                     carry.  */
                  goto subtract;
              }
            /* q_star = q*,
               c1 = (r[j+n]*beta+r[j+n-1]) - q* * b[n-1] (>=0, <beta).  */
            {
              mp_twolimb_t c2 = /* c1*beta+r[j+n-2] */
                ((mp_twolimb_t) c1 << GMP_LIMB_BITS) | r_ptr[j + b_len - 2];
              mp_twolimb_t c3 = /* b[n-2] * q* */
                (mp_twolimb_t) b_2msd * (mp_twolimb_t) q_star;
              /* While c2 < c3, increase c2 and decrease c3.
                 Consider c3-c2.  While it is > 0, decrease it by
                 b[n-1]*beta+b[n-2].  Because of b[n-1]*beta+b[n-2] >= beta^2/2
                 this can happen only twice.  */
              if (c3 > c2)
                {
                  q_star = q_star - 1; /* q* := q* - 1 */
                  if (c3 - c2 > b_msdd)
                    q_star = q_star - 1; /* q* := q* - 1 */
                }
            }
            if (q_star > 0)
              subtract:
              {
                /* Subtract r := r - b * q* * beta^j.  */
                mp_limb_t cr;
                {
                  const mp_limb_t *sourceptr = b_ptr;
                  mp_limb_t *destptr = r_ptr + j;
                  mp_twolimb_t carry = 0;
                  size_t count;
                  for (count = b_len; count > 0; count--)
                    {
                      /* Here 0 <= carry <= q*.  */
                      carry =
                        carry
                        + (mp_twolimb_t) q_star * (mp_twolimb_t) *sourceptr++
                        + (mp_limb_t) ~(*destptr);
                      /* Here 0 <= carry <= beta*q* + beta-1.  */
                      *destptr++ = ~(mp_limb_t) carry;
                      carry = carry >> GMP_LIMB_BITS; /* <= q* */
                    }
                  cr = (mp_limb_t) carry;
                }
                /* Subtract cr from r_ptr[j + b_len], then forget about
                   r_ptr[j + b_len].  */
                if (cr > r_ptr[j + b_len])
                  {
                    /* Subtraction gave a carry.  */
                    q_star = q_star - 1; /* q* := q* - 1 */
                    /* Add b back.  */
                    {
                      const mp_limb_t *sourceptr = b_ptr;
                      mp_limb_t *destptr = r_ptr + j;
                      mp_limb_t carry = 0;
                      size_t count;
                      for (count = b_len; count > 0; count--)
                        {
                          mp_limb_t source1 = *sourceptr++;
                          mp_limb_t source2 = *destptr;
                          *destptr++ = source1 + source2 + carry;
                          carry =
                            (carry
                             ? source1 >= (mp_limb_t) ~source2
                             : source1 > (mp_limb_t) ~source2);
                        }
                    }
                    /* Forget about the carry and about r[j+n].  */
                  }
              }
            /* q* is determined.  Store it as q[j].  */
            q_ptr[j] = q_star;
            if (j == 0)
              break;
            j--;
          }
      }
      r_len = b_len;
      /* Normalise q.  */
      if (q_ptr[q_len - 1] == 0)
        q_len--;
# if 0 /* Not needed here, since we need r only to compare it with b/2, and
          b is shifted left by s bits.  */
      /* Shift r right by s bits.  */
      if (s > 0)
        {
          mp_limb_t ptr = r_ptr + r_len;
          mp_twolimb_t accu = 0;
          size_t count;
          for (count = r_len; count > 0; count--)
            {
              accu = (mp_twolimb_t) (mp_limb_t) accu << GMP_LIMB_BITS;
              accu += (mp_twolimb_t) *--ptr << (GMP_LIMB_BITS - s);
              *ptr = (mp_limb_t) (accu >> GMP_LIMB_BITS);
            }
        }
# endif
      /* Normalise r.  */
      while (r_len > 0 && r_ptr[r_len - 1] == 0)
        r_len--;
    }
  /* Compare r << 1 with b.  */
  if (r_len > b_len)
    goto increment_q;
  {
    size_t i;
    for (i = b_len;;)
      {
        mp_limb_t r_i =
          (i <= r_len && i > 0 ? r_ptr[i - 1] >> (GMP_LIMB_BITS - 1) : 0)
          | (i < r_len ? r_ptr[i] << 1 : 0);
        mp_limb_t b_i = (i < b_len ? b_ptr[i] : 0);
        if (r_i > b_i)
          goto increment_q;
        if (r_i < b_i)
          goto keep_q;
        if (i == 0)
          break;
        i--;
      }
  }
  if (q_len > 0 && ((q_ptr[0] & 1) != 0))
    /* q is odd.  */
    increment_q:
    {
      size_t i;
      for (i = 0; i < q_len; i++)
        if (++(q_ptr[i]) != 0)
          goto keep_q;
      q_ptr[q_len++] = 1;
    }
  keep_q:
  if (tmp_roomptr != NULL)
    free (tmp_roomptr);
  q->limbs = q_ptr;
  q->nlimbs = q_len;
  return roomptr;
}

The product performs operations on a memory buffer, but it reads from or writes to a memory location outside the buffer's intended boundary. This may result in read or write operations on unexpected memory locations that could be linked to other variables, data structures, or internal program data.