--- /dev/null
+/* -*-c-*-
+ *
+ * Calculating %$n$%th roots
+ *
+ * (c) 2020 Straylight/Edgeware
+ */
+
+/*----- Licensing notice --------------------------------------------------*
+ *
+ * This file is part of Catacomb.
+ *
+ * Catacomb is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU Library General Public License as published
+ * by the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * Catacomb is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Library General Public License for more details.
+ *
+ * You should have received a copy of the GNU Library General Public
+ * License along with Catacomb. If not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+ * USA.
+ */
+
+/*----- Header files ------------------------------------------------------*/
+
+#include "mp.h"
+#include "mpint.h"
+#include "mplimits.h"
+#include "primeiter.h"
+
+/*----- Main code ---------------------------------------------------------*/
+
+/* --- @mp_nthrt@ --- *
+ *
+ * Arguments: @mp *d@ = fake destination
+ * @mp *a@ = an integer
+ * @mp *n@ = a strictly positive integer
+ * @int *exectp_out@ = set nonzero if an exact solution is found
+ *
+ * Returns: The integer %$\bigl\lfloor \sqrt[n]{a} \bigr\rfloor$%.
+ *
+ * Use: Return an approximation to the %$n$%th root of %$a$%.
+ * Specifically, it returns the largest integer %$x$% such that
+ * %$x^n \le a$%. If %$x^n = a$% then @*exactp_out@ is set
+ * nonzero; otherwise it is set zero. (If @exactp_out@ is null
+ * then this information is discarded.)
+ *
+ * The exponent %$n$% must be strictly positive: it's not clear
+ * to me what the right answer is for %$n \le 0$%. If %$a$% is
+ * negative then %$n$% must be odd; otherwise there is no real
+ * solution.
+ */
+
+mp *mp_nthrt(mp *d, mp *a, mp *n, int *exactp_out)
+{
+ /* We want to find %$x$% such that %$x^n = a$%. Newton--Raphson says we
+ * start with a guess %$x_i$%, and refine it to a better guess %$x_{i+1}$%
+ * which is where the tangent to the curve %$y = x^n - a$% at %$x = x_i$%
+ * meets the %$x$%-axis. The tangent meets the curve at
+ * %($x_i, x_i^n - a)$%, and has gradient %$n x_i^{n-1}$%, and hence meets
+ * the %$x$%-axis at %$x_{i+1} = x_i - (x_i^n - a)/(n x_i^{n-1})$%.
+ *
+ * We're working with plain integers here, and we actually want
+ * %$\lfloor x \rfloor$%. We can check that we're close enough because
+ * %$x^n$% is monotonic for positive %$x$%: if %$x_i^n > a$% then we're too
+ * large; if %$(xi + i)^n \le a$% then we're too small; otherwise we're
+ * done.
+ */
+
+ mp *ai = MP_NEW, *bi = MP_NEW, *xi = MP_NEW,
+ *nm1 = MP_NEW, *delta = MP_NEW;
+ int cmp;
+ unsigned f = 0;
+#define f_neg 1u
+
+ /* Firstly, deal with a zero or negative %$xa%. */
+ if (!MP_NEGP(a))
+ MP_COPY(a);
+ else {
+ assert(MP_ODDP(n));
+ f |= f_neg;
+ a = mp_neg(MP_NEW, a);
+ }
+
+ /* Secondly, reject %$n \le 0$%.
+ *
+ * Clearly %$x^0 = 1$% for nonzero %$x$%, so if %$a$% is zero then we could
+ * return anything, but maybe %$1$% for concreteness; but what do we return
+ * for %$a \ne 1$%? For %$n < 0$% there's just no hope.
+ */
+ assert(MP_POSP(n));
+
+ /* Pick a starting point. This is rather important to get right. In
+ * particular, for large %$n$%, if we our initial guess too small, then the
+ * next iteration is a wild overestimate and it takes a long time to
+ * converge back.
+ */
+ nm1 = mp_sub(nm1, n, MP_ONE);
+ xi = mp_fromulong(xi, mp_bits(a));
+ xi = mp_add(xi, xi, nm1); mp_div(&xi, 0, xi, n);
+ assert(MP_CMP(xi, <=, MP_ULONG_MAX));
+ xi = mp_setbit(xi, MP_ZERO, mp_toulong(xi));
+
+ /* The main iteration. */
+ for (;;) {
+ ai = mp_exp(ai, xi, n);
+ bi = mp_add(bi, xi, MP_ONE); bi = mp_exp(bi, bi, n);
+ cmp = mp_cmp(ai, a);
+ if (cmp <= 0 && MP_CMP(a, <, bi)) break;
+ ai = mp_sub(ai, ai, a);
+ bi = mp_exp(bi, xi, nm1); bi = mp_mul(bi, bi, n);
+ mp_div(&delta, 0, ai, bi);
+ if (MP_ZEROP(delta) && cmp > 0) xi = mp_sub(xi, xi, MP_ONE);
+ else xi = mp_sub(xi, xi, delta);
+ }
+
+ /* Final sorting out of negative inputs.
+ *
+ * If the input %$a$% is not an exact %$n$%th root, then we must round
+ * %%\emph{up}%% so that, after negation, we implement the correct floor
+ * behaviour.
+ */
+ if (f&f_neg) {
+ if (cmp) xi = mp_add(xi, xi, MP_ONE);
+ xi = mp_neg(xi, xi);
+ }
+
+ /* Free up the temporary things. */
+ MP_DROP(a); MP_DROP(ai); MP_DROP(bi);
+ MP_DROP(nm1); MP_DROP(delta); mp_drop(d);
+
+ /* Done. */
+ if (exactp_out) *exactp_out = (cmp == 0);
+ return (xi);
+
+#undef f_neg
+}
+
+/* --- @mp_perfect_power_p@ --- *
+ *
+ * Arguments: @mp **x@ = where to write the base
+ * @mp **n@ = where to write the exponent
+ * @mp *a@ = an integer
+ *
+ * Returns: Nonzero if %$a$% is a perfect power.
+ *
+ * Use: Returns whether an integer %$a$% is a perfect power, i.e.,
+ * whether it can be written in the form %$a = x^n$% where
+ * %$|x| > 1$% and %$n > 1$% are integers. If this is possible,
+ * then (a) store %$x$% and the largest such %$n$% in @*x@ and
+ * @*n@, and return nonzero; otherwise, store %$x = a$% and
+ * %$n = 1$% and return zero. (Either @x@ or @n@, or both, may
+ * be null to discard these outputs.)
+ *
+ * Note that %$-1$%, %$0$% and %$1$% are not considered perfect
+ * powers by this definition. (The exponent is not well-defined
+ * in these cases, but it seemed better to implement a function
+ * which worked for all integers.) Note also that %$-4$% is not
+ * a perfect power since it has no real square root.
+ */
+
+int mp_perfect_power_p(mp **x, mp **n, mp *a)
+{
+ mp *r = MP_ONE, *p = MP_NEW, *t = MP_NEW;
+ int exactp, rc = 0;
+ primeiter pi;
+ unsigned f = 0;
+#define f_neg 1u
+
+ MP_COPY(a);
+ if (MP_ZEROP(a)) goto done;
+ primeiter_create(&pi, MP_NEGP(a) ? MP_THREE : MP_TWO);
+ if (MP_NEGP(a)) { a = mp_neg(a, a); f |= f_neg; }
+ p = primeiter_next(&pi, p);
+ for (;;) {
+ t = mp_nthrt(t, a, p, &exactp);
+ if (MP_EQ(t, MP_ONE))
+ break;
+ else if (!exactp) {
+ if (MP_EQ(t, MP_ONE)) break;
+ p = primeiter_next(&pi, p);
+ } else {
+ r = mp_mul(r, r, p);
+ MP_DROP(a); a = t; t = MP_NEW;
+ rc = 1;
+ }
+ }
+ primeiter_destroy(&pi);
+done:
+ 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; }
+ mp_drop(a); mp_drop(p); mp_drop(r); mp_drop(t);
+ return (rc);
+
+#undef f_neg
+}
+
+/*----- Test rig ----------------------------------------------------------*/
+
+#ifdef TEST_RIG
+
+#include <mLib/testrig.h>
+
+static int vrf_nthrt(dstr *v)
+{
+ mp *a = *(mp **)v[0].buf;
+ mp *n = *(mp **)v[1].buf;
+ mp *r0 = *(mp **)v[2].buf;
+ int exactp0 = *(int *)v[3].buf, exactp;
+ mp *r = mp_nthrt(MP_NEW, a, n, &exactp);
+ int ok = 1;
+
+ if (!MP_EQ(r, r0) || exactp != exactp0) {
+ ok = 0;
+ fputs("\n*** nthrt failed", stderr);
+ fputs("\n*** a = ", stderr); mp_writefile(a, stderr, 10);
+ fputs("\n*** n = ", stderr); mp_writefile(n, stderr, 10);
+ fputs("\n*** r calc = ", stderr); mp_writefile(r, stderr, 10);
+ fprintf(stderr, "\n*** ex calc = %d", exactp);
+ fputs("\n*** r exp = ", stderr); mp_writefile(r0, stderr, 10);
+ fprintf(stderr, "\n*** ex exp = %d", exactp0);
+ fputc('\n', stderr);
+ }
+
+ mp_drop(a); mp_drop(n); mp_drop(r); mp_drop(r0);
+ assert(mparena_count(MPARENA_GLOBAL) == 0);
+
+ return (ok);
+}
+
+static int vrf_ppp(dstr *v)
+{
+ mp *a = *(mp **)v[0].buf;
+ int ret0 = *(int *)v[1].buf;
+ mp *x0 = *(mp **)v[2].buf;
+ mp *n0 = *(mp **)v[3].buf;
+ mp *x = MP_NEW, *n = MP_NEW;
+ int ret = mp_perfect_power_p(&x, &n, a);
+ int ok = 1;
+
+ if (ret != ret0 || !MP_EQ(x, x0) || !MP_EQ(n, n0)) {
+ ok = 0;
+ fputs("\n*** perfect_power_p failed", stderr);
+ fputs("\n*** a = ", stderr); mp_writefile(a, stderr, 10);
+ fprintf(stderr, "\n*** r calc = %d", ret);
+ fputs("\n*** x calc = ", stderr); mp_writefile(x, stderr, 10);
+ fputs("\n*** n calc = ", stderr); mp_writefile(n, stderr, 10);
+ fprintf(stderr, "\n*** r exp = %d", ret0);
+ fputs("\n*** x exp = ", stderr); mp_writefile(x0, stderr, 10);
+ fputs("\n*** n exp = ", stderr); mp_writefile(n0, stderr, 10);
+ fputc('\n', stderr);
+ }
+
+ mp_drop(a); mp_drop(x); mp_drop(n); mp_drop(x0); mp_drop(n0);
+ assert(mparena_count(MPARENA_GLOBAL) == 0);
+
+ return (ok);
+}
+
+static test_chunk tests[] = {
+ { "nthrt", vrf_nthrt, { &type_mp, &type_mp, &type_mp, &type_int, 0 } },
+ { "perfect-power-p", vrf_ppp, { &type_mp, &type_int, &type_mp, &type_mp, 0 } },
+ { 0, 0, { 0 } },
+};
+
+int main(int argc, char *argv[])
+{
+ sub_init();
+ test_run(argc, argv, tests, SRCDIR "/t/mp");
+ return (0);
+}
+
+#endif
+
+/*----- That's all, folks -------------------------------------------------*/
extern mp *mp_sqrt(mp */*d*/, mp */*a*/);
+/* --- @mp_nthrt@ --- *
+ *
+ * Arguments: @mp *d@ = fake destination
+ * @mp *a@ = an integer
+ * @mp *n@ = a strictly positive integer
+ * @int *exectp_out@ = set nonzero if an exact solution is found
+ *
+ * Returns: The integer %$\bigl\lfloor \sqrt[n]{a} \bigr\rfloor$%.
+ *
+ * Use: Return an approximation to the %$n$%th root of %$a$%.
+ * Specifically, it returns the largest integer %$x$% such that
+ * %$x^n \le a$%. If %$x^n = a$% then @*exactp_out@ is set
+ * nonzero; otherwise it is set zero. (If @exactp_out@ is null
+ * then this information is discarded.)
+ *
+ * The exponent %$n$% must be strictly positive: it's not clear
+ * to me what the right answer is for %$n \le 0$%. If %$a$% is
+ * negative then %$n$% must be odd; otherwise there is no real
+ * solution.
+ */
+
+extern mp *mp_nthrt(mp */*d*/, mp */*a*/, mp */*n*/, int */*exactp_out*/);
+
+/* --- @mp_perfect_power_p@ --- *
+ *
+ * Arguments: @mp **x@ = where to write the base
+ * @mp **n@ = where to write the exponent
+ * @mp *a@ = an integer
+ *
+ * Returns: Nonzero if %$a$% is a perfect power.
+ *
+ * Use: Returns whether an integer %$a$% is a perfect power, i.e.,
+ * whether it can be written in the form %$a = x^n$% where
+ * %$|x| > 1$% and %$n > 1$% are integers. If this is possible,
+ * then (a) store %$x$% and the largest such %$n$% in @*x@ and
+ * @*n@, and return nonzero; otherwise, store %$x = a$% and
+ * %$n = 1$% and return zero. (Either @x@ or @n@, or both, may
+ * be null to discard these outputs.)
+ *
+ * Note that %$-1$%, %$0$% and %$1$% are not considered perfect
+ * powers by this definition. (The exponent is not well-defined
+ * in these cases, but it seemed better to implement a function
+ * which worked for all integers.) Note also that %$-4$% is not
+ * a perfect power since it has no real square root.
+ */
+
+extern int mp_perfect_power_p(mp **/*x*/, mp **/*n*/, mp */*a*/);
+
/* --- @mp_gcd@ --- *
*
* Arguments: @mp **gcd, **xx, **yy@ = where to write the results
~mdw / catacomb / commitdiff