Алгоритмы НОД для арифметики произвольной точности

Я полностью застрял в этом вопросе, поэтому я ищу любую помощь.

Я думаю, что все знают об основных алгоритмах вычисления НОД, таких как двоичный или евклидов НОД. Не проблема реализовать такой метод для вычисления двух чисел с одинарной точностью. На самом деле, это всего лишь пара ударов.

Мне нужно реализовать этот метод (на языке C) для чисел с множественной точностью (более 10 ^ 5 бит). Доступно несколько библиотек GNU (GNU MP, MPFR, MPIR), и у них есть средства для определения чисел с множественной точностью и выполнения над ними действий. Это выглядит как одно число с многократной точностью, хранящееся в памяти как пара частей с одинарной точностью, также называемых «конечностями».

В них реализованы некоторые методы нахождения gcd(a, b), но их, на самом деле, сложно использовать для моих нужд. Двоичный метод вычисления НОД используется только в том случае, если a и b содержат ровно две конечности. Метод HGCD используется, когда min(a,b) содержит более (т.е. 630) конечностей и т. д. Мне трудно понять, как любой из этих методов можно расширить для использования с любой длиной a и b. Я также обнаружил, что разные версии библиотек GNU содержат разные версии и методы алгоритмов GCD.

Вопрос: я хочу выяснить, можно ли заставить двоичный алгоритм НОД работать с целыми числами многократной точности любой длины в терминах "конечностей", и если это возможно - получить какую-либо помощь или идеи как это реализовать на C. Есть ли у кого-нибудь идеи или части кода, как это реализовать?

Я хотел бы рассмотреть любые советы или любое другое решение для решения этой проблемы.

Вот часть бинарного GCD-метода GNU MP для (a = b = 2 конечностей), если кто-нибудь взглянет:

/* Use binary algorithm to compute G <-- GCD (U, V) for usize, vsize == 2.
   Both U and V must be odd. */
static inline mp_size_t
gcd_2 (mp_ptr gp, mp_srcptr up, mp_srcptr vp)
{
  printf("gcd_2 invoked\n");
  mp_limb_t u0, u1, v0, v1;
  mp_size_t gn;

  u0 = up[0];
  u1 = up[1];
  v0 = vp[0];
  v1 = vp[1];

  ASSERT (u0 & 1);
  ASSERT (v0 & 1);

  /* Check for u0 != v0 needed to ensure that argument to
   * count_trailing_zeros is non-zero. */
  while (u1 != v1 && u0 != v0)
    {
      unsigned long int r;
      if (u1 > v1)
  {
    sub_ddmmss (u1, u0, u1, u0, v1, v0);
    count_trailing_zeros (r, u0);
    u0 = ((u1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (u0 >> r);
    u1 >>= r;
  }
      else  /* u1 < v1.  */
  {
    sub_ddmmss (v1, v0, v1, v0, u1, u0);
    count_trailing_zeros (r, v0);
    v0 = ((v1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (v0 >> r);
    v1 >>= r;
  }
    }

  gp[0] = u0, gp[1] = u1, gn = 1 + (u1 != 0);

  /* If U == V == GCD, done.  Otherwise, compute GCD (V, |U - V|).  */
  if (u1 == v1 && u0 == v0)
    return gn;

  v0 = (u0 == v0) ? ((u1 > v1) ? u1-v1 : v1-u1) : ((u0 > v0) ? u0-v0 : v0-u0);
  gp[0] = mpn_gcd_1 (gp, gn, v0);

  return 1;
}

CodePaste вышеперечисленного.