Mercurial > dropbear

comparison libtommath/bn_mp_div.c @ 1692:1051e4eea25a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update LibTomMath to 1.2.0 (#84) * update C files * update other files * update headers * update makefiles * remove mp_set/get_double() * use ltm 1.2.0 API * update ltm_desc * use bundled tommath if system-tommath is too old * XMALLOC etc. were changed to MP_MALLOC etc.
author Steffen Jaeckel <s@jaeckel.eu>
date Tue, 26 May 2020 17:36:47 +0200
parents f52919ffd3b1
children
comparison
equal deleted inserted replaced
1691:2d3745d58843 1692:1051e4eea25a
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 /* SPDX-License-Identifier: Unlicense */
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 * SPDX-License-Identifier: Unlicense
13 */
14 5
15 #ifdef BN_MP_DIV_SMALL 6 #ifdef BN_MP_DIV_SMALL
16 7
17 /* slower bit-bang division... also smaller */ 8 /* slower bit-bang division... also smaller */
18 int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d) 9 mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
19 { 10 {
20 mp_int ta, tb, tq, q; 11 mp_int ta, tb, tq, q;
21 int res, n, n2; 12 int n, n2;
13 mp_err err;
22 14
23 /* is divisor zero ? */ 15 /* is divisor zero ? */
24 if (mp_iszero(b) == MP_YES) { 16 if (MP_IS_ZERO(b)) {
25 return MP_VAL; 17 return MP_VAL;
26 } 18 }
27 19
28 /* if a < b then q=0, r = a */ 20 /* if a < b then q=0, r = a */
29 if (mp_cmp_mag(a, b) == MP_LT) { 21 if (mp_cmp_mag(a, b) == MP_LT) {
30 if (d != NULL) { 22 if (d != NULL) {
31 res = mp_copy(a, d); 23 err = mp_copy(a, d);
32 } else { 24 } else {
33 res = MP_OKAY; 25 err = MP_OKAY;
34 } 26 }
35 if (c != NULL) { 27 if (c != NULL) {
36 mp_zero(c); 28 mp_zero(c);
37 } 29 }
38 return res; 30 return err;
39 } 31 }
40 32
41 /* init our temps */ 33 /* init our temps */
42 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) { 34 if ((err = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) {
43 return res; 35 return err;
44 } 36 }
45 37
46 38
47 mp_set(&tq, 1uL); 39 mp_set(&tq, 1uL);
48 n = mp_count_bits(a) - mp_count_bits(b); 40 n = mp_count_bits(a) - mp_count_bits(b);
49 if (((res = mp_abs(a, &ta)) != MP_OKAY) || 41 if ((err = mp_abs(a, &ta)) != MP_OKAY) goto LBL_ERR;
50 ((res = mp_abs(b, &tb)) != MP_OKAY) || 42 if ((err = mp_abs(b, &tb)) != MP_OKAY) goto LBL_ERR;
51 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || 43 if ((err = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) goto LBL_ERR;
52 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { 44 if ((err = mp_mul_2d(&tq, n, &tq)) != MP_OKAY) goto LBL_ERR;
53 goto LBL_ERR;
54 }
55 45
56 while (n-- >= 0) { 46 while (n-- >= 0) {
57 if (mp_cmp(&tb, &ta) != MP_GT) { 47 if (mp_cmp(&tb, &ta) != MP_GT) {
58 if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) || 48 if ((err = mp_sub(&ta, &tb, &ta)) != MP_OKAY) goto LBL_ERR;
59 ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) { 49 if ((err = mp_add(&q, &tq, &q)) != MP_OKAY) goto LBL_ERR;
60 goto LBL_ERR; 50 }
61 } 51 if ((err = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) goto LBL_ERR;
62 } 52 if ((err = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY) goto LBL_ERR;
63 if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
64 ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
65 goto LBL_ERR;
66 }
67 } 53 }
68 54
69 /* now q == quotient and ta == remainder */ 55 /* now q == quotient and ta == remainder */
70 n = a->sign; 56 n = a->sign;
71 n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 57 n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
72 if (c != NULL) { 58 if (c != NULL) {
73 mp_exch(c, &q); 59 mp_exch(c, &q);
74 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; 60 c->sign = MP_IS_ZERO(c) ? MP_ZPOS : n2;
75 } 61 }
76 if (d != NULL) { 62 if (d != NULL) {
77 mp_exch(d, &ta); 63 mp_exch(d, &ta);
78 d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n; 64 d->sign = MP_IS_ZERO(d) ? MP_ZPOS : n;
79 } 65 }
80 LBL_ERR: 66 LBL_ERR:
81 mp_clear_multi(&ta, &tb, &tq, &q, NULL); 67 mp_clear_multi(&ta, &tb, &tq, &q, NULL);
82 return res; 68 return err;
83 } 69 }
84 70
85 #else 71 #else
86 72
87 /* integer signed division. 73 /* integer signed division.
95 * case that y has fewer than three digits, etc.. 81 * case that y has fewer than three digits, etc..
96 * 82 *
97 * The overall algorithm is as described as 83 * The overall algorithm is as described as
98 * 14.20 from HAC but fixed to treat these cases. 84 * 14.20 from HAC but fixed to treat these cases.
99 */ 85 */
100 int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d) 86 mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
101 { 87 {
102 mp_int q, x, y, t1, t2; 88 mp_int q, x, y, t1, t2;
103 int res, n, t, i, norm, neg; 89 int n, t, i, norm;
90 mp_sign neg;
91 mp_err err;
104 92
105 /* is divisor zero ? */ 93 /* is divisor zero ? */
106 if (mp_iszero(b) == MP_YES) { 94 if (MP_IS_ZERO(b)) {
107 return MP_VAL; 95 return MP_VAL;
108 } 96 }
109 97
110 /* if a < b then q=0, r = a */ 98 /* if a < b then q=0, r = a */
111 if (mp_cmp_mag(a, b) == MP_LT) { 99 if (mp_cmp_mag(a, b) == MP_LT) {
112 if (d != NULL) { 100 if (d != NULL) {
113 res = mp_copy(a, d); 101 err = mp_copy(a, d);
114 } else { 102 } else {
115 res = MP_OKAY; 103 err = MP_OKAY;
116 } 104 }
117 if (c != NULL) { 105 if (c != NULL) {
118 mp_zero(c); 106 mp_zero(c);
119 } 107 }
120 return res; 108 return err;
121 } 109 }
122 110
123 if ((res = mp_init_size(&q, a->used + 2)) != MP_OKAY) { 111 if ((err = mp_init_size(&q, a->used + 2)) != MP_OKAY) {
124 return res; 112 return err;
125 } 113 }
126 q.used = a->used + 2; 114 q.used = a->used + 2;
127 115
128 if ((res = mp_init(&t1)) != MP_OKAY) { 116 if ((err = mp_init(&t1)) != MP_OKAY) goto LBL_Q;
129 goto LBL_Q; 117
130 } 118 if ((err = mp_init(&t2)) != MP_OKAY) goto LBL_T1;
131 119
132 if ((res = mp_init(&t2)) != MP_OKAY) { 120 if ((err = mp_init_copy(&x, a)) != MP_OKAY) goto LBL_T2;
133 goto LBL_T1; 121
134 } 122 if ((err = mp_init_copy(&y, b)) != MP_OKAY) goto LBL_X;
135
136 if ((res = mp_init_copy(&x, a)) != MP_OKAY) {
137 goto LBL_T2;
138 }
139
140 if ((res = mp_init_copy(&y, b)) != MP_OKAY) {
141 goto LBL_X;
142 }
143 123
144 /* fix the sign */ 124 /* fix the sign */
145 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 125 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
146 x.sign = y.sign = MP_ZPOS; 126 x.sign = y.sign = MP_ZPOS;
147 127
148 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */ 128 /* normalize both x and y, ensure that y >= b/2, [b == 2**MP_DIGIT_BIT] */
149 norm = mp_count_bits(&y) % DIGIT_BIT; 129 norm = mp_count_bits(&y) % MP_DIGIT_BIT;
150 if (norm < (DIGIT_BIT - 1)) { 130 if (norm < (MP_DIGIT_BIT - 1)) {
151 norm = (DIGIT_BIT - 1) - norm; 131 norm = (MP_DIGIT_BIT - 1) - norm;
152 if ((res = mp_mul_2d(&x, norm, &x)) != MP_OKAY) { 132 if ((err = mp_mul_2d(&x, norm, &x)) != MP_OKAY) goto LBL_Y;
153 goto LBL_Y; 133 if ((err = mp_mul_2d(&y, norm, &y)) != MP_OKAY) goto LBL_Y;
154 }
155 if ((res = mp_mul_2d(&y, norm, &y)) != MP_OKAY) {
156 goto LBL_Y;
157 }
158 } else { 134 } else {
159 norm = 0; 135 norm = 0;
160 } 136 }
161 137
162 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */ 138 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
163 n = x.used - 1; 139 n = x.used - 1;
164 t = y.used - 1; 140 t = y.used - 1;
165 141
166 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */ 142 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
167 if ((res = mp_lshd(&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */ 143 /* y = y*b**{n-t} */
168 goto LBL_Y; 144 if ((err = mp_lshd(&y, n - t)) != MP_OKAY) goto LBL_Y;
169 }
170 145
171 while (mp_cmp(&x, &y) != MP_LT) { 146 while (mp_cmp(&x, &y) != MP_LT) {
172 ++(q.dp[n - t]); 147 ++(q.dp[n - t]);
173 if ((res = mp_sub(&x, &y, &x)) != MP_OKAY) { 148 if ((err = mp_sub(&x, &y, &x)) != MP_OKAY) goto LBL_Y;
174 goto LBL_Y;
175 }
176 } 149 }
177 150
178 /* reset y by shifting it back down */ 151 /* reset y by shifting it back down */
179 mp_rshd(&y, n - t); 152 mp_rshd(&y, n - t);
180 153
185 } 158 }
186 159
187 /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 160 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
188 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */ 161 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
189 if (x.dp[i] == y.dp[t]) { 162 if (x.dp[i] == y.dp[t]) {
190 q.dp[(i - t) - 1] = ((mp_digit)1 << (mp_digit)DIGIT_BIT) - (mp_digit)1; 163 q.dp[(i - t) - 1] = ((mp_digit)1 << (mp_digit)MP_DIGIT_BIT) - (mp_digit)1;
191 } else { 164 } else {
192 mp_word tmp; 165 mp_word tmp;
193 tmp = (mp_word)x.dp[i] << (mp_word)DIGIT_BIT; 166 tmp = (mp_word)x.dp[i] << (mp_word)MP_DIGIT_BIT;
194 tmp |= (mp_word)x.dp[i - 1]; 167 tmp |= (mp_word)x.dp[i - 1];
195 tmp /= (mp_word)y.dp[t]; 168 tmp /= (mp_word)y.dp[t];
196 if (tmp > (mp_word)MP_MASK) { 169 if (tmp > (mp_word)MP_MASK) {
197 tmp = MP_MASK; 170 tmp = MP_MASK;
198 } 171 }
211 /* find left hand */ 184 /* find left hand */
212 mp_zero(&t1); 185 mp_zero(&t1);
213 t1.dp[0] = ((t - 1) < 0) ? 0u : y.dp[t - 1]; 186 t1.dp[0] = ((t - 1) < 0) ? 0u : y.dp[t - 1];
214 t1.dp[1] = y.dp[t]; 187 t1.dp[1] = y.dp[t];
215 t1.used = 2; 188 t1.used = 2;
216 if ((res = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) { 189 if ((err = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) goto LBL_Y;
217 goto LBL_Y;
218 }
219 190
220 /* find right hand */ 191 /* find right hand */
221 t2.dp[0] = ((i - 2) < 0) ? 0u : x.dp[i - 2]; 192 t2.dp[0] = ((i - 2) < 0) ? 0u : x.dp[i - 2];
222 t2.dp[1] = ((i - 1) < 0) ? 0u : x.dp[i - 1]; 193 t2.dp[1] = x.dp[i - 1]; /* i >= 1 always holds */
223 t2.dp[2] = x.dp[i]; 194 t2.dp[2] = x.dp[i];
224 t2.used = 3; 195 t2.used = 3;
225 } while (mp_cmp_mag(&t1, &t2) == MP_GT); 196 } while (mp_cmp_mag(&t1, &t2) == MP_GT);
226 197
227 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */ 198 /* 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) { 199 if ((err = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) goto LBL_Y;
229 goto LBL_Y; 200
230 } 201 if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) goto LBL_Y;
231 202
232 if ((res = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) { 203 if ((err = mp_sub(&x, &t1, &x)) != MP_OKAY) goto LBL_Y;
233 goto LBL_Y;
234 }
235
236 if ((res = mp_sub(&x, &t1, &x)) != MP_OKAY) {
237 goto LBL_Y;
238 }
239 204
240 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */ 205 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
241 if (x.sign == MP_NEG) { 206 if (x.sign == MP_NEG) {
242 if ((res = mp_copy(&y, &t1)) != MP_OKAY) { 207 if ((err = mp_copy(&y, &t1)) != MP_OKAY) goto LBL_Y;
243 goto LBL_Y; 208 if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) goto LBL_Y;
244 } 209 if ((err = mp_add(&x, &t1, &x)) != MP_OKAY) goto LBL_Y;
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 210
252 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & MP_MASK; 211 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & MP_MASK;
253 } 212 }
254 } 213 }
255 214
265 mp_exch(&q, c); 224 mp_exch(&q, c);
266 c->sign = neg; 225 c->sign = neg;
267 } 226 }
268 227
269 if (d != NULL) { 228 if (d != NULL) {
270 if ((res = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY) { 229 if ((err = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY) goto LBL_Y;
271 goto LBL_Y;
272 }
273 mp_exch(&x, d); 230 mp_exch(&x, d);
274 } 231 }
275 232
276 res = MP_OKAY; 233 err = MP_OKAY;
277 234
278 LBL_Y: 235 LBL_Y:
279 mp_clear(&y); 236 mp_clear(&y);
280 LBL_X: 237 LBL_X:
281 mp_clear(&x); 238 mp_clear(&x);
283 mp_clear(&t2); 240 mp_clear(&t2);
284 LBL_T1: 241 LBL_T1:
285 mp_clear(&t1); 242 mp_clear(&t1);
286 LBL_Q: 243 LBL_Q:
287 mp_clear(&q); 244 mp_clear(&q);
288 return res; 245 return err;
289 } 246 }
290 247
291 #endif 248 #endif
292 249
293 #endif 250 #endif
294
295 /* ref: HEAD -> master, tag: v1.1.0 */
296 /* git commit: 08549ad6bc8b0cede0b357a9c341c5c6473a9c55 */
297 /* commit time: 2019-01-28 20:32:32 +0100 */