---

euclid.c

解説を見る。

#include "multiprec.h"
#include <string.h>

void
ymp_euclid(mp_ref_t g, mp_cref_t a, mp_cref_t b)
{
  struct {
    size_t used;
    digit_t *digits;
  } r0, r1;

  r0.used = a->used; r1.used = b->used;
  YMP_TEMP_ALLOCATE(r0.digits, digit_t, r0.used);
  YMP_TEMP_ALLOCATE(r1.digits, digit_t, r1.used);
  memcpy(r0.digits, a->digits, r0.used*sizeof(digit_t));
  memcpy(r1.digits, b->digits, r1.used*sizeof(digit_t));

  while (r1.used !=0 && (r1.used != 1 || r1.digits[0] != 0))
    {
      register size_t tmp_used;
      register digit_t *tmp_digits;
      if (r0.used >= r1.used)
        {
          r0.used = ymp_divmodabs(NULL, r0.digits, 
                                  r0.digits, r1.digits, r0.used, r1.used);
        }
      tmp_used = r0.used; r0.used = r1.used; r1.used = tmp_used;
      tmp_digits = r0.digits; r0.digits = r1.digits; r1.digits = tmp_digits;
    }
  ymp_assign_array(g, r0.used, r0.digits);
}


void
ymp_euclid_ex(mp_ref_t g, mp_ref_t s, mp_ref_t t, mp_cref_t a, mp_cref_t b)
{
  struct multiprec r0, r1, s0, s1, t0, t1, q, tmp;

  ymp_alloca_initialize(&r0, a->used);
  ymp_alloca_initialize(&r1, b->used);
  ymp_alloca_initialize(&q, MAX(a->used, b->used));
  ymp_assign(&r0, a); ymp_assign(&r1, b);
  
  ymp_initialize_by_digit(&s0, 1); ymp_initialize_by_digit(&s1, 0);
  ymp_initialize_by_digit(&t0, 0); ymp_initialize_by_digit(&t1, 1);
  ymp_initialize(&tmp);

  while (!ymp_is_zero(&r1))
    {
      ymp_divmod(&q, &r0, &r0, &r1);
      ymp_mul(&tmp, &q, &s1);
      ymp_sub(&s0, &s0, &tmp);
      ymp_mul(&tmp, &q, &t1);
      ymp_sub(&t0, &t0, &tmp);

      ymp_swap(&r0, &r1);
      ymp_swap(&s0, &s1);
      ymp_swap(&t0, &t1);
    }

  if (g) ymp_assign(g, &r0);
  if (s) ymp_assign(s, &s0);
  if (t) ymp_assign(t, &t0);
  ymp_destroy(&s0); ymp_destroy(&s1);
  ymp_destroy(&t0); ymp_destroy(&t1);
  ymp_destroy(&tmp);
}


int
ymp_modinv(mp_ref_t result, mp_cref_t self, mp_cref_t modulus)
{
  struct multiprec r0, r1, s0, s1, q, tmp;
  register size_t len = MAX(self->used, modulus->used);

  ymp_alloca_initialize(&r0, self->used);          /* |r0| <= |self|        */
  ymp_alloca_initialize(&r1, modulus->used);       /* |r1| <= |modulus|     */
  ymp_alloca_initialize(&s0, modulus->used+len+1); /* |s0| <= |r0|*|q|+|s1| */
  ymp_alloca_initialize(&s1, modulus->used+len+1); /* |s1| <= |r1|*|q|+|s0| */
  ymp_alloca_initialize(&q, len);                  /* |q|  <= max(|r0|,|r1|)*/
  ymp_alloca_initialize(&tmp, modulus->used+len);  /* |tmp| <= |r_i|*|q|    */

  ymp_assign(&r0, self); ymp_assign(&r1, modulus);
  ymp_assign_digit(&s0, 1); ymp_assign_digit(&s1, 0);

  while (!ymp_is_zero(&r1))
    {
      ymp_divmod(&q, &r0, &r0, &r1);
      ymp_mul(&tmp, &s1, &q);
      ymp_sub(&s0, &s0, &tmp);
      ymp_mod(&s0, &s0, modulus);

      ymp_swap(&r0, &r1);
      ymp_swap(&s0, &s1);
    }
  
  if (r0.used == 1 && r0.digits[0] == 1)
    {
      ymp_assign(result, &s0);
      if (r0.sign) result->sign = !result->sign;
      return 0;
    }
  else
    {
      return 1;
    }
}


int
ymp_binary_sunzi(mp_ref_t result, 
                 const mp_cref_t remainders[2], const mp_cref_t moduli[2])
{
  /*
   * H. L. Garnerの高速解法
   * c_{0,1} := m_0^{-1} (mod m_1)
   *
   * v_0 := u_0 mod m_0
   * v_1 := (u_1 - u_0) * c_{0,1} mod m_1
   *
   * u = v_1 * m_0 + v_0
   */

  struct multiprec m_0_inv; /* m_0^{-1} (mod m_1) */
  struct multiprec v[2];
  register size_t len;

  ymp_alloca_initialize(&m_0_inv, moduli[1]->used);
  if (ymp_modinv(&m_0_inv, moduli[0], moduli[1])) return 1;

  ymp_alloca_initialize(&v[0], moduli[0]->used);
  ymp_mod(&v[0], remainders[0], moduli[0]);

  len = MAX(remainders[1]->used, moduli[0]->used + moduli[1]->used);
  ymp_alloca_initialize(&v[1], len);

  ymp_sub(&v[1], remainders[1], &v[0]);
  ymp_modmul(&v[1], &v[1], &m_0_inv, moduli[1]);
  ymp_mul(&v[1], &v[1], moduli[0]);

  ymp_add(result, &v[1], &v[0]);
  return 0;
}

---