Mercurial > dropbear
annotate libtommath/bn_mp_n_root_ex.c @ 1604:99ca393afc56
It turns out you can't have a single-quote in an #error
| author | Matt Johnston <matt@ucc.asn.au> |
|---|---|
| date | Thu, 08 Mar 2018 22:37:54 +0800 |
| parents | 8bba51a55704 |
| children | f52919ffd3b1 |
| rev | line source |
|---|---|
| 1436 | 1 #include <tommath_private.h> |
| 2 #ifdef BN_MP_N_ROOT_EX_C | |
| 3 /* LibTomMath, multiple-precision integer library -- Tom St Denis | |
| 4 * | |
| 5 * LibTomMath is a library that provides multiple-precision | |
| 6 * integer arithmetic as well as number theoretic functionality. | |
| 7 * | |
| 8 * The library was designed directly after the MPI library by | |
| 9 * Michael Fromberger but has been written from scratch with | |
| 10 * additional optimizations in place. | |
| 11 * | |
| 12 * The library is free for all purposes without any express | |
| 13 * guarantee it works. | |
| 14 * | |
| 15 * Tom St Denis, [email protected], http://libtom.org | |
| 16 */ | |
| 17 | |
| 18 /* find the n'th root of an integer | |
| 19 * | |
| 20 * Result found such that (c)**b <= a and (c+1)**b > a | |
| 21 * | |
| 22 * This algorithm uses Newton's approximation | |
| 23 * x[i+1] = x[i] - f(x[i])/f'(x[i]) | |
| 24 * which will find the root in log(N) time where | |
| 25 * each step involves a fair bit. This is not meant to | |
| 26 * find huge roots [square and cube, etc]. | |
| 27 */ | |
| 28 int mp_n_root_ex (mp_int * a, mp_digit b, mp_int * c, int fast) | |
| 29 { | |
| 30 mp_int t1, t2, t3; | |
| 31 int res, neg; | |
| 32 | |
| 33 /* input must be positive if b is even */ | |
| 34 if (((b & 1) == 0) && (a->sign == MP_NEG)) { | |
| 35 return MP_VAL; | |
| 36 } | |
| 37 | |
| 38 if ((res = mp_init (&t1)) != MP_OKAY) { | |
| 39 return res; | |
| 40 } | |
| 41 | |
| 42 if ((res = mp_init (&t2)) != MP_OKAY) { | |
| 43 goto LBL_T1; | |
| 44 } | |
| 45 | |
| 46 if ((res = mp_init (&t3)) != MP_OKAY) { | |
| 47 goto LBL_T2; | |
| 48 } | |
| 49 | |
| 50 /* if a is negative fudge the sign but keep track */ | |
| 51 neg = a->sign; | |
| 52 a->sign = MP_ZPOS; | |
| 53 | |
| 54 /* t2 = 2 */ | |
| 55 mp_set (&t2, 2); | |
| 56 | |
| 57 do { | |
| 58 /* t1 = t2 */ | |
| 59 if ((res = mp_copy (&t2, &t1)) != MP_OKAY) { | |
| 60 goto LBL_T3; | |
| 61 } | |
| 62 | |
| 63 /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */ | |
| 64 | |
| 65 /* t3 = t1**(b-1) */ | |
| 66 if ((res = mp_expt_d_ex (&t1, b - 1, &t3, fast)) != MP_OKAY) { | |
| 67 goto LBL_T3; | |
| 68 } | |
| 69 | |
| 70 /* numerator */ | |
| 71 /* t2 = t1**b */ | |
| 72 if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) { | |
| 73 goto LBL_T3; | |
| 74 } | |
| 75 | |
| 76 /* t2 = t1**b - a */ | |
| 77 if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) { | |
| 78 goto LBL_T3; | |
| 79 } | |
| 80 | |
| 81 /* denominator */ | |
| 82 /* t3 = t1**(b-1) * b */ | |
| 83 if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) { | |
| 84 goto LBL_T3; | |
| 85 } | |
| 86 | |
| 87 /* t3 = (t1**b - a)/(b * t1**(b-1)) */ | |
| 88 if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) { | |
| 89 goto LBL_T3; | |
| 90 } | |
| 91 | |
| 92 if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) { | |
| 93 goto LBL_T3; | |
| 94 } | |
| 95 } while (mp_cmp (&t1, &t2) != MP_EQ); | |
| 96 | |
| 97 /* result can be off by a few so check */ | |
| 98 for (;;) { | |
| 99 if ((res = mp_expt_d_ex (&t1, b, &t2, fast)) != MP_OKAY) { | |
| 100 goto LBL_T3; | |
| 101 } | |
| 102 | |
| 103 if (mp_cmp (&t2, a) == MP_GT) { | |
| 104 if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) { | |
| 105 goto LBL_T3; | |
| 106 } | |
| 107 } else { | |
| 108 break; | |
| 109 } | |
| 110 } | |
| 111 | |
| 112 /* reset the sign of a first */ | |
| 113 a->sign = neg; | |
| 114 | |
| 115 /* set the result */ | |
| 116 mp_exch (&t1, c); | |
| 117 | |
| 118 /* set the sign of the result */ | |
| 119 c->sign = neg; | |
| 120 | |
| 121 res = MP_OKAY; | |
| 122 | |
| 123 LBL_T3:mp_clear (&t3); | |
| 124 LBL_T2:mp_clear (&t2); | |
| 125 LBL_T1:mp_clear (&t1); | |
| 126 return res; | |
| 127 } | |
| 128 #endif | |
| 129 | |
|
1470
8bba51a55704
Update to libtommath v1.0.1
Matt Johnston <matt@ucc.asn.au>
parents:
1436
diff
changeset
|
130 /* ref: $Format:%D$ */ |
|
8bba51a55704
Update to libtommath v1.0.1
Matt Johnston <matt@ucc.asn.au>
parents:
1436
diff
changeset
|
131 /* git commit: $Format:%H$ */ |
|
8bba51a55704
Update to libtommath v1.0.1
Matt Johnston <matt@ucc.asn.au>
parents:
1436
diff
changeset
|
132 /* commit time: $Format:%ai$ */ |