/src/dropbear/libtommath/bn_mp_div.c

Source (jump to first uncovered line)
#include "tommath_private.h"
#ifdef BN_MP_DIV_C
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
/* SPDX-License-Identifier: Unlicense */

#ifdef BN_MP_DIV_SMALL

/* slower bit-bang division... also smaller */
mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
{
   mp_int ta, tb, tq, q;
   int     n, n2;
   mp_err err;

   /* is divisor zero ? */
   if (MP_IS_ZERO(b)) {
      return MP_VAL;
   }

   /* if a < b then q=0, r = a */
   if (mp_cmp_mag(a, b) == MP_LT) {
      if (d != NULL) {
         err = mp_copy(a, d);
      } else {
         err = MP_OKAY;
      }
      if (c != NULL) {
         mp_zero(c);
      }
      return err;
   }

   /* init our temps */
   if ((err = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) {
      return err;
   }


   mp_set(&tq, 1uL);
   n = mp_count_bits(a) - mp_count_bits(b);
   if ((err = mp_abs(a, &ta)) != MP_OKAY)                         goto LBL_ERR;
   if ((err = mp_abs(b, &tb)) != MP_OKAY)                         goto LBL_ERR;
   if ((err = mp_mul_2d(&tb, n, &tb)) != MP_OKAY)                 goto LBL_ERR;
   if ((err = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)                 goto LBL_ERR;

   while (n-- >= 0) {
      if (mp_cmp(&tb, &ta) != MP_GT) {
         if ((err = mp_sub(&ta, &tb, &ta)) != MP_OKAY)            goto LBL_ERR;
         if ((err = mp_add(&q, &tq, &q)) != MP_OKAY)              goto LBL_ERR;
      }
      if ((err = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY)        goto LBL_ERR;
      if ((err = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)        goto LBL_ERR;
   }

   /* now q == quotient and ta == remainder */
   n  = a->sign;
   n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
   if (c != NULL) {
      mp_exch(c, &q);
      c->sign  = MP_IS_ZERO(c) ? MP_ZPOS : n2;
   }
   if (d != NULL) {
      mp_exch(d, &ta);
      d->sign = MP_IS_ZERO(d) ? MP_ZPOS : n;
   }
LBL_ERR:
   mp_clear_multi(&ta, &tb, &tq, &q, NULL);
   return err;
}

#else

/* integer signed division.
 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
 * HAC pp.598 Algorithm 14.20
 *
 * Note that the description in HAC is horribly
 * incomplete.  For example, it doesn't consider
 * the case where digits are removed from 'x' in
 * the inner loop.  It also doesn't consider the
 * case that y has fewer than three digits, etc..
 *
 * The overall algorithm is as described as
 * 14.20 from HAC but fixed to treat these cases.
*/
mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
{
   mp_int  q, x, y, t1, t2;
   int     n, t, i, norm;
   mp_sign neg;
   mp_err  err;

   /* is divisor zero ? */
   if (MP_IS_ZERO(b)) {
      return MP_VAL;
   }

   /* if a < b then q=0, r = a */
   if (mp_cmp_mag(a, b) == MP_LT) {
      if (d != NULL) {
         err = mp_copy(a, d);
      } else {
         err = MP_OKAY;
      }
      if (c != NULL) {
         mp_zero(c);
      }
      return err;
   }

   if ((err = mp_init_size(&q, a->used + 2)) != MP_OKAY) {
      return err;
   }
   q.used = a->used + 2;

   if ((err = mp_init(&t1)) != MP_OKAY)                           goto LBL_Q;

   if ((err = mp_init(&t2)) != MP_OKAY)                           goto LBL_T1;

   if ((err = mp_init_copy(&x, a)) != MP_OKAY)                    goto LBL_T2;

   if ((err = mp_init_copy(&y, b)) != MP_OKAY)                    goto LBL_X;

   /* fix the sign */
   neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
   x.sign = y.sign = MP_ZPOS;

   /* normalize both x and y, ensure that y >= b/2, [b == 2**MP_DIGIT_BIT] */
   norm = mp_count_bits(&y) % MP_DIGIT_BIT;
   if (norm < (MP_DIGIT_BIT - 1)) {
      norm = (MP_DIGIT_BIT - 1) - norm;
      if ((err = mp_mul_2d(&x, norm, &x)) != MP_OKAY)             goto LBL_Y;
      if ((err = mp_mul_2d(&y, norm, &y)) != MP_OKAY)             goto LBL_Y;
   } else {
      norm = 0;
   }

   /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
   n = x.used - 1;
   t = y.used - 1;

   /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
   /* y = y*b**{n-t} */
   if ((err = mp_lshd(&y, n - t)) != MP_OKAY)                     goto LBL_Y;

   while (mp_cmp(&x, &y) != MP_LT) {
      ++(q.dp[n - t]);
      if ((err = mp_sub(&x, &y, &x)) != MP_OKAY)                  goto LBL_Y;
   }

   /* reset y by shifting it back down */
   mp_rshd(&y, n - t);

   /* step 3. for i from n down to (t + 1) */
   for (i = n; i >= (t + 1); i--) {
      if (i > x.used) {
         continue;
      }

      /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
       * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
      if (x.dp[i] == y.dp[t]) {
         q.dp[(i - t) - 1] = ((mp_digit)1 << (mp_digit)MP_DIGIT_BIT) - (mp_digit)1;
      } else {
         mp_word tmp;
         tmp = (mp_word)x.dp[i] << (mp_word)MP_DIGIT_BIT;
         tmp |= (mp_word)x.dp[i - 1];
         tmp /= (mp_word)y.dp[t];
         if (tmp > (mp_word)MP_MASK) {
            tmp = MP_MASK;
         }
         q.dp[(i - t) - 1] = (mp_digit)(tmp & (mp_word)MP_MASK);
      }

      /* while (q{i-t-1} * (yt * b + y{t-1})) >
               xi * b**2 + xi-1 * b + xi-2

         do q{i-t-1} -= 1;
      */
      q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1uL) & (mp_digit)MP_MASK;
      do {
         q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & (mp_digit)MP_MASK;

         /* find left hand */
         mp_zero(&t1);
         t1.dp[0] = ((t - 1) < 0) ? 0u : y.dp[t - 1];
         t1.dp[1] = y.dp[t];
         t1.used = 2;
         if ((err = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) goto LBL_Y;

         /* find right hand */
         t2.dp[0] = ((i - 2) < 0) ? 0u : x.dp[i - 2];
         t2.dp[1] = x.dp[i - 1]; /* i >= 1 always holds */
         t2.dp[2] = x.dp[i];
         t2.used = 3;
      } while (mp_cmp_mag(&t1, &t2) == MP_GT);

      /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
      if ((err = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) goto LBL_Y;

      if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY)           goto LBL_Y;

      if ((err = mp_sub(&x, &t1, &x)) != MP_OKAY)                 goto LBL_Y;

      /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
      if (x.sign == MP_NEG) {
         if ((err = mp_copy(&y, &t1)) != MP_OKAY)                 goto LBL_Y;
         if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY)        goto LBL_Y;
         if ((err = mp_add(&x, &t1, &x)) != MP_OKAY)              goto LBL_Y;

         q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & MP_MASK;
      }
   }

   /* now q is the quotient and x is the remainder
    * [which we have to normalize]
    */

   /* get sign before writing to c */
   x.sign = (x.used == 0) ? MP_ZPOS : a->sign;

   if (c != NULL) {
      mp_clamp(&q);
      mp_exch(&q, c);
      c->sign = neg;
   }

   if (d != NULL) {
      if ((err = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY)       goto LBL_Y;
      mp_exch(&x, d);
   }

   err = MP_OKAY;

LBL_Y:
   mp_clear(&y);
LBL_X:
   mp_clear(&x);
LBL_T2:
   mp_clear(&t2);
LBL_T1:
   mp_clear(&t1);
LBL_Q:
   mp_clear(&q);
   return err;
}

#endif

#endif

Coverage Report

Created: 2023-03-20 06:28

Line	Count	Source (jump to first uncovered line)
1		#include "tommath_private.h"
2		#ifdef BN_MP_DIV_C
3		/* LibTomMath, multiple-precision integer library -- Tom St Denis */
4		/* SPDX-License-Identifier: Unlicense */
5
6		#ifdef BN_MP_DIV_SMALL
7
8		/* slower bit-bang division... also smaller */
9		mp_err mp_div(const mp_int a, const mp_int b, mp_int c, mp_int d)
10		{
11		mp_int ta, tb, tq, q;
12		int n, n2;
13		mp_err err;
14
15		/* is divisor zero ? */
16		if (MP_IS_ZERO(b)) {
17		return MP_VAL;
18		}
19
20		/* if a < b then q=0, r = a */
21		if (mp_cmp_mag(a, b) == MP_LT) {
22		if (d != NULL) {
23		err = mp_copy(a, d);
24		} else {
25		err = MP_OKAY;
26		}
27		if (c != NULL) {
28		mp_zero(c);
29		}
30		return err;
31		}
32
33		/* init our temps */
34		if ((err = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) {
35		return err;
36		}
37
38
39		mp_set(&tq, 1uL);
40		n = mp_count_bits(a) - mp_count_bits(b);
41		if ((err = mp_abs(a, &ta)) != MP_OKAY) goto LBL_ERR;
42		if ((err = mp_abs(b, &tb)) != MP_OKAY) goto LBL_ERR;
43		if ((err = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) goto LBL_ERR;
44		if ((err = mp_mul_2d(&tq, n, &tq)) != MP_OKAY) goto LBL_ERR;
45
46		while (n-- >= 0) {
47		if (mp_cmp(&tb, &ta) != MP_GT) {
48		if ((err = mp_sub(&ta, &tb, &ta)) != MP_OKAY) goto LBL_ERR;
49		if ((err = mp_add(&q, &tq, &q)) != MP_OKAY) goto LBL_ERR;
50		}
51		if ((err = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) goto LBL_ERR;
52		if ((err = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY) goto LBL_ERR;
53		}
54
55		/* now q == quotient and ta == remainder */
56		n = a->sign;
57		n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
58		if (c != NULL) {
59		mp_exch(c, &q);
60		c->sign = MP_IS_ZERO(c) ? MP_ZPOS : n2;
61		}
62		if (d != NULL) {
63		mp_exch(d, &ta);
64		d->sign = MP_IS_ZERO(d) ? MP_ZPOS : n;
65		}
66		LBL_ERR:
67		mp_clear_multi(&ta, &tb, &tq, &q, NULL);
68		return err;
69		}
70
71		#else
72
73		/* integer signed division.
74		* c*b + d == a [e.g. a/b, c=quotient, d=remainder]
75		* HAC pp.598 Algorithm 14.20
76		*
77		* Note that the description in HAC is horribly
78		* incomplete. For example, it doesn't consider
79		* the case where digits are removed from 'x' in
80		* the inner loop. It also doesn't consider the
81		* case that y has fewer than three digits, etc..
82		*
83		* The overall algorithm is as described as
84		* 14.20 from HAC but fixed to treat these cases.
85		*/
86		mp_err mp_div(const mp_int a, const mp_int b, mp_int c, mp_int d)
87	4	{
88	4	mp_int q, x, y, t1, t2;
89	4	int n, t, i, norm;
90	4	mp_sign neg;
91	4	mp_err err;
92
93		/* is divisor zero ? */
94	4	if (MP_IS_ZERO(b)) {
95	0	return MP_VAL;
96	0	}
97
98		/* if a < b then q=0, r = a */
99	4	if (mp_cmp_mag(a, b) == MP_LT) {
100	0	if (d != NULL) {
101	0	err = mp_copy(a, d);
102	0	} else {
103	0	err = MP_OKAY;
104	0	}
105	0	if (c != NULL) {
106	0	mp_zero(c);
107	0	}
108	0	return err;
109	0	}
110
111	4	if ((err = mp_init_size(&q, a->used + 2)) != MP_OKAY) {
112	0	return err;
113	0	}
114	4	q.used = a->used + 2;
115
116	4	if ((err = mp_init(&t1)) != MP_OKAY) goto LBL_Q;
117
118	4	if ((err = mp_init(&t2)) != MP_OKAY) goto LBL_T1;
119
120	4	if ((err = mp_init_copy(&x, a)) != MP_OKAY) goto LBL_T2;
121
122	4	if ((err = mp_init_copy(&y, b)) != MP_OKAY) goto LBL_X;
123
124		/* fix the sign */
125	4	neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
126	4	x.sign = y.sign = MP_ZPOS;
127
128		/* normalize both x and y, ensure that y >= b/2, [b == 2*MP_DIGIT_BIT] /
129	4	norm = mp_count_bits(&y) % MP_DIGIT_BIT;
130	4	if (norm < (MP_DIGIT_BIT - 1)) {
131	4	norm = (MP_DIGIT_BIT - 1) - norm;
132	4	if ((err = mp_mul_2d(&x, norm, &x)) != MP_OKAY) goto LBL_Y;
133	4	if ((err = mp_mul_2d(&y, norm, &y)) != MP_OKAY) goto LBL_Y;
134	4	} else {
135	0	norm = 0;
136	0	}
137
138		/* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
139	4	n = x.used - 1;
140	4	t = y.used - 1;
141
142		/* while (x >= ybn-t) do { q[n-t] += 1; x -= yb*{n-t} } /
143		/* y = yb{n-t} /
144	4	if ((err = mp_lshd(&y, n - t)) != MP_OKAY) goto LBL_Y;
145
146	4	while (mp_cmp(&x, &y) != MP_LT) {
147	0	++(q.dp[n - t]);
148	0	if ((err = mp_sub(&x, &y, &x)) != MP_OKAY) goto LBL_Y;
149	0	}
150
151		/* reset y by shifting it back down */
152	4	mp_rshd(&y, n - t);
153
154		/* step 3. for i from n down to (t + 1) */
155	24	for (i = n; i >= (t + 1); i--) {
156	20	if (i > x.used) {
157	0	continue;
158	0	}
159
160		/* step 3.1 if xi == yt then set q{i-t-1} to b-1,
161		* otherwise set q{i-t-1} to (xib + x{i-1})/yt /
162	20	if (x.dp[i] == y.dp[t]) {
163	0	q.dp[(i - t) - 1] = ((mp_digit)1 << (mp_digit)MP_DIGIT_BIT) - (mp_digit)1;
164	20	} else {
165	20	mp_word tmp;
166	20	tmp = (mp_word)x.dp[i] << (mp_word)MP_DIGIT_BIT;
167	20	tmp \|= (mp_word)x.dp[i - 1];
168	20	tmp /= (mp_word)y.dp[t];
169	20	if (tmp > (mp_word)MP_MASK) {
170	0	tmp = MP_MASK;
171	0	}
172	20	q.dp[(i - t) - 1] = (mp_digit)(tmp & (mp_word)MP_MASK);
173	20	}
174
175		/* while (q{i-t-1} * (yt * b + y{t-1})) >
176		xi * b*2 + xi-1 b + xi-2
177
178		do q{i-t-1} -= 1;
179		*/
180	20	q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1uL) & (mp_digit)MP_MASK;
181	20	do {
182	20	q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & (mp_digit)MP_MASK;
183
184		/* find left hand */
185	20	mp_zero(&t1);
186	20	t1.dp[0] = ((t - 1) < 0) ? 0u : y.dp[t - 1];
187	20	t1.dp[1] = y.dp[t];
188	20	t1.used = 2;
189	20	if ((err = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) goto LBL_Y;
190
191		/* find right hand */
192	20	t2.dp[0] = ((i - 2) < 0) ? 0u : x.dp[i - 2];
193	20	t2.dp[1] = x.dp[i - 1]; /* i >= 1 always holds */
194	20	t2.dp[2] = x.dp[i];
195	20	t2.used = 3;
196	20	} while (mp_cmp_mag(&t1, &t2) == MP_GT);
197
198		/* step 3.3 x = x - q{i-t-1} * y * b*{i-t-1} /
199	20	if ((err = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) goto LBL_Y;
200
201	20	if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) goto LBL_Y;
202
203	20	if ((err = mp_sub(&x, &t1, &x)) != MP_OKAY) goto LBL_Y;
204
205		/* if x < 0 then { x = x + yb{i-t-1}; q{i-t-1} -= 1; } /
206	20	if (x.sign == MP_NEG) {
207	0	if ((err = mp_copy(&y, &t1)) != MP_OKAY) goto LBL_Y;
208	0	if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) goto LBL_Y;
209	0	if ((err = mp_add(&x, &t1, &x)) != MP_OKAY) goto LBL_Y;
210
211	0	q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & MP_MASK;
212	0	}
213	20	}
214
215		/* now q is the quotient and x is the remainder
216		* [which we have to normalize]
217		*/
218
219		/* get sign before writing to c */
220	4	x.sign = (x.used == 0) ? MP_ZPOS : a->sign;
221
222	4	if (c != NULL) {
223	0	mp_clamp(&q);
224	0	mp_exch(&q, c);
225	0	c->sign = neg;
226	0	}
227
228	4	if (d != NULL) {
229	4	if ((err = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY) goto LBL_Y;
230	4	mp_exch(&x, d);
231	4	}
232
233	4	err = MP_OKAY;
234
235	4	LBL_Y:
236	4	mp_clear(&y);
237	4	LBL_X:
238	4	mp_clear(&x);
239	4	LBL_T2:
240	4	mp_clear(&t2);
241	4	LBL_T1:
242	4	mp_clear(&t1);
243	4	LBL_Q:
244	4	mp_clear(&q);
245	4	return err;
246	4	}
247
248		#endif
249
250		#endif