From: Mark Wooding <mdw@distorted.org.uk>
Date: Wed, 14 Oct 2020 02:03:41 +0000 (+0100)
Subject: math/mp-nthrt.c: Add commentary for `mp_perfect_power_p'.
X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/commitdiff_plain/36109084f51f185bb40e28248e07d66ae91346ba

math/mp-nthrt.c: Add commentary for `mp_perfect_power_p'.

This is quite simple, really, but it doesn't hurt to explain what's
going on.
---

diff --git a/math/mp-nthrt.c b/math/mp-nthrt.c
index 77e6ebb0..b253d2fe 100644
--- a/math/mp-nthrt.c
+++ b/math/mp-nthrt.c
@@ -171,28 +171,65 @@ int mp_perfect_power_p(mp **x, mp **n, mp *a)
   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