40 ivec
find(
const bvec &invector)
42 it_assert(invector.size() > 0,
"find(): vector cannot be empty");
43 ivec temp(invector.size());
45 for (
int i = 0;i < invector.size();i++) {
46 if (invector(i) ==
bin(1)) {
51 temp.set_size(pos,
true);
57 #define CREATE_SET_FUNS(typef,typem,name,value) \
58 typef name(int size) \
65 typem name(int rows, int cols) \
67 typem t(rows, cols); \
72 #define CREATE_EYE_FUN(type,name,zero,one) \
73 type name(int size) { \
76 for (int i=0; i<size; i++) \
81 CREATE_SET_FUNS(vec, mat,
ones, 1.0)
83 CREATE_SET_FUNS(ivec, imat,
ones_i, 1)
84 CREATE_SET_FUNS(cvec, cmat,
ones_c, std::complex<
double>(1.0))
86 CREATE_SET_FUNS(vec, mat,
zeros, 0.0)
87 CREATE_SET_FUNS(bvec, bmat,
zeros_b, bin(0))
88 CREATE_SET_FUNS(ivec, imat,
zeros_i, 0)
89 CREATE_SET_FUNS(cvec, cmat,
zeros_c, std::complex<
double>(0.0))
91 CREATE_EYE_FUN(mat,
eye, 0.0, 1.0)
92 CREATE_EYE_FUN(bmat,
eye_b, bin(0), bin(1))
93 CREATE_EYE_FUN(imat,
eye_i, 0, 1)
94 CREATE_EYE_FUN(cmat,
eye_c, std::complex<
double>(0.0), std::complex<
double>(1.0))
97 void eye(
int size, Mat<T> &m)
99 m.set_size(size, size,
false);
101 for (
int i = size - 1; i >= 0; i--)
125 double step = (to - from) /
double(points - 1);
127 for (i = 0; i < points; i++)
128 output(i) = from + i * step;
135 it_assert(K > 0,
"zigzag_space:() K must be positive");
139 for (
int k = 0; k < K; k++) {
141 for (
int i = 1; i <
length(Nn); i += 2) {
150 for (
int i = 0; i < n; i++) {
159 it_assert(size > 0,
"hadamard(): size is not a power of 2");
160 int logsize =
ceil_i(::
log2(static_cast<double>(size)));
161 it_assert(
pow2i(logsize) == size,
"hadamard(): size is not a power of 2");
166 for (
int i = 0; i < logsize; ++i) {
168 for (
int k = 0; k <
pow2; ++k) {
169 for (
int l = 0; l <
pow2; ++l) {
171 H(k + pow2, l) = H(k, l);
172 H(k, l + pow2) = H(k, l);
173 H(k + pow2, l + pow2) = (-1) * H(k, l);
182 int quadratic_residue;
189 for (i = 0; i < (p - 1) / 2; i++) {
190 quadratic_residue = ((i + 1) * (i + 1)) % p;
192 for (j = 0; j < p; j++) {
193 out(j, (j + quadratic_residue) % p) = 1;
198 for (i = 0; i < p; i++) {
211 out.set_submatrix(0, 0, 1, n - 1, 1);
212 out.set_submatrix(1, n - 1, 0, 0, 1);
224 cvec c_conj =
conj(c);
225 for (
int i = 1; i < s; ++i) {
226 for (
int j = 0; j < s - i; ++j) {
227 output(i + j, j) = c_conj(i);
231 for (
int j = 0; j < s; ++j) {
232 for (
int i = 0; i < s - j; ++i) {
233 output(i, i + j) = c(j);
246 plane1 < dim && plane2 < dim && plane1 != plane2,
247 "Invalid arguments to rotation_matrix()");
249 m.set_size(dim, dim,
false);
251 for (
int i = 0; i < dim; i++)
254 m(plane1, plane1) = c;
255 m(plane1, plane2) = -s;
256 m(plane2, plane1) = s;
257 m(plane2, plane2) = c;
262 void house(
const vec &x, vec &v,
double &beta)
273 sigma =
sum(
sqr(x(1, n - 1)));
282 v(0) = -sigma / (x(0) + mu);
283 beta = 2 *
sqr(v(0)) / (sigma +
sqr(v(0)));
288 void givens(
double a,
double b,
double &c,
double &s)
297 if (fabs(b) > fabs(a)) {
323 if (fabs(b) > fabs(a)) {
361 if (fabs(b) > fabs(a)) {
386 template void eye(
int, mat &);
388 template void eye(
int, bmat &);
390 template void eye(
int, imat &);
392 template void eye(
int, cmat &);