1 /* $NetBSD: bn_mp_div.c,v 1.1.1.2 2014/04/24 12:45:31 pettai Exp $ */
5 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7 * LibTomMath is a library that provides multiple-precision
8 * integer arithmetic as well as number theoretic functionality.
10 * The library was designed directly after the MPI library by
11 * Michael Fromberger but has been written from scratch with
12 * additional optimizations in place.
14 * The library is free for all purposes without any express
17 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
20 #ifdef BN_MP_DIV_SMALL
22 /* slower bit-bang division... also smaller */
23 int mp_div(mp_int
* a
, mp_int
* b
, mp_int
* c
, mp_int
* d
)
28 /* is divisor zero ? */
29 if (mp_iszero (b
) == 1) {
33 /* if a < b then q=0, r = a */
34 if (mp_cmp_mag (a
, b
) == MP_LT
) {
47 if ((res
= mp_init_multi(&ta
, &tb
, &tq
, &q
, NULL
) != MP_OKAY
)) {
53 n
= mp_count_bits(a
) - mp_count_bits(b
);
54 if (((res
= mp_abs(a
, &ta
)) != MP_OKAY
) ||
55 ((res
= mp_abs(b
, &tb
)) != MP_OKAY
) ||
56 ((res
= mp_mul_2d(&tb
, n
, &tb
)) != MP_OKAY
) ||
57 ((res
= mp_mul_2d(&tq
, n
, &tq
)) != MP_OKAY
)) {
62 if (mp_cmp(&tb
, &ta
) != MP_GT
) {
63 if (((res
= mp_sub(&ta
, &tb
, &ta
)) != MP_OKAY
) ||
64 ((res
= mp_add(&q
, &tq
, &q
)) != MP_OKAY
)) {
68 if (((res
= mp_div_2d(&tb
, 1, &tb
, NULL
)) != MP_OKAY
) ||
69 ((res
= mp_div_2d(&tq
, 1, &tq
, NULL
)) != MP_OKAY
)) {
74 /* now q == quotient and ta == remainder */
76 n2
= (a
->sign
== b
->sign
? MP_ZPOS
: MP_NEG
);
79 c
->sign
= (mp_iszero(c
) == MP_YES
) ? MP_ZPOS
: n2
;
83 d
->sign
= (mp_iszero(d
) == MP_YES
) ? MP_ZPOS
: n
;
86 mp_clear_multi(&ta
, &tb
, &tq
, &q
, NULL
);
92 /* integer signed division.
93 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
94 * HAC pp.598 Algorithm 14.20
96 * Note that the description in HAC is horribly
97 * incomplete. For example, it doesn't consider
98 * the case where digits are removed from 'x' in
99 * the inner loop. It also doesn't consider the
100 * case that y has fewer than three digits, etc..
102 * The overall algorithm is as described as
103 * 14.20 from HAC but fixed to treat these cases.
105 int mp_div (mp_int
* a
, mp_int
* b
, mp_int
* c
, mp_int
* d
)
107 mp_int q
, x
, y
, t1
, t2
;
108 int res
, n
, t
, i
, norm
, neg
;
110 /* is divisor zero ? */
111 if (mp_iszero (b
) == 1) {
115 /* if a < b then q=0, r = a */
116 if (mp_cmp_mag (a
, b
) == MP_LT
) {
118 res
= mp_copy (a
, d
);
128 if ((res
= mp_init_size (&q
, a
->used
+ 2)) != MP_OKAY
) {
131 q
.used
= a
->used
+ 2;
133 if ((res
= mp_init (&t1
)) != MP_OKAY
) {
137 if ((res
= mp_init (&t2
)) != MP_OKAY
) {
141 if ((res
= mp_init_copy (&x
, a
)) != MP_OKAY
) {
145 if ((res
= mp_init_copy (&y
, b
)) != MP_OKAY
) {
150 neg
= (a
->sign
== b
->sign
) ? MP_ZPOS
: MP_NEG
;
151 x
.sign
= y
.sign
= MP_ZPOS
;
153 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
154 norm
= mp_count_bits(&y
) % DIGIT_BIT
;
155 if (norm
< (int)(DIGIT_BIT
-1)) {
156 norm
= (DIGIT_BIT
-1) - norm
;
157 if ((res
= mp_mul_2d (&x
, norm
, &x
)) != MP_OKAY
) {
160 if ((res
= mp_mul_2d (&y
, norm
, &y
)) != MP_OKAY
) {
167 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
171 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
172 if ((res
= mp_lshd (&y
, n
- t
)) != MP_OKAY
) { /* y = y*b**{n-t} */
176 while (mp_cmp (&x
, &y
) != MP_LT
) {
178 if ((res
= mp_sub (&x
, &y
, &x
)) != MP_OKAY
) {
183 /* reset y by shifting it back down */
186 /* step 3. for i from n down to (t + 1) */
187 for (i
= n
; i
>= (t
+ 1); i
--) {
192 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
193 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
194 if (x
.dp
[i
] == y
.dp
[t
]) {
195 q
.dp
[i
- t
- 1] = ((((mp_digit
)1) << DIGIT_BIT
) - 1);
198 tmp
= ((mp_word
) x
.dp
[i
]) << ((mp_word
) DIGIT_BIT
);
199 tmp
|= ((mp_word
) x
.dp
[i
- 1]);
200 tmp
/= ((mp_word
) y
.dp
[t
]);
201 if (tmp
> (mp_word
) MP_MASK
)
203 q
.dp
[i
- t
- 1] = (mp_digit
) (tmp
& (mp_word
) (MP_MASK
));
206 /* while (q{i-t-1} * (yt * b + y{t-1})) >
207 xi * b**2 + xi-1 * b + xi-2
211 q
.dp
[i
- t
- 1] = (q
.dp
[i
- t
- 1] + 1) & MP_MASK
;
213 q
.dp
[i
- t
- 1] = (q
.dp
[i
- t
- 1] - 1) & MP_MASK
;
217 t1
.dp
[0] = (t
- 1 < 0) ? 0 : y
.dp
[t
- 1];
220 if ((res
= mp_mul_d (&t1
, q
.dp
[i
- t
- 1], &t1
)) != MP_OKAY
) {
224 /* find right hand */
225 t2
.dp
[0] = (i
- 2 < 0) ? 0 : x
.dp
[i
- 2];
226 t2
.dp
[1] = (i
- 1 < 0) ? 0 : x
.dp
[i
- 1];
229 } while (mp_cmp_mag(&t1
, &t2
) == MP_GT
);
231 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
232 if ((res
= mp_mul_d (&y
, q
.dp
[i
- t
- 1], &t1
)) != MP_OKAY
) {
236 if ((res
= mp_lshd (&t1
, i
- t
- 1)) != MP_OKAY
) {
240 if ((res
= mp_sub (&x
, &t1
, &x
)) != MP_OKAY
) {
244 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
245 if (x
.sign
== MP_NEG
) {
246 if ((res
= mp_copy (&y
, &t1
)) != MP_OKAY
) {
249 if ((res
= mp_lshd (&t1
, i
- t
- 1)) != MP_OKAY
) {
252 if ((res
= mp_add (&x
, &t1
, &x
)) != MP_OKAY
) {
256 q
.dp
[i
- t
- 1] = (q
.dp
[i
- t
- 1] - 1UL) & MP_MASK
;
260 /* now q is the quotient and x is the remainder
261 * [which we have to normalize]
264 /* get sign before writing to c */
265 x
.sign
= x
.used
== 0 ? MP_ZPOS
: a
->sign
;
274 mp_div_2d (&x
, norm
, &x
, NULL
);
282 LBL_T2
:mp_clear (&t2
);
283 LBL_T1
:mp_clear (&t1
);
292 /* Source: /cvs/libtom/libtommath/bn_mp_div.c,v */
294 /* Date: 2006/12/28 01:25:13 */