---

perfect_power.c

解説を見る。

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

inline static int
ymp_odd_root(mp_ref_t root, mp_cref_t source, size_t pow)
{
/* 
 * 2-adic Newton iteration. `Modern computer algebra' p.264-268参照
 */
  size_t clog2, k, two_exp;
  struct multiprec g, s, tmp;

  /* clog2 := */
  clog2 = ymp_abs_clog2(source->digits, source->used);
  if (clog2 == 0)
    {
      if (root) ymp_assign_digit(root, 1);
      return 1;
    }

  /* k := 2**{n*k} > <digits> なる最小のk */
  k = clog2 / pow + (clog2 % pow != 0);

  ymp_alloca_initialize(&g, k * pow + 1);
  ymp_alloca_initialize(&s, source->len + 1);
  ymp_alloca_initialize(&tmp, k * pow + 1);
  ymp_assign_digit(&g, 1);
  ymp_assign_digit(&s, 1);

  for (two_exp=2; two_exp<k; two_exp<<=1)
    {
      /* tmp := f(g_i)s_i = (g_i**pow - source)) * s_i  (mod 2**two_exp) */
      ymp_modpow_2exp_z(&tmp, &g, two_exp, pow);
      ymp_sub(&tmp, &tmp, source); ymp_mod_2exp(&tmp, &tmp, two_exp);
      ymp_modmul_2exp(&tmp, &tmp, &s, two_exp);

      /* g_{i+1} := g_i - tmp */
      ymp_sub(&g, &g, &tmp); ymp_mod_2exp(&g, &g, two_exp);

      /* tmp := f'(g_{i+1})* s_i**2 
       *      = pow * g_{i+1}**(pow-1) * s_i**2 (mod 2**two_exp)
       */
      ymp_modpow_2exp_z(&tmp, &g, two_exp, pow-1);
      ymp_mul_digit(&tmp, &tmp, pow); ymp_mod_2exp(&tmp, &tmp, two_exp);
      ymp_modmul_2exp(&tmp, &tmp, &s, two_exp); 
      ymp_modmul_2exp(&tmp, &tmp, &s, two_exp);

      /* s_{i+1} := 2*s_i - tmp */
      ymp_add(&s, &s, &s); ymp_mod_2exp(&s, &s, two_exp);
      ymp_sub(&s, &s, &tmp); ymp_mod_2exp(&s, &s, two_exp);
    }

  /* g = g_i - f(g_i)s_i  ( mod 2**k )*/
  ymp_modpow_2exp_z(&tmp, &g, k, pow);
  ymp_sub(&tmp, &tmp, source); ymp_mod_2exp(&tmp, &tmp, k);
  ymp_modmul_2exp(&tmp, &tmp, &s, k);
  ymp_sub(&g, &g, &tmp); ymp_mod_2exp(&g, &g, k);

  /* if (g < 0) g += 2**k */
  if (g.sign)
    {
      ymp_assign_2exp(&tmp, k);
      ymp_add(&g, &g, &tmp);
    }


  ymp_pow(&tmp, &g, pow);
  if (ymp_compare(&tmp, source) == 0)
    {
      if (root) ymp_assign(root, &g);
      return 1;
    }
  else
    {
      return 0;
    }
}


inline static int
ymp_nsqrt_mod2exp(mp_ref_t root, mp_cref_t source, size_t e_mod)
{
  /*
   * MATHEMATICS OF COMPUTATION Volume 67, Number 223, 
   * July 1998, Pages 1253-1283
   *
   * Daniel J. Bernstein
   * "DETECTING PERFECT POWERS IN ESSENTIALLY LINEAR TIME": Section 21
   *
   * Algorithm S2 and Lemma 21.4
   */
  
  size_t n_digits;
  size_t e;
  struct multiprec tmp, r;

  if ((source->digits[0] & 0x03) != 1) return 0;
  if (e_mod == 1) { ymp_assign_digit(root, 1); return 1;}
  if ((source->digits[0] & 0x07) != 1) return 0;
  if (e_mod == 2) { ymp_assign_digit(root, 1); return 1;}

  n_digits = e_mod / DIGIT_BIT + (e_mod % DIGIT_BIT != 0);
  ymp_alloca_initialize(&r, n_digits+2);
  ymp_alloca_initialize(&tmp, n_digits+2);

  ymp_assign_digit(&r, 1);
  for (e = 3; e < e_mod; e = 2*e-1)
    {
      /* tmp := (3r - source*r**3) mod 2**(e+1) */
      ymp_modmul_2exp(&tmp, source, &r, e+1);
      ymp_modmul_2exp(&tmp, &tmp, &r, e+1);
      ymp_sub_digit(&tmp, &tmp, 3);
      ymp_neg(&tmp, &tmp);
      ymp_modmul_2exp(&tmp, &tmp, &r, e+1);

      /* r := tmp/2 mod 2**e */
      ymp_div_2exp(&r, &tmp, 1);
    }

  /* tmp := (3r - source*r**3) mod 2**(e_mod+1) */
  ymp_modmul_2exp(&tmp, source, &r, e_mod+1);
  ymp_modmul_2exp(&tmp, &tmp, &r, e_mod+1);
  ymp_sub_digit(&tmp, &tmp, 3);
  ymp_neg(&tmp, &tmp);
  ymp_modmul_2exp(&tmp, &tmp, &r, e_mod+1);

  /* r := tmp/2 mod 2**e */
  ymp_div_2exp(root, &tmp, 1);
  return 1;
}



