|
2
|
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 } |
