1 // SPDX-License-Identifier: GPL-2.0-or-later
2 /* mpihelp-div.c - MPI helper functions
3 * Copyright (C) 1994, 1996 Free Software Foundation, Inc.
4 * Copyright (C) 1998, 1999 Free Software Foundation, Inc.
5 *
6 * This file is part of GnuPG.
7 *
8 * Note: This code is heavily based on the GNU MP Library.
9 * Actually it's the same code with only minor changes in the
10 * way the data is stored; this is to support the abstraction
11 * of an optional secure memory allocation which may be used
12 * to avoid revealing of sensitive data due to paging etc.
13 * The GNU MP Library itself is published under the LGPL;
14 * however I decided to publish this code under the plain GPL.
15 */
16
17 #include "mpi-internal.h"
18 #include "longlong.h"
19
20 #ifndef UMUL_TIME
21 #define UMUL_TIME 1
22 #endif
23 #ifndef UDIV_TIME
24 #define UDIV_TIME UMUL_TIME
25 #endif
26
27
28 mpi_limb_t
mpihelp_mod_1(mpi_ptr_t dividend_ptr,mpi_size_t dividend_size,mpi_limb_t divisor_limb)29 mpihelp_mod_1(mpi_ptr_t dividend_ptr, mpi_size_t dividend_size,
30 mpi_limb_t divisor_limb)
31 {
32 mpi_size_t i;
33 mpi_limb_t n1, n0, r;
34 mpi_limb_t dummy __maybe_unused;
35
36 /* Botch: Should this be handled at all? Rely on callers? */
37 if (!dividend_size)
38 return 0;
39
40 /* If multiplication is much faster than division, and the
41 * dividend is large, pre-invert the divisor, and use
42 * only multiplications in the inner loop.
43 *
44 * This test should be read:
45 * Does it ever help to use udiv_qrnnd_preinv?
46 * && Does what we save compensate for the inversion overhead?
47 */
48 if (UDIV_TIME > (2 * UMUL_TIME + 6)
49 && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME) {
50 int normalization_steps;
51
52 normalization_steps = count_leading_zeros(divisor_limb);
53 if (normalization_steps) {
54 mpi_limb_t divisor_limb_inverted;
55
56 divisor_limb <<= normalization_steps;
57
58 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
59 * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
60 * most significant bit (with weight 2**N) implicit.
61 *
62 * Special case for DIVISOR_LIMB == 100...000.
63 */
64 if (!(divisor_limb << 1))
65 divisor_limb_inverted = ~(mpi_limb_t)0;
66 else
67 udiv_qrnnd(divisor_limb_inverted, dummy,
68 -divisor_limb, 0, divisor_limb);
69
70 n1 = dividend_ptr[dividend_size - 1];
71 r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
72
73 /* Possible optimization:
74 * if (r == 0
75 * && divisor_limb > ((n1 << normalization_steps)
76 * | (dividend_ptr[dividend_size - 2] >> ...)))
77 * ...one division less...
78 */
79 for (i = dividend_size - 2; i >= 0; i--) {
80 n0 = dividend_ptr[i];
81 UDIV_QRNND_PREINV(dummy, r, r,
82 ((n1 << normalization_steps)
83 | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))),
84 divisor_limb, divisor_limb_inverted);
85 n1 = n0;
86 }
87 UDIV_QRNND_PREINV(dummy, r, r,
88 n1 << normalization_steps,
89 divisor_limb, divisor_limb_inverted);
90 return r >> normalization_steps;
91 } else {
92 mpi_limb_t divisor_limb_inverted;
93
94 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
95 * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
96 * most significant bit (with weight 2**N) implicit.
97 *
98 * Special case for DIVISOR_LIMB == 100...000.
99 */
100 if (!(divisor_limb << 1))
101 divisor_limb_inverted = ~(mpi_limb_t)0;
102 else
103 udiv_qrnnd(divisor_limb_inverted, dummy,
104 -divisor_limb, 0, divisor_limb);
105
106 i = dividend_size - 1;
107 r = dividend_ptr[i];
108
109 if (r >= divisor_limb)
110 r = 0;
111 else
112 i--;
113
114 for ( ; i >= 0; i--) {
115 n0 = dividend_ptr[i];
116 UDIV_QRNND_PREINV(dummy, r, r,
117 n0, divisor_limb, divisor_limb_inverted);
118 }
119 return r;
120 }
121 } else {
122 if (UDIV_NEEDS_NORMALIZATION) {
123 int normalization_steps;
124
125 normalization_steps = count_leading_zeros(divisor_limb);
126 if (normalization_steps) {
127 divisor_limb <<= normalization_steps;
128
129 n1 = dividend_ptr[dividend_size - 1];
130 r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
131
132 /* Possible optimization:
133 * if (r == 0
134 * && divisor_limb > ((n1 << normalization_steps)
135 * | (dividend_ptr[dividend_size - 2] >> ...)))
136 * ...one division less...
137 */
138 for (i = dividend_size - 2; i >= 0; i--) {
139 n0 = dividend_ptr[i];
140 udiv_qrnnd(dummy, r, r,
141 ((n1 << normalization_steps)
142 | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))),
143 divisor_limb);
144 n1 = n0;
145 }
146 udiv_qrnnd(dummy, r, r,
147 n1 << normalization_steps,
148 divisor_limb);
149 return r >> normalization_steps;
150 }
151 }
152 /* No normalization needed, either because udiv_qrnnd doesn't require
153 * it, or because DIVISOR_LIMB is already normalized.
154 */
155 i = dividend_size - 1;
156 r = dividend_ptr[i];
157
158 if (r >= divisor_limb)
159 r = 0;
160 else
161 i--;
162
163 for (; i >= 0; i--) {
164 n0 = dividend_ptr[i];
165 udiv_qrnnd(dummy, r, r, n0, divisor_limb);
166 }
167 return r;
168 }
169 }
170
171 /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
172 * the NSIZE-DSIZE least significant quotient limbs at QP
173 * and the DSIZE long remainder at NP. If QEXTRA_LIMBS is
174 * non-zero, generate that many fraction bits and append them after the
175 * other quotient limbs.
176 * Return the most significant limb of the quotient, this is always 0 or 1.
177 *
178 * Preconditions:
179 * 0. NSIZE >= DSIZE.
180 * 1. The most significant bit of the divisor must be set.
181 * 2. QP must either not overlap with the input operands at all, or
182 * QP + DSIZE >= NP must hold true. (This means that it's
183 * possible to put the quotient in the high part of NUM, right after the
184 * remainder in NUM.
185 * 3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero.
186 */
187
188 mpi_limb_t
mpihelp_divrem(mpi_ptr_t qp,mpi_size_t qextra_limbs,mpi_ptr_t np,mpi_size_t nsize,mpi_ptr_t dp,mpi_size_t dsize)189 mpihelp_divrem(mpi_ptr_t qp, mpi_size_t qextra_limbs,
190 mpi_ptr_t np, mpi_size_t nsize, mpi_ptr_t dp, mpi_size_t dsize)
191 {
192 mpi_limb_t most_significant_q_limb = 0;
193
194 switch (dsize) {
195 case 0:
196 /* We are asked to divide by zero, so go ahead and do it! (To make
197 the compiler not remove this statement, return the value.) */
198 /*
199 * existing clients of this function have been modified
200 * not to call it with dsize == 0, so this should not happen
201 */
202 return 1 / dsize;
203
204 case 1:
205 {
206 mpi_size_t i;
207 mpi_limb_t n1;
208 mpi_limb_t d;
209
210 d = dp[0];
211 n1 = np[nsize - 1];
212
213 if (n1 >= d) {
214 n1 -= d;
215 most_significant_q_limb = 1;
216 }
217
218 qp += qextra_limbs;
219 for (i = nsize - 2; i >= 0; i--)
220 udiv_qrnnd(qp[i], n1, n1, np[i], d);
221 qp -= qextra_limbs;
222
223 for (i = qextra_limbs - 1; i >= 0; i--)
224 udiv_qrnnd(qp[i], n1, n1, 0, d);
225
226 np[0] = n1;
227 }
228 break;
229
230 case 2:
231 {
232 mpi_size_t i;
233 mpi_limb_t n1, n0, n2;
234 mpi_limb_t d1, d0;
235
236 np += nsize - 2;
237 d1 = dp[1];
238 d0 = dp[0];
239 n1 = np[1];
240 n0 = np[0];
241
242 if (n1 >= d1 && (n1 > d1 || n0 >= d0)) {
243 sub_ddmmss(n1, n0, n1, n0, d1, d0);
244 most_significant_q_limb = 1;
245 }
246
247 for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--) {
248 mpi_limb_t q;
249 mpi_limb_t r;
250
251 if (i >= qextra_limbs)
252 np--;
253 else
254 np[0] = 0;
255
256 if (n1 == d1) {
257 /* Q should be either 111..111 or 111..110. Need special
258 * treatment of this rare case as normal division would
259 * give overflow. */
260 q = ~(mpi_limb_t) 0;
261
262 r = n0 + d1;
263 if (r < d1) { /* Carry in the addition? */
264 add_ssaaaa(n1, n0, r - d0,
265 np[0], 0, d0);
266 qp[i] = q;
267 continue;
268 }
269 n1 = d0 - (d0 != 0 ? 1 : 0);
270 n0 = -d0;
271 } else {
272 udiv_qrnnd(q, r, n1, n0, d1);
273 umul_ppmm(n1, n0, d0, q);
274 }
275
276 n2 = np[0];
277 q_test:
278 if (n1 > r || (n1 == r && n0 > n2)) {
279 /* The estimated Q was too large. */
280 q--;
281 sub_ddmmss(n1, n0, n1, n0, 0, d0);
282 r += d1;
283 if (r >= d1) /* If not carry, test Q again. */
284 goto q_test;
285 }
286
287 qp[i] = q;
288 sub_ddmmss(n1, n0, r, n2, n1, n0);
289 }
290 np[1] = n1;
291 np[0] = n0;
292 }
293 break;
294
295 default:
296 {
297 mpi_size_t i;
298 mpi_limb_t dX, d1, n0;
299
300 np += nsize - dsize;
301 dX = dp[dsize - 1];
302 d1 = dp[dsize - 2];
303 n0 = np[dsize - 1];
304
305 if (n0 >= dX) {
306 if (n0 > dX
307 || mpihelp_cmp(np, dp, dsize - 1) >= 0) {
308 mpihelp_sub_n(np, np, dp, dsize);
309 n0 = np[dsize - 1];
310 most_significant_q_limb = 1;
311 }
312 }
313
314 for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--) {
315 mpi_limb_t q;
316 mpi_limb_t n1, n2;
317 mpi_limb_t cy_limb;
318
319 if (i >= qextra_limbs) {
320 np--;
321 n2 = np[dsize];
322 } else {
323 n2 = np[dsize - 1];
324 MPN_COPY_DECR(np + 1, np, dsize - 1);
325 np[0] = 0;
326 }
327
328 if (n0 == dX) {
329 /* This might over-estimate q, but it's probably not worth
330 * the extra code here to find out. */
331 q = ~(mpi_limb_t) 0;
332 } else {
333 mpi_limb_t r;
334
335 udiv_qrnnd(q, r, n0, np[dsize - 1], dX);
336 umul_ppmm(n1, n0, d1, q);
337
338 while (n1 > r
339 || (n1 == r
340 && n0 > np[dsize - 2])) {
341 q--;
342 r += dX;
343 if (r < dX) /* I.e. "carry in previous addition?" */
344 break;
345 n1 -= n0 < d1;
346 n0 -= d1;
347 }
348 }
349
350 /* Possible optimization: We already have (q * n0) and (1 * n1)
351 * after the calculation of q. Taking advantage of that, we
352 * could make this loop make two iterations less. */
353 cy_limb = mpihelp_submul_1(np, dp, dsize, q);
354
355 if (n2 != cy_limb) {
356 mpihelp_add_n(np, np, dp, dsize);
357 q--;
358 }
359
360 qp[i] = q;
361 n0 = np[dsize - 1];
362 }
363 }
364 }
365
366 return most_significant_q_limb;
367 }
368
369 /****************
370 * Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
371 * Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
372 * Return the single-limb remainder.
373 * There are no constraints on the value of the divisor.
374 *
375 * QUOT_PTR and DIVIDEND_PTR might point to the same limb.
376 */
377
378 mpi_limb_t
mpihelp_divmod_1(mpi_ptr_t quot_ptr,mpi_ptr_t dividend_ptr,mpi_size_t dividend_size,mpi_limb_t divisor_limb)379 mpihelp_divmod_1(mpi_ptr_t quot_ptr,
380 mpi_ptr_t dividend_ptr, mpi_size_t dividend_size,
381 mpi_limb_t divisor_limb)
382 {
383 mpi_size_t i;
384 mpi_limb_t n1, n0, r;
385 mpi_limb_t dummy __maybe_unused;
386
387 if (!dividend_size)
388 return 0;
389
390 /* If multiplication is much faster than division, and the
391 * dividend is large, pre-invert the divisor, and use
392 * only multiplications in the inner loop.
393 *
394 * This test should be read:
395 * Does it ever help to use udiv_qrnnd_preinv?
396 * && Does what we save compensate for the inversion overhead?
397 */
398 if (UDIV_TIME > (2 * UMUL_TIME + 6)
399 && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME) {
400 int normalization_steps;
401
402 normalization_steps = count_leading_zeros(divisor_limb);
403 if (normalization_steps) {
404 mpi_limb_t divisor_limb_inverted;
405
406 divisor_limb <<= normalization_steps;
407
408 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
409 * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
410 * most significant bit (with weight 2**N) implicit.
411 */
412 /* Special case for DIVISOR_LIMB == 100...000. */
413 if (!(divisor_limb << 1))
414 divisor_limb_inverted = ~(mpi_limb_t)0;
415 else
416 udiv_qrnnd(divisor_limb_inverted, dummy,
417 -divisor_limb, 0, divisor_limb);
418
419 n1 = dividend_ptr[dividend_size - 1];
420 r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
421
422 /* Possible optimization:
423 * if (r == 0
424 * && divisor_limb > ((n1 << normalization_steps)
425 * | (dividend_ptr[dividend_size - 2] >> ...)))
426 * ...one division less...
427 */
428 for (i = dividend_size - 2; i >= 0; i--) {
429 n0 = dividend_ptr[i];
430 UDIV_QRNND_PREINV(quot_ptr[i + 1], r, r,
431 ((n1 << normalization_steps)
432 | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))),
433 divisor_limb, divisor_limb_inverted);
434 n1 = n0;
435 }
436 UDIV_QRNND_PREINV(quot_ptr[0], r, r,
437 n1 << normalization_steps,
438 divisor_limb, divisor_limb_inverted);
439 return r >> normalization_steps;
440 } else {
441 mpi_limb_t divisor_limb_inverted;
442
443 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
444 * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
445 * most significant bit (with weight 2**N) implicit.
446 */
447 /* Special case for DIVISOR_LIMB == 100...000. */
448 if (!(divisor_limb << 1))
449 divisor_limb_inverted = ~(mpi_limb_t) 0;
450 else
451 udiv_qrnnd(divisor_limb_inverted, dummy,
452 -divisor_limb, 0, divisor_limb);
453
454 i = dividend_size - 1;
455 r = dividend_ptr[i];
456
457 if (r >= divisor_limb)
458 r = 0;
459 else
460 quot_ptr[i--] = 0;
461
462 for ( ; i >= 0; i--) {
463 n0 = dividend_ptr[i];
464 UDIV_QRNND_PREINV(quot_ptr[i], r, r,
465 n0, divisor_limb, divisor_limb_inverted);
466 }
467 return r;
468 }
469 } else {
470 if (UDIV_NEEDS_NORMALIZATION) {
471 int normalization_steps;
472
473 normalization_steps = count_leading_zeros(divisor_limb);
474 if (normalization_steps) {
475 divisor_limb <<= normalization_steps;
476
477 n1 = dividend_ptr[dividend_size - 1];
478 r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
479
480 /* Possible optimization:
481 * if (r == 0
482 * && divisor_limb > ((n1 << normalization_steps)
483 * | (dividend_ptr[dividend_size - 2] >> ...)))
484 * ...one division less...
485 */
486 for (i = dividend_size - 2; i >= 0; i--) {
487 n0 = dividend_ptr[i];
488 udiv_qrnnd(quot_ptr[i + 1], r, r,
489 ((n1 << normalization_steps)
490 | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))),
491 divisor_limb);
492 n1 = n0;
493 }
494 udiv_qrnnd(quot_ptr[0], r, r,
495 n1 << normalization_steps,
496 divisor_limb);
497 return r >> normalization_steps;
498 }
499 }
500 /* No normalization needed, either because udiv_qrnnd doesn't require
501 * it, or because DIVISOR_LIMB is already normalized.
502 */
503 i = dividend_size - 1;
504 r = dividend_ptr[i];
505
506 if (r >= divisor_limb)
507 r = 0;
508 else
509 quot_ptr[i--] = 0;
510
511 for (; i >= 0; i--) {
512 n0 = dividend_ptr[i];
513 udiv_qrnnd(quot_ptr[i], r, r, n0, divisor_limb);
514 }
515 return r;
516 }
517 }
518