bn_mp_n_root_ex.c 2.93 KB
 Dominique Dumont committed Jan 06, 2019 1 ``````#include "tommath_private.h" `````` Dominique Dumont committed Feb 29, 2016 2 3 4 5 6 7 8 9 10 11 ``````#ifdef BN_MP_N_ROOT_EX_C /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * Michael Fromberger but has been written from scratch with * additional optimizations in place. * `````` Dominique Dumont committed Jan 06, 2019 12 `````` * SPDX-License-Identifier: Unlicense `````` Dominique Dumont committed Feb 29, 2016 13 14 15 16 17 18 19 20 21 22 23 24 `````` */ /* find the n'th root of an integer * * Result found such that (c)**b <= a and (c+1)**b > a * * 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]. */ `````` Dominique Dumont committed Jan 06, 2019 25 ``````int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast) `````` Dominique Dumont committed Feb 29, 2016 26 ``````{ `````` Dominique Dumont committed Jan 06, 2019 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 `````` mp_int t1, t2, t3, a_; int res; /* input must be positive if b is even */ if (((b & 1u) == 0u) && (a->sign == MP_NEG)) { return MP_VAL; } if ((res = mp_init(&t1)) != MP_OKAY) { return res; } if ((res = mp_init(&t2)) != MP_OKAY) { goto LBL_T1; } if ((res = mp_init(&t3)) != MP_OKAY) { goto LBL_T2; } /* if a is negative fudge the sign but keep track */ a_ = *a; a_.sign = MP_ZPOS; /* t2 = 2 */ mp_set(&t2, 2uL); do { /* t1 = t2 */ if ((res = mp_copy(&t2, &t1)) != MP_OKAY) { `````` Dominique Dumont committed Feb 29, 2016 57 58 59 `````` goto LBL_T3; } `````` Dominique Dumont committed Jan 06, 2019 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 `````` /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */ /* t3 = t1**(b-1) */ 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) { goto LBL_T3; } /* t2 = t1**b - a */ if ((res = mp_sub(&t2, &a_, &t2)) != MP_OKAY) { goto LBL_T3; } /* denominator */ /* t3 = t1**(b-1) * b */ if ((res = mp_mul_d(&t3, b, &t3)) != MP_OKAY) { goto LBL_T3; } /* t3 = (t1**b - a)/(b * t1**(b-1)) */ if ((res = mp_div(&t2, &t3, &t3, NULL)) != MP_OKAY) { goto LBL_T3; } if ((res = mp_sub(&t1, &t3, &t2)) != MP_OKAY) { goto LBL_T3; } } while (mp_cmp(&t1, &t2) != MP_EQ); /* result can be off by a few so check */ 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; } } else { break; } } `````` Dominique Dumont committed Feb 29, 2016 108 `````` `````` Dominique Dumont committed Jan 06, 2019 109 110 `````` /* set the result */ mp_exch(&t1, c); `````` Dominique Dumont committed Feb 29, 2016 111 `````` `````` Dominique Dumont committed Jan 06, 2019 112 113 `````` /* set the sign of the result */ c->sign = a->sign; `````` Dominique Dumont committed Feb 29, 2016 114 `````` `````` Dominique Dumont committed Jan 06, 2019 115 `````` res = MP_OKAY; `````` Dominique Dumont committed Feb 29, 2016 116 `````` `````` Dominique Dumont committed Jan 06, 2019 117 118 119 120 121 122 123 ``````LBL_T3: mp_clear(&t3); LBL_T2: mp_clear(&t2); LBL_T1: mp_clear(&t1); return res; `````` Dominique Dumont committed Feb 29, 2016 124 125 126 ``````} #endif `````` Dominique Dumont committed Jan 29, 2019 127 128 129 ``````/* ref: HEAD -> master, tag: v1.1.0 */ /* git commit: 08549ad6bc8b0cede0b357a9c341c5c6473a9c55 */ /* commit time: 2019-01-28 20:32:32 +0100 */``````