chiark - git - mdw - catacomb/blob - math/mp-nthrt.c

   1 /* -*-c-*-
   2  *
   3  * Calculating %$n$%th roots
   4  *
   5  * (c) 2020 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software: you can redistribute it and/or modify it
  13  * under the terms of the GNU Library General Public License as published
  14  * by the Free Software Foundation; either version 2 of the License, or
  15  * (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful, but
  18  * WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb.  If not, write to the Free Software
  24  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
  25  * USA.
  26  */
  27
  28 /*----- Header files ------------------------------------------------------*/
  29
  30 #include "mp.h"
  31 #include "mpint.h"
  32 #include "mplimits.h"
  33 #include "primeiter.h"
  34
  35 /*----- Main code ---------------------------------------------------------*/
  36
  37 /* --- @mp_nthrt@ --- *
  38  *
  39  * Arguments:   @mp *d@ = fake destination
  40  *              @mp *a@ = an integer
  41  *              @mp *n@ = a strictly positive integer
  42  *              @int *exectp_out@ = set nonzero if an exact solution is found
  43  *
  44  * Returns:     The integer %$\bigl\lfloor \sqrt[n]{a} \bigr\rfloor$%.
  45  *
  46  * Use:         Return an approximation to the %$n$%th root of %$a$%.
  47  *              Specifically, it returns the largest integer %$x$% such that
  48  *              %$x^n \le a$%.  If %$x^n = a$% then @*exactp_out@ is set
  49  *              nonzero; otherwise it is set zero.  (If @exactp_out@ is null
  50  *              then this information is discarded.)
  51  *
  52  *              The exponent %$n$% must be strictly positive: it's not clear
  53  *              to me what the right answer is for %$n \le 0$%.  If %$a$% is
  54  *              negative then %$n$% must be odd; otherwise there is no real
  55  *              solution.
  56  */
  57
  58 mp *mp_nthrt(mp *d, mp *a, mp *n, int *exactp_out)
  59 {
  60   /* We want to find %$x$% such that %$x^n = a$%.  Newton--Raphson says we
  61    * start with a guess %$x_i$%, and refine it to a better guess %$x_{i+1}$%
  62    * which is where the tangent to the curve %$y = x^n - a$% at %$x = x_i$%
  63    * meets the %$x$%-axis.  The tangent meets the curve at
  64    * %($x_i, x_i^n - a)$%, and has gradient %$n x_i^{n-1}$%, and hence meets
  65    * the %$x$%-axis at %$x_{i+1} = x_i - (x_i^n - a)/(n x_i^{n-1})$%.
  66    *
  67    * We're working with plain integers here, and we actually want
  68    * %$\lfloor x \rfloor$%.  We can check that we're close enough because
  69    * %$x^n$% is monotonic for positive %$x$%: if %$x_i^n > a$% then we're too
  70    * large; if %$(xi + i)^n \le a$% then we're too small; otherwise we're
  71    * done.
  72    */
  73
  74   mp *ai = MP_NEW, *bi = MP_NEW, *xi = MP_NEW,
  75     *nm1 = MP_NEW, *delta = MP_NEW;
  76   int cmp;
  77   unsigned f = 0;
  78 #define f_neg 1u
  79
  80   /* Firstly, deal with a zero or negative %$xa%. */
  81   if (!MP_NEGP(a))
  82     MP_COPY(a);
  83   else {
  84     assert(MP_ODDP(n));
  85     f |= f_neg;
  86     a = mp_neg(MP_NEW, a);
  87   }
  88
  89   /* Secondly, reject %$n \le 0$%.
  90    *
  91    * Clearly %$x^0 = 1$% for nonzero %$x$%, so if %$a$% is zero then we could
  92    * return anything, but maybe %$1$% for concreteness; but what do we return
  93    * for %$a \ne 1$%?  For %$n < 0$% there's just no hope.
  94    */
  95   assert(MP_POSP(n));
  96
  97   /* Pick a starting point.  This is rather important to get right.  In
  98    * particular, for large %$n$%, if our initial guess is too small, then the
  99    * next iteration is a wild overestimate and it takes a long time to
 100    * converge back.
 101    */
 102   nm1 = mp_sub(nm1, n, MP_ONE);
 103   xi = mp_fromulong(xi, mp_bits(a));
 104   xi = mp_add(xi, xi, nm1); mp_div(&xi, 0, xi, n);
 105   assert(MP_CMP(xi, <=, MP_ULONG_MAX));
 106   xi = mp_setbit(xi, MP_ZERO, mp_toulong(xi));
 107
 108   /* The main iteration. */
 109   for (;;) {
 110     ai = mp_exp(ai, xi, n);
 111     bi = mp_add(bi, xi, MP_ONE); bi = mp_exp(bi, bi, n);
 112     cmp = mp_cmp(ai, a);
 113     if (cmp <= 0 && MP_CMP(a, <, bi)) break;
 114     ai = mp_sub(ai, ai, a);
 115     bi = mp_exp(bi, xi, nm1); bi = mp_mul(bi, bi, n);
 116     mp_div(&delta, 0, ai, bi);
 117     if (MP_ZEROP(delta) && cmp > 0) xi = mp_sub(xi, xi, MP_ONE);
 118     else xi = mp_sub(xi, xi, delta);
 119   }
 120
 121   /* Final sorting out of negative inputs.
 122    *
 123    * If the input %$a$% is not an exact %$n$%th root, then we must round
 124    * %%\emph{up}%% so that, after negation, we implement the correct floor
 125    * behaviour.
 126    */
 127   if (f&f_neg) {
 128     if (cmp) xi = mp_add(xi, xi, MP_ONE);
 129     xi = mp_neg(xi, xi);
 130   }
 131
 132   /* Free up the temporary things. */
 133   MP_DROP(a); MP_DROP(ai); MP_DROP(bi);
 134   MP_DROP(nm1); MP_DROP(delta); mp_drop(d);
 135
 136   /* Done. */
 137   if (exactp_out) *exactp_out = (cmp == 0);
 138   return (xi);
 139
 140 #undef f_neg
 141 }
 142
 143 /* --- @mp_perfect_power_p@ --- *
 144  *
 145  * Arguments:   @mp **x@ = where to write the base
 146  *              @mp **n@ = where to write the exponent
 147  *              @mp *a@ = an integer
 148  *
 149  * Returns:     Nonzero if %$a$% is a perfect power.
 150  *
 151  * Use:         Returns whether an integer %$a$% is a perfect power, i.e.,
 152  *              whether it can be written in the form %$a = x^n$% where
 153  *              %$|x| > 1$% and %$n > 1$% are integers.  If this is possible,
 154  *              then (a) store %$x$% and the largest such %$n$% in @*x@ and
 155  *              @*n@, and return nonzero; otherwise, store %$x = a$% and
 156  *              %$n = 1$% and return zero.  (Either @x@ or @n@, or both, may
 157  *              be null to discard these outputs.)
 158  *
 159  *              Note that %$-1$%, %$0$% and %$1$% are not considered perfect
 160  *              powers by this definition.  (The exponent is not well-defined
 161  *              in these cases, but it seemed better to implement a function
 162  *              which worked for all integers.)  Note also that %$-4$% is not
 163  *              a perfect power since it has no real square root.
 164  */
 165
 166 int mp_perfect_power_p(mp **x, mp **n, mp *a)
 167 {
 168   mp *r = MP_ONE, *p = MP_NEW, *t = MP_NEW;
 169   int exactp, rc = 0;
 170   primeiter pi;
 171   unsigned f = 0;
 172 #define f_neg 1u
 173
 174   /* Initial setup.  We'll modify @a@ as we go, so we have to destroy it when
 175    * we're done; therefore we should take a copy now.
 176    */
 177   MP_COPY(a);
 178
 179   /* According to the definition above, @0@ is %%\emph{not}%% a perfect
 180    * power.  Dispose of this special case early.
 181    */
 182   if (MP_ZEROP(a)) goto done;
 183
 184   /* Start by initializing a prime iterator.  If %$a$% is negative, then it
 185    * can't be an even power so start iterating at 3; otherwise we must start
 186    * with 2.
 187    */
 188   primeiter_create(&pi, MP_NEGP(a) ? MP_THREE : MP_TWO);
 189   if (MP_NEGP(a)) { a = mp_neg(a, a); f |= f_neg; }
 190
 191   /* Prime the pump for the main loop. */
 192   p = primeiter_next(&pi, p);
 193
 194   /* Here we go. */
 195   for (;;) {
 196
 197     /* See if %$a$% is a perfect %$p$%th power. */
 198     t = mp_nthrt(t, a, p, &exactp);
 199     if (MP_EQ(t, MP_ONE))
 200       /* No, and moreover %$2^p > a$%.  We won't get anywhere by examining
 201        * larger primes, so stop here.
 202        */
 203       break;
 204
 205     else if (!exactp)
 206       /* No.  Oh, well: let's try the next prime, then. */
 207
 208       p = primeiter_next(&pi, p);
 209     else {
 210       /* Yes!  We've found that %$a = a'^p$% for some %$a'$%.  We know
 211        * (induction) that %$a'$% is not a perfect %$q$%th power for any prime
 212        * %$q < p$%, so we can write %$a = x^n$% if and only if we can write
 213        * %$a' = x'^{n'}$%, in which case we have %$a = a'^p = x'^{p n'}$%.
 214        */
 215       r = mp_mul(r, r, p);
 216       MP_DROP(a); a = t; t = MP_NEW;
 217       rc = 1;
 218     }
 219   }
 220
 221   /* We're all done. */
 222   primeiter_destroy(&pi);
 223
 224 done:
 225   /* Return the resulting decomposition. */
 226   if (x) { mp_drop(*x); *x = f&f_neg ? mp_neg(a, a) : a; a = 0; }
 227   if (n) { mp_drop(*n); *n = r; r = 0; }
 228
 229   /* Clean everything up. */
 230   mp_drop(a); mp_drop(p); mp_drop(r); mp_drop(t);
 231
 232   /* All done. */
 233   return (rc);
 234
 235 #undef f_neg
 236 }
 237
 238 /*----- Test rig ----------------------------------------------------------*/
 239
 240 #ifdef TEST_RIG
 241
 242 #include <mLib/testrig.h>
 243
 244 static int vrf_nthrt(dstr *v)
 245 {
 246   mp *a = *(mp **)v[0].buf;
 247   mp *n = *(mp **)v[1].buf;
 248   mp *r0 = *(mp **)v[2].buf;
 249   int exactp0 = *(int *)v[3].buf, exactp;
 250   mp *r = mp_nthrt(MP_NEW, a, n, &exactp);
 251   int ok = 1;
 252
 253   if (!MP_EQ(r, r0) || exactp != exactp0) {
 254     ok = 0;
 255     fputs("\n*** nthrt failed", stderr);
 256     fputs("\n***       a = ", stderr); mp_writefile(a, stderr, 10);
 257     fputs("\n***       n = ", stderr); mp_writefile(n, stderr, 10);
 258     fputs("\n***  r calc = ", stderr); mp_writefile(r, stderr, 10);
 259     fprintf(stderr, "\n*** ex calc = %d", exactp);
 260     fputs("\n***   r exp = ", stderr); mp_writefile(r0, stderr, 10);
 261     fprintf(stderr, "\n***  ex exp = %d", exactp0);
 262     fputc('\n', stderr);
 263   }
 264
 265   mp_drop(a); mp_drop(n); mp_drop(r); mp_drop(r0);
 266   assert(mparena_count(MPARENA_GLOBAL) == 0);
 267
 268   return (ok);
 269 }
 270
 271 static int vrf_ppp(dstr *v)
 272 {
 273   mp *a = *(mp **)v[0].buf;
 274   int ret0 = *(int *)v[1].buf;
 275   mp *x0 = *(mp **)v[2].buf;
 276   mp *n0 = *(mp **)v[3].buf;
 277   mp *x = MP_NEW, *n = MP_NEW;
 278   int ret = mp_perfect_power_p(&x, &n, a);
 279   int ok = 1;
 280
 281   if (ret != ret0 || !MP_EQ(x, x0) || !MP_EQ(n, n0)) {
 282     ok = 0;
 283     fputs("\n*** perfect_power_p failed", stderr);
 284     fputs("\n***       a = ", stderr); mp_writefile(a, stderr, 10);
 285     fprintf(stderr, "\n***  r calc = %d", ret);
 286     fputs("\n***  x calc = ", stderr); mp_writefile(x, stderr, 10);
 287     fputs("\n***  n calc = ", stderr); mp_writefile(n, stderr, 10);
 288     fprintf(stderr, "\n***   r exp = %d", ret0);
 289     fputs("\n***   x exp = ", stderr); mp_writefile(x0, stderr, 10);
 290     fputs("\n***   n exp = ", stderr); mp_writefile(n0, stderr, 10);
 291     fputc('\n', stderr);
 292   }
 293
 294   mp_drop(a); mp_drop(x); mp_drop(n); mp_drop(x0); mp_drop(n0);
 295   assert(mparena_count(MPARENA_GLOBAL) == 0);
 296
 297   return (ok);
 298 }
 299
 300 static test_chunk tests[] = {
 301   { "nthrt", vrf_nthrt, { &type_mp, &type_mp, &type_mp, &type_int, 0 } },
 302   { "perfect-power-p", vrf_ppp, { &type_mp, &type_int, &type_mp, &type_mp, 0 } },
 303   { 0, 0, { 0 } },
 304 };
 305
 306 int main(int argc, char *argv[])
 307 {
 308   sub_init();
 309   test_run(argc, argv, tests, SRCDIR "/t/mp");
 310   return (0);
 311 }
 312
 313 #endif
 314
 315 /*----- That's all, folks -------------------------------------------------*/