Mercurial > dropbear

comparison bn_mp_div.c @ 2:86e0b50a9b58 libtommath-orig ltm-0.30-orig

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ltm 0.30 orig import
author Matt Johnston <matt@ucc.asn.au>
date Mon, 31 May 2004 18:25:22 +0000
parents
children d29b64170cf0
comparison
equal deleted inserted replaced
-1:000000000000 2:86e0b50a9b58
1 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2 *
3 * LibTomMath is a library that provides multiple-precision
4 * integer arithmetic as well as number theoretic functionality.
5 *
6 * The library was designed directly after the MPI library by
7 * Michael Fromberger but has been written from scratch with
8 * additional optimizations in place.
9 *
10 * The library is free for all purposes without any express
11 * guarantee it works.
12 *
13 * Tom St Denis, [email protected], http://math.libtomcrypt.org
14 */
15 #include <tommath.h>
16
17 /* integer signed division.
18 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
19 * HAC pp.598 Algorithm 14.20
20 *
21 * Note that the description in HAC is horribly
22 * incomplete. For example, it doesn't consider
23 * the case where digits are removed from 'x' in
24 * the inner loop. It also doesn't consider the
25 * case that y has fewer than three digits, etc..
26 *
27 * The overall algorithm is as described as
28 * 14.20 from HAC but fixed to treat these cases.
29 */
30 int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
31 {
32 mp_int q, x, y, t1, t2;
33 int res, n, t, i, norm, neg;
34
35 /* is divisor zero ? */
36 if (mp_iszero (b) == 1) {
37 return MP_VAL;
38 }
39
40 /* if a < b then q=0, r = a */
41 if (mp_cmp_mag (a, b) == MP_LT) {
42 if (d != NULL) {
43 res = mp_copy (a, d);
44 } else {
45 res = MP_OKAY;
46 }
47 if (c != NULL) {
48 mp_zero (c);
49 }
50 return res;
51 }
52
53 if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) {
54 return res;
55 }
56 q.used = a->used + 2;
57
58 if ((res = mp_init (&t1)) != MP_OKAY) {
59 goto __Q;
60 }
61
62 if ((res = mp_init (&t2)) != MP_OKAY) {
63 goto __T1;
64 }
65
66 if ((res = mp_init_copy (&x, a)) != MP_OKAY) {
67 goto __T2;
68 }
69
70 if ((res = mp_init_copy (&y, b)) != MP_OKAY) {
71 goto __X;
72 }
73
74 /* fix the sign */
75 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
76 x.sign = y.sign = MP_ZPOS;
77
78 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
79 norm = mp_count_bits(&y) % DIGIT_BIT;
80 if (norm < (int)(DIGIT_BIT-1)) {
81 norm = (DIGIT_BIT-1) - norm;
82 if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
83 goto __Y;
84 }
85 if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) {
86 goto __Y;
87 }
88 } else {
89 norm = 0;
90 }
91
92 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
93 n = x.used - 1;
94 t = y.used - 1;
95
96 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
97 if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
98 goto __Y;
99 }
100
101 while (mp_cmp (&x, &y) != MP_LT) {
102 ++(q.dp[n - t]);
103 if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {
104 goto __Y;
105 }
106 }
107
108 /* reset y by shifting it back down */
109 mp_rshd (&y, n - t);
110
111 /* step 3. for i from n down to (t + 1) */
112 for (i = n; i >= (t + 1); i--) {
113 if (i > x.used) {
114 continue;
115 }
116
117 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
118 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
119 if (x.dp[i] == y.dp[t]) {
120 q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
121 } else {
122 mp_word tmp;
123 tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
124 tmp |= ((mp_word) x.dp[i - 1]);
125 tmp /= ((mp_word) y.dp[t]);
126 if (tmp > (mp_word) MP_MASK)
127 tmp = MP_MASK;
128 q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
129 }
130
131 /* while (q{i-t-1} * (yt * b + y{t-1})) >
132 xi * b**2 + xi-1 * b + xi-2
133
134 do q{i-t-1} -= 1;
135 */
136 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
137 do {
138 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
139
140 /* find left hand */
141 mp_zero (&t1);
142 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
143 t1.dp[1] = y.dp[t];
144 t1.used = 2;
145 if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
146 goto __Y;
147 }
148
149 /* find right hand */
150 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
151 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
152 t2.dp[2] = x.dp[i];
153 t2.used = 3;
154 } while (mp_cmp_mag(&t1, &t2) == MP_GT);
155
156 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
157 if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
158 goto __Y;
159 }
160
161 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
162 goto __Y;
163 }
164
165 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
166 goto __Y;
167 }
168
169 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
170 if (x.sign == MP_NEG) {
171 if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
172 goto __Y;
173 }
174 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
175 goto __Y;
176 }
177 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
178 goto __Y;
179 }
180
181 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
182 }
183 }
184
185 /* now q is the quotient and x is the remainder
186 * [which we have to normalize]
187 */
188
189 /* get sign before writing to c */
190 x.sign = a->sign;
191
192 if (c != NULL) {
193 mp_clamp (&q);
194 mp_exch (&q, c);
195 c->sign = neg;
196 }
197
198 if (d != NULL) {
199 mp_div_2d (&x, norm, &x, NULL);
200 mp_exch (&x, d);
201 }
202
203 res = MP_OKAY;
204
205 __Y:mp_clear (&y);
206 __X:mp_clear (&x);
207 __T2:mp_clear (&t2);
208 __T1:mp_clear (&t1);
209 __Q:mp_clear (&q);
210 return res;
211 }