From 5fd49349f6fa4717b81603f67bf243ae97babcc8 Mon Sep 17 00:00:00 2001
From: h_vogt <h_vogt>
Date: Sun, 19 Dec 2010 11:05:03 +0000
Subject: [PATCH] update to RTS noise

---
 ChangeLog                  |   5 ++
 src/main.c                 |   1 -
 src/maths/misc/Makefile.am |   3 +-
 src/maths/misc/randnumb.c  |  10 ++-
 src/maths/misc/rnorrexp.c  | 133 -------------------------------------
 visualc/vngspice.vcproj    |  12 ++--
 6 files changed, 19 insertions(+), 145 deletions(-)
 delete mode 100644 src/maths/misc/rnorrexp.c

diff --git a/ChangeLog b/ChangeLog
index a4c1039dc..da4798ff4 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2010-12-19 Holger Vogt
+ * simplify exp variate generation for RTS noise:
+   rnorrexp.c removed, 
+   update to randnumb.c, main.c, /maths/misc/makefile.am, visualc/vngspice.vcproj
+
 2010-12-18 Dietmar Warning
  * devices/bsim*/*ask.c: correct multiplier m for certain currents, conductances, charges 
    and capacitances if they are requested by plotting and printing, e.g. @m1[gds]
diff --git a/src/main.c b/src/main.c
index 15463a521..d0c43df91 100644
--- a/src/main.c
+++ b/src/main.c
@@ -849,7 +849,6 @@ main(int argc, char **argv)
     cp_program = ft_sim->simulator;
 
     srand(getpid());
-    zigset(getpid());
     TausSeed();
 
     /* --- Process command line options --- */
diff --git a/src/maths/misc/Makefile.am b/src/maths/misc/Makefile.am
index 588669247..b5ff770fc 100644
--- a/src/maths/misc/Makefile.am
+++ b/src/maths/misc/Makefile.am
@@ -15,8 +15,7 @@ libmathmisc_la_SOURCES = \
 	scalb.c		\
 	norm.h		\
 	norm.c		\
-	randnumb.c  \
-	rnorrexp.c
+	randnumb.c
 
 
 EXTRA_DIST = test_accuracy.c test_erfc.c
diff --git a/src/maths/misc/randnumb.c b/src/maths/misc/randnumb.c
index 991ac1f00..9e06b66a8 100644
--- a/src/maths/misc/randnumb.c
+++ b/src/maths/misc/randnumb.c
@@ -65,6 +65,7 @@ double CombLCGTaus(void);
 float  CombLCGTaus2(void);
 unsigned int CombLCGTausInt(void);
 unsigned int CombLCGTausInt2(void);
+double exprand(double);
 
 void checkseed(void);
 double drand(void);
@@ -84,7 +85,6 @@ void checkseed(void)
    if (cp_getvar("rndseed", CP_NUM, &newseed)) {
       if ((newseed > 0) && (oldseed != newseed)) {
          srand(newseed); //srandom(newseed);
-         zigset(newseed);
          TausSeed();
          oldseed = newseed;
          printf("Seed value for random number generator is set to %d\n", newseed);
@@ -241,3 +241,11 @@ void rgauss(double* py1, double* py2)
 }	
 
 
+/* return an exponentially distributed random number */
+double exprand(double mean)
+{
+    double expval;
+    checkseed();
+    expval = -log(CombLCGTaus()) * mean;
+    return expval;
+}
diff --git a/src/maths/misc/rnorrexp.c b/src/maths/misc/rnorrexp.c
deleted file mode 100644
index 0d7dd56e6..000000000
--- a/src/maths/misc/rnorrexp.c
+++ /dev/null
@@ -1,133 +0,0 @@
-/* The ziggurat method for RNOR and REXP
-Combine the code below with the main program in which you want
-normal or exponential variates.   Then use of RNOR in any expression
-will provide a standard normal variate with mean zero, variance 1,
-while use of REXP in any expression will provide an exponential variate
-with density exp(-x),x>0.
-Before using RNOR or REXP in your main, insert a command such as
-zigset(86947731 );
-with your own choice of seed value>0, rather than 86947731.
-(If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
-For details of the method, see Marsaglia and Tsang, "The ziggurat method
-for generating random variables", Journ. Statistical Software.
-
-http://www.jstatsoft.org/v05/i08/supp/1
-download 2010-12-18
-*/
-
-#include <math.h>
-static unsigned long jz,jsr=123456789;
-
-#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
-#define UNI (.5f + (signed) SHR3*.2328306e-9f)
-#define IUNI SHR3
-
-static long hz;
-static unsigned long iz, kn[128], ke[256];
-static float wn[128],fn[128], we[256],fe[256];
-
-#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
-#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())
-
-/* prototypes */
-double exprand(double);
-float nfix(void);
-float efix(void);
-void zigset(unsigned long jsrseed);
-
-/* nfix() generates variates from the residue when rejection in RNOR occurs. */
-
-float nfix(void)
-{
-const float r = 3.442620f;     /* The start of the right tail */
-static float x, y;
- for(;;)
-  {  x=hz*wn[iz];      /* iz==0, handles the base strip */
-     if(iz==0)
-       { do{ x=-log(UNI)*0.2904764; y=-log(UNI);}	/* .2904764 is 1/r */
-        while(y+y<x*x);
-        return (hz>0)? r+x : -r-x;
-       }
-                         /* iz>0, handle the wedges of other strips */
-      if( fn[iz]+UNI*(fn[iz-1]-fn[iz]) < exp(-.5*x*x) ) return x;
-
-     /* initiate, try to exit for(;;) for loop*/
-      hz=SHR3;
-      iz=hz&127;
-      if(fabs(hz)<kn[iz]) return (hz*wn[iz]);
-  }
-
-}
-
-/* efix() generates variates from the residue when rejection in REXP occurs. */
-float efix(void)
-{ float x;
- for(;;)
-  {  if(iz==0) return (7.69711f-log(UNI));          /* iz==0 */
-     x=jz*we[iz]; if( fe[iz]+UNI*(fe[iz-1]-fe[iz]) < exp(-x) ) return (x);
-
-      /* initiate, try to exit for(;;) loop */
-   jz=SHR3;
-   iz=(jz&255);
-   if(jz<ke[iz]) return (jz*we[iz]);
-  }
-}
-
-/* return an exponential random number */
-double exprand(double mean)
-{
-    return (double)REXP * mean;
-}
-
-/*--------This procedure sets the seed and creates the tables------*/
-
-void zigset(unsigned long jsrseed)
-{  const double m1 = 2147483648.0, m2 = 4294967296.;
-   double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
-   double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
-   int i;
-   jsr^=jsrseed;
-
-/* Set up tables for RNOR */
-   q=vn/exp(-.5*dn*dn);
-   kn[0]=(dn/q)*m1;
-   kn[1]=0;
-
-   wn[0]=q/m1;
-   wn[127]=dn/m1;
-
-   fn[0]=1.;
-   fn[127]=exp(-.5*dn*dn);
-
-    for(i=126;i>=1;i--)
-    {dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
-     kn[i+1]=(dn/tn)*m1;
-     tn=dn;
-     fn[i]=exp(-.5*dn*dn);
-     wn[i]=dn/m1;
-    }
-
-/* Set up tables for REXP */
-    q = ve/exp(-de);
-    ke[0]=(de/q)*m2;
-    ke[1]=0;
-
-    we[0]=q/m2;
-    we[255]=de/m2;
-
-    fe[0]=1.;
-    fe[255]=exp(-de);
-
-   for(i=254;i>=1;i--)
-  {de=-log(ve/de+exp(-de));
-   ke[i+1]= (de/te)*m2;
-   te=de;
-   fe[i]=exp(-de);
-   we[i]=de/m2;
-  }
-}
-
-
-
-
-
diff --git a/visualc/vngspice.vcproj b/visualc/vngspice.vcproj
index 348a7a7f1..479b20f2d 100644
--- a/visualc/vngspice.vcproj
+++ b/visualc/vngspice.vcproj
@@ -1891,11 +1891,11 @@
 				>
 			</File>
 			<File
-				RelativePath="..\src\include\grid.h"
+				RelativePath="..\src\frontend\plotting\grid.h"
 				>
 			</File>
 			<File
-				RelativePath="..\src\frontend\plotting\grid.h"
+				RelativePath="..\src\include\grid.h"
 				>
 			</File>
 			<File
@@ -2003,11 +2003,11 @@
 				>
 			</File>
 			<File
-				RelativePath="..\src\spicelib\parser\inp.h"
+				RelativePath="..\src\frontend\inp.h"
 				>
 			</File>
 			<File
-				RelativePath="..\src\frontend\inp.h"
+				RelativePath="..\src\spicelib\parser\inp.h"
 				>
 			</File>
 			<File
@@ -8130,10 +8130,6 @@
 				RelativePath="..\src\spicelib\devices\res\restemp.c"
 				>
 			</File>
-			<File
-				RelativePath="..\src\maths\misc\rnorrexp.c"
-				>
-			</File>
 			<File
 				RelativePath="..\src\frontend\runcoms.c"
 				>