Mercurial > dropbear

comparison libtommath/bn_mp_div.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 1051e4eea25a
comparison
equal deleted inserted replaced
1654:cc0fc5131c5c 1655:f52919ffd3b1
1 #include <tommath_private.h> 1 #include "tommath_private.h"
2 #ifdef BN_MP_DIV_C 2 #ifdef BN_MP_DIV_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 #ifdef BN_MP_DIV_SMALL 15 #ifdef BN_MP_DIV_SMALL
19 16
20 /* slower bit-bang division... also smaller */ 17 /* slower bit-bang division... also smaller */
21 int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d) 18 int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
22 { 19 {
23 mp_int ta, tb, tq, q; 20 mp_int ta, tb, tq, q;
24 int res, n, n2; 21 int res, n, n2;
25 22
26 /* is divisor zero ? */ 23 /* is divisor zero ? */
27 if (mp_iszero (b) == MP_YES) { 24 if (mp_iszero(b) == MP_YES) {
28 return MP_VAL; 25 return MP_VAL;
29 } 26 }
30 27
31 /* if a < b then q=0, r = a */ 28 /* if a < b then q=0, r = a */
32 if (mp_cmp_mag (a, b) == MP_LT) { 29 if (mp_cmp_mag(a, b) == MP_LT) {
33 if (d != NULL) { 30 if (d != NULL) {
34 res = mp_copy (a, d); 31 res = mp_copy(a, d);
35 } else { 32 } else {
36 res = MP_OKAY; 33 res = MP_OKAY;
37 } 34 }
38 if (c != NULL) { 35 if (c != NULL) {
39 mp_zero (c); 36 mp_zero(c);
40 } 37 }
41 return res; 38 return res;
42 } 39 }
43 40
44 /* init our temps */ 41 /* init our temps */
45 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) { 42 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) {
46 return res; 43 return res;
47 } 44 }
48 45
49 46
50 mp_set(&tq, 1); 47 mp_set(&tq, 1uL);
51 n = mp_count_bits(a) - mp_count_bits(b); 48 n = mp_count_bits(a) - mp_count_bits(b);
52 if (((res = mp_abs(a, &ta)) != MP_OKAY) || 49 if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
53 ((res = mp_abs(b, &tb)) != MP_OKAY) || 50 ((res = mp_abs(b, &tb)) != MP_OKAY) ||
54 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || 51 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
55 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { 52 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
56 goto LBL_ERR; 53 goto LBL_ERR;
57 } 54 }
58 55
59 while (n-- >= 0) { 56 while (n-- >= 0) {
60 if (mp_cmp(&tb, &ta) != MP_GT) { 57 if (mp_cmp(&tb, &ta) != MP_GT) {
61 if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) || 58 if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
62 ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) { 59 ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
63 goto LBL_ERR; 60 goto LBL_ERR;
64 } 61 }
65 } 62 }
66 if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) || 63 if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
67 ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) { 64 ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
68 goto LBL_ERR; 65 goto LBL_ERR;
69 } 66 }
70 } 67 }
71 68
72 /* now q == quotient and ta == remainder */ 69 /* now q == quotient and ta == remainder */
73 n = a->sign; 70 n = a->sign;
74 n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 71 n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
75 if (c != NULL) { 72 if (c != NULL) {
76 mp_exch(c, &q); 73 mp_exch(c, &q);
77 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; 74 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
78 } 75 }
79 if (d != NULL) { 76 if (d != NULL) {
80 mp_exch(d, &ta); 77 mp_exch(d, &ta);
81 d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n; 78 d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
82 } 79 }
83 LBL_ERR: 80 LBL_ERR:
84 mp_clear_multi(&ta, &tb, &tq, &q, NULL); 81 mp_clear_multi(&ta, &tb, &tq, &q, NULL);
85 return res; 82 return res;
86 } 83 }
87 84
98 * case that y has fewer than three digits, etc.. 95 * case that y has fewer than three digits, etc..
99 * 96 *
100 * The overall algorithm is as described as 97 * The overall algorithm is as described as
101 * 14.20 from HAC but fixed to treat these cases. 98 * 14.20 from HAC but fixed to treat these cases.
102 */ 99 */
103 int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d) 100 int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
104 { 101 {
105 mp_int q, x, y, t1, t2; 102 mp_int q, x, y, t1, t2;
106 int res, n, t, i, norm, neg; 103 int res, n, t, i, norm, neg;
107 104
108 /* is divisor zero ? */ 105 /* is divisor zero ? */
109 if (mp_iszero (b) == MP_YES) { 106 if (mp_iszero(b) == MP_YES) {
110 return MP_VAL; 107 return MP_VAL;
111 } 108 }
112 109
113 /* if a < b then q=0, r = a */ 110 /* if a < b then q=0, r = a */
114 if (mp_cmp_mag (a, b) == MP_LT) { 111 if (mp_cmp_mag(a, b) == MP_LT) {
115 if (d != NULL) { 112 if (d != NULL) {
116 res = mp_copy (a, d); 113 res = mp_copy(a, d);
117 } else { 114 } else {
118 res = MP_OKAY; 115 res = MP_OKAY;
119 } 116 }
120 if (c != NULL) { 117 if (c != NULL) {
121 mp_zero (c); 118 mp_zero(c);
122 } 119 }
123 return res; 120 return res;
124 } 121 }
125 122
126 if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) { 123 if ((res = mp_init_size(&q, a->used + 2)) != MP_OKAY) {
127 return res; 124 return res;
128 } 125 }
129 q.used = a->used + 2; 126 q.used = a->used + 2;
130 127
131 if ((res = mp_init (&t1)) != MP_OKAY) { 128 if ((res = mp_init(&t1)) != MP_OKAY) {
132 goto LBL_Q; 129 goto LBL_Q;
133 } 130 }
134 131
135 if ((res = mp_init (&t2)) != MP_OKAY) { 132 if ((res = mp_init(&t2)) != MP_OKAY) {
136 goto LBL_T1; 133 goto LBL_T1;
137 } 134 }
138 135
139 if ((res = mp_init_copy (&x, a)) != MP_OKAY) { 136 if ((res = mp_init_copy(&x, a)) != MP_OKAY) {
140 goto LBL_T2; 137 goto LBL_T2;
141 } 138 }
142 139
143 if ((res = mp_init_copy (&y, b)) != MP_OKAY) { 140 if ((res = mp_init_copy(&y, b)) != MP_OKAY) {
144 goto LBL_X; 141 goto LBL_X;
145 } 142 }
146 143
147 /* fix the sign */ 144 /* fix the sign */
148 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 145 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
149 x.sign = y.sign = MP_ZPOS; 146 x.sign = y.sign = MP_ZPOS;
150 147
151 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */ 148 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
152 norm = mp_count_bits(&y) % DIGIT_BIT; 149 norm = mp_count_bits(&y) % DIGIT_BIT;
153 if (norm < (int)(DIGIT_BIT-1)) { 150 if (norm < (DIGIT_BIT - 1)) {
154 norm = (DIGIT_BIT-1) - norm; 151 norm = (DIGIT_BIT - 1) - norm;
155 if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) { 152 if ((res = mp_mul_2d(&x, norm, &x)) != MP_OKAY) {
156 goto LBL_Y; 153 goto LBL_Y;
157 } 154 }
158 if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) { 155 if ((res = mp_mul_2d(&y, norm, &y)) != MP_OKAY) {
159 goto LBL_Y; 156 goto LBL_Y;
160 } 157 }
161 } else { 158 } else {
162 norm = 0; 159 norm = 0;
163 } 160 }
164 161
165 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */ 162 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
166 n = x.used - 1; 163 n = x.used - 1;
167 t = y.used - 1; 164 t = y.used - 1;
168 165
169 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */ 166 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
170 if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */ 167 if ((res = mp_lshd(&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
171 goto LBL_Y;
172 }
173
174 while (mp_cmp (&x, &y) != MP_LT) {
175 ++(q.dp[n - t]);
176 if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {
177 goto LBL_Y; 168 goto LBL_Y;
178 } 169 }
179 } 170
180 171 while (mp_cmp(&x, &y) != MP_LT) {
181 /* reset y by shifting it back down */ 172 ++(q.dp[n - t]);
182 mp_rshd (&y, n - t); 173 if ((res = mp_sub(&x, &y, &x)) != MP_OKAY) {
183 174 goto LBL_Y;
184 /* step 3. for i from n down to (t + 1) */ 175 }
185 for (i = n; i >= (t + 1); i--) { 176 }
186 if (i > x.used) { 177
187 continue; 178 /* reset y by shifting it back down */
188 } 179 mp_rshd(&y, n - t);
189 180
190 /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 181 /* step 3. for i from n down to (t + 1) */
191 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */ 182 for (i = n; i >= (t + 1); i--) {
192 if (x.dp[i] == y.dp[t]) { 183 if (i > x.used) {
193 q.dp[(i - t) - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1); 184 continue;
194 } else { 185 }
195 mp_word tmp; 186
196 tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT); 187 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
197 tmp |= ((mp_word) x.dp[i - 1]); 188 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
198 tmp /= ((mp_word) y.dp[t]); 189 if (x.dp[i] == y.dp[t]) {
199 if (tmp > (mp_word) MP_MASK) { 190 q.dp[(i - t) - 1] = ((mp_digit)1 << (mp_digit)DIGIT_BIT) - (mp_digit)1;
200 tmp = MP_MASK; 191 } else {
201 } 192 mp_word tmp;
202 q.dp[(i - t) - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK)); 193 tmp = (mp_word)x.dp[i] << (mp_word)DIGIT_BIT;
203 } 194 tmp |= (mp_word)x.dp[i - 1];
204 195 tmp /= (mp_word)y.dp[t];
205 /* while (q{i-t-1} * (yt * b + y{t-1})) > 196 if (tmp > (mp_word)MP_MASK) {
206 xi * b**2 + xi-1 * b + xi-2 197 tmp = MP_MASK;
207 198 }
208 do q{i-t-1} -= 1; 199 q.dp[(i - t) - 1] = (mp_digit)(tmp & (mp_word)MP_MASK);
200 }
201
202 /* while (q{i-t-1} * (yt * b + y{t-1})) >
203 xi * b**2 + xi-1 * b + xi-2
204
205 do q{i-t-1} -= 1;
206 */
207 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1uL) & (mp_digit)MP_MASK;
208 do {
209 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & (mp_digit)MP_MASK;
210
211 /* find left hand */
212 mp_zero(&t1);
213 t1.dp[0] = ((t - 1) < 0) ? 0u : y.dp[t - 1];
214 t1.dp[1] = y.dp[t];
215 t1.used = 2;
216 if ((res = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
217 goto LBL_Y;
218 }
219
220 /* find right hand */
221 t2.dp[0] = ((i - 2) < 0) ? 0u : x.dp[i - 2];
222 t2.dp[1] = ((i - 1) < 0) ? 0u : x.dp[i - 1];
223 t2.dp[2] = x.dp[i];
224 t2.used = 3;
225 } while (mp_cmp_mag(&t1, &t2) == MP_GT);
226
227 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
228 if ((res = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
229 goto LBL_Y;
230 }
231
232 if ((res = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) {
233 goto LBL_Y;
234 }
235
236 if ((res = mp_sub(&x, &t1, &x)) != MP_OKAY) {
237 goto LBL_Y;
238 }
239
240 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
241 if (x.sign == MP_NEG) {
242 if ((res = mp_copy(&y, &t1)) != MP_OKAY) {
243 goto LBL_Y;
244 }
245 if ((res = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) {
246 goto LBL_Y;
247 }
248 if ((res = mp_add(&x, &t1, &x)) != MP_OKAY) {
249 goto LBL_Y;
250 }
251
252 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & MP_MASK;
253 }
254 }
255
256 /* now q is the quotient and x is the remainder
257 * [which we have to normalize]
209 */ 258 */
210 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1) & MP_MASK; 259
211 do { 260 /* get sign before writing to c */
212 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1) & MP_MASK; 261 x.sign = (x.used == 0) ? MP_ZPOS : a->sign;
213 262
214 /* find left hand */ 263 if (c != NULL) {
215 mp_zero (&t1); 264 mp_clamp(&q);
216 t1.dp[0] = ((t - 1) < 0) ? 0 : y.dp[t - 1]; 265 mp_exch(&q, c);
217 t1.dp[1] = y.dp[t]; 266 c->sign = neg;
218 t1.used = 2; 267 }
219 if ((res = mp_mul_d (&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) { 268
220 goto LBL_Y; 269 if (d != NULL) {
221 } 270 if ((res = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY) {
222 271 goto LBL_Y;
223 /* find right hand */ 272 }
224 t2.dp[0] = ((i - 2) < 0) ? 0 : x.dp[i - 2]; 273 mp_exch(&x, d);
225 t2.dp[1] = ((i - 1) < 0) ? 0 : x.dp[i - 1]; 274 }
226 t2.dp[2] = x.dp[i]; 275
227 t2.used = 3; 276 res = MP_OKAY;
228 } while (mp_cmp_mag(&t1, &t2) == MP_GT); 277
229 278 LBL_Y:
230 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */ 279 mp_clear(&y);
231 if ((res = mp_mul_d (&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) { 280 LBL_X:
232 goto LBL_Y; 281 mp_clear(&x);
233 } 282 LBL_T2:
234 283 mp_clear(&t2);
235 if ((res = mp_lshd (&t1, (i - t) - 1)) != MP_OKAY) { 284 LBL_T1:
236 goto LBL_Y; 285 mp_clear(&t1);
237 } 286 LBL_Q:
238 287 mp_clear(&q);
239 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) { 288 return res;
240 goto LBL_Y;
241 }
242
243 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
244 if (x.sign == MP_NEG) {
245 if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
246 goto LBL_Y;
247 }
248 if ((res = mp_lshd (&t1, (i - t) - 1)) != MP_OKAY) {
249 goto LBL_Y;
250 }
251 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
252 goto LBL_Y;
253 }
254
255 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1UL) & MP_MASK;
256 }
257 }
258
259 /* now q is the quotient and x is the remainder
260 * [which we have to normalize]
261 */
262
263 /* get sign before writing to c */
264 x.sign = (x.used == 0) ? MP_ZPOS : a->sign;
265
266 if (c != NULL) {
267 mp_clamp (&q);
268 mp_exch (&q, c);
269 c->sign = neg;
270 }
271
272 if (d != NULL) {
273 if ((res = mp_div_2d (&x, norm, &x, NULL)) != MP_OKAY) {
274 goto LBL_Y;
275 }
276 mp_exch (&x, d);
277 }
278
279 res = MP_OKAY;
280
281 LBL_Y:mp_clear (&y);
282 LBL_X:mp_clear (&x);
283 LBL_T2:mp_clear (&t2);
284 LBL_T1:mp_clear (&t1);
285 LBL_Q:mp_clear (&q);
286 return res;
287 } 289 }
288 290
289 #endif 291 #endif
290 292
291 #endif 293 #endif
292 294
293 /* ref: $Format:%D$ */ 295 /* ref: HEAD -> master, tag: v1.1.0 */
294 /* git commit: $Format:%H$ */ 296 /* git commit: 08549ad6bc8b0cede0b357a9c341c5c6473a9c55 */
295 /* commit time: $Format:%ai$ */ 297 /* commit time: 2019-01-28 20:32:32 +0100 */