inline static int
ymp_odd_sqrt(mp_ref_t root, mp_cref_t source)
{
  size_t n_digits, n_bits;
  size_t clog2, b;
  struct multiprec r, tmp;

  clog2 = ymp_abs_clog2(source->digits, source->used);
  b = clog2/2 + clog2%2;

  n_digits = b / DIGIT_BIT; n_bits = b % DIGIT_BIT;
  ymp_alloca_initialize(&r, n_digits + 1);

  if (!ymp_nsqrt_mod2exp(&r, source, b)) return 0;

  ymp_modmul_2exp(&r, &r, source, b);
  
  ymp_alloca_initialize(&tmp, MAX(2*(n_digits+1), 2*(r.used+1)));

  if (r.sign)
    {
      ymp_assign_2exp(&tmp, b);
      ymp_add(&r, &tmp, &r);
    }
  
  ymp_mul(&tmp, &r, &r);
  if (ymp_compare(source, &tmp) == 0)
    {
      if (root) ymp_assign(root, &r);
      return 1;
    }

  ymp_assign_2exp(&tmp, b);
  ymp_sub(&tmp, &tmp, &r);
  ymp_mul(&tmp, &tmp, &tmp);
  if (ymp_compare(source, &tmp) == 0)
    {
      if (root) ymp_assign(root, &r);
      return 1;
    }

  return 0;
}


int
ymp_proot(mp_ref_t result, mp_cref_t self, size_t pow)
{
  struct multiprec tmp;
  size_t n_square;            /* powの素因数分解の中の2の個数 */
  size_t e2_digits, e2_bits;  /* 下位の0の並んだ部分の、digit数, bit数 */
  register digit_t ld;



  /* 正数以外は除く */
  if (self->sign || ymp_is_zero(self)) return 0;



  /* e2 := {2**e が <self> を割る最大の e} */
  for (e2_digits=0; e2_digits<self->used; ++e2_digits)
    {
      if (self->digits[e2_digits] != 0) break;
    }
  e2_bits = 0;
  for (ld = self->digits[e2_digits]; (ld & 0x1) == 0; ld>>=1)
    {
      ++e2_bits;
    }

  /* 素因数分解の中の2の冪指数をチェック */
  if ((e2_digits * DIGIT_BIT + e2_bits) % pow != 0) return 0;

  /* 以下、2の冪は考えなくて良い */


  n_square = 0;
  while ( (pow & 0x01) == 0 ) { pow >>= 1; ++n_square; }


  ymp_alloca_initialize(&tmp, self->used - e2_digits);

  ymp_div_2exp(&tmp, self, e2_digits * DIGIT_BIT + e2_bits);
  if (!ymp_odd_root(&tmp, &tmp, pow))
    {
      return 0;
    }


  /* 以下、冪指数のうち2冪の部分のみ調べる */
  while (n_square-- > 0)
    {
      if (!ymp_odd_sqrt(&tmp, &tmp)) return 0;
    }
  
  ymp_assign(result, &tmp);
  return 1;
}



int
ymp_is_perfect_power(mp_cref_t self)
{
  struct multiprec tmp;
  size_t e2_digits, e2_bits, e2;  /* 下位の0の並んだ部分の、digit数, bit数 */
  size_t pow1, pow2, max_pow;
  register digit_t ld;

  /* 正数以外は除く */
  if (self->sign || ymp_is_zero(self)) return 0;

  /* e2 := {2**e が <self> を割る最大の e} */
  for (e2_digits=0; e2_digits<self->used; ++e2_digits)
    {
      if (self->digits[e2_digits] != 0) break;
    }
  e2_bits = 0;
  for (ld = self->digits[e2_digits]; (ld & 0x1) == 0; ld>>=1)
    {
      ++e2_bits;
    }
  e2 = e2_digits * DIGIT_BIT + e2_bits;

  
  ymp_alloca_initialize(&tmp, self->used - e2_digits);
  ymp_div_2exp(&tmp, self, e2);
  
  /* 完全平方か? */
  if (e2 % 2 == 0 && ymp_odd_sqrt(NULL, &tmp)) return 1;

  /* 完全立方か? */
  if (e2 % 3 == 0 && ymp_odd_root(NULL, &tmp, 3)) return 1;

  /* 完全奇素数乗か ? */
  max_pow = ymp_abs_clog2(tmp.digits, tmp.used);
  for (pow1=1, pow2=5; 1; pow1 += 6, pow2 += 6)
    {
      if (pow1 > max_pow) break;
      if (e2 % pow1 == 0 && ymp_odd_root(NULL, &tmp, pow1)) return 1;
      if (pow2 > max_pow) break;
      if (e2 % pow2 == 0 && ymp_odd_root(NULL, &tmp, pow2)) return 1;
    }

  return 0;
}

---