1 files changed, 95 insertions, 0 deletions

diff --git a/src/mt/random.cc b/src/mt/random.cc
new file mode 100644
index 0000000..87a4c8c
--- /dev/null
+++ b/src/mt/random.cc
@@ -0,0 +1,95 @@
+/*************************************************************************
+ *
+ * HTCd - Copyright (C) 1998-2006 Henrik Rydberg
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#include <mt/random.h>
+#include <mt/mlock.h>
+#include <time.h>
+
+//////////////////////////////////////////////////////
+//
+// This file implements some routines from NumRec
+//
+
+const int NTAB=32; // NOTE: change in header if changing this
+
+const int IA=16807;
+const int IM=2147483647;
+const float AM=1/float(IM);
+const int IQ=127773;
+const int IR=2836;
+const int NDIV=1+(IM-1)/NTAB;
+const float EPS=1.2e-7;
+const float RNMX=1-EPS;
+
+static int idum,iy,iv[32],flag;
+static float extra;
+static int inited;
+static mutex_t mutex;
+
+unsigned long RandomInit(unsigned long seed)
+{
+  inited=1;
+  idum=seed?seed:1;
+  for(int j=NTAB+7;j>=0;j--) {
+    int k=idum/IQ;
+    idum=IA*(idum-k*IQ)-IR*k;
+    if(idum<0) idum+=IM;
+    if(j<NTAB) iv[j]=idum;
+  }
+  iy=iv[0];
+  return seed;
+}
+
+float Uniform()
+{
+  MLOCK(mutex);
+  if(!inited) RandomInit(time(0));
+  int k=idum/IQ;
+  idum=IA*(idum-k*IQ)-IR*k;
+  if(idum<0) idum+=IM;
+  int j=iy/NDIV;
+  iy=iv[j];
+  iv[j]=idum;
+  float temp=AM*iy; if(temp>RNMX) temp=RNMX;
+  return temp;
+}
+
+float Exponential()
+{
+  return -log(Uniform());
+}
+
+float Gaussian()
+{
+  MLOCK(mutex);
+  float v1,v2,rsq;
+  if(flag) {
+    flag=0;
+    return extra;
+  } else {
+    do {
+      v1=2*Uniform()-1;
+      v2=2*Uniform()-1;
+      rsq=v1*v1+v2*v2;
+    } while(rsq>=1||rsq<=0);
+    float fac=sqrt(-2*log(rsq)/rsq);
+    flag=1; extra=v1*fac;
+    return v2*fac;
+  }
+}