annotate clustalomega/clustal-omega-1.0.2/src/squid/sre_random.c @ 1:bc707542e5de

Uploaded
author clustalomega
date Thu, 21 Jul 2011 13:35:08 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
1 /* sre_random.c
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
2 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
3 * Portable random number generator, and sampling routines.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
4 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
5 * SRE, Tue Oct 1 15:24:11 2002 [St. Louis]
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
6 * CVS $Id: sre_random.c,v 1.1 2002/10/09 14:26:09 eddy Exp)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
7 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
8
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
9 #include <stdio.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
10 #include <stdlib.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
11 #include <math.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
12 #include "sre_random.h"
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
13
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
14 static int sre_randseed = 42; /* default seed for sre_random() */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
15
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
16 /* Function: sre_random()
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
17 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
18 * Purpose: Return a uniform deviate x, 0.0 <= x < 1.0.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
19 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
20 * sre_randseed is a static variable, set
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
21 * by sre_srandom(). When it is non-zero,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
22 * we re-seed.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
23 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
24 * Implements L'Ecuyer's algorithm for combining output
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
25 * of two linear congruential generators, plus a Bays-Durham
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
26 * shuffle. This is essentially ran2() from Numerical Recipes,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
27 * sans their nonhelpful Rand/McNally-esque code obfuscation.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
28 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
29 * Overflow errors are avoided by Schrage's algorithm:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
30 * az % m = a(z%q) - r(z/q) (+m if <0)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
31 * where q=m/a, r=m%a
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
32 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
33 * Requires that long int's have at least 32 bits.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
34 * This function uses statics and is NOT THREADSAFE.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
35 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
36 * Reference: Press et al. Numerical Recipes in C, 1992.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
37 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
38 * Reliable and portable, but slow. Benchmarks on wrasse,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
39 * using Linux gcc and Linux glibc rand() (see randspeed, in Testsuite):
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
40 * sre_random(): 0.5 usec/call
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
41 * rand(): 0.2 usec/call
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
42 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
43 double
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
44 sre_random(void)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
45 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
46 static long rnd1; /* random number from LCG1 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
47 static long rnd2; /* random number from LCG2 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
48 static long rnd; /* random number we return */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
49 static long tbl[64]; /* table for Bays/Durham shuffle */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
50 long x,y;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
51 int i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
52
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
53 /* Magic numbers a1,m1, a2,m2 from L'Ecuyer, for 2 LCGs.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
54 * q,r derive from them (q=m/a, r=m%a) and are needed for Schrage's algorithm.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
55 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
56 long a1 = 40014;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
57 long m1 = 2147483563;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
58 long q1 = 53668;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
59 long r1 = 12211;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
60
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
61 long a2 = 40692;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
62 long m2 = 2147483399;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
63 long q2 = 52774;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
64 long r2 = 3791;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
65
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
66 if (sre_randseed > 0)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
67 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
68 rnd1 = sre_randseed;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
69 rnd2 = sre_randseed;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
70 /* Fill the table for Bays/Durham */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
71 for (i = 0; i < 64; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
72 x = a1*(rnd1%q1); /* LCG1 in action... */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
73 y = r1*(rnd1/q1);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
74 rnd1 = x-y;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
75 if (rnd1 < 0) rnd1 += m1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
76
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
77 x = a2*(rnd2%q2); /* LCG2 in action... */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
78 y = r2*(rnd2/q2);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
79 rnd2 = x-y;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
80 if (rnd2 < 0) rnd2 += m2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
81
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
82 tbl[i] = rnd1-rnd2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
83 if (tbl[i] < 0) tbl[i] += m1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
84 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
85 sre_randseed = 0; /* drop the flag. */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
86 }/* end of initialization*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
87
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
88
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
89 x = a1*(rnd1%q1); /* LCG1 in action... */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
90 y = r1*(rnd1/q1);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
91 rnd1 = x-y;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
92 if (rnd1 < 0) rnd1 += m1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
93
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
94 x = a2*(rnd2%q2); /* LCG2 in action... */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
95 y = r2*(rnd2/q2);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
96 rnd2 = x-y;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
97 if (rnd2 < 0) rnd2 += m2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
98
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
99 /* Choose our random number from the table... */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
100 i = (int) (((double) rnd / (double) m1) * 64.);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
101 rnd = tbl[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
102 /* and replace with a new number by L'Ecuyer. */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
103 tbl[i] = rnd1-rnd2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
104 if (tbl[i] < 0) tbl[i] += m1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
105
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
106 return ((double) rnd / (double) m1);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
107 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
108
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
109 /* Function: sre_srandom()
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
110 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
111 * Purpose: Initialize with a random seed. Seed must be
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
112 * >= 0 to work; we silently enforce this.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
113 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
114 void
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
115 sre_srandom(int seed)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
116 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
117 if (seed < 0) seed = -1 * seed;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
118 if (seed == 0) seed = 42;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
119 sre_randseed = seed;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
120 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
121
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
122 /* Function: sre_random_positive()
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
123 * Date: SRE, Wed Apr 17 13:34:32 2002 [St. Louis]
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
124 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
125 * Purpose: Assure 0 < x < 1 (positive uniform deviate)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
126 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
127 double
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
128 sre_random_positive(void)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
129 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
130 double x;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
131 do { x = sre_random(); } while (x == 0.0);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
132 return x;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
133 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
134
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
135 /* Function: ExponentialRandom()
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
136 * Date: SRE, Mon Sep 6 21:24:29 1999 [St. Louis]
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
137 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
138 * Purpose: Pick an exponentially distributed random variable
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
139 * 0 > x >= infinity
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
140 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
141 * Args: (void)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
142 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
143 * Returns: x
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
144 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
145 double
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
146 ExponentialRandom(void)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
147 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
148 double x;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
149
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
150 do x = sre_random(); while (x == 0.0);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
151 return -log(x);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
152 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
153
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
154 /* Function: Gaussrandom()
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
155 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
156 * Pick a Gaussian-distributed random variable
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
157 * with some mean and standard deviation, and
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
158 * return it.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
159 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
160 * Based on RANLIB.c public domain implementation.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
161 * Thanks to the authors, Barry W. Brown and James Lovato,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
162 * University of Texas, M.D. Anderson Cancer Center, Houston TX.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
163 * Their implementation is from Ahrens and Dieter, "Extensions
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
164 * of Forsythe's method for random sampling from the normal
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
165 * distribution", Math. Comput. 27:927-937 (1973).
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
166 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
167 * Impenetrability of the code is to be blamed on its FORTRAN/f2c lineage.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
168 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
169 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
170 double
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
171 Gaussrandom(double mean, double stddev)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
172 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
173 static double a[32] = {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
174 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904, 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
175 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
176 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
177 1.862732,2.153875
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
178 };
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
179 static double d[31] = {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
180 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
181 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
182 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
183 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
184 };
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
185 static double t[31] = {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
186 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
187 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
188 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
189 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
190 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
191 };
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
192 static double h[31] = {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
193 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
194 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
195 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
196 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
197 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
198 };
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
199 static long i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
200 static double snorm,u,s,ustar,aa,w,y,tt;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
201
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
202 u = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
203 s = 0.0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
204 if(u > 0.5) s = 1.0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
205 u += (u-s);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
206 u = 32.0*u;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
207 i = (long) (u);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
208 if(i == 32) i = 31;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
209 if(i == 0) goto S100;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
210 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
211 * START CENTER
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
212 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
213 ustar = u-(double)i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
214 aa = *(a+i-1);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
215 S40:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
216 if(ustar <= *(t+i-1)) goto S60;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
217 w = (ustar-*(t+i-1))**(h+i-1);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
218 S50:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
219 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
220 * EXIT (BOTH CASES)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
221 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
222 y = aa+w;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
223 snorm = y;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
224 if(s == 1.0) snorm = -y;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
225 return (stddev*snorm + mean);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
226 S60:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
227 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
228 * CENTER CONTINUED
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
229 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
230 u = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
231 w = u*(*(a+i)-aa);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
232 tt = (0.5*w+aa)*w;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
233 goto S80;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
234 S70:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
235 tt = u;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
236 ustar = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
237 S80:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
238 if(ustar > tt) goto S50;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
239 u = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
240 if(ustar >= u) goto S70;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
241 ustar = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
242 goto S40;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
243 S100:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
244 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
245 * START TAIL
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
246 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
247 i = 6;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
248 aa = *(a+31);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
249 goto S120;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
250 S110:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
251 aa += *(d+i-1);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
252 i += 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
253 S120:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
254 u += u;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
255 if(u < 1.0) goto S110;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
256 u -= 1.0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
257 S140:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
258 w = u**(d+i-1);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
259 tt = (0.5*w+aa)*w;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
260 goto S160;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
261 S150:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
262 tt = u;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
263 S160:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
264 ustar = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
265 if(ustar > tt) goto S50;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
266 u = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
267 if(ustar >= u) goto S150;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
268 u = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
269 goto S140;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
270 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
271
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
272
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
273 /* Functions: DChoose(), FChoose()
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
274 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
275 * Purpose: Make a random choice from a normalized distribution.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
276 * DChoose() is for double-precision vectors;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
277 * FChoose() is for single-precision float vectors.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
278 * Returns the number of the choice.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
279 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
280 int
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
281 DChoose(double *p, int N)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
282 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
283 double roll; /* random fraction */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
284 double sum; /* integrated prob */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
285 int i; /* counter over the probs */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
286
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
287 roll = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
288 sum = 0.0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
289 for (i = 0; i < N; i++)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
290 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
291 sum += p[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
292 if (roll < sum) return i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
293 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
294 return (int) (sre_random() * N); /* bulletproof */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
295 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
296 int
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
297 FChoose(float *p, int N)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
298 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
299 float roll; /* random fraction */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
300 float sum; /* integrated prob */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
301 int i; /* counter over the probs */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
302
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
303 roll = sre_random();
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
304 sum = 0.0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
305 for (i = 0; i < N; i++)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
306 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
307 sum += p[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
308 if (roll < sum) return i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
309 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
310 return (int) (sre_random() * N); /* bulletproof */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
311 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
312
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
313