| 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 | MP_COPY(a); |
| 175 | if (MP_ZEROP(a)) goto done; |
| 176 | primeiter_create(&pi, MP_NEGP(a) ? MP_THREE : MP_TWO); |
| 177 | if (MP_NEGP(a)) { a = mp_neg(a, a); f |= f_neg; } |
| 178 | p = primeiter_next(&pi, p); |
| 179 | for (;;) { |
| 180 | t = mp_nthrt(t, a, p, &exactp); |
| 181 | if (MP_EQ(t, MP_ONE)) |
| 182 | break; |
| 183 | else if (!exactp) |
| 184 | p = primeiter_next(&pi, p); |
| 185 | else { |
| 186 | r = mp_mul(r, r, p); |
| 187 | MP_DROP(a); a = t; t = MP_NEW; |
| 188 | rc = 1; |
| 189 | } |
| 190 | } |
| 191 | primeiter_destroy(&pi); |
| 192 | done: |
| 193 | if (x) { mp_drop(*x); *x = f&f_neg ? mp_neg(a, a) : a; a = 0; } |
| 194 | if (n) { mp_drop(*n); *n = r; r = 0; } |
| 195 | mp_drop(a); mp_drop(p); mp_drop(r); mp_drop(t); |
| 196 | return (rc); |
| 197 | |
| 198 | #undef f_neg |
| 199 | } |
| 200 | |
| 201 | /*----- Test rig ----------------------------------------------------------*/ |
| 202 | |
| 203 | #ifdef TEST_RIG |
| 204 | |
| 205 | #include <mLib/testrig.h> |
| 206 | |
| 207 | static int vrf_nthrt(dstr *v) |
| 208 | { |
| 209 | mp *a = *(mp **)v[0].buf; |
| 210 | mp *n = *(mp **)v[1].buf; |
| 211 | mp *r0 = *(mp **)v[2].buf; |
| 212 | int exactp0 = *(int *)v[3].buf, exactp; |
| 213 | mp *r = mp_nthrt(MP_NEW, a, n, &exactp); |
| 214 | int ok = 1; |
| 215 | |
| 216 | if (!MP_EQ(r, r0) || exactp != exactp0) { |
| 217 | ok = 0; |
| 218 | fputs("\n*** nthrt failed", stderr); |
| 219 | fputs("\n*** a = ", stderr); mp_writefile(a, stderr, 10); |
| 220 | fputs("\n*** n = ", stderr); mp_writefile(n, stderr, 10); |
| 221 | fputs("\n*** r calc = ", stderr); mp_writefile(r, stderr, 10); |
| 222 | fprintf(stderr, "\n*** ex calc = %d", exactp); |
| 223 | fputs("\n*** r exp = ", stderr); mp_writefile(r0, stderr, 10); |
| 224 | fprintf(stderr, "\n*** ex exp = %d", exactp0); |
| 225 | fputc('\n', stderr); |
| 226 | } |
| 227 | |
| 228 | mp_drop(a); mp_drop(n); mp_drop(r); mp_drop(r0); |
| 229 | assert(mparena_count(MPARENA_GLOBAL) == 0); |
| 230 | |
| 231 | return (ok); |
| 232 | } |
| 233 | |
| 234 | static int vrf_ppp(dstr *v) |
| 235 | { |
| 236 | mp *a = *(mp **)v[0].buf; |
| 237 | int ret0 = *(int *)v[1].buf; |
| 238 | mp *x0 = *(mp **)v[2].buf; |
| 239 | mp *n0 = *(mp **)v[3].buf; |
| 240 | mp *x = MP_NEW, *n = MP_NEW; |
| 241 | int ret = mp_perfect_power_p(&x, &n, a); |
| 242 | int ok = 1; |
| 243 | |
| 244 | if (ret != ret0 || !MP_EQ(x, x0) || !MP_EQ(n, n0)) { |
| 245 | ok = 0; |
| 246 | fputs("\n*** perfect_power_p failed", stderr); |
| 247 | fputs("\n*** a = ", stderr); mp_writefile(a, stderr, 10); |
| 248 | fprintf(stderr, "\n*** r calc = %d", ret); |
| 249 | fputs("\n*** x calc = ", stderr); mp_writefile(x, stderr, 10); |
| 250 | fputs("\n*** n calc = ", stderr); mp_writefile(n, stderr, 10); |
| 251 | fprintf(stderr, "\n*** r exp = %d", ret0); |
| 252 | fputs("\n*** x exp = ", stderr); mp_writefile(x0, stderr, 10); |
| 253 | fputs("\n*** n exp = ", stderr); mp_writefile(n0, stderr, 10); |
| 254 | fputc('\n', stderr); |
| 255 | } |
| 256 | |
| 257 | mp_drop(a); mp_drop(x); mp_drop(n); mp_drop(x0); mp_drop(n0); |
| 258 | assert(mparena_count(MPARENA_GLOBAL) == 0); |
| 259 | |
| 260 | return (ok); |
| 261 | } |
| 262 | |
| 263 | static test_chunk tests[] = { |
| 264 | { "nthrt", vrf_nthrt, { &type_mp, &type_mp, &type_mp, &type_int, 0 } }, |
| 265 | { "perfect-power-p", vrf_ppp, { &type_mp, &type_int, &type_mp, &type_mp, 0 } }, |
| 266 | { 0, 0, { 0 } }, |
| 267 | }; |
| 268 | |
| 269 | int main(int argc, char *argv[]) |
| 270 | { |
| 271 | sub_init(); |
| 272 | test_run(argc, argv, tests, SRCDIR "/t/mp"); |
| 273 | return (0); |
| 274 | } |
| 275 | |
| 276 | #endif |
| 277 | |
| 278 | /*----- That's all, folks -------------------------------------------------*/ |