botan.math.bigint.divide source code

1 /**
2 * Division
3 * 
4 * Copyright:
5 * (C) 1999-2007 Jack Lloyd
6 * (C) 2014-2015 Etienne Cimon
7 *
8 * License:
9 * Botan is released under the Simplified BSD License (see LICENSE.md)
10 */
11 module botan.math.bigint.divide;
12 
13 import botan.constants;
14 import botan.math.bigint.bigint;
15 import botan.math.mp.mp_core;
16 import botan.constants;
17 
18 /**
19 * BigInt Division
20 * Params:
21 *  x = an integer
22 *  y_arg = a non-zero integer
23 *  q = will be set to x / y
24 *  r = will be set to x % y
25 */
26 void divide()(auto const ref BigInt x, auto const ref BigInt y_arg, ref BigInt q, ref BigInt r)
27 {
28     /*
29     * Solve x = q * y + r
30     */
31     if (y_arg.isZero())
32         throw new BigInt.DivideByZero();
33     
34     BigInt y = y_arg.dup;
35     const size_t y_words = y.sigWords();
36     
37     r = x.dup;
38     q = 0;
39     
40     r.setSign(BigInt.Positive);
41     y.setSign(BigInt.Positive);
42     
43     int compare = r.cmp(y);
44     
45     if (compare == 0)
46     {
47         q = 1;
48         r = 0;
49     }
50     else if (compare > 0)
51     {
52         size_t shifts = 0;
53         word y_top = y.wordAt(y.sigWords()-1);
54         while (y_top < MP_WORD_TOP_BIT) { y_top <<= 1; ++shifts; }
55         y <<= shifts;
56         r <<= shifts;
57         
58         const size_t n = r.sigWords() - 1, t = y_words - 1;
59         
60         if (n < t)
61             throw new InternalError("BigInt division word sizes");
62         
63         q.growTo(n - t + 1);
64         
65         word* q_words = q.mutablePtr();
66         
67         if (n <= t)
68         {
69             while (r > y) { r -= y; ++q; }
70             r >>= shifts;
71             signFixup(x, y_arg, q, r);
72             return;
73         }
74         
75         BigInt temp = y << (MP_WORD_BITS * (n-t));
76         
77         while (r >= temp) { r -= temp; q_words[n-t] += 1; }
78         
79         for (size_t j = n; j != t; --j)
80         {
81             const word x_j0  = r.wordAt(j);
82             const word x_j1 = r.wordAt(j-1);
83             const word y_t  = y.wordAt(t);
84             
85             if (x_j0 == y_t)
86                 q_words[j-t-1] = MP_WORD_MAX;
87             else
88                 q_words[j-t-1] = bigint_divop(x_j0, x_j1, y_t);
89             
90             while (divisionCheck(q_words[j-t-1], y_t, y.wordAt(t-1), x_j0, x_j1, r.wordAt(j-2)))
91             {
92                 q_words[j-t-1] -= 1;
93             }
94 			auto y_1 = (y * q_words[j-t-1]);
95 			auto j_t_1 = (j-t-1);
96             r -= y_1 << (MP_WORD_BITS * j_t_1);
97             
98             if (r.isNegative())
99             {
100                 r += y << (MP_WORD_BITS * (j-t-1));
101                 q_words[j-t-1] -= 1;
102             }
103         }
104         r >>= shifts;
105     }
106     
107     signFixup(x, y_arg, q, r);
108 }
109 
110 private:
111 /*
112 * Handle signed operands, if necessary
113 */
114 void signFixup(const ref BigInt x, const ref BigInt y, ref BigInt q, ref BigInt r)
115 {
116     if (x.sign() == BigInt.Negative)
117     {
118         q.flipSign();
119         if (r.isNonzero()) { --q; r = y.abs() - r; }
120     }
121     if (y.sign() == BigInt.Negative)
122         q.flipSign();
123 }
124 
125 bool divisionCheck(word q, word y2, word y1, word x3, word x2, word x1)
126 {
127     // Compute (y3,y2,y1) = (y2,y1) * q
128     
129     word y3 = 0;
130     y1 = word_madd2(q, y1, &y3);
131     y2 = word_madd2(q, y2, &y3);
132 
133     // Return (y3,y2,y1) >? (x3,x2,x1)
134     
135     if (y3 > x3) return true;
136     if (y3 < x3) return false;
137     
138     if (y2 > x2) return true;
139     if (y2 < x2) return false;
140     
141     if (y1 > x1) return true;
142     if (y1 < x1) return false;
143     
144     return false;
145 }