1 /*$T src/random.c GC 1.140 10/29/11 13:34:15 */
2
3 /*$I0
4
5 This file is part of CWStudio.
6
7 Copyright 2008-2011 Lukasz Komsta, SP8QED
8
9 CWStudio is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13
14 CWStudio is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with CWStudio. If not, see <http://www.gnu.org/licenses/>.
21
22 */
23 #include "cwgen.h"
24
25 /*
26 =======================================================================================================================
27 Generate standard normal pseudorandom numbers using Box-Muller transform of uniform deviates.
28 =======================================================================================================================
29 */
cw_rand_norm(long int length,unsigned int seed)30 floating *cw_rand_norm(long int length, unsigned int seed)
31 {
32 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
33 int *uniform;
34 long int i;
35 floating *normal, u1, u2, r;
36 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
37
38 if((uniform = cw_rand_int(1000, length + length, seed)) == NULL) return(NULL);
39 if((normal = cw_malloc(length * sizeof(floating))) == NULL) return(NULL);
40 for(i = 0; i < length; i++) {
41 u1 = ((floating) (uniform[2 * i] + 1)) / 1000.0;
42 u2 = ((floating) (uniform[2 * i + 1] + 1)) / 1000.0;
43 r = cw_sqrt(-2 * cw_log(u1));
44 normal[i] = r * cw_cos(6.283185 * u2);
45 }
46
47 cw_free(uniform);
48 return(normal);
49 }
50
51 /*
52 =======================================================================================================================
53 Generate random integers of given shape. Positive shape value increases probability of higher values, negative
54 increase the lower ones.
55 =======================================================================================================================
56 */
cw_rand_shaped(int range,int shape,long int length,unsigned int seed)57 int *cw_rand_shaped(int range, int shape, long int length, unsigned int seed)
58 {
59 /*~~~~~~~~~~~~~~~~~*/
60 int *pvector;
61 long int i;
62 floating *normals;
63 /*~~~~~~~~~~~~~~~~~*/
64
65 /* If shape parameter is equal to zero, the result is equal to cw_rand_int */
66 if(shape == 0) return(cw_rand_int(range, length, seed));
67
68 /* Memory allocation */
69 if((pvector = (int *) cw_malloc(length * sizeof(int))) == NULL) return(NULL);
70
71 /* Generate equivalent length of random numbers */
72 if((normals = cw_rand_norm(length, seed)) == NULL) return(NULL);
73
74 /* Transform to integers */
75 for(i = 0; i < length; i++) {
76 normals[i] = cw_fabs(normals[i]) / (floating) (abs(shape) / 3.0) * (floating) range;
77 pvector[i] = (((int) normals[i]) % range);
78 }
79
80 /* If negative shape, revert the result */
81 if(shape > 0)
82 for(i = 0; i < length; i++) pvector[i] = range - 1 - pvector[i];
83
84 return(pvector);
85 }
86
87 /*
88 =======================================================================================================================
89 Generate random integers of given range, generated using given seed. Algorithm "Multiply-With-Carry" of G.
90 Marsaglia
91 =======================================================================================================================
92 */
cw_rand_int(int range,long int length,unsigned int seed)93 int *cw_rand_int(int range, long int length, unsigned int seed)
94 {
95 /*~~~~~~~~~~~~~~~~~~~~~*/
96 int *pvector;
97 long int i;
98 unsigned int seed2;
99 /*~~~~~~~~~~~~~~~~~~~~~*/
100
101 if((pvector = (int *) cw_malloc(length * sizeof(int))) == NULL) return(NULL);
102 seed2 = seed;
103 for(i = 0; i < length; i++) {
104 seed = 36969 * (seed & 65535) + (seed >> 16);
105 seed2 = 18000 * (seed2 & 65535) + (seed2 >> 16);
106 *(pvector + i) = ((seed << 16) + seed2) % range;
107 }
108
109 return(pvector);
110 }
111
112 /*
113 =======================================================================================================================
114 Generate autocorrelated uniform random floating variables using given seed. Used to simulate a drift (QSB, detune).
115 =======================================================================================================================
116 */
cw_rand_corr(int length,floating corr,unsigned int seed)117 floating *cw_rand_corr(int length, floating corr, unsigned int seed)
118 {
119 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
120 int *runif;
121 int i;
122 floating max, min, x, u, w, *pvector;
123 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
124
125 max = 0;
126 min = 1;
127 if((pvector = (floating *) cw_malloc(length * sizeof(floating))) == NULL) return(NULL);
128 if((runif = cw_rand_int(1000, length, seed)) == NULL) return(NULL);
129 u = 0.9999;
130 for(i = 0; i < length; i++) {
131 w = (floating) * (runif + i) / 1000.0 * corr;
132 x = u + w;
133 if(x <= corr)
134 *(pvector + i) = x * x / (2 * corr);
135 else if((x > corr) && (x <= 1))
136 *(pvector + i) = (2 * x - corr) / 2;
137 else
138 *(pvector + i) = 1 - ((1 + corr - x) * (1 + corr - x)) / (2 * corr);
139 u = *(pvector + i);
140 if(u < min) min = u;
141 if(u > max) max = u;
142 }
143
144 cw_free(runif);
145 for(i = 0; i < length; i++) {
146 *(pvector + i) -= min;
147 *(pvector + i) /= max - min;
148 }
149
150 return(pvector);
151 }
152
153 /*
154 =======================================================================================================================
155 Generate random groups of given charset, started by VVV
156 =======================================================================================================================
157 */
cw_rand_groups(int ngroup,int shape,const char * charset,unsigned int seed)158 char *cw_rand_groups(int ngroup, int shape, const char *charset, unsigned int seed)
159 {
160 /*~~~~~~~~~~~~~~~~~~~~~~~~~~*/
161 int length, charsetlength;
162 char *pgroups;
163 int *randoms, i;
164 /*~~~~~~~~~~~~~~~~~~~~~~~~~~*/
165
166 length = (ngroup + 1) * 6 * sizeof(char);
167 charsetlength = strlen(charset);
168 if((pgroups = (char *) cw_malloc(length + 1)) == NULL) return(NULL);
169 *(pgroups + length) = '\x00';
170 if((randoms = cw_rand_shaped(charsetlength, shape, length, seed)) == NULL) return(NULL);
171 for(i = 0; i < length; i++) strncpy(pgroups + i, &charset[*(randoms + i)], 1);
172 strncpy(pgroups, "vvv =\n", 6);
173 for(i = 2; i < ngroup + 2; i++) {
174 if((i - 1) % 5)
175 *(pgroups + (6 * i) - 1) = ' ';
176 else
177 *(pgroups + (6 * i) - 1) = '\n';
178 }
179
180 cw_free(randoms);
181 return(pgroups);
182 }
183
184 /*
185 =======================================================================================================================
186 Generate group of random words, started by VVV
187 =======================================================================================================================
188 */
cw_rand_words(int nwords,int shape,int wordset,unsigned int seed)189 char *cw_rand_words(int nwords, int shape, int wordset, unsigned int seed)
190 {
191 /*~~~~~~~~~~~~~~~~*/
192 int length;
193 char *pwords;
194 int *randoms, i;
195 /*~~~~~~~~~~~~~~~~*/
196
197 if((randoms = cw_rand_shaped(wordset, shape, nwords, seed)) == NULL) return(NULL);
198 length = 6;
199 for(i = 0; i < nwords; i++) length += strlen(cw_words[randoms[i]]) + 2;
200 if((pwords = (char *) cw_malloc(length + 1)) == NULL) return(NULL);
201 strcpy(pwords, "vvv =\n");
202 for(i = 0; i < nwords; i++) {
203 strcat(pwords, cw_words[*(randoms + i)]);
204 if((i > 0) && ((i % 5) == 0))
205 strcat(pwords, "\n");
206 else
207 strcat(pwords, " ");
208 }
209
210 cw_free(randoms);
211 return(pwords);
212 }
213
214 /*
215 =======================================================================================================================
216 Generate group of random calls, started by VVV
217 =======================================================================================================================
218 */
cw_rand_calls(int ncalls,int shape,unsigned int seed)219 char *cw_rand_calls(int ncalls, int shape, unsigned int seed)
220 {
221 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
222 int length;
223 char *pcalls;
224 int *randoms, i, ind1, ind2;
225 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
226
227 if((randoms = cw_rand_shaped(22720, shape, ncalls, seed)) == NULL) return(NULL);
228 length = 6;
229 for(i = 0; i < ncalls; i++) {
230 ind1 = (int) (randoms[i] / 22720);
231 ind2 = randoms[i] % 22720;
232 length += strlen(cw_calls[ind1][ind2]) + 2;
233 }
234
235 if((pcalls = (char *) cw_malloc(length + 1)) == NULL) return(NULL);
236 strcpy(pcalls, "vvv =\n");
237 for(i = 0; i < ncalls; i++) {
238 ind1 = (int) (randoms[i] / 22720);
239 ind2 = randoms[i] % 22720;
240 strcat(pcalls, cw_calls[ind1][ind2]);
241 if((i > 0) && ((i % 5) == 0))
242 strcat(pcalls, "\n");
243 else
244 strcat(pcalls, " ");
245 }
246
247 cw_free(randoms);
248 return(pcalls);
249 }
250