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
|
#include "orthogonal.h"
void vector_replace(q_t n, double *v1, const double *v2) {
// writes vector v2 of length n to memory located at v1
for (q_t i = 0; i < n; i++) {
v1[i] = v2[i];
}
}
void vector_add(q_t n, double *v1, const double *v2) {
// adds vector v2 of length n to vector v1
for (q_t i = 0; i < n; i++) {
v1[i] += v2[i];
}
}
void vector_subtract(q_t n, double *v1, const double *v2) {
// subtracts vector v2 of length n from vector v1
for (q_t i = 0; i < n; i++) {
v1[i] -= v2[i];
}
}
double *vector_rotate(q_t n, const double *rot, const double *vec) {
// multiplies n by n rotation matrix rot to vector vec
double *rot_vec = (double *)malloc(n * sizeof(double));
double prod = 0.0;
for (q_t i = 0; i < n; i++) {
prod += rot[i] * vec[i];
}
for (q_t i = 0; i < n; i++) {
rot_vec[i] = vec[i] - 2 * prod * rot[i];
}
return rot_vec;
}
double *vector_rotate_inverse(q_t n, const double *rot, const double *vec) {
double *rot_vec = (double *)calloc(n, sizeof(double));
for (q_t i = 0; i < n; i++) {
for (q_t j = 0; j < n; j++) {
rot_vec[i] += rot[n * j + i] * vec[j];
}
}
return rot_vec;
}
double vector_dot(q_t n, const double *v1, const double *v2) {
double dot = 0;
for (q_t i = 0; i < n; i++) {
dot += v1[i] * v2[i];
}
return dot;
}
double *orthogonal_rotate(q_t n, const double *r, const double *m) {
double *mul = (double *)calloc(n * n, sizeof(double));
for (q_t i = 0; i < n; i++) {
double akOki = 0;
for (q_t k = 0; k < n; k++) {
akOki += r[k] * m[n * k + i];
}
for (q_t j = 0; j < n; j++) {
mul[n * j + i] = m[n * j + i] - 2 * akOki * r[j];
}
}
return mul;
}
double *gen_rot(gsl_rng *r, q_t n) {
double *v = (double *)malloc(n * sizeof(double));
double v2 = 0;
for (q_t i = 0; i < n; i++) {
v[i] = gsl_ran_ugaussian(r);
v2 += v[i] * v[i];
}
double magv = sqrt(v2);
for (q_t i = 0; i < n; i++) {
v[i] /= magv;
}
return v;
}
|
