3 * Calculating %$n$%th roots
5 * (c) 2020 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
33 #include "primeiter.h"
35 /*----- Main code ---------------------------------------------------------*/
37 /* --- @mp_nthrt@ --- *
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
44 * Returns: The integer %$\bigl\lfloor \sqrt[n]{a} \bigr\rfloor$%.
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.)
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
58 mp
*mp_nthrt(mp
*d
, mp
*a
, mp
*n
, int *exactp_out
)
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})$%.
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
74 mp
*ai
= MP_NEW
, *bi
= MP_NEW
, *xi
= MP_NEW
,
75 *nm1
= MP_NEW
, *delta
= MP_NEW
;
80 /* Firstly, deal with a zero or negative %$xa%. */
86 a
= mp_neg(MP_NEW
, a
);
89 /* Secondly, reject %$n \le 0$%.
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.
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
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
));
108 /* The main iteration. */
110 ai
= mp_exp(ai
, xi
, n
);
111 bi
= mp_add(bi
, xi
, MP_ONE
); bi
= mp_exp(bi
, bi
, n
);
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
);
121 /* Final sorting out of negative inputs.
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
128 if (cmp
) xi
= mp_add(xi
, xi
, MP_ONE
);
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
);
137 if (exactp_out
) *exactp_out
= (cmp
== 0);
143 /* --- @mp_perfect_power_p@ --- *
145 * Arguments: @mp **x@ = where to write the base
146 * @mp **n@ = where to write the exponent
147 * @mp *a@ = an integer
149 * Returns: Nonzero if %$a$% is a perfect power.
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.)
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.
166 int mp_perfect_power_p(mp
**x
, mp
**n
, mp
*a
)
168 mp
*r
= MP_ONE
, *p
= MP_NEW
, *t
= MP_NEW
;
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.
179 /* According to the definition above, @0@ is %%\emph{not}%% a perfect
180 * power. Dispose of this special case early.
182 if (MP_ZEROP(a
)) goto done
;
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
188 primeiter_create(&pi
, MP_NEGP(a
) ? MP_THREE
: MP_TWO
);
189 if (MP_NEGP(a
)) { a
= mp_neg(a
, a
); f
|= f_neg
; }
191 /* Prime the pump for the main loop. */
192 p
= primeiter_next(&pi
, p
);
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.
206 /* No. Oh, well: let's try the next prime, then. */
208 p
= primeiter_next(&pi
, p
);
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'}$%.
216 MP_DROP(a
); a
= t
; t
= MP_NEW
;
221 /* We're all done. */
222 primeiter_destroy(&pi
);
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; }
229 /* Clean everything up. */
230 mp_drop(a
); mp_drop(p
); mp_drop(r
); mp_drop(t
);
238 /*----- Test rig ----------------------------------------------------------*/
242 #include <mLib/testrig.h>
244 static int vrf_nthrt(dstr
*v
)
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
);
253 if (!MP_EQ(r
, r0
) || exactp
!= exactp0
) {
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
);
265 mp_drop(a
); mp_drop(n
); mp_drop(r
); mp_drop(r0
);
266 assert(mparena_count(MPARENA_GLOBAL
) == 0);
271 static int vrf_ppp(dstr
*v
)
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
);
281 if (ret
!= ret0
|| !MP_EQ(x
, x0
) || !MP_EQ(n
, n0
)) {
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);
294 mp_drop(a
); mp_drop(x
); mp_drop(n
); mp_drop(x0
); mp_drop(n0
);
295 assert(mparena_count(MPARENA_GLOBAL
) == 0);
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 } },
306 int main(int argc
, char *argv
[])
309 test_run(argc
, argv
, tests
, SRCDIR
"/t/mp");
315 /*----- That's all, folks -------------------------------------------------*/