Commit | Line | Data |
---|---|---|
8a693e0d MW |
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 | |
487795c6 | 98 | * particular, for large %$n$%, if our initial guess is too small, then the |
8a693e0d MW |
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 | ||
36109084 MW |
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 | */ | |
8a693e0d | 177 | MP_COPY(a); |
36109084 MW |
178 | |
179 | /* According to the definition above, @0@ is %%\emph{not}%% a perfect | |
180 | * power. Dispose of this special case early. | |
181 | */ | |
8a693e0d | 182 | if (MP_ZEROP(a)) goto done; |
36109084 MW |
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 | */ | |
8a693e0d MW |
188 | primeiter_create(&pi, MP_NEGP(a) ? MP_THREE : MP_TWO); |
189 | if (MP_NEGP(a)) { a = mp_neg(a, a); f |= f_neg; } | |
36109084 MW |
190 | |
191 | /* Prime the pump for the main loop. */ | |
8a693e0d | 192 | p = primeiter_next(&pi, p); |
36109084 MW |
193 | |
194 | /* Here we go. */ | |
8a693e0d | 195 | for (;;) { |
36109084 MW |
196 | |
197 | /* See if %$a$% is a perfect %$p$%th power. */ | |
8a693e0d MW |
198 | t = mp_nthrt(t, a, p, &exactp); |
199 | if (MP_EQ(t, MP_ONE)) | |
36109084 MW |
200 | /* No, and moreover %$2^p > a$%. We won't get anywhere by examining |
201 | * larger primes, so stop here. | |
202 | */ | |
8a693e0d | 203 | break; |
36109084 | 204 | |
30b7a66d | 205 | else if (!exactp) |
36109084 MW |
206 | /* No. Oh, well: let's try the next prime, then. */ |
207 | ||
8a693e0d | 208 | p = primeiter_next(&pi, p); |
30b7a66d | 209 | else { |
36109084 MW |
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 | */ | |
8a693e0d MW |
215 | r = mp_mul(r, r, p); |
216 | MP_DROP(a); a = t; t = MP_NEW; | |
217 | rc = 1; | |
218 | } | |
219 | } | |
36109084 MW |
220 | |
221 | /* We're all done. */ | |
8a693e0d | 222 | primeiter_destroy(&pi); |
36109084 | 223 | |
8a693e0d | 224 | done: |
36109084 | 225 | /* Return the resulting decomposition. */ |
8a693e0d MW |
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; } | |
36109084 MW |
228 | |
229 | /* Clean everything up. */ | |
8a693e0d | 230 | mp_drop(a); mp_drop(p); mp_drop(r); mp_drop(t); |
36109084 MW |
231 | |
232 | /* All done. */ | |
8a693e0d MW |
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 -------------------------------------------------*/ |