unsigned f = 0;
#define f_neg 1u
+ /* Initial setup. We'll modify @a@ as we go, so we have to destroy it when
+ * we're done; therefore we should take a copy now.
+ */
MP_COPY(a);
+
+ /* According to the definition above, @0@ is %%\emph{not}%% a perfect
+ * power. Dispose of this special case early.
+ */
if (MP_ZEROP(a)) goto done;
+
+ /* Start by initializing a prime iterator. If %$a$% is negative, then it
+ * can't be an even power so start iterating at 3; otherwise we must start
+ * with 2.
+ */
primeiter_create(&pi, MP_NEGP(a) ? MP_THREE : MP_TWO);
if (MP_NEGP(a)) { a = mp_neg(a, a); f |= f_neg; }
+
+ /* Prime the pump for the main loop. */
p = primeiter_next(&pi, p);
+
+ /* Here we go. */
for (;;) {
+
+ /* See if %$a$% is a perfect %$p$%th power. */
t = mp_nthrt(t, a, p, &exactp);
if (MP_EQ(t, MP_ONE))
+ /* No, and moreover %$2^p > a$%. We won't get anywhere by examining
+ * larger primes, so stop here.
+ */
break;
+
else if (!exactp)
+ /* No. Oh, well: let's try the next prime, then. */
+
p = primeiter_next(&pi, p);
else {
+ /* Yes! We've found that %$a = a'^p$% for some %$a'$%. We know
+ * (induction) that %$a'$% is not a perfect %$q$%th power for any prime
+ * %$q < p$%, so we can write %$a = x^n$% if and only if we can write
+ * %$a' = x'^{n'}$%, in which case we have %$a = a'^p = x'^{p n'}$%.
+ */
r = mp_mul(r, r, p);
MP_DROP(a); a = t; t = MP_NEW;
rc = 1;
}
}
+
+ /* We're all done. */
primeiter_destroy(&pi);
+
done:
+ /* Return the resulting decomposition. */
if (x) { mp_drop(*x); *x = f&f_neg ? mp_neg(a, a) : a; a = 0; }
if (n) { mp_drop(*n); *n = r; r = 0; }
+
+ /* Clean everything up. */
mp_drop(a); mp_drop(p); mp_drop(r); mp_drop(t);
+
+ /* All done. */
return (rc);
#undef f_neg