@@ -35,161 +35,15 @@ namespace boost {
35
35
*/
36
36
template <typename Backend>
37
37
class montgomery_params : virtual public base_params<Backend> {
38
- protected:
39
- template <typename Number>
40
- inline void initialize_montgomery_params (const Number& p) {
41
- this ->initialize_base_params (p);
42
- find_const_variables (p);
43
- }
44
-
45
- inline void initialize_montgomery_params (const montgomery_params<Backend>& p) {
46
- this ->initialize_base_params (p);
47
- find_const_variables (p);
48
- }
49
-
50
- /*
51
- * Compute -input^-1 mod 2^limb_bits. Throws an exception if input
52
- * is even. If input is odd, then input and 2^n are relatively prime
53
- * and an inverse exists.
54
- */
55
- limb_type monty_inverse (limb_type a) {
56
- assert (a % 2 == 1 );
57
-
58
- limb_type b = 1 ;
59
- limb_type r = 0 ;
60
-
61
- for (size_t i = 0 ; i != sizeof (limb_type) * CHAR_BIT; ++i) {
62
- const limb_type bi = b % 2 ;
63
- r >>= 1 ;
64
- r += bi << (sizeof (limb_type) * CHAR_BIT - 1 );
65
-
66
- b -= a * bi;
67
- b >>= 1 ;
68
- }
69
-
70
- // Now invert in addition space
71
- r = (~static_cast <limb_type>(0 ) - r) + 1 ;
72
-
73
- return r;
74
- }
75
-
76
- template <typename T>
77
- void find_const_variables (const T& pp) {
78
- Backend p = pp;
79
- if (p <= 0 || !(p % 2 )) {
80
- return ;
81
- }
82
-
83
- m_p_words = this ->m_mod .size ();
84
-
85
- m_p_dash = monty_inverse (this ->m_mod .limbs ()[0 ]);
86
-
87
- Backend r;
88
-
89
- boost::multiprecision::default_ops::eval_bit_set (r, m_p_words * sizeof (limb_type) * CHAR_BIT);
90
-
91
- m_r2 = r * r;
92
- barrett_params<Backend> barrettParams (this ->m_mod );
93
- barrettParams.barrett_reduce (m_r2);
94
- }
95
-
96
38
public:
97
39
montgomery_params () : base_params<Backend>() {
98
40
}
99
41
100
42
template <typename Number>
101
43
explicit montgomery_params (const Number& p) : base_params<Backend>(p) {
102
- initialize_montgomery_params (p);
103
- }
104
-
105
- inline const Backend& r2 () const {
106
- return m_r2;
107
- }
108
-
109
- inline limb_type p_dash () const {
110
- return m_p_dash;
111
- }
112
-
113
- inline size_t p_words () const {
114
- return m_p_words;
115
- }
116
-
117
- template <class V >
118
- montgomery_params& operator =(const V& v) {
119
- initialize_montgomery_params (v);
120
- return *this ;
121
- }
122
-
123
- inline void montgomery_reduce (Backend& result) const {
124
- using boost::multiprecision::default_ops::eval_multiply_add;
125
- using boost::multiprecision::default_ops::eval_right_shift;
126
- using boost::multiprecision::default_ops::eval_add;
127
-
128
- typedef cpp_int_modular_backend<sizeof (limb_type) * CHAR_BIT * 3 >
129
- cpp_three_int_backend;
130
-
131
- const size_t p_size = m_p_words;
132
- const limb_type p_dash = m_p_dash;
133
- const size_t z_size = 2 * (p_words () + 1 );
134
-
135
- boost::container::vector<limb_type> z (
136
- result.size (), 0 ); // container::vector<limb_type, alloc> z(result.size(), 0);
137
- for (size_t i = 0 ; i < result.size (); ++i) {
138
- z[i] = result.limbs ()[i];
139
- }
140
-
141
- if (result.size () < z_size) {
142
- result.resize (z_size, z_size);
143
- z.resize (z_size, 0 );
144
- }
145
-
146
- cpp_three_int_backend w (z[0 ]);
147
-
148
- result.limbs ()[0 ] = w.limbs ()[0 ] * p_dash;
149
-
150
- eval_multiply_add (w, result.limbs ()[0 ], this ->m_mod .limbs ()[0 ]);
151
- eval_right_shift (w, sizeof (limb_type) * CHAR_BIT);
152
-
153
- for (size_t i = 1 ; i != p_size; ++i) {
154
- for (size_t j = 0 ; j < i; ++j) {
155
- eval_multiply_add (w, result.limbs ()[j], this ->m_mod .limbs ()[i - j]);
156
- }
157
-
158
- eval_add (w, z[i]);
159
-
160
- result.limbs ()[i] = w.limbs ()[0 ] * p_dash;
161
-
162
- eval_multiply_add (w, result.limbs ()[i], this ->m_mod .limbs ()[0 ]);
163
-
164
- eval_right_shift (w, sizeof (limb_type) * CHAR_BIT);
165
- }
166
-
167
- for (size_t i = 0 ; i != p_size; ++i) {
168
- for (size_t j = i + 1 ; j != p_size; ++j) {
169
- eval_multiply_add (w, result.limbs ()[j], this ->m_mod .limbs ()[p_size + i - j]);
170
- }
171
-
172
- eval_add (w, z[p_size + i]);
173
-
174
- result.limbs ()[i] = w.limbs ()[0 ];
175
-
176
- eval_right_shift (w, sizeof (limb_type) * CHAR_BIT);
177
- }
178
-
179
- eval_add (w, z[z_size - 1 ]);
180
-
181
- result.limbs ()[p_size] = w.limbs ()[0 ];
182
- result.limbs ()[p_size + 1 ] = w.limbs ()[1 ];
183
-
184
- // TODO(mart we cannot resize any more.
185
- if (result.size () != p_size + 1 ) {
186
- result.resize (p_size + 1 , p_size + 1 );
187
- }
188
44
}
189
45
190
46
protected:
191
- Backend m_r2;
192
- limb_type m_p_dash;
193
47
size_t m_p_words;
194
48
195
49
};
0 commit comments