1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
|
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<title>~/ntl-11.4.2/doc/mat_lzz_p.cpp.html</title>
<meta name="Generator" content="Vim/8.0">
<meta name="plugin-version" content="vim7.4_v2">
<meta name="syntax" content="cpp">
<meta name="settings" content="use_css,pre_wrap,no_foldcolumn,expand_tabs,prevent_copy=">
<meta name="colorscheme" content="macvim">
<style type="text/css">
<!--
pre { white-space: pre-wrap; font-family: monospace; color: #000000; background-color: #ffffff; }
body { font-family: monospace; color: #000000; background-color: #ffffff; }
* { font-size: 1em; }
.String { color: #4a708b; }
.PreProc { color: #1874cd; }
.Statement { color: #b03060; font-weight: bold; }
.Comment { color: #0000ee; font-style: italic; }
.Type { color: #008b00; font-weight: bold; }
.Boolean { color: #cd0000; }
-->
</style>
<script type='text/javascript'>
<!--
-->
</script>
</head>
<body>
<pre id='vimCodeElement'>
<span class="Comment">/*</span><span class="Comment">*************************************************************************\</span>
<span class="Comment">MODULE: mat_zz_p</span>
<span class="Comment">SUMMARY:</span>
<span class="Comment">Defines the class mat_zz_p.</span>
<span class="Comment">Note that the modulus p need not be a prime, except as indicated below.</span>
<span class="Comment">IMPLEMENTATION NOTES: </span>
<span class="Comment">Starting with NTL version 9.7.0 (and 9.7.1), many of the routines here have</span>
<span class="Comment">been optimized to take better advantage of specific hardware features available</span>
<span class="Comment">on 64-bit Intel CPU's. Currently, the mul, inv, determinant, solve, gauss,</span>
<span class="Comment">kernel, and image routines are fastest for p up to 23-bits long (assuming the</span>
<span class="Comment">CPU supports AVX instructions). After that, performance degrades in three</span>
<span class="Comment">stages: stage 1: up to 28-bits; stage 2: up to 31-bits; stage 3: 32-bits and</span>
<span class="Comment">up. </span>
<span class="Comment">For primes up to 23-bits, AVX floating point instructions are used. After</span>
<span class="Comment">that, ordinary integer arithmetic is used. In a future version, I may exploit</span>
<span class="Comment">AVX2 integer instructions to get better stage 2 performance. And in the more</span>
<span class="Comment">distant future, AVX512 instructions will be used, when they become available.</span>
<span class="Comment">On older Intel machines, or non-Intel machines that have "long long" support,</span>
<span class="Comment">one still gets optimizations corresponding to the three stages above. On</span>
<span class="Comment">32-bit machines, one still gets three stages, just with smaller crossover</span>
<span class="Comment">points.</span>
<span class="Comment">\*************************************************************************</span><span class="Comment">*/</span>
<span class="PreProc">#include </span><span class="String"><NTL/matrix.h></span>
<span class="PreProc">#include </span><span class="String">"vec_vec_zz_p.h"</span>
<span class="Type">typedef</span> Mat<zz_p> mat_zz_p; <span class="Comment">// backward compatibility</span>
<span class="Type">void</span> add(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> mat_zz_p& B);
<span class="Comment">// X = A + B</span>
<span class="Type">void</span> sub(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> mat_zz_p& B);
<span class="Comment">// X = A - B</span>
<span class="Type">void</span> mul(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> mat_zz_p& B);
<span class="Comment">// X = A * B</span>
<span class="Type">void</span> mul(vec_zz_p& x, <span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> vec_zz_p& b);
<span class="Comment">// x = A * b</span>
<span class="Type">void</span> mul(vec_zz_p& x, <span class="Type">const</span> vec_zz_p& a, <span class="Type">const</span> mat_zz_p& B);
<span class="Comment">// x = a * B</span>
<span class="Type">void</span> mul(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, zz_p b);
<span class="Type">void</span> mul(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">long</span> b);
<span class="Comment">// X = A * b</span>
<span class="Type">void</span> mul(mat_zz_p& X, zz_p a, <span class="Type">const</span> mat_zz_p& B);
<span class="Type">void</span> mul(mat_zz_p& X, <span class="Type">long</span> a, <span class="Type">const</span> mat_zz_p& B);
<span class="Comment">// X = a * B</span>
<span class="Type">void</span> transpose(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A);
mat_zz_p transpose(<span class="Type">const</span> mat_zz_p& A);
<span class="Comment">// X = transpose of A</span>
<span class="Type">void</span> determinant(zz_p& d, <span class="Type">const</span> mat_zz_p& A);
zz_p determinant(<span class="Type">const</span> mat_zz_p& a);
<span class="Comment">// d = determinant(A)</span>
<span class="Type">void</span> solve(zz_p& d, vec_zz_p& x, <span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> vec_zz_p& b);
<span class="Comment">// A is an n x n matrix, b is a length n vector. Computes d = determinant(A).</span>
<span class="Comment">// If d != 0, solves x*A = b (so x and b are treated as a row vectors).</span>
<span class="Type">void</span> solve(zz_p& d, <span class="Type">const</span> mat_zz_p& A, vec_zz_p& x, <span class="Type">const</span> vec_zz_p& b);
<span class="Comment">// A is an n x n matrix, b is a length n vector. Computes d = determinant(A).</span>
<span class="Comment">// If d != 0, solves A*x = b (so x and b are treated as a column vectors).</span>
<span class="Type">void</span> inv(zz_p& d, mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A);
<span class="Comment">// A is an n x n matrix. Computes d = determinant(A). If d != 0,</span>
<span class="Comment">// computes X = A^{-1}.</span>
<span class="Type">void</span> inv(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A);
mat_zz_p inv(<span class="Type">const</span> mat_zz_p& A);
<span class="Comment">// X = A^{-1}; error is raised if A is singular</span>
<span class="Type">void</span> power(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> ZZ& e);
mat_zz_p power(<span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> ZZ& e);
<span class="Type">void</span> power(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">long</span> e);
mat_zz_p power(<span class="Type">const</span> mat_zz_p& A, <span class="Type">long</span> e);
<span class="Comment">// X = A^e; e may be negative (in which case A must be nonsingular).</span>
<span class="Comment">// NOTE: the routines determinant, solve, inv, and power (with negative</span>
<span class="Comment">// exponent) all require that the modulus p is prime: during elimination, if a</span>
<span class="Comment">// non-zero pivot element does not have an inverse, and error is raised. The</span>
<span class="Comment">// following "relaxed" versions of these routines will also work with prime</span>
<span class="Comment">// powers, if the optional parameter relax is true (which is the default).</span>
<span class="Comment">// However, note that in these relaxed routines, if a computed determinant</span>
<span class="Comment">// value is zero, this may not be the true determinant: all that you can assume</span>
<span class="Comment">// is that the true determinant is not invertible mod p. If the parameter</span>
<span class="Comment">// relax==false, then these routines behave identically to their "unrelaxed"</span>
<span class="Comment">// counterparts.</span>
<span class="Type">void</span> relaxed_determinant(zz_p& d, <span class="Type">const</span> mat_zz_p& A, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
zz_p relaxed_determinant(<span class="Type">const</span> mat_zz_p& a, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
<span class="Type">void</span> relaxed_solve(zz_p& d, vec_zz_p& x, <span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> vec_zz_p& b, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
<span class="Type">void</span> relaxed_solve(zz_p& d, <span class="Type">const</span> mat_zz_p& A, vec_zz_p& x, <span class="Type">const</span> vec_zz_p& b, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
<span class="Type">void</span> relaxed_inv(zz_p& d, mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
<span class="Type">void</span> relaxed_inv(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
mat_zz_p relaxed_inv(<span class="Type">const</span> mat_zz_p& A, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
<span class="Type">void</span> relaxed_power(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> ZZ& e, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
mat_zz_p relaxed_power(<span class="Type">const</span> mat_zz_p& A, <span class="Type">const</span> ZZ& e, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
<span class="Type">void</span> relaxed_power(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A, <span class="Type">long</span> e, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
mat_zz_p relaxed_power(<span class="Type">const</span> mat_zz_p& A, <span class="Type">long</span> e, <span class="Type">bool</span> relax=<span class="Boolean">true</span>);
<span class="Type">void</span> sqr(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A);
mat_zz_p sqr(<span class="Type">const</span> mat_zz_p& A);
<span class="Comment">// X = A*A </span>
<span class="Type">void</span> ident(mat_zz_p& X, <span class="Type">long</span> n);
mat_zz_p ident_mat_zz_p(<span class="Type">long</span> n);
<span class="Comment">// X = n x n identity matrix</span>
<span class="Type">long</span> IsIdent(<span class="Type">const</span> mat_zz_p& A, <span class="Type">long</span> n);
<span class="Comment">// test if A is the n x n identity matrix</span>
<span class="Type">void</span> diag(mat_zz_p& X, <span class="Type">long</span> n, zz_p d);
mat_zz_p diag(<span class="Type">long</span> n, zz_p d);
<span class="Comment">// X = n x n diagonal matrix with d on diagonal</span>
<span class="Type">long</span> IsDiag(<span class="Type">const</span> mat_zz_p& A, <span class="Type">long</span> n, zz_p d);
<span class="Comment">// test if X is an n x n diagonal matrix with d on diagonal</span>
<span class="Type">void</span> random(mat_zz_p& x, <span class="Type">long</span> n, <span class="Type">long</span> m); <span class="Comment">// x = random n x m matrix</span>
mat_zz_p random_mat_zz_p(<span class="Type">long</span> n, <span class="Type">long</span> m);
<span class="Type">long</span> gauss(mat_zz_p& M);
<span class="Type">long</span> gauss(mat_zz_p& M, <span class="Type">long</span> w);
<span class="Comment">// Performs unitary row operations so as to bring M into row echelon</span>
<span class="Comment">// form. If the optional argument w is supplied, stops when first w</span>
<span class="Comment">// columns are in echelon form. The return value is the rank (or the</span>
<span class="Comment">// rank of the first w columns).</span>
<span class="Type">void</span> image(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A);
<span class="Comment">// The rows of X are computed as basis of A's row space. X is is row</span>
<span class="Comment">// echelon form</span>
<span class="Type">void</span> kernel(mat_zz_p& X, <span class="Type">const</span> mat_zz_p& A);
<span class="Comment">// Computes a basis for the kernel of the map x -> x*A. where x is a</span>
<span class="Comment">// row vector.</span>
<span class="Comment">// NOTE: the gauss, image, and kernel routines all require that</span>
<span class="Comment">// the modulus p is prime. </span>
<span class="Comment">// miscellaneous:</span>
<span class="Type">void</span> clear(mat_zz_p& a);
<span class="Comment">// x = 0 (dimension unchanged)</span>
<span class="Type">long</span> IsZero(<span class="Type">const</span> mat_zz_p& a);
<span class="Comment">// test if a is the zero matrix (any dimension)</span>
<span class="Comment">// operator notation:</span>
mat_zz_p <span class="Statement">operator</span>+(<span class="Type">const</span> mat_zz_p& a, <span class="Type">const</span> mat_zz_p& b);
mat_zz_p <span class="Statement">operator</span>-(<span class="Type">const</span> mat_zz_p& a, <span class="Type">const</span> mat_zz_p& b);
mat_zz_p <span class="Statement">operator</span>*(<span class="Type">const</span> mat_zz_p& a, <span class="Type">const</span> mat_zz_p& b);
mat_zz_p <span class="Statement">operator</span>-(<span class="Type">const</span> mat_zz_p& a);
<span class="Comment">// matrix/scalar multiplication:</span>
mat_zz_p <span class="Statement">operator</span>*(<span class="Type">const</span> mat_zz_p& a, zz_p b);
mat_zz_p <span class="Statement">operator</span>*(<span class="Type">const</span> mat_zz_p& a, <span class="Type">long</span> b);
mat_zz_p <span class="Statement">operator</span>*(zz_p a, <span class="Type">const</span> mat_zz_p& b);
mat_zz_p <span class="Statement">operator</span>*(<span class="Type">long</span> a, <span class="Type">const</span> mat_zz_p& b);
<span class="Comment">// matrix/vector multiplication:</span>
vec_zz_p <span class="Statement">operator</span>*(<span class="Type">const</span> mat_zz_p& a, <span class="Type">const</span> vec_zz_p& b);
vec_zz_p <span class="Statement">operator</span>*(<span class="Type">const</span> vec_zz_p& a, <span class="Type">const</span> mat_zz_p& b);
<span class="Comment">// assignment operator notation:</span>
mat_zz_p& <span class="Statement">operator</span>+=(mat_zz_p& x, <span class="Type">const</span> mat_zz_p& a);
mat_zz_p& <span class="Statement">operator</span>-=(mat_zz_p& x, <span class="Type">const</span> mat_zz_p& a);
mat_zz_p& <span class="Statement">operator</span>*=(mat_zz_p& x, <span class="Type">const</span> mat_zz_p& a);
mat_zz_p& <span class="Statement">operator</span>*=(mat_zz_p& x, zz_p a);
mat_zz_p& <span class="Statement">operator</span>*=(mat_zz_p& x, <span class="Type">long</span> a);
vec_zz_p& <span class="Statement">operator</span>*=(vec_zz_p& x, <span class="Type">const</span> mat_zz_p& a);
</pre>
</body>
</html>
<!-- vim: set foldmethod=manual : -->
|
