Mercurial > dropbear

comparison libtommath/bn_mp_n_root_ex.c @ 1655:f52919ffd3b1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update ltm to 1.1.0 and enable FIPS 186.4 compliant key-generation (#79) * make key-generation compliant to FIPS 186.4 * fix includes in tommath_class.h * update fuzzcorpus instead of error-out * fixup fuzzing make-targets * update Makefile.in * apply necessary patches to ltm sources * clean-up not required ltm files * update to vanilla ltm 1.1.0 this already only contains the required files * remove set/get double
author Steffen Jaeckel <s_jaeckel@gmx.de>
date Mon, 16 Sep 2019 15:50:38 +0200
parents 8bba51a55704
children
comparison
equal deleted inserted replaced
1654:cc0fc5131c5c 1655:f52919ffd3b1
1 #include <tommath_private.h> 1 #include "tommath_private.h"
2 #ifdef BN_MP_N_ROOT_EX_C 2 #ifdef BN_MP_N_ROOT_EX_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis 3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4 * 4 *
5 * LibTomMath is a library that provides multiple-precision 5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality. 6 * integer arithmetic as well as number theoretic functionality.
7 * 7 *
8 * The library was designed directly after the MPI library by 8 * The library was designed directly after the MPI library by
9 * Michael Fromberger but has been written from scratch with 9 * Michael Fromberger but has been written from scratch with
10 * additional optimizations in place. 10 * additional optimizations in place.
11 * 11 *
12 * The library is free for all purposes without any express 12 * SPDX-License-Identifier: Unlicense
13 * guarantee it works.
14 *
15 * Tom St Denis, [email protected], http://libtom.org
16 */ 13 */
17 14
18 /* find the n'th root of an integer 15 /* find the n'th root of an integer
19 * 16 *
20 * Result found such that (c)**b <= a and (c+1)**b > a 17 * Result found such that (c)**b <= a and (c+1)**b > a
23 * x[i+1] = x[i] - f(x[i])/f'(x[i]) 20 * x[i+1] = x[i] - f(x[i])/f'(x[i])
24 * which will find the root in log(N) time where 21 * which will find the root in log(N) time where
25 * each step involves a fair bit. This is not meant to 22 * each step involves a fair bit. This is not meant to
26 * find huge roots [square and cube, etc]. 23 * find huge roots [square and cube, etc].
27 */ 24 */
28 int mp_n_root_ex (mp_int * a, mp_digit b, mp_int * c, int fast) 25 int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
29 { 26 {
30 mp_int t1, t2, t3; 27 mp_int t1, t2, t3, a_;
31 int res, neg; 28 int res;
32 29
33 /* input must be positive if b is even */ 30 /* input must be positive if b is even */
34 if (((b & 1) == 0) && (a->sign == MP_NEG)) { 31 if (((b & 1u) == 0u) && (a->sign == MP_NEG)) {
35 return MP_VAL; 32 return MP_VAL;
36 } 33 }
37 34
38 if ((res = mp_init (&t1)) != MP_OKAY) { 35 if ((res = mp_init(&t1)) != MP_OKAY) {
39 return res; 36 return res;
40 } 37 }
41 38
42 if ((res = mp_init (&t2)) != MP_OKAY) { 39 if ((res = mp_init(&t2)) != MP_OKAY) {
43 goto LBL_T1; 40 goto LBL_T1;
44 } 41 }
45 42
46 if ((res = mp_init (&t3)) != MP_OKAY) { 43 if ((res = mp_init(&t3)) != MP_OKAY) {
47 goto LBL_T2; 44 goto LBL_T2;
48 } 45 }
49 46
50 /* if a is negative fudge the sign but keep track */ 47 /* if a is negative fudge the sign but keep track */
51 neg = a->sign; 48 a_ = *a;
52 a->sign = MP_ZPOS; 49 a_.sign = MP_ZPOS;
53 50
54 /* t2 = 2 */ 51 /* t2 = 2 */
55 mp_set (&t2, 2); 52 mp_set(&t2, 2uL);
56 53
57 do { 54 do {
58 /* t1 = t2 */ 55 /* t1 = t2 */
59 if ((res = mp_copy (&t2, &t1)) != MP_OKAY) { 56 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; 57 goto LBL_T3;
106 } 58 }
107 } else {
108 break;
109 }
110 }
111 59
112 /* reset the sign of a first */ 60 /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
113 a->sign = neg;
114 61
115 /* set the result */ 62 /* t3 = t1**(b-1) */
116 mp_exch (&t1, c); 63 if ((res = mp_expt_d_ex(&t1, b - 1u, &t3, fast)) != MP_OKAY) {
64 goto LBL_T3;
65 }
117 66
118 /* set the sign of the result */ 67 /* numerator */
119 c->sign = neg; 68 /* t2 = t1**b */
69 if ((res = mp_mul(&t3, &t1, &t2)) != MP_OKAY) {
70 goto LBL_T3;
71 }
120 72
121 res = MP_OKAY; 73 /* t2 = t1**b - a */
74 if ((res = mp_sub(&t2, &a_, &t2)) != MP_OKAY) {
75 goto LBL_T3;
76 }
122 77
123 LBL_T3:mp_clear (&t3); 78 /* denominator */
124 LBL_T2:mp_clear (&t2); 79 /* t3 = t1**(b-1) * b */
125 LBL_T1:mp_clear (&t1); 80 if ((res = mp_mul_d(&t3, b, &t3)) != MP_OKAY) {
126 return res; 81 goto LBL_T3;
82 }
83
84 /* t3 = (t1**b - a)/(b * t1**(b-1)) */
85 if ((res = mp_div(&t2, &t3, &t3, NULL)) != MP_OKAY) {
86 goto LBL_T3;
87 }
88
89 if ((res = mp_sub(&t1, &t3, &t2)) != MP_OKAY) {
90 goto LBL_T3;
91 }
92 } while (mp_cmp(&t1, &t2) != MP_EQ);
93
94 /* result can be off by a few so check */
95 for (;;) {
96 if ((res = mp_expt_d_ex(&t1, b, &t2, fast)) != MP_OKAY) {
97 goto LBL_T3;
98 }
99
100 if (mp_cmp(&t2, &a_) == MP_GT) {
101 if ((res = mp_sub_d(&t1, 1uL, &t1)) != MP_OKAY) {
102 goto LBL_T3;
103 }
104 } else {
105 break;
106 }
107 }
108
109 /* set the result */
110 mp_exch(&t1, c);
111
112 /* set the sign of the result */
113 c->sign = a->sign;
114
115 res = MP_OKAY;
116
117 LBL_T3:
118 mp_clear(&t3);
119 LBL_T2:
120 mp_clear(&t2);
121 LBL_T1:
122 mp_clear(&t1);
123 return res;
127 } 124 }
128 #endif 125 #endif
129 126
130 /* ref: $Format:%D$ */ 127 /* ref: HEAD -> master, tag: v1.1.0 */
131 /* git commit: $Format:%H$ */ 128 /* git commit: 08549ad6bc8b0cede0b357a9c341c5c6473a9c55 */
132 /* commit time: $Format:%ai$ */ 129 /* commit time: 2019-01-28 20:32:32 +0100 */