-
Notifications
You must be signed in to change notification settings - Fork 0
/
random-unit-vectors.c
148 lines (113 loc) · 3.5 KB
/
random-unit-vectors.c
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
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define M_PI acos(-1.0)
double get_rand();
void box_draw(double *v, int d);
void polar_3d(double *v, int d);
void equal_area_proj(double *v, int d);
void discard_corners(double *v, int d);
void generate_and_write(int d, int count, const char* path, void(*generate)(double *v, int d));
int main(void)
{
srand((unsigned int)time(NULL));
mkdir("output");
mkdir("output/data");
generate_and_write(3, 1000, "output/data/box-draw-1000.dat", box_draw);
generate_and_write(3, 10000, "output/data/box-draw-10000.dat", box_draw);
generate_and_write(3, 1000000, "output/data/box-draw-1000000.dat", box_draw);
generate_and_write(3, 2000000, "output/data/box-draw-2000000.dat", box_draw);
generate_and_write(3, 2000000, "output/data/polar-3d-2000000.dat", polar_3d);
generate_and_write(3, 2000000, "output/data/equal-area-proj-2000000.dat", equal_area_proj);
generate_and_write(3, 2000000, "output/data/discard-corners-2000000.dat", discard_corners);
return 0;
}
double get_rand()
{
/* There are better random number generators out there
* e.g., MT19937, but in the interests of not having
* any dependencies just use the stdlib one here, it's
* good enough. */
return (double)rand()/RAND_MAX;
}
void box_draw(double *v, int d)
{
/* Draw from a d-cube */
for (int i = 0; i < d; i++)
v[i] = 2 * get_rand() - 1;
/* Determine the length */
double length = 0;
for (int i = 0; i < d; i++)
length += v[i] * v[i];
length = sqrt(length);
/* Apply normalisation */
for (int i = 0; i < d; i++)
v[i] /= length;
}
void polar_3d(double *v, int d)
{
/* This only works in d == 3 (although can be extended) */
(void)d;
/* Theta is aziumthutal angle, phi is polar angle */
double theta = 2 * M_PI * get_rand();
double phi = M_PI * get_rand();
/* Note, y and z coordinates switched from the
* presentation in the post in order to to rotate
* the sphere giving a front view instead of a
* top view. The front view is more enlightening. */
v[0] = sin(theta) * cos(phi);
v[1] = cos(theta);
v[2] = sin(theta) * sin(phi);
}
void equal_area_proj(double *v, int d)
{
/* This only works in d == 3 (for now) */
(void)d;
/* Theta is aziumuthal angle */
double theta = 2 * M_PI * get_rand();
double z = 2 * get_rand() - 1;
/* From the inverse Lambert equal area projection */
v[0] = sqrt((1 - z) / 2) * cos(theta);
v[1] = sqrt((1 - z) / 2) * sin(theta);
v[2] = z;
}
void discard_corners(double *v, int d)
{
double length;
/* Do until we get a point inside the sphere */
do
{
for (int i = 0; i < d; i++)
v[i] = 2 * get_rand() - 1;
/* Determine l^2 */
length = 0;
for (int i = 0; i < d; i++)
length += v[i] * v[i];
} while (length > 1);
/* For a unit sphere r = r^2, so when comparing the
* length in the do-loop we really only need to compare
* l^2, save time by taking the square root outside. */
length = sqrt(length);
/* Apply normalisation */
for (int i = 0; i < d; i++)
v[i] /= length;
}
void generate_and_write(int d, int count, const char* path, void(*generate)(double *v, int d))
{
printf("Generaring %s...", path);
/* Watch out, this isn't supported by MSVC, depsite being valid C99. */
double vec[d];
FILE* fp = fopen(path, "w");
for (int k = 0; k < count; ++k)
{
/* Generate a vector using the passed generator */
generate(vec, d);
/* Write it to a file */
for (int i = 0; i < d; ++i)
fprintf(fp, "%lf ", vec[i]);
fprintf(fp, "\n");
}
fclose(fp);
printf(" Done.\n");
}