2 * This file is part of TISEAN
4 * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
6 * TISEAN is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
11 * TISEAN is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
16 * You should have received a copy of the GNU General Public License
17 * along with TISEAN; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20 /*Author: Rainer Hegger. Last modified: Feb 12, 2006 */
22 Sep 5, 2004 : add extern definition of check_alloc
23 Feb 12, 2006: add was_set to avoid multiple initialisations
42 #define M_PI 3.1415926535897932385E0
45 extern void check_alloc(void*);
47 static unsigned long *rnd_array,rnd69,*rnd1279,factor;
48 static unsigned long *nexti,rndtime,rndtime1,rndtime2,rndtime3,*next1279;
49 static unsigned long t1279,t1279_1,t1279_2,t1279_3;
50 static double lo_limit;
51 static char rnd_init_was_set=0;
53 void rnd_init(unsigned long iseed)
56 unsigned long z,index;
58 if (rnd_init_was_set == 1)
63 if (sizeof(long) == 8) {
65 factor=factor*factor*factor*13;
69 lo_limit=(double)ULONG_MAX;
71 check_alloc(rnd_array=(unsigned long *)malloc(9689*sizeof(unsigned long)));
72 check_alloc(nexti=(unsigned long *)malloc(9689*sizeof(long)));
73 check_alloc(rnd1279=(unsigned long *)malloc(1279*sizeof(unsigned long)));
74 check_alloc(next1279=(unsigned long *)malloc(1279*sizeof(long)));
76 rnd_array[0]=rnd1279[0]=iseed;
79 nexti[0]=next1279[0]=1;
81 for (i=1;i<9689;i++) {
82 rnd_array[i]=factor*rnd_array[i-1]+1;
86 for (i=1;i<1279;i++) {
87 rnd1279[i]=factor*rnd1279[i-1]+1;
90 nexti[9688]=next1279[1278]=0;
92 for (i=1;i<2000;i++) {
94 z=rnd1279[((index>>10)%1279)];
97 rnd1279[((index>>10)%1279)] += z;
100 nexti[9688]=next1279[1278]=0;
111 unsigned long rnd_long(void)
113 rndtime=nexti[rndtime];
114 rndtime1=nexti[rndtime1];
115 rndtime2=nexti[rndtime2];
116 rndtime3=nexti[rndtime3];
117 rnd_array[rndtime] ^= rnd_array[rndtime1]
118 ^rnd_array[rndtime2]^rnd_array[rndtime3];
120 return rnd_array[rndtime];
123 unsigned long rnd_1279(void)
125 t1279=next1279[t1279];
126 t1279_1=next1279[t1279_1];
127 t1279_2=next1279[t1279_2];
128 t1279_3=next1279[t1279_3];
130 rnd1279[t1279] += (rnd1279[t1279_1] + rnd1279[t1279_2]
132 return rnd1279[t1279];
135 unsigned long rnd69069(void)
137 return (rnd69=rnd69*factor+1);
140 double gaussian(double sigma)
142 static unsigned long gausscount=0;
146 if (!(gausscount++ & 0x1)) {
147 phi=2.0*M_PI*(double)rnd_1279()/lo_limit;
148 u=(double)rnd_1279()/lo_limit;
149 r=sqrt(-2.0*sigma*sigma*log(u));