Mercurial > dropbear
annotate bn_mp_n_root.c @ 200:c5c969ed76f3 libtommath
propagate from branch 'au.asn.ucc.matt.ltm-orig' (head 7fa10cba9535de3461cedb14b877c24858826204)
to branch 'au.asn.ucc.matt.dropbear.ltm' (head fc26f60de0370ab0a281fa41a2d13fb17c9d90a8)
| author | Matt Johnston <matt@ucc.asn.au> |
|---|---|
| date | Wed, 11 May 2005 16:15:27 +0000 |
| parents | d8254fc979e9 |
| children |
| rev | line source |
|---|---|
| 142 | 1 #include <tommath.h> |
| 2 #ifdef BN_MP_N_ROOT_C | |
| 2 | 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://math.libtomcrypt.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 (mp_int * a, mp_digit b, mp_int * c) | |
| 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
43 goto LBL_T1; |
| 2 | 44 } |
| 45 | |
| 46 if ((res = mp_init (&t3)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
47 goto LBL_T2; |
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
60 goto LBL_T3; |
| 2 | 61 } |
| 62 | |
| 63 /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */ | |
| 64 | |
| 65 /* t3 = t1**(b-1) */ | |
| 66 if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
67 goto LBL_T3; |
| 2 | 68 } |
| 69 | |
| 70 /* numerator */ | |
| 71 /* t2 = t1**b */ | |
| 72 if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
73 goto LBL_T3; |
| 2 | 74 } |
| 75 | |
| 76 /* t2 = t1**b - a */ | |
| 77 if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
78 goto LBL_T3; |
| 2 | 79 } |
| 80 | |
| 81 /* denominator */ | |
| 82 /* t3 = t1**(b-1) * b */ | |
| 83 if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
84 goto LBL_T3; |
| 2 | 85 } |
| 86 | |
| 87 /* t3 = (t1**b - a)/(b * t1**(b-1)) */ | |
| 88 if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
89 goto LBL_T3; |
| 2 | 90 } |
| 91 | |
| 92 if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
93 goto LBL_T3; |
| 2 | 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 (&t1, b, &t2)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
100 goto LBL_T3; |
| 2 | 101 } |
| 102 | |
| 103 if (mp_cmp (&t2, a) == MP_GT) { | |
| 104 if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
105 goto LBL_T3; |
| 2 | 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 | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
123 LBL_T3:mp_clear (&t3); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
124 LBL_T2:mp_clear (&t2); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
125 LBL_T1:mp_clear (&t1); |
| 2 | 126 return res; |
| 127 } | |
| 142 | 128 #endif |