website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_c/routines/rand.c

   1 /*
   2  *   This file is part of TISEAN
   3  *
   4  *   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   5  *
   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.
  10  *
  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.
  15  *
  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
  19  */
  20 /*Author: Rainer Hegger. Last modified: Feb 12, 2006 */
  21 /* Changes:
  22    Sep 5, 2004 : add extern definition of check_alloc
  23    Feb 12, 2006: add was_set to avoid multiple initialisations
  24 */
  25
  26 #define __RANDOM
  27
  28 #ifndef _STDLIB_H
  29 #include <stdlib.h>
  30 #endif
  31
  32 #ifndef _LIMITS_H
  33 #include <limits.h>
  34 #endif
  35
  36 #ifndef _MATH_H
  37 #include <math.h>
  38 #endif
  39
  40
  41 #ifndef M_PI
  42 #define M_PI 3.1415926535897932385E0
  43 #endif
  44
  45 extern void check_alloc(void*);
  46
  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;
  52
  53 void rnd_init(unsigned long iseed)
  54 {
  55   int i;
  56   unsigned long z,index;
  57
  58   if (rnd_init_was_set == 1)
  59     return ;
  60
  61   rnd_init_was_set=1;
  62
  63   if (sizeof(long) == 8) {
  64     factor=13*13*13*13;
  65     factor=factor*factor*factor*13;
  66   }
  67   else
  68     factor=69069;
  69   lo_limit=(double)ULONG_MAX;
  70
  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)));
  75
  76   rnd_array[0]=rnd1279[0]=iseed;
  77   rnd69=iseed;
  78   index=iseed;
  79   nexti[0]=next1279[0]=1;
  80
  81   for (i=1;i<9689;i++) {
  82     rnd_array[i]=factor*rnd_array[i-1]+1;
  83     nexti[i]=i+1;
  84   }
  85
  86   for (i=1;i<1279;i++) {
  87     rnd1279[i]=factor*rnd1279[i-1]+1;
  88     next1279[i]=i+1;
  89   }
  90   nexti[9688]=next1279[1278]=0;
  91
  92   for (i=1;i<2000;i++) {
  93     index=factor*index+1;
  94     z=rnd1279[((index>>10)%1279)];
  95     z=(z<<10)+(z>>10);
  96     index=factor*index+1;
  97     rnd1279[((index>>10)%1279)] += z;
  98   }
  99
 100   nexti[9688]=next1279[1278]=0;
 101   rndtime=9688;
 102   rndtime1=9688-157;
 103   rndtime2=9688-314;
 104   rndtime3=9688-471;
 105   t1279=1278;
 106   t1279_1=1278-216;
 107   t1279_2=1278-299;
 108   t1279_3=1278-598;
 109 }
 110
 111 unsigned long rnd_long(void)
 112 {
 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];
 119
 120   return rnd_array[rndtime];
 121 }
 122
 123 unsigned long rnd_1279(void)
 124 {
 125   t1279=next1279[t1279];
 126   t1279_1=next1279[t1279_1];
 127   t1279_2=next1279[t1279_2];
 128   t1279_3=next1279[t1279_3];
 129
 130   rnd1279[t1279] += (rnd1279[t1279_1] + rnd1279[t1279_2]
 131                      + rnd1279[t1279_3]);
 132   return rnd1279[t1279];
 133 }
 134
 135 unsigned long rnd69069(void)
 136 {
 137   return (rnd69=rnd69*factor+1);
 138 }
 139
 140 double gaussian(double sigma)
 141 {
 142   static unsigned long gausscount=0;
 143   double x,r,u,phi;
 144   static double y;
 145
 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));
 150     x=r*cos(phi);
 151     y=r*sin(phi);
 152
 153     return x;
 154   }
 155   else
 156     return y;
 157 }