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 we our initial guess 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
;
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
);
180 t
= mp_nthrt(t
, a
, p
, &exactp
);
181 if (MP_EQ(t
, MP_ONE
))
184 if (MP_EQ(t
, MP_ONE
)) break;
185 p
= primeiter_next(&pi
, p
);
188 MP_DROP(a
); a
= t
; t
= MP_NEW
;
192 primeiter_destroy(&pi
);
194 if (x
) { mp_drop(*x
); *x
= f
&f_neg ?
mp_neg(a
, a
) : a
; a
= 0; }
195 if (n
) { mp_drop(*n
); *n
= r
; r
= 0; }
196 mp_drop(a
); mp_drop(p
); mp_drop(r
); mp_drop(t
);
202 /*----- Test rig ----------------------------------------------------------*/
206 #include <mLib/testrig.h>
208 static int vrf_nthrt(dstr
*v
)
210 mp
*a
= *(mp
**)v
[0].buf
;
211 mp
*n
= *(mp
**)v
[1].buf
;
212 mp
*r0
= *(mp
**)v
[2].buf
;
213 int exactp0
= *(int *)v
[3].buf
, exactp
;
214 mp
*r
= mp_nthrt(MP_NEW
, a
, n
, &exactp
);
217 if (!MP_EQ(r
, r0
) || exactp
!= exactp0
) {
219 fputs("\n*** nthrt failed", stderr
);
220 fputs("\n*** a = ", stderr
); mp_writefile(a
, stderr
, 10);
221 fputs("\n*** n = ", stderr
); mp_writefile(n
, stderr
, 10);
222 fputs("\n*** r calc = ", stderr
); mp_writefile(r
, stderr
, 10);
223 fprintf(stderr
, "\n*** ex calc = %d", exactp
);
224 fputs("\n*** r exp = ", stderr
); mp_writefile(r0
, stderr
, 10);
225 fprintf(stderr
, "\n*** ex exp = %d", exactp0
);
229 mp_drop(a
); mp_drop(n
); mp_drop(r
); mp_drop(r0
);
230 assert(mparena_count(MPARENA_GLOBAL
) == 0);
235 static int vrf_ppp(dstr
*v
)
237 mp
*a
= *(mp
**)v
[0].buf
;
238 int ret0
= *(int *)v
[1].buf
;
239 mp
*x0
= *(mp
**)v
[2].buf
;
240 mp
*n0
= *(mp
**)v
[3].buf
;
241 mp
*x
= MP_NEW
, *n
= MP_NEW
;
242 int ret
= mp_perfect_power_p(&x
, &n
, a
);
245 if (ret
!= ret0
|| !MP_EQ(x
, x0
) || !MP_EQ(n
, n0
)) {
247 fputs("\n*** perfect_power_p failed", stderr
);
248 fputs("\n*** a = ", stderr
); mp_writefile(a
, stderr
, 10);
249 fprintf(stderr
, "\n*** r calc = %d", ret
);
250 fputs("\n*** x calc = ", stderr
); mp_writefile(x
, stderr
, 10);
251 fputs("\n*** n calc = ", stderr
); mp_writefile(n
, stderr
, 10);
252 fprintf(stderr
, "\n*** r exp = %d", ret0
);
253 fputs("\n*** x exp = ", stderr
); mp_writefile(x0
, stderr
, 10);
254 fputs("\n*** n exp = ", stderr
); mp_writefile(n0
, stderr
, 10);
258 mp_drop(a
); mp_drop(x
); mp_drop(n
); mp_drop(x0
); mp_drop(n0
);
259 assert(mparena_count(MPARENA_GLOBAL
) == 0);
264 static test_chunk tests
[] = {
265 { "nthrt", vrf_nthrt
, { &type_mp
, &type_mp
, &type_mp
, &type_int
, 0 } },
266 { "perfect-power-p", vrf_ppp
, { &type_mp
, &type_int
, &type_mp
, &type_mp
, 0 } },
270 int main(int argc
, char *argv
[])
273 test_run(argc
, argv
, tests
, SRCDIR
"/t/mp");
279 /*----- That's all, folks -------------------------------------------------*/