~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/math/mathUtil_rand.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathUtil_rand.c,v $
 37  * Date    :  $Date: 2007/12/07 01:08:01 $
 38  * Version :  $Revision: 1.9 $
 39  * CVS Id  :  $Id: mathUtil_rand.c,v 1.9 2007/12/07 01:08:01 paul Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Random number generation (various)
 44  *
 45  *
 46  *********
 47 */
 48 /** 
 49  *  @file
 50  *  @brief  Various Random number generators.   
 51  *
 52  *
 53  * 
 54 */
 55 
 56 #include "mathUtil_rand.h"
 57 
 58 #if HAVE_CONFIG_H
 59 #include <config.h>
 60 #endif
 61 
 62 
 63 #include <stdio.h>
 64 #include <math.h>
 65 #if STDC_HEADERS
 66 #include <stdlib.h>
 67 #include <stddef.h>
 68 #else
 69 #if HAVE_STDLIB_H
 70 #include <stdlib.h>
 71 #endif                          /* HAVE_STDLIB_H */
 72 #endif
 73 #if TIME_WITH_SYS_TIME
 74 #include <sys/time.h>
 75 #include <time.h>
 76 #else
 77 #if HAVE_SYS_TIME_H
 78 #include <sys/time.h>
 79 #else
 80 #include <time.h>
 81 #endif
 82 #endif                          /* TIME_WITH_SYS_TIME */
 83 
 84 
 85 #include <tina/sys/sysDef.h>
 86 #include <tina/sys/sysPro.h>
 87 #include <tina/math/mathDef.h>
 88 #include <tina/math/mathPro.h>
 89 
 90 /* ipoole */
 91 #ifdef _PCC
 92 long            random()
 93 {
 94         long            result;
 95         result = rand() << 15;
 96         //leaves bit 14 unset...
 97                 result += rand();
 98         return result;
 99 }
100 
101 #endif
102 
103 #define OUR_RAND_MAX (((unsigned long) 1<<31)-1)
104 
105 #define RAND_STEP (1.0/(OUR_RAND_MAX+1.0))
106 #define RAND_ITMAX      100
107 #define RAND_EPS        3.0e-7
108 
109 double          rand_1(void)
110 {
111         return (random() * RAND_STEP);
112 }
113 
114 int             rand_bit(void)
115 {
116         return (random() & 01);
117 }
118 
119 int             rand_int(int a, int b)
120 {
121         return (a + random() % (b - a));
122 }
123 
124 double          rand_unif(double x, double y)
125 {
126         return (x + (y - x) * random() * RAND_STEP);
127 }
128 
129 static double   rand_normal_1(void)
130 {
131         static Bool     saved = false;          /* static data! */
132         static double   v1, v2;                         /* static data! */
133         double          r, k;
134 
135         if (saved == false)
136         {
137                 do
138                 {
139                         v1 = 2.0 * rand_1() - 1.0;
140                         v2 = 2.0 * rand_1() - 1.0;
141                         r = v1 * v1 + v2 * v2;
142                 } while (r >= 1.0 && r != 0.0);
143                 k = sqrt(-2.0 * log(r) / r);
144                 v1 *= k;
145                 v2 *= k;
146                 saved = true;
147                 return (v1);
148         } else
149         {
150                 saved = false;
151                 return (v2);
152         }
153 }
154 
155 double          rand_normal(double mu, double sigma)
156 {
157         return (mu + sigma * rand_normal_1());
158 }
159 
160 double rand_poisson(float xm)
161 {
162    static float sq, alxm, g, oldm=(-1.0);
163    float em, t, y;
164 
165    if (xm < 12.0)
166    {
167       if (xm != oldm)
168       {
169          oldm=xm;
170          g = exp(-xm);
171       }
172       em = -1;
173       t = 1.0;
174       do
175       {
176          em++;
177          t *= rand_unif(0.0, 1.0);
178       } while (t > g);
179    }
180    else
181    {
182       if (xm != oldm)
183       {
184          oldm = xm;
185          sq = sqrt(2.0*xm);
186          alxm=log(xm);
187          g=xm*alxm-gammln(xm+1.0);
188       }
189       do
190       {
191          do
192          {
193             y = tan(PI*rand_unif(0.0, 1.0));
194             em = sq*y+xm;
195          } while (em < 0.0);
196          em = floor(em);
197          t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
198       } while(rand_unif(0.0, 1.0) > t);
199    }
200    return (double)em;
201 }
202 
203 double          gammln(double x)
204 {
205         double          tmp, sum;
206         static double   cof[6] =                                        /* static data! */
207         {76.18009173, -86.50532033, 24.01409822,
208         -1.231739516, 0.120858003e-2, -0.536382e-5};
209         int             j;
210 
211         x -= 1.0;
212         tmp = x + 5.5;
213         tmp -= (x + 0.5) * log(tmp);
214         sum = 1.0;
215         for (j = 0; j <= 5; j++)
216         {
217                 x += 1.0;
218                 sum += cof[j] / x;
219         }
220         return (-tmp + log(2.50662827465 * sum));
221 }
222 
223 static void     gser(double *gamser, double a, double x, double *gln)
224 {
225         int             n;
226         double          sum, del, ap;
227 
228         *gln = gammln(a);
229         if (x <= 0.0)
230         {
231                 if (x < 0.0)
232                         return;
233                 *gamser = 0.0;
234                 return;
235         } else
236         {
237                 ap = a;
238                 del = sum = 1.0 / a;
239                 for (n = 1; n <= RAND_ITMAX; n++)
240                 {
241                         ap += 1.0;
242                         del *= x / ap;
243                         sum += del;
244                         if (fabs(del) < fabs(sum) * RAND_EPS)
245                         {
246                                 *gamser = sum * exp(-x + a * log(x) - (*gln));
247                                 return;
248                         }
249                 }
250                 message("a too large, RAND_ITMAX too small in routine GSER");
251                 return;
252         }
253 }
254 
255 static void     gcf(double *gammcf, double a, double x, double *gln)
256 {
257         int             n;
258         double          gold = 0.0, g, fac = 1.0, b1 = 1.0;
259         double          b0 = 0.0, anf, ana, an, a1, a0 = 1.0;
260 
261         *gln = gammln(a);
262         a1 = x;
263         for (n = 1; n <= RAND_ITMAX; n++)
264         {
265                 an = (double) n;
266                 ana = an - a;
267                 a0 = (a1 + a0 * ana) * fac;
268                 b0 = (b1 + b0 * ana) * fac;
269                 anf = an * fac;
270                 a1 = x * a0 + anf * a1;
271                 b1 = x * b0 + anf * b1;
272                 if (a1)
273                 {
274                         fac = 1.0 / a1;
275                         g = b1 * fac;
276                         if (fabs((g - gold) / g) < RAND_EPS)
277                         {
278                                 *gammcf = exp(-x + a * log(x) - (*gln)) * g;
279                                 return;
280                         }
281                         gold = g;
282                 }
283         }
284 }
285 
286 double          gammp(double a, double x)
287 {
288         double          gamser, gammcf, gln;
289 
290         if (x < 0.0 || a <= 0.0)
291                 return (0.0);
292         if (x < (a + 1.0))
293         {
294                 gser(&gamser, a, x, &gln);
295                 return (gamser);
296         } else
297         {
298                 gcf(&gammcf, a, x, &gln);
299                 return (1.0 - gammcf);
300         }
301 }
302 
303 double          gammq(double a, double x)
304 {
305         double          gamser, gammcf, gln;
306 
307         if (x < 0.0 || a <= 0.0)
308                 return (0.0);
309         if (x < (a + 1.0))
310         {
311                 gser(&gamser, a, x, &gln);
312                 return 1.0 - gamser;
313         } else
314         {
315                 gcf(&gammcf, a, x, &gln);
316                 return (gammcf);
317         }
318 }
319 
320 double          chisq(double x, int n)
321 {
322         if (x == 0.0)
323                 return (1.0);
324         if (n == 0)
325                 return (0.0);
326         return (gammq(n * 0.5, x * 0.5));
327 }
328 
329 
330 
331 /* Seed random with microsecond part of time */
332 void            rand_time_seed(void)
333 {
334         time_t          tm;
335 
336         (void) time(&tm);
337         srand((long) tm);
338 }
339 
340 
341 #undef OUR_RAND_MAX
342 #undef RAND_STEP
343 #undef RAND_ITMAX
344 #undef RAND_EPS
345 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.