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 import std.algorithm : max;
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(const(BigInt)* x, const(BigInt)* y_arg, BigInt* q, 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             if (y_t == 0)
85                 throw new InvalidArgument("Division by zero due to too large bn");
86             
87             if (x_j0 == y_t)
88                 q_words[j-t-1] = MP_WORD_MAX;
89             else
90                 q_words[j-t-1] = bigint_divop(x_j0, x_j1, y_t);
91             
92             while (divisionCheck(q_words[j-t-1], y_t, y.wordAt(t-1), x_j0, x_j1, r.wordAt(j-2)))
93             {
94                 q_words[j-t-1] -= 1;
95             }
96 			auto y_1 = (y * q_words[j-t-1]);
97 			auto j_t_1 = (j-t-1);
98             y_1 <<= (MP_WORD_BITS * j_t_1);
99             *r -= y_1;
100             
101             if (r.isNegative())
102             {
103                 y <<= (MP_WORD_BITS * (j-t-1));
104                 *r += y;
105                 q_words[j-t-1] -= 1;
106             }
107         }
108         *r >>= shifts;
109     }
110     
111     signFixup(x, y_arg, q, r);
112 }
113 
114 private:
115 /*
116 * Handle signed operands, if necessary
117 */
118 void signFixup(const(BigInt)* x, const(BigInt)* y, BigInt* q, BigInt* r)
119 {
120     if (x.sign() == BigInt.Negative)
121     {
122         q.flipSign();
123         if (r.isNonzero()) { --*q; *r = y.abs() - *r; }
124     }
125     if (y.sign() == BigInt.Negative)
126         q.flipSign();
127 }
128 
129 bool divisionCheck(word q, word y2, word y1, word x3, word x2, word x1)
130 {
131     // Compute (y3,y2,y1) = (y2,y1) * q
132     
133     word y3 = 0;
134     y1 = word_madd2(q, y1, &y3);
135     y2 = word_madd2(q, y2, &y3);
136 
137     // Return (y3,y2,y1) >? (x3,x2,x1)
138     
139     if (y3 > x3) return true;
140     if (y3 < x3) return false;
141     
142     if (y2 > x2) return true;
143     if (y2 < x2) return false;
144     
145     if (y1 > x1) return true;
146     if (y1 < x1) return false;
147     
148     return false;
149 }