diff --git a/bn_mp_n_root_ex.c b/bn_mp_n_root_ex.c index ebc08ba..4efa7aa 100644 --- a/bn_mp_n_root_ex.c +++ b/bn_mp_n_root_ex.c @@ -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$ */ diff --git a/callgraph.txt b/callgraph.txt index 2a9aed0..790d841 100644 --- a/callgraph.txt +++ b/callgraph.txt @@ -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 diff --git a/demo/test.c b/demo/test.c index d720802..3cbc8d2 100644 --- a/demo/test.c +++ b/demo/test.c @@ -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), diff --git a/doc/bn.tex b/doc/bn.tex index 2213b2d..c83b2e3 100644 --- a/doc/bn.tex +++ b/doc/bn.tex @@ -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. diff --git a/tommath_class.h b/tommath_class.h index 96b2ba0..d529477 100644 --- a/tommath_class.h +++ b/tommath_class.h @@ -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