AuroraMiddleware/libtommath
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

changed seed to make nth-root usable

Browse Source
This commit is contained in:
czurnieden 2019-04-01 03:41:26 +02:00
parent 432f995ff3
commit 984d3ff679
5 changed files with 395 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
bn_mp_n_root_ex.c
View File

@ -19,13 +19,13 @@
* This algorithm uses Newton's approximation
* x[i+1] = x[i] - f(x[i])/f'(x[i])
* which will find the root in log(N) time where
* each step involves a fair bit. This is not meant to
* find huge roots [square and cube, etc].
* each step involves a fair bit.
*/
int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
{
mp_int t1, t2, t3, a_;
int res;
int res, cmp;
int ilog2;
/* input must be positive if b is even */
if (((b & 1u) == 0u) && (a->sign == MP_NEG)) {
@ -48,9 +48,49 @@ int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
a_ = *a;
a_.sign = MP_ZPOS;
/* t2 = 2 */
mp_set(&t2, 2uL);
/* Compute seed: 2^(log_2(n)/b + 2)*/
ilog2 = mp_count_bits(a);
/*
GCC and clang do not understand the sizeof tests and complain,
icc (the Intel compiler) seems to understand, at least it doesn't complain.
2 of 3 say these macros are necessary, so there they are.
*/
#if ( !(defined MP_8BIT) && !(defined MP_16BIT) )
/*
The type of mp_digit might be larger than an int.
If "b" is larger than INT_MAX it is also larger than
log_2(n) because the bit-length of the "n" is measured
with an int and hence the root is always < 2 (two).
*/
if (sizeof(mp_digit) >= sizeof(int)) {
if (b > (mp_digit)(INT_MAX/2)) {
mp_set(c, 1uL);
c->sign = a->sign;
res = MP_OKAY;
goto LBL_T3;
}
}
#endif
/* "b" is smaller than INT_MAX, we can cast safely */
if (ilog2 < (int)b) {
mp_set(c, 1uL);
c->sign = a->sign;
res = MP_OKAY;
goto LBL_T3;
}
ilog2 = ilog2 / ((int)b);
if (ilog2 == 0) {
mp_set(c, 1uL);
c->sign = a->sign;
res = MP_OKAY;
goto LBL_T3;
}
/* Start value must be larger than root */
ilog2 += 2;
if ((res = mp_2expt(&t2,ilog2)) != MP_OKAY) {
goto LBL_T3;
}
do {
/* t1 = t2 */
if ((res = mp_copy(&t2, &t1)) != MP_OKAY) {
@ -63,7 +103,6 @@ int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
if ((res = mp_expt_d_ex(&t1, b - 1u, &t3, fast)) != MP_OKAY) {
goto LBL_T3;
}
/* numerator */
/* t2 = t1**b */
if ((res = mp_mul(&t3, &t1, &t2)) != MP_OKAY) {
@ -89,14 +128,39 @@ int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
if ((res = mp_sub(&t1, &t3, &t2)) != MP_OKAY) {
goto LBL_T3;
}
/*
Number of rounds is at most log_2(root). If it is more it
got stuck, so break out of the loop and do the rest manually.
*/
if (ilog2-- == 0) {
break;
}
} while (mp_cmp(&t1, &t2) != MP_EQ);
/* result can be off by a few so check */
/* Loop beneath can overshoot by one if found root is smaller than actual root */
for (;;) {
if ((res = mp_expt_d_ex(&t1, b, &t2, fast)) != MP_OKAY) {
goto LBL_T3;
}
cmp = mp_cmp(&t2, &a_);
if (cmp == MP_EQ) {
res = MP_OKAY;
goto LBL_T3;
}
if (cmp == MP_LT) {
if ((res = mp_add_d(&t1, 1uL, &t1)) != MP_OKAY) {
goto LBL_T3;
}
} else {
break;
}
}
/* correct overshoot from above or from recurrence */
for (;;) {
if ((res = mp_expt_d_ex(&t1, b, &t2, fast)) != MP_OKAY) {
goto LBL_T3;
}
if (mp_cmp(&t2, &a_) == MP_GT) {
if ((res = mp_sub_d(&t1, 1uL, &t1)) != MP_OKAY) {
goto LBL_T3;
@ -123,7 +187,6 @@ LBL_T1:
return res;
}
#endif
/* ref: $Format:%D$ */
/* git commit: $Format:%H$ */
/* commit time: $Format:%ai$ */

101
callgraph.txt
View File

@ -2811,8 +2811,12 @@ BN_MP_IS_SQUARE_C
| +--->BN_MP_N_ROOT_C
| | +--->BN_MP_N_ROOT_EX_C
| | | +--->BN_MP_INIT_C
| | | +--->BN_MP_COUNT_BITS_C
| | | +--->BN_MP_SET_C
| | | | +--->BN_MP_ZERO_C
| | | +--->BN_MP_2EXPT_C
| | | | +--->BN_MP_ZERO_C
| | | | +--->BN_MP_GROW_C
| | | +--->BN_MP_COPY_C
| | | | +--->BN_MP_GROW_C
| | | +--->BN_MP_EXPT_D_EX_C
@ -3039,7 +3043,6 @@ BN_MP_IS_SQUARE_C
| | | | +--->BN_MP_ZERO_C
| | | | +--->BN_MP_INIT_MULTI_C
| | | | | +--->BN_MP_CLEAR_C
| | | | +--->BN_MP_COUNT_BITS_C
| | | | +--->BN_MP_ABS_C
| | | | +--->BN_MP_MUL_2D_C
| | | | | +--->BN_MP_GROW_C
@ -3073,10 +3076,13 @@ BN_MP_IS_SQUARE_C
| | | | +--->BN_MP_CLEAR_C
| | | +--->BN_MP_CMP_C
| | | | +--->BN_MP_CMP_MAG_C
| | | +--->BN_MP_ADD_D_C
| | | | +--->BN_MP_GROW_C
| | | | +--->BN_MP_SUB_D_C
| | | | | +--->BN_MP_CLAMP_C
| | | | +--->BN_MP_CLAMP_C
| | | +--->BN_MP_SUB_D_C
| | | | +--->BN_MP_GROW_C
| | | | +--->BN_MP_ADD_D_C
| | | | | +--->BN_MP_CLAMP_C
| | | | +--->BN_MP_CLAMP_C
| | | +--->BN_MP_EXCH_C
| | | +--->BN_MP_CLEAR_C
@ -4055,8 +4061,12 @@ BN_MP_NEG_C
BN_MP_N_ROOT_C
+--->BN_MP_N_ROOT_EX_C
| +--->BN_MP_INIT_C
| +--->BN_MP_COUNT_BITS_C
| +--->BN_MP_SET_C
| | +--->BN_MP_ZERO_C
| +--->BN_MP_2EXPT_C
| | +--->BN_MP_ZERO_C
| | +--->BN_MP_GROW_C
| +--->BN_MP_COPY_C
| | +--->BN_MP_GROW_C
| +--->BN_MP_EXPT_D_EX_C
@ -4283,7 +4293,6 @@ BN_MP_N_ROOT_C
| | +--->BN_MP_ZERO_C
| | +--->BN_MP_INIT_MULTI_C
| | | +--->BN_MP_CLEAR_C
| | +--->BN_MP_COUNT_BITS_C
| | +--->BN_MP_ABS_C
| | +--->BN_MP_MUL_2D_C
| | | +--->BN_MP_GROW_C
@ -4317,10 +4326,13 @@ BN_MP_N_ROOT_C
| | +--->BN_MP_CLEAR_C
| +--->BN_MP_CMP_C
| | +--->BN_MP_CMP_MAG_C
| +--->BN_MP_ADD_D_C
| | +--->BN_MP_GROW_C
| | +--->BN_MP_SUB_D_C
| | | +--->BN_MP_CLAMP_C
| | +--->BN_MP_CLAMP_C
| +--->BN_MP_SUB_D_C
| | +--->BN_MP_GROW_C
| | +--->BN_MP_ADD_D_C
| | | +--->BN_MP_CLAMP_C
| | +--->BN_MP_CLAMP_C
| +--->BN_MP_EXCH_C
| +--->BN_MP_CLEAR_C
@ -4328,8 +4340,12 @@ BN_MP_N_ROOT_C
BN_MP_N_ROOT_EX_C
+--->BN_MP_INIT_C
+--->BN_MP_COUNT_BITS_C
+--->BN_MP_SET_C
| +--->BN_MP_ZERO_C
+--->BN_MP_2EXPT_C
| +--->BN_MP_ZERO_C
| +--->BN_MP_GROW_C
+--->BN_MP_COPY_C
| +--->BN_MP_GROW_C
+--->BN_MP_EXPT_D_EX_C
@ -4556,7 +4572,6 @@ BN_MP_N_ROOT_EX_C
| +--->BN_MP_ZERO_C
| +--->BN_MP_INIT_MULTI_C
| | +--->BN_MP_CLEAR_C
| +--->BN_MP_COUNT_BITS_C
| +--->BN_MP_ABS_C
| +--->BN_MP_MUL_2D_C
| | +--->BN_MP_GROW_C
@ -4590,10 +4605,13 @@ BN_MP_N_ROOT_EX_C
| +--->BN_MP_CLEAR_C
+--->BN_MP_CMP_C
| +--->BN_MP_CMP_MAG_C
+--->BN_MP_ADD_D_C
| +--->BN_MP_GROW_C
| +--->BN_MP_SUB_D_C
| | +--->BN_MP_CLAMP_C
| +--->BN_MP_CLAMP_C
+--->BN_MP_SUB_D_C
| +--->BN_MP_GROW_C
| +--->BN_MP_ADD_D_C
| | +--->BN_MP_CLAMP_C
| +--->BN_MP_CLAMP_C
+--->BN_MP_EXCH_C
+--->BN_MP_CLEAR_C
@ -5725,8 +5743,12 @@ BN_MP_PRIME_FROBENIUS_UNDERWOOD_C
| | | +--->BN_MP_N_ROOT_C
| | | | +--->BN_MP_N_ROOT_EX_C
| | | | | +--->BN_MP_INIT_C
| | | | | +--->BN_MP_COUNT_BITS_C
| | | | | +--->BN_MP_SET_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | +--->BN_MP_2EXPT_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_COPY_C
| | | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_EXPT_D_EX_C
@ -5953,7 +5975,6 @@ BN_MP_PRIME_FROBENIUS_UNDERWOOD_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | | +--->BN_MP_INIT_MULTI_C
| | | | | | | +--->BN_MP_CLEAR_C
| | | | | | +--->BN_MP_COUNT_BITS_C
| | | | | | +--->BN_MP_ABS_C
| | | | | | +--->BN_MP_MUL_2D_C
| | | | | | | +--->BN_MP_GROW_C
@ -5987,10 +6008,13 @@ BN_MP_PRIME_FROBENIUS_UNDERWOOD_C
| | | | | | +--->BN_MP_CLEAR_C
| | | | | +--->BN_MP_CMP_C
| | | | | | +--->BN_MP_CMP_MAG_C
| | | | | +--->BN_MP_ADD_D_C
| | | | | | +--->BN_MP_GROW_C
| | | | | | +--->BN_MP_SUB_D_C
| | | | | | | +--->BN_MP_CLAMP_C
| | | | | | +--->BN_MP_CLAMP_C
| | | | | +--->BN_MP_SUB_D_C
| | | | | | +--->BN_MP_GROW_C
| | | | | | +--->BN_MP_ADD_D_C
| | | | | | | +--->BN_MP_CLAMP_C
| | | | | | +--->BN_MP_CLAMP_C
| | | | | +--->BN_MP_EXCH_C
| | | | | +--->BN_MP_CLEAR_C
@ -8035,8 +8059,12 @@ BN_MP_PRIME_IS_PRIME_C
| | +--->BN_MP_N_ROOT_C
| | | +--->BN_MP_N_ROOT_EX_C
| | | | +--->BN_MP_INIT_C
| | | | +--->BN_MP_COUNT_BITS_C
| | | | +--->BN_MP_SET_C
| | | | | +--->BN_MP_ZERO_C
| | | | +--->BN_MP_2EXPT_C
| | | | | +--->BN_MP_ZERO_C
| | | | | +--->BN_MP_GROW_C
| | | | +--->BN_MP_COPY_C
| | | | | +--->BN_MP_GROW_C
| | | | +--->BN_MP_EXPT_D_EX_C
@ -8263,7 +8291,6 @@ BN_MP_PRIME_IS_PRIME_C
| | | | | +--->BN_MP_ZERO_C
| | | | | +--->BN_MP_INIT_MULTI_C
| | | | | | +--->BN_MP_CLEAR_C
| | | | | +--->BN_MP_COUNT_BITS_C
| | | | | +--->BN_MP_ABS_C
| | | | | +--->BN_MP_MUL_2D_C
| | | | | | +--->BN_MP_GROW_C
@ -8297,10 +8324,13 @@ BN_MP_PRIME_IS_PRIME_C
| | | | | +--->BN_MP_CLEAR_C
| | | | +--->BN_MP_CMP_C
| | | | | +--->BN_MP_CMP_MAG_C
| | | | +--->BN_MP_ADD_D_C
| | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_SUB_D_C
| | | | | | +--->BN_MP_CLAMP_C
| | | | | +--->BN_MP_CLAMP_C
| | | | +--->BN_MP_SUB_D_C
| | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_ADD_D_C
| | | | | | +--->BN_MP_CLAMP_C
| | | | | +--->BN_MP_CLAMP_C
| | | | +--->BN_MP_EXCH_C
| | | | +--->BN_MP_CLEAR_C
@ -11420,6 +11450,10 @@ BN_MP_PRIME_NEXT_PRIME_C
| | +--->BN_MP_SQRT_C
| | | +--->BN_MP_N_ROOT_C
| | | | +--->BN_MP_N_ROOT_EX_C
| | | | | +--->BN_MP_COUNT_BITS_C
| | | | | +--->BN_MP_2EXPT_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_COPY_C
| | | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_EXPT_D_EX_C
@ -11646,7 +11680,6 @@ BN_MP_PRIME_NEXT_PRIME_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | | +--->BN_MP_INIT_MULTI_C
| | | | | | | +--->BN_MP_CLEAR_C
| | | | | | +--->BN_MP_COUNT_BITS_C
| | | | | | +--->BN_MP_ABS_C
| | | | | | +--->BN_MP_MUL_2D_C
| | | | | | | +--->BN_MP_GROW_C
@ -13616,8 +13649,12 @@ BN_MP_PRIME_RANDOM_EX_C
| | | +--->BN_MP_N_ROOT_C
| | | | +--->BN_MP_N_ROOT_EX_C
| | | | | +--->BN_MP_INIT_C
| | | | | +--->BN_MP_COUNT_BITS_C
| | | | | +--->BN_MP_SET_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | +--->BN_MP_2EXPT_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_COPY_C
| | | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_EXPT_D_EX_C
@ -13844,7 +13881,6 @@ BN_MP_PRIME_RANDOM_EX_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | | +--->BN_MP_INIT_MULTI_C
| | | | | | | +--->BN_MP_CLEAR_C
| | | | | | +--->BN_MP_COUNT_BITS_C
| | | | | | +--->BN_MP_ABS_C
| | | | | | +--->BN_MP_MUL_2D_C
| | | | | | | +--->BN_MP_GROW_C
@ -13878,10 +13914,13 @@ BN_MP_PRIME_RANDOM_EX_C
| | | | | | +--->BN_MP_CLEAR_C
| | | | | +--->BN_MP_CMP_C
| | | | | | +--->BN_MP_CMP_MAG_C
| | | | | +--->BN_MP_ADD_D_C
| | | | | | +--->BN_MP_GROW_C
| | | | | | +--->BN_MP_SUB_D_C
| | | | | | | +--->BN_MP_CLAMP_C
| | | | | | +--->BN_MP_CLAMP_C
| | | | | +--->BN_MP_SUB_D_C
| | | | | | +--->BN_MP_GROW_C
| | | | | | +--->BN_MP_ADD_D_C
| | | | | | | +--->BN_MP_CLAMP_C
| | | | | | +--->BN_MP_CLAMP_C
| | | | | +--->BN_MP_EXCH_C
| | | | | +--->BN_MP_CLEAR_C
@ -15916,8 +15955,12 @@ BN_MP_PRIME_STRONG_LUCAS_SELFRIDGE_C
| | | +--->BN_MP_N_ROOT_C
| | | | +--->BN_MP_N_ROOT_EX_C
| | | | | +--->BN_MP_INIT_C
| | | | | +--->BN_MP_COUNT_BITS_C
| | | | | +--->BN_MP_SET_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | +--->BN_MP_2EXPT_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_COPY_C
| | | | | | +--->BN_MP_GROW_C
| | | | | +--->BN_MP_EXPT_D_EX_C
@ -16144,7 +16187,6 @@ BN_MP_PRIME_STRONG_LUCAS_SELFRIDGE_C
| | | | | | +--->BN_MP_ZERO_C
| | | | | | +--->BN_MP_INIT_MULTI_C
| | | | | | | +--->BN_MP_CLEAR_C
| | | | | | +--->BN_MP_COUNT_BITS_C
| | | | | | +--->BN_MP_ABS_C
| | | | | | +--->BN_MP_MUL_2D_C
| | | | | | | +--->BN_MP_GROW_C
@ -16178,10 +16220,13 @@ BN_MP_PRIME_STRONG_LUCAS_SELFRIDGE_C
| | | | | | +--->BN_MP_CLEAR_C
| | | | | +--->BN_MP_CMP_C
| | | | | | +--->BN_MP_CMP_MAG_C
| | | | | +--->BN_MP_ADD_D_C
| | | | | | +--->BN_MP_GROW_C
| | | | | | +--->BN_MP_SUB_D_C
| | | | | | | +--->BN_MP_CLAMP_C
| | | | | | +--->BN_MP_CLAMP_C
| | | | | +--->BN_MP_SUB_D_C
| | | | | | +--->BN_MP_GROW_C
| | | | | | +--->BN_MP_ADD_D_C
| | | | | | | +--->BN_MP_CLAMP_C
| | | | | | +--->BN_MP_CLAMP_C
| | | | | +--->BN_MP_EXCH_C
| | | | | +--->BN_MP_CLEAR_C
@ -20056,8 +20101,12 @@ BN_MP_SQRT_C
+--->BN_MP_N_ROOT_C
| +--->BN_MP_N_ROOT_EX_C
| | +--->BN_MP_INIT_C
| | +--->BN_MP_COUNT_BITS_C
| | +--->BN_MP_SET_C
| | | +--->BN_MP_ZERO_C
| | +--->BN_MP_2EXPT_C
| | | +--->BN_MP_ZERO_C
| | | +--->BN_MP_GROW_C
| | +--->BN_MP_COPY_C
| | | +--->BN_MP_GROW_C
| | +--->BN_MP_EXPT_D_EX_C
@ -20284,7 +20333,6 @@ BN_MP_SQRT_C
| | | +--->BN_MP_ZERO_C
| | | +--->BN_MP_INIT_MULTI_C
| | | | +--->BN_MP_CLEAR_C
| | | +--->BN_MP_COUNT_BITS_C
| | | +--->BN_MP_ABS_C
| | | +--->BN_MP_MUL_2D_C
| | | | +--->BN_MP_GROW_C
@ -20318,10 +20366,13 @@ BN_MP_SQRT_C
| | | +--->BN_MP_CLEAR_C
| | +--->BN_MP_CMP_C
| | | +--->BN_MP_CMP_MAG_C
| | +--->BN_MP_ADD_D_C
| | | +--->BN_MP_GROW_C
| | | +--->BN_MP_SUB_D_C
| | | | +--->BN_MP_CLAMP_C
| | | +--->BN_MP_CLAMP_C
| | +--->BN_MP_SUB_D_C
| | | +--->BN_MP_GROW_C
| | | +--->BN_MP_ADD_D_C
| | | | +--->BN_MP_CLAMP_C
| | | +--->BN_MP_CLAMP_C
| | +--->BN_MP_EXCH_C
| | +--->BN_MP_CLEAR_C

243
demo/test.c
View File

@ -1362,6 +1362,248 @@ LTM_ERR:
return EXIT_FAILURE;
}
/*
Cannot test mp_exp(_d) without mp_n_root and vice versa.
So one of the two has to be tested from scratch.
Numbers generated by
for i in {1..10}
do
seed=$(head -c 10000 /dev/urandom | tr -dc '[:digit:]' | head -c 120);
echo $seed;
convertbase $seed 10 64;
done
(The program "convertbase" uses libtommath's to/from_radix functions)
Roots were precalculated with Pari/GP
default(realprecision,1000);
for(n=3,100,r = floor(a^(1/n));printf("\"" r "\", "))
All numbers as strings to simplifiy things, especially for the
low-mp branch.
*/
static int test_mp_n_root(void)
{
mp_int a, c, r;
int e;
int i, j;
const char *input[] = {
"4n9cbk886QtLQmofprid3l2Q0GD8Yv979Lh8BdZkFE8g2pDUUSMBET/+M/YFyVZ3mBp",
"5NlgzHhmIX05O5YoW5yW5reAlVNtRAlIcN2dfoATnNdc1Cw5lHZUTwNthmK6/ZLKfY6",
"3gweiHDX+ji5utraSe46IJX+uuh7iggs63xIpMP5MriU4Np+LpHI5are8RzS9pKh9xP",
"5QOJUSKMrfe7LkeyJOlupS8h7bjT+TXmZkDzOjZtfj7mdA7cbg0lRX3CuafhjIrpK8S",
"4HtYFldVkyVbrlg/s7kmaA7j45PvLQm+1bbn6ehgP8tVoBmGbv2yDQI1iQQze4AlHyN",
"3bwCUx79NAR7c68OPSp5ZabhZ9aBEr7rWNTO2oMY7zhbbbw7p6shSMxqE9K9nrTNucf",
"4j5RGb78TfuYSzrXn0z6tiAoWiRI81hGY3el9AEa9S+gN4x/AmzotHT2Hvj6lyBpE7q",
"4lwg30SXqZhEHNsl5LIXdyu7UNt0VTWebP3m7+WUL+hsnFW9xJe7UnzYngZsvWh14IE",
"1+tcqFeRuGqjRADRoRUJ8gL4UUSFQVrVVoV6JpwVcKsuBq5G0pABn0dLcQQQMViiVRj",
"hXwxuFySNSFcmbrs/coz4FUAaUYaOEt+l4V5V8vY71KyBvQPxRq/6lsSrG2FHvWDax"
};
/* roots 3-100 of the above */
const char *root[10][100] = {
{
"9163694094944489658600517465135586130944",
"936597377180979771960755204040", "948947857956884030956907",
"95727185767390496595", "133844854039712620", "967779611885360",
"20926191452627", "974139547476", "79203891950", "9784027073",
"1667309744", "365848129", "98268452", "31109156", "11275351",
"4574515", "2040800", "986985", "511525", "281431", "163096",
"98914", "62437", "40832", "27556", "19127", "13614", "9913",
"7367", "5577", "4294", "3357", "2662", "2138", "1738", "1428",
"1185", "993", "839", "715", "613", "530", "461", "403", "355",
"314", "279", "249", "224", "202", "182", "166", "151", "138",
"126", "116", "107", "99", "92", "85", "79", "74", "69", "65", "61",
"57", "54", "51", "48", "46", "43", "41", "39", "37", "36", "34",
"32", "31", "30", "28", "27", "26", "25", "24", "23", "23", "22",
"21", "20", "20", "19", "18", "18", "17", "17", "16", "16", "15"
}, {
"9534798256755061606359588498764080011382",
"964902943621813525741417593772", "971822399862464674540423",
"97646291566833512831", "136141536090599560", "982294733581430",
"21204945933335", "985810529393", "80066084985", "9881613813",
"1682654547", "368973625", "99051783", "31341581", "11354620",
"4604882", "2053633", "992879", "514434", "282959", "163942",
"99406", "62736", "41020", "27678", "19208", "13670", "9952",
"7395", "5598", "4310", "3369", "2671", "2145", "1744", "1433",
"1189", "996", "842", "717", "615", "531", "462", "404", "356",
"315", "280", "250", "224", "202", "183", "166", "151", "138",
"127", "116", "107", "99", "92", "85", "80", "74", "70", "65", "61",
"58", "54", "51", "48", "46", "43", "41", "39", "37", "36", "34",
"32", "31", "30", "29", "27", "26", "25", "24", "23", "23", "22",
"21", "20", "20", "19", "18", "18", "17", "17", "16", "16", "15"
}, {
"8398539113202579297642815367509019445624",
"877309458945432597462853440936", "900579899458998599215071",
"91643543761699761637", "128935656335800903", "936647990947203",
"20326748623514", "948988882684", "77342677787", "9573063447",
"1634096832", "359076114", "96569670", "30604705", "11103188",
"4508519", "2012897", "974160", "505193", "278105", "161251",
"97842", "61788", "40423", "27291", "18949", "13492", "9826",
"7305", "5532", "4260", "3332", "2642", "2123", "1726", "1418",
"1177", "986", "834", "710", "610", "527", "458", "401", "353",
"312", "278", "248", "223", "201", "181", "165", "150", "137",
"126", "116", "107", "99", "91", "85", "79", "74", "69", "65", "61",
"57", "54", "51", "48", "46", "43", "41", "39", "37", "35", "34",
"32", "31", "30", "28", "27", "26", "25", "24", "23", "22", "22",
"21", "20", "20", "19", "18", "18", "17", "17", "16", "16", "15"
}, {
"9559098494021810340217797724866627755195",
"966746709063325235560830083787", "973307706084821682248292",
"97770642291138756434", "136290128605981259", "983232784778520",
"21222944848922", "986563584410", "80121684894", "9887903837",
"1683643206", "369174929", "99102220", "31356542", "11359721",
"4606836", "2054458", "993259", "514621", "283057", "163997",
"99437", "62755", "41032", "27686", "19213", "13674", "9955",
"7397", "5599", "4311", "3370", "2672", "2146", "1744", "1433",
"1189", "996", "842", "717", "615", "532", "462", "404", "356",
"315", "280", "250", "224", "202", "183", "166", "151", "138",
"127", "116", "107", "99", "92", "86", "80", "74", "70", "65", "61",
"58", "54", "51", "48", "46", "43", "41", "39", "37", "36", "34",
"32", "31", "30", "29", "27", "26", "25", "24", "23", "23", "22",
"21", "20", "20", "19", "18", "18", "17", "17", "16", "16", "15"
}, {
"8839202025813295923132694443541993309220",
"911611499784863252820288596270", "928640961450376817534853",
"94017030509441723821", "131792686685970629", "954783483196511",
"20676214073400", "963660189823", "78428929840", "9696237956",
"1653495486", "363032624", "97562430", "30899570", "11203842",
"4547110", "2029216", "981661", "508897", "280051", "162331",
"98469", "62168", "40663", "27446", "19053", "13563", "9877",
"7341", "5558", "4280", "3347", "2654", "2132", "1733", "1424",
"1182", "990", "837", "713", "612", "529", "460", "402", "354",
"313", "279", "249", "223", "201", "182", "165", "150", "138",
"126", "116", "107", "99", "92", "85", "79", "74", "69", "65", "61",
"57", "54", "51", "48", "46", "43", "41", "39", "37", "36", "34",
"32", "31", "30", "28", "27", "26", "25", "24", "23", "23", "22",
"21", "20", "20", "19", "18", "18", "17", "17", "16", "16", "15"
}, {
"8338442683973420410660145045849076963795",
"872596990706967613912664152945", "896707843885562730147307",
"91315073695274540969", "128539440806486007", "934129001105825",
"20278149285734", "946946589774", "77191347471", "9555892093",
"1631391010", "358523975", "96431070", "30563524", "11089126",
"4503126", "2010616", "973111", "504675", "277833", "161100",
"97754", "61734", "40390", "27269", "18934", "13482", "9819",
"7300", "5528", "4257", "3330", "2641", "2122", "1725", "1417",
"1177", "986", "833", "710", "609", "527", "458", "401", "353",
"312", "278", "248", "222", "200", "181", "165", "150", "137",
"126", "116", "107", "99", "91", "85", "79", "74", "69", "65", "61",
"57", "54", "51", "48", "46", "43", "41", "39", "37", "35", "34",
"32", "31", "30", "28", "27", "26", "25", "24", "23", "22", "22",
"21", "20", "20", "19", "18", "18", "17", "17", "16", "16", "15"
}, {
"9122818552483814953977703257848970704164",
"933462289569511464780529972314", "946405863353935713909178",
"95513446972056321834", "133588658082928446",
"966158521967027", "20895030642048", "972833934108",
"79107381638", "9773098125", "1665590516", "365497822",
"98180628", "31083090", "11266459", "4571108", "2039360",
"986323", "511198", "281260", "163001", "98858",
"62404", "40811", "27543", "19117", "13608", "9908",
"7363", "5575", "4292", "3356", "2661", "2138",
"1737", "1428", "1185", "993", "839", "714", "613",
"530", "461", "403", "355", "314", "279", "249",
"224", "202", "182", "165", "151", "138", "126",
"116", "107", "99", "92", "85", "79", "74", "69",
"65", "61", "57", "54", "51", "48", "46", "43",
"41", "39", "37", "36", "34", "32", "31", "30",
"28", "27", "26", "25", "24", "23", "23", "22",
"21", "20", "20", "19", "18", "18", "17", "17",
"16", "16", "15"
}, {
"9151329724083804100369546479681933027521",
"935649419557299174433860420387", "948179413831316112751907",
"95662582675170358900", "133767426788182384",
"967289728859610", "20916775466497", "973745045600",
"79174731802", "9780725058", "1666790321", "365742295",
"98241919", "31101281", "11272665", "4573486", "2040365",
"986785", "511426", "281380", "163067", "98897",
"62427", "40826", "27552", "19124", "13612", "9911",
"7366", "5576", "4294", "3357", "2662", "2138",
"1738", "1428", "1185", "993", "839", "715", "613",
"530", "461", "403", "355", "314", "279", "249",
"224", "202", "182", "165", "151", "138", "126",
"116", "107", "99", "92", "85", "79", "74", "69",
"65", "61", "57", "54", "51", "48", "46", "43",
"41", "39", "37", "36", "34", "32", "31", "30",
"28", "27", "26", "25", "24", "23", "23", "22",
"21", "20", "20", "19", "18", "18", "17", "17",
"16", "16", "15"
}, {
"6839396355168045468586008471269923213531",
"752078770083218822016981965090", "796178899357307807726034",
"82700643015444840424", "118072966296549115",
"867224751770392", "18981881485802", "892288574037",
"73130030771", "9093989389", "1558462688", "343617470",
"92683740", "29448679", "10708016", "4356820", "1948676",
"944610", "490587", "270425", "156989", "95362",
"60284", "39477", "26675", "18536", "13208", "9627",
"7161", "5426", "4181", "3272", "2596", "2087",
"1697", "1395", "1159", "971", "821", "700", "601",
"520", "452", "396", "348", "308", "274", "245",
"220", "198", "179", "163", "148", "136", "124",
"114", "106", "98", "91", "84", "78", "73", "68",
"64", "60", "57", "53", "50", "48", "45", "43",
"41", "39", "37", "35", "34", "32", "31", "29",
"28", "27", "26", "25", "24", "23", "22", "22",
"21", "20", "19", "19", "18", "18", "17", "17",
"16", "16", "15"
}, {
"4788090721380022347683138981782307670424",
"575601315594614059890185238256", "642831903229558719812840",
"69196031110028430211", "101340693763170691",
"758683936560287", "16854690815260", "801767985909",
"66353290503", "8318415180", "1435359033", "318340531",
"86304307", "27544217", "10054988", "4105446", "1841996",
"895414", "466223", "257591", "149855", "91205",
"57758", "37886", "25639", "17842", "12730", "9290",
"6918", "5248", "4048", "3170", "2518", "2026",
"1649", "1357", "1128", "946", "800", "682", "586",
"507", "441", "387", "341", "302", "268", "240",
"215", "194", "176", "160", "146", "133", "122",
"112", "104", "96", "89", "83", "77", "72", "67",
"63", "59", "56", "53", "50", "47", "45", "42",
"40", "38", "36", "35", "33", "32", "30", "29",
"28", "27", "26", "25", "24", "23", "22", "21",
"21", "20", "19", "19", "18", "17", "17", "16",
"16", "15", "15"
}
};
if ((e = mp_init_multi(&a, &c, &r, NULL)) != MP_OKAY) {
return EXIT_FAILURE;
}
#ifdef MP_8BIT
for (i = 0; i < 1; i++) {
#else
for (i = 0; i < 10; i++) {
#endif
mp_read_radix(&a, input[i], 64);
#ifdef MP_8BIT
for (j = 3; j < 10; j++) {
#else
for (j = 3; j < 100; j++) {
#endif
mp_n_root(&a, (mp_digit) j, &c);
mp_read_radix(&r, root[i][j-3], 10);
if (mp_cmp(&r, &c) != MP_EQ) {
fprintf(stderr, "mp_n_root failed at input #%d, root #%d\n", i, j);
goto LTM_ERR;
}
}
}
mp_clear_multi(&a, &c, &r, NULL);
return EXIT_SUCCESS;
LTM_ERR:
mp_clear_multi(&a, &c, &r, NULL);
return EXIT_FAILURE;
}
int unit_tests(void)
{
static const struct {
@ -1387,6 +1629,7 @@ int unit_tests(void)
T(mp_read_radix),
T(mp_reduce_2k),
T(mp_reduce_2k_l),
T(mp_n_root),
T(mp_set_double),
T(mp_sqrt),
T(mp_sqrtmod_prime),

11
doc/bn.tex
View File

@ -1830,18 +1830,11 @@ It calculates $c = a \mod 2^b$.
\begin{alltt}
int mp_n_root (mp_int * a, mp_digit b, mp_int * c)
\end{alltt}
This computes $c = a^{1/b}$ such that $c^b \le a$ and $(c+1)^b > a$. The implementation of this function is not
ideal for values of $b$ greater than three. It will work but become very slow. So unless you are working with very small
numbers (less than 1000 bits) I'd avoid $b > 3$ situations. Will return a positive root only for even roots and return
This computes $c = a^{1/b}$ such that $c^b \le a$ and $(c+1)^b > a$. Will return a positive root only for even roots and return
a root with the sign of the input for odd roots. For example, performing $4^{1/2}$ will return $2$ whereas $(-8)^{1/3}$
will return $-2$.
This algorithm uses the ``Newton Approximation'' method and will converge on the correct root fairly quickly. Since
the algorithm requires raising $a$ to the power of $b$ it is not ideal to attempt to find roots for large
values of $b$. If particularly large roots are required then a factor method could be used instead. For example,
$a^{1/16}$ is equivalent to $\left (a^{1/4} \right)^{1/4}$ or simply
$\left ( \left ( \left ( a^{1/2} \right )^{1/2} \right )^{1/2} \right )^{1/2}$
This algorithm uses the ``Newton Approximation'' method and will converge on the correct root fairly quickly.
The square root $c = a^{1/2}$ (with the same conditions $c^2 \le a$ and $(c+1)^2 > a$) is implemented with a faster algorithm.

3
tommath_class.h
View File

@ -669,7 +669,9 @@
#if defined(BN_MP_N_ROOT_EX_C)
# define BN_MP_INIT_C
# define BN_MP_COUNT_BITS_C
# define BN_MP_SET_C
# define BN_MP_2EXPT_C
# define BN_MP_COPY_C
# define BN_MP_EXPT_D_EX_C
# define BN_MP_MUL_C
@ -677,6 +679,7 @@
# define BN_MP_MUL_D_C
# define BN_MP_DIV_C
# define BN_MP_CMP_C
# define BN_MP_ADD_D_C
# define BN_MP_SUB_D_C
# define BN_MP_EXCH_C
# define BN_MP_CLEAR_C