From 8a693e0d881d17cabea7051347c495800d80a316 Mon Sep 17 00:00:00 2001 From: Mark Wooding Date: Wed, 22 Jul 2020 23:44:38 +0100 Subject: [PATCH] math/mp-nthrt.c: Implement nth-root, and perfect-power detection. --- math/Makefile.am | 2 + math/mp-nthrt.c | 279 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ math/mp.h | 48 ++++++++++ math/t/mp | 33 +++++++ 4 files changed, 362 insertions(+) create mode 100644 math/mp-nthrt.c diff --git a/math/Makefile.am b/math/Makefile.am index 37d88d1a..8c53940c 100644 --- a/math/Makefile.am +++ b/math/Makefile.am @@ -100,6 +100,8 @@ libmath_la_SOURCES += mp-modsqrt.c TESTS += mp-modsqrt.t$(EXEEXT) libmath_la_SOURCES += mp-sqrt.c TESTS += mp-sqrt.t$(EXEEXT) +libmath_la_SOURCES += mp-nthrt.c +TESTS += mp-nthrt.t$(EXEEXT) libmath_la_SOURCES += mp-test.c EXTRA_DIST += t/mp diff --git a/math/mp-nthrt.c b/math/mp-nthrt.c new file mode 100644 index 00000000..a46590af --- /dev/null +++ b/math/mp-nthrt.c @@ -0,0 +1,279 @@ +/* -*-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 + +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 -------------------------------------------------*/ diff --git a/math/mp.h b/math/mp.h index 0434c282..23e71e36 100644 --- a/math/mp.h +++ b/math/mp.h @@ -923,6 +923,54 @@ extern mp *mp_leastcongruent(mp */*d*/, mp */*b*/, mp */*r*/, mp */*m*/); 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 diff --git a/math/t/mp b/math/t/mp index b50ea8dd..8d0942d7 100644 --- a/math/t/mp +++ b/math/t/mp @@ -115,6 +115,39 @@ sqrt { 14565040310136678240 3816417208; } +nthrt { + 0 27 0 1; + 1 1 1 1; + 99 2 9 0; + 100 2 10 1; + 101 2 10 0; + + 2432442434617858985744608211960343276041892998697958012044143077567778256072696563501333460622383626492631158845093813667916645390906408185968436731121086804986194010729874783817632607960227495980162127756247771205609001938726 37 1234566 0; + 2432442434617858985744608211960343276041892998697958012044143077567778256072696563501333460622383626492631158845093813667916645390906408185968436731121086804986194010729874783817632607960227495980162127756247771205609001938727 37 1234567 1; + 2432442434617858985744608211960343276041892998697958012044143077567778256072696563501333460622383626492631158845093813667916645390906408185968436731121086804986194010729874783817632607960227495980162127756247771205609001938728 37 1234567 0; + + -26 3 -3 0; + -27 3 -3 1; + -28 3 -4 0; +} + +perfect-power-p { + 0 0 0 1; + 1 0 1 1; + -1 0 -1 1; + + -4 0 -4 1; + -8 1 -2 3; + -64 1 -4 3; + + 80 0 80 1; + 81 1 3 4; + 82 0 82 1; + + 2432442434617858985744608211960343276041892998697958012044143077567778256072696563501333460622383626492631158845093813667916645390906408185968436731121086804986194010729874783817632607960227495980162127756247771205609001938727 1 1234567 37; + 42467986438630307821661186973460619303572935864570185492440237295438188325895624954701633272891586038903309915601221633039963682793644006385615600911359153716020597273608200491915551536581527267184634993651215467730190125236224 1 12 210; +} + gcd { # --- Simple tests --- -- 2.11.0