Eliminated random_lite library

20 files changed, 211 insertions, 4412 deletions

diff --git a/scipy/base/function_base.py b/scipy/base/function_base.py
index a1ed19449..4e88fdba1 100644
--- a/scipy/base/function_base.py
+++ b/scipy/base/function_base.py
@@ -1,4 +1,5 @@
 import types
+import math, operator
 import numeric as _nx
 from numeric import ones, zeros, arange, concatenate, array, asarray, empty
 from umath import pi, multiply, add, arctan2, maximum, minimum, frompyfunc, \
@@ -11,12 +12,55 @@ from _compiled_base import digitize, bincount, _insert
 
 __all__ = ['round','logspace','linspace','fix','mod',
            'select','piecewise','trim_zeros','alen','amax', 'amin', 'ptp',
-           'copy',
+           'copy', 'iterable', 'base_repr', 'binary_repr', 
            'prod','cumprod', 'diff','gradient','angle','unwrap','sort_complex',
            'disp','unique','extract','insert','nansum','nanmax','nanargmax',
            'nanargmin','nanmin', 'vectorize','asarray_chkfinite',
            'average','histogram','bincount','digitize']
 
+_lkup = {'0':'000',
+         '1':'001',
+         '2':'010',
+         '3':'011',
+         '4':'100',
+         '5':'101',
+         '6':'110',
+         '7':'111',
+         'L':''}
+
+def binary_repr(num):
+    """Return the binary representation of the input number as a string.
+
+    This is abuut 25x faster than using base_repr with base 2.
+    """
+    ostr = oct(num)
+    bin = ''
+    for ch in ostr[1:]:
+        bin += _lkup[ch]
+    ind = 0
+    while bin[ind] == '0':
+        ind += 1
+    return bin[ind:]
+
+def base_repr (number, base=2, padding=0):
+    """Return the representation of a number in any given base.
+    """
+    chars = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
+
+    lnb = math.log(base)
+    res = padding*chars[0]
+    if number == 0:
+        return res + chars[0]
+    exponent = int (math.log (number)/lnb)
+    while(exponent >= 0):
+        term = long(base)**exponent
+        lead_digit = int(number / term)
+        res += chars[lead_digit]
+        number -= term*lead_digit
+        exponent -= 1
+    return res
+#end Fernando's utilities
+
 
 
 def logspace(start,stop,num=50,endpoint=1):
@@ -53,7 +97,7 @@ def linspace(start,stop,num=50,endpoint=1,retstep=0):
         return y
 
 def iterable(y):
-    try: len(y)
+    try: iter(y)
     except: return 0
     return 1
 
diff --git a/scipy/base/index_tricks.py b/scipy/base/index_tricks.py
index cf914f0d8..98b872f89 100644
--- a/scipy/base/index_tricks.py
+++ b/scipy/base/index_tricks.py
@@ -2,6 +2,7 @@
 
 import types
 import numeric as _nx
+from numeric import asarray
 __all__ = ['mgrid','ogrid','r_', 'row', 'c_', 'col', 'index_exp',
            'ix_']
 
diff --git a/scipy/base/matrix.py b/scipy/base/matrix.py
index 3710dc4da..34a4673c8 100644
--- a/scipy/base/matrix.py
+++ b/scipy/base/matrix.py
@@ -1,6 +1,7 @@
 
 import numeric as N
 from numeric import ArrayType, concatenate
+from function_base import binary_repr
 import types
 import string as str_
 
@@ -41,25 +42,6 @@ def _convert_from_string(data):
         newdata.append(newrow)
     return newdata
 
-_lkup = {'0':'000',
-         '1':'001',
-         '2':'010',
-         '3':'011',
-         '4':'100',
-         '5':'101',
-         '6':'110',
-         '7':'111'}
-
-def _binary(num):
-    ostr = oct(num)
-    bin = ''
-    for ch in ostr[1:]:
-        bin += _lkup[ch]
-    ind = 0
-    while bin[ind] == '0':
-        ind += 1
-    return bin[ind:]
-
 
 class matrix(N.ndarray):
     __array_priority__ = 10.0
@@ -163,7 +145,7 @@ class matrix(N.ndarray):
                 return result
             # binary decomposition to reduce the number of Matrix
             #  Multiplies for other > 3.
-            beta = _binary(other)
+            beta = binary_repr(other)
             t = len(beta)
             Z,q = x.copy(),0
             while beta[t-q-1] == '0':
diff --git a/scipy/base/numeric.py b/scipy/base/numeric.py
index 07e06620d..fb4a511dd 100644
--- a/scipy/base/numeric.py
+++ b/scipy/base/numeric.py
@@ -19,13 +19,7 @@ frombuffer = multiarray.frombuffer
 where = multiarray.where
 concatenate = multiarray.concatenate
 fastCopyAndTranspose = multiarray._fastCopyAndTranspose
-#def where(condition, x=None, y=None):
-#    """where(condition,x,y) is shaped like condition and has elements of x and
-#    y where condition is respectively true or false.
-#    """
-#    if (x is None) or (y is None):
-#        return nonzero(condition)
-#    return choose(not_equal(condition, 0), (y, x))
+register_dtype = multiarray.register_dtype
 
 def asarray(a, dtype=None):
     """asarray(a,dtype=None) returns a as a NumPy array.  Unlike array(),
@@ -50,7 +44,29 @@ def ensure_array(a, dtype=None):
         except AttributeError:
             a = array(a,dtype,copy=1)  # copy irrelevant
             dtype = None
-            
+
+def isfortran(a):
+    flags = a.flags
+    return flags['FORTRAN'] and a.ndim > 1
+
+# from Fernando Perez's IPython
+def zeros_like(a):
+    """Return an array of zeros of the shape and typecode of a.
+
+    If you don't explicitly need the array to be zeroed, you should instead
+    use empty_like(), which is faster as it only allocates memory."""
+
+    return zeros(a.shape,a.dtype,a.flags['FORTRAN'] and a.ndim > 1)
+
+def empty_like(a):
+    """Return an empty (uninitialized) array of the shape and typecode of a.
+
+    Note that this does NOT initialize the returned array.  If you require
+    your array to be initialized, you should use zeros_like().
+
+    """
+    return empty(a.shape,a.dtype,a.flags['FORTRAN'] and a.ndim > 1)
+# end Fernando's utilities
 
 _mode_from_name_dict = {'v': 0,
                         's' : 1,
@@ -205,14 +221,17 @@ def indices(dimensions, dtype=intp):
         lst.append( add.accumulate(tmp, i, )-1 )
     return array(lst)
 
-def fromfunction(function, dimensions):
+def fromfunction(function, dimensions, **kwargs):
     """fromfunction(function, dimensions) returns an array constructed by
     calling function on a tuple of number grids.  The function should
     accept as many arguments as there are dimensions which is a list of
     numbers indicating the length of the desired output for each axis.
+
+    The function can also accept keyword arguments which will be
+    passed in as well. 
     """
     args = indices(dimensions)
-    return function(*args)
+    return function(*args,**kwargs)
     
 
 from cPickle import load, loads
diff --git a/scipy/base/polynomial.py b/scipy/base/polynomial.py
index b21e30238..d4ce7080f 100644
--- a/scipy/base/polynomial.py
+++ b/scipy/base/polynomial.py
@@ -1,21 +1,43 @@
 ## Automatically adapted for scipy Sep 19, 2005 by convertcode.py
 
 import numeric
+_nx = numeric
 from numeric import *
 from scimath import *
 
+
 from type_check import isscalar
-from twodim_base import diag
+from twodim_base import diag, vander
 from shape_base import hstack, atleast_1d
 from function_base import trim_zeros, sort_complex
+eigvals = None
+lstsq = None
 
-__all__ = ['poly','roots','polyint','polyder','polyadd','polysub','polymul',
-           'polydiv','polyval','poly1d']
- 
-def get_eigval_func():
+def get_linalg_funcs():
+    global eigvals, lstsq
     import scipy.linalg
     eigvals = scipy.linalg.eigvals
-    return eigvals
+    lstsq = scipy.linalg.lstsq
+    return
+
+def _eigvals(arg):
+    try:
+        return eigvals(arg)
+    except TypeError:
+        get_linalg_funcs()
+        return eigvals(arg)
+
+def _lstsq(*args):
+    try:
+        return lstsq(*args)
+    except TypeError:
+        get_linalg_funcs()
+        return lstsq(*args)
+    
+
+__all__ = ['poly','roots','polyint','polyder','polyadd','polysub','polymul',
+           'polydiv','polyval','poly1d','poly1d','polyfit']
+ 
 
 def poly(seq_of_zeros):
     """ Return a sequence representing a polynomial given a sequence of roots.
@@ -31,8 +53,7 @@ def poly(seq_of_zeros):
     seq_of_zeros = atleast_1d(seq_of_zeros)    
     sh = shape(seq_of_zeros)
     if len(sh) == 2 and sh[0] == sh[1]:
-        eig = get_eigval_func()
-        seq_of_zeros=eig(seq_of_zeros)
+        seq_of_zeros=_eigvals(seq_of_zeros)
     elif len(sh) ==1:
         pass
     else:
@@ -64,7 +85,6 @@ def roots(p):
         p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
     """
     # If input is scalar, this makes it an array
-    eig = get_eigval_func()
     p = atleast_1d(p)
     if len(p.shape) != 1:
         raise ValueError,"Input must be a rank-1 array."
@@ -87,7 +107,7 @@ def roots(p):
         # build companion matrix and find its eigenvalues (the roots)
         A = diag(ones((N-2,),p.dtypechar),-1)
         A[0,:] = -p[1:] / p[0]
-        roots = eig(A)
+        roots = _eigvals(A)
     else:
         return array([])
 
@@ -145,6 +165,52 @@ def polyder(p,m=1):
             val = poly1d(val)
         return val
 
+def polyfit(x,y,N):
+    """
+
+    Do a best fit polynomial of order N of y to x.  Return value is a
+    vector of polynomial coefficients [pk ... p1 p0].  Eg, for N=2
+
+      p2*x0^2 +  p1*x0 + p0 = y1
+      p2*x1^2 +  p1*x1 + p0 = y1
+      p2*x2^2 +  p1*x2 + p0 = y2
+      .....
+      p2*xk^2 +  p1*xk + p0 = yk
+
+
+    Method: if X is a the Vandermonde Matrix computed from x (see
+    http://mathworld.wolfram.com/VandermondeMatrix.html), then the
+    polynomial least squares solution is given by the 'p' in
+
+      X*p = y
+
+    where X is a len(x) x N+1 matrix, p is a N+1 length vector, and y
+    is a len(x) x 1 vector
+
+    This equation can be solved as
+
+      p = (XT*X)^-1 * XT * y
+
+    where XT is the transpose of X and -1 denotes the inverse.
+
+    For more info, see
+    http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html,
+    but note that the k's and n's in the superscripts and subscripts
+    on that page.  The linear algebra is correct, however.
+
+    See also polyval
+
+    """
+    x = asarray(x)+0.
+    y = asarray(y)+0.
+    y = reshape(y, (len(y),1))
+    X = vander(x, N+1)
+    c,resids,rank,s = _lstsq(X,y)
+    c.shape = (N+1,)
+    return c
+
+
+
 def polyval(p,x):
     """Evaluate the polynomial p at x.  If x is a polynomial then composition.
 
@@ -156,6 +222,9 @@ def polyval(p,x):
       x can be a sequence and p(x) will be returned for all elements of x.
       or x can be another polynomial and the composite polynomial p(x) will be
       returned.
+
+      Notice:  This can produce inaccurate results for polynomials with significant
+      variability. Use carefully.
     """
     p = asarray(p)
     if isinstance(x,poly1d):
@@ -168,7 +237,7 @@ def polyval(p,x):
     return y
 
 def polyadd(a1,a2):
-    """Adds two polynomials represented as lists
+    """Adds two polynomials represented as sequences
     """
     truepoly = (isinstance(a1,poly1d) or isinstance(a2,poly1d))
     a1,a2 = map(atleast_1d,(a1,a2))
@@ -186,7 +255,7 @@ def polyadd(a1,a2):
     return val
 
 def polysub(a1,a2):
-    """Subtracts two polynomials represented as lists
+    """Subtracts two polynomials represented as sequences
     """
     truepoly = (isinstance(a1,poly1d) or isinstance(a2,poly1d))
     a1,a2 = map(atleast_1d,(a1,a2))
@@ -205,7 +274,7 @@ def polysub(a1,a2):
 
 
 def polymul(a1,a2):
-    """Multiplies two polynomials represented as lists.
+    """Multiplies two polynomials represented as sequences.
     """
     truepoly = (isinstance(a1,poly1d) or isinstance(a2,poly1d))
     val = _nx.convolve(a1,a2)
@@ -213,19 +282,9 @@ def polymul(a1,a2):
         val = poly1d(val)
     return val
 
-def polydiv(a1,a2):
-    """Computes q and r polynomials so that a1(s) = q(s)*a2(s) + r(s)
-    """
-    truepoly = (isinstance(a1,poly1d) or isinstance(a2,poly1d))
-    q, r = deconvolve(a1,a2)
-    while _nx.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
-        r = r[1:]
-    if truepoly:
-        q, r = map(poly1d,(q,r))
-    return q, r
 
 def deconvolve(signal, divisor):
-    """Deconvolves divisor out of signal.
+    """Deconvolves divisor out of signal.  Requires scipy.signal library
     """
     try:
         import scipy.signal
@@ -245,6 +304,18 @@ def deconvolve(signal, divisor):
         rem = num - _nx.convolve(den,quot,mode=2)
     return quot, rem
 
+def polydiv(a1,a2):
+    """Computes q and r polynomials so that a1(s) = q(s)*a2(s) + r(s)
+    """
+    truepoly = (isinstance(a1,poly1d) or isinstance(a2,poly1d))
+    q, r = deconvolve(a1,a2)
+    while _nx.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
+        r = r[1:]
+    if truepoly:
+        q, r = map(poly1d,(q,r))
+    return q, r
+
+
 import re
 _poly_mat = re.compile(r"[*][*]([0-9]*)")
 def _raise_power(astr, wrap=70):
@@ -429,7 +500,7 @@ class poly1d:
     def __getattr__(self, key):
         if key in ['r','roots']:
             return roots(self.coeffs)
-        elif key in ['S1','coef','coefficients']:
+        elif key in ['c','coef','coefficients']:
             return self.coeffs
         elif key in ['o']:
             return self.order
diff --git a/scipy/base/scimath.py b/scipy/base/scimath.py
index 99fe112d5..18d75b2df 100644
--- a/scipy/base/scimath.py
+++ b/scipy/base/scimath.py
@@ -16,6 +16,8 @@ from type_check import isreal
 __all__.extend([key for key in dir(_nx.umath) \
                 if key[0]!='_' and key not in __all__])
 
+_ln2 = log(2.0)
+
 def _tocomplex(arr):
     if arr.dtypechar in ['f', 'h', 'B', 'b','H']:
         return arr.astype('F')
@@ -57,7 +59,7 @@ def log2(x):
     """ Take log base 2 of x.
     """
     x = _fix_real_lt_zero(x)
-    return log(x)/log(2)
+    return log(x)/_ln2
 
 def power(x, p):
     x = _fix_real_lt_zero(x)
diff --git a/scipy/base/src/arrayobject.c b/scipy/base/src/arrayobject.c
index 12e2c9208..e1f3a5576 100644
--- a/scipy/base/src/arrayobject.c
+++ b/scipy/base/src/arrayobject.c
@@ -3849,11 +3849,12 @@ array_descr_get(PyArrayObject *self)
 
 	/* hand this off to the typeobject */
 	/* or give default */
-
-	res = PyObject_GetAttrString((PyObject *)self->descr->typeobj, 
-				      "__array_descr__");
-	if (res) return res;
-	PyErr_Clear();
+	if (PyArray_ISUSERDEF(self)) {
+		res = PyObject_GetAttrString((PyObject *)self->descr->typeobj, 
+					     "__array_descr__");
+		if (res) return res;
+		PyErr_Clear();
+	}
 	/* get default */
 	dobj = PyTuple_New(2);
 	if (dobj == NULL) return NULL;
diff --git a/scipy/base/src/multiarraymodule.c b/scipy/base/src/multiarraymodule.c
index 0d4326b9f..4ad5d03f4 100644
--- a/scipy/base/src/multiarraymodule.c
+++ b/scipy/base/src/multiarraymodule.c
@@ -1117,7 +1117,7 @@ PyArray_Concatenate(PyObject *op, int axis)
 		return NULL;
 	}
 
-	if ((0 < axis) && (axis < MAX_DIMS)) 
+	if ((axis < 0) || ((0 < axis) && (axis < MAX_DIMS)))
 		return _swap_and_concat(op, axis, n);
 	
 	ret = NULL;
diff --git a/scipy/base/twodim_base.py b/scipy/base/twodim_base.py
index 1ac7b43e7..f8cc367c8 100644
--- a/scipy/base/twodim_base.py
+++ b/scipy/base/twodim_base.py
@@ -2,7 +2,8 @@
 
 """
 
-__all__ = ['diag','eye','fliplr','flipud','rot90','tri','triu','tril']
+__all__ = ['diag','eye','fliplr','flipud','rot90','tri','triu','tril',
+           'vander']
 
 from numeric import *
 import sys
@@ -40,19 +41,19 @@ def rot90(m, k=1):
     elif k == 2: return fliplr(flipud(m))
     else: return fliplr(transpose(m))  # k==3
     
-def eye(N, M=None, k=0, typecode='d'):
+def eye(N, M=None, k=0, dtype='d'):
     """ eye returns a N-by-M 2-d array where the  k-th diagonal is all ones, 
         and everything else is zeros.
     """
     if M is None: M = N
     if type(M) == type('d'): 
-        typecode = M
+        dtype = M
         M = N
     m = equal(subtract.outer(arange(N), arange(M)),-k)
-    if typecode is None:
+    if dtype is None:
         return m
     else:
-        return m.astype(typecode)
+        return m.astype(dtype)
 
 def diag(v, k=0):
     """ returns the k-th diagonal if v is a array or returns a array 
@@ -63,9 +64,9 @@ def diag(v, k=0):
     if len(s)==1:
         n = s[0]+abs(k)
         if k > 0:
-            v = concatenate((zeros(k, v.typecode()),v))
+            v = concatenate((zeros(k, v.dtype),v))
         elif k < 0:
-            v = concatenate((v,zeros(-k, v.typecode())))
+            v = concatenate((v,zeros(-k, v.dtype)))
         return eye(n, k=k)*v
     elif len(s)==2:
         v = add.reduce(eye(s[0], s[1], k=k)*v)
@@ -110,4 +111,21 @@ def triu(m, k=0):
 
 
 
+# from matplotlib
+def vander(x, N=None):
+    """
+    X = vander(x,N=None)
+
+    The Vandermonde matrix of vector x.  The i-th column of X is the
+    the i-th power of x.  N is the maximum power to compute; if N is
+    None it defaults to len(x).
+
+    """
+    if N is None: N=len(x)
+    X = ones( (len(x),N), x.dtypechar)
+    for i in range(N-1):
+        X[:,i] = x**(N-i-1)
+    return X
+
+
 
diff --git a/scipy/corelib/random_lite/com.c b/scipy/corelib/random_lite/com.c
deleted file mode 100644
index d68b91f38..000000000
--- a/scipy/corelib/random_lite/com.c
+++ /dev/null
@@ -1,391 +0,0 @@
-#include "ranlib.h"
-#include "Python.h"
-#include <stdio.h>
-/* #include <stdlib.h> */
-void advnst(long k)
-/*
-**********************************************************************
-     void advnst(long k)
-               ADV-a-N-ce ST-ate
-     Advances the state  of  the current  generator  by 2^K values  and
-     resets the initial seed to that value.
-     This is  a  transcription from   Pascal to  Fortran    of  routine
-     Advance_State from the paper
-     L'Ecuyer, P. and  Cote, S. "Implementing  a  Random Number Package
-     with  Splitting   Facilities."  ACM  Transactions  on Mathematical
-     Software, 17:98-111 (1991)
-                              Arguments
-     k -> The generator is advanced by2^K values
-**********************************************************************
-*/
-{
-#define numg 32L
-extern void gsrgs(long getset,long *qvalue);
-extern void gscgn(long getset,long *g);
-extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
-static long g,i,ib1,ib2;
-static long qrgnin;
-/*
-     Abort unless random number generator initialized
-*/
-    gsrgs(0L,&qrgnin);
-    if(qrgnin) goto S10;
-    fputs(" ADVNST called before random generator initialized - ABORT\n",stderr);
-    PyErr_SetString (PyExc_RuntimeError, "Described above.");
-    return;
-S10:
-    gscgn(0L,&g);
-    ib1 = Xa1;
-    ib2 = Xa2;
-    for(i=1; i<=k; i++) {
-        ib1 = mltmod(ib1,ib1,Xm1);
-        if (PyErr_Occurred ()) return;
-        ib2 = mltmod(ib2,ib2,Xm2);
-        if (PyErr_Occurred ()) return;
-    }
-    ib1 = mltmod(ib1,*(Xcg1+g-1),Xm1);
-    if (PyErr_Occurred ()) return;
-    ib2 = mltmod(ib2,*(Xcg2+g-1),Xm2);
-    if (PyErr_Occurred ()) return;
-    setsd(ib1,ib2);
-/*
-     NOW, IB1 = A1**K AND IB2 = A2**K
-*/
-#undef numg
-}
-void getsd(long *iseed1,long *iseed2)
-/*
-**********************************************************************
-     void getsd(long *iseed1,long *iseed2)
-               GET SeeD
-     Returns the value of two integer seeds of the current generator
-     This  is   a  transcription from  Pascal   to  Fortran  of routine
-     Get_State from the paper
-     L'Ecuyer, P. and  Cote,  S. "Implementing a Random Number  Package
-     with   Splitting Facilities."  ACM  Transactions   on Mathematical
-     Software, 17:98-111 (1991)
-                              Arguments
-     iseed1 <- First integer seed of generator G
-     iseed2 <- Second integer seed of generator G
-**********************************************************************
-*/
-{
-#define numg 32L
-extern void gsrgs(long getset,long *qvalue);
-extern void gscgn(long getset,long *g);
-extern long Xcg1[],Xcg2[];
-static long g;
-static long qrgnin;
-/*
-     Abort unless random number generator initialized
-*/
-    gsrgs(0L,&qrgnin);
-    if(qrgnin) goto S10;
-    fprintf(stderr,"%s\n",
-      " GETSD called before random number generator  initialized -- abort!\n");
-    PyErr_SetString (PyExc_RuntimeError, "Described above.");
-    return;
-S10:
-    gscgn(0L,&g);
-    *iseed1 = *(Xcg1+g-1);
-    *iseed2 = *(Xcg2+g-1);
-#undef numg
-}
-long ignlgi(void)
-/*
-**********************************************************************
-     long ignlgi(void)
-               GeNerate LarGe Integer
-     Returns a random integer following a uniform distribution over
-     (1, 2147483562) using the current generator.
-     This is a transcription from Pascal to Fortran of routine
-     Random from the paper
-     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
-     with Splitting Facilities." ACM Transactions on Mathematical
-     Software, 17:98-111 (1991)
-**********************************************************************
-*/
-{
-#define numg 32L
-extern void gsrgs(long getset,long *qvalue);
-extern void gssst(long getset,long *qset);
-extern void gscgn(long getset,long *g);
-extern void inrgcm(void);
-extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
-extern long Xqanti[];
-static long ignlgi,curntg,k,s1,s2,z;
-static long qqssd,qrgnin;
-/*
-     IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO.
-     IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO
-     THIS ROUTINE  2) A CALL TO SETALL.
-*/
-    gsrgs(0L,&qrgnin);
-    if(!qrgnin) inrgcm();
-    gssst(0,&qqssd);
-    if(!qqssd) setall(1234567890L,123456789L);
-/*
-     Get Current Generator
-*/
-    gscgn(0L,&curntg);
-    s1 = *(Xcg1+curntg-1);
-    s2 = *(Xcg2+curntg-1);
-    k = s1/53668L;
-    s1 = Xa1*(s1-k*53668L)-k*12211;
-    if(s1 < 0) s1 += Xm1;
-    k = s2/52774L;
-    s2 = Xa2*(s2-k*52774L)-k*3791;
-    if(s2 < 0) s2 += Xm2;
-    *(Xcg1+curntg-1) = s1;
-    *(Xcg2+curntg-1) = s2;
-    z = s1-s2;
-    if(z < 1) z += (Xm1-1);
-    if(*(Xqanti+curntg-1)) z = Xm1-z;
-    ignlgi = z;
-    return ignlgi;
-#undef numg
-}
-void initgn(long isdtyp)
-/*
-**********************************************************************
-     void initgn(long isdtyp)
-          INIT-ialize current G-e-N-erator
-     Reinitializes the state of the current generator
-     This is a transcription from Pascal to Fortran of routine
-     Init_Generator from the paper
-     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
-     with Splitting Facilities." ACM Transactions on Mathematical
-     Software, 17:98-111 (1991)
-                              Arguments
-     isdtyp -> The state to which the generator is to be set
-          isdtyp = -1  => sets the seeds to their initial value
-          isdtyp =  0  => sets the seeds to the first value of
-                          the current block
-          isdtyp =  1  => sets the seeds to the first value of
-                          the next block
-**********************************************************************
-*/
-{
-#define numg 32L
-extern void gsrgs(long getset,long *qvalue);
-extern void gscgn(long getset,long *g);
-extern long Xm1,Xm2,Xa1w,Xa2w,Xig1[],Xig2[],Xlg1[],Xlg2[],Xcg1[],Xcg2[];
-static long g;
-static long qrgnin;
-/*
-     Abort unless random number generator initialized
-*/
-    gsrgs(0L,&qrgnin);
-    if(qrgnin) goto S10;
-    fprintf(stderr,"%s\n",
-      " INITGN called before random number generator  initialized -- abort!");
-    PyErr_SetString (PyExc_RuntimeError, "Described above.");
-    return;
-S10:
-    gscgn(0L,&g);
-    if(-1 != isdtyp) goto S20;
-    *(Xlg1+g-1) = *(Xig1+g-1);
-    *(Xlg2+g-1) = *(Xig2+g-1);
-    goto S50;
-S20:
-    if(0 != isdtyp) goto S30;
-    goto S50;
-S30:
-/*
-     do nothing
-*/
-    if(1 != isdtyp) goto S40;
-    *(Xlg1+g-1) = mltmod(Xa1w,*(Xlg1+g-1),Xm1);
-    if (PyErr_Occurred ()) return;
-    *(Xlg2+g-1) = mltmod(Xa2w,*(Xlg2+g-1),Xm2);
-    if (PyErr_Occurred ()) return;
-    goto S50;
-S40:
-    fprintf(stderr,"%s\n","isdtyp not in range in INITGN");
-    PyErr_SetString (PyExc_ValueError, "Described above.");
-    return;
-S50:
-    *(Xcg1+g-1) = *(Xlg1+g-1);
-    *(Xcg2+g-1) = *(Xlg2+g-1);
-#undef numg
-}
-void inrgcm(void)
-/*
-**********************************************************************
-     void inrgcm(void)
-          INitialize Random number Generator CoMmon
-                              Function
-     Initializes common area  for random number  generator.  This saves
-     the  nuisance  of  a  BLOCK DATA  routine  and the  difficulty  of
-     assuring that the routine is loaded with the other routines.
-**********************************************************************
-*/
-{
-#define numg 32L
-extern void gsrgs(long getset,long *qvalue);
-extern long Xm1,Xm2,Xa1,Xa2,Xa1w,Xa2w,Xa1vw,Xa2vw;
-extern long Xqanti[];
-static long T1;
-static long i;
-/*
-     V=20;                            W=30;
-     A1W = MOD(A1**(2**W),M1)         A2W = MOD(A2**(2**W),M2)
-     A1VW = MOD(A1**(2**(V+W)),M1)    A2VW = MOD(A2**(2**(V+W)),M2)
-   If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
-    An efficient way to precompute a**(2*j) MOD m is to start with
-    a and square it j times modulo m using the function MLTMOD.
-*/
-    Xm1 = 2147483563L;
-    Xm2 = 2147483399L;
-    Xa1 = 40014L;
-    Xa2 = 40692L;
-    Xa1w = 1033780774L;
-    Xa2w = 1494757890L;
-    Xa1vw = 2082007225L;
-    Xa2vw = 784306273L;
-    for(i=0; i<numg; i++) *(Xqanti+i) = 0;
-    T1 = 1;
-/*
-     Tell the world that common has been initialized
-*/
-    gsrgs(1L,&T1);
-#undef numg
-}
-void setall(long iseed1,long iseed2)
-/*
-**********************************************************************
-     void setall(long iseed1,long iseed2)
-               SET ALL random number generators
-     Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
-     initial seeds of the other generators are set accordingly, and
-     all generators states are set to these seeds.
-     This is a transcription from Pascal to Fortran of routine
-     Set_Initial_Seed from the paper
-     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
-     with Splitting Facilities." ACM Transactions on Mathematical
-     Software, 17:98-111 (1991)
-                              Arguments
-     iseed1 -> First of two integer seeds
-     iseed2 -> Second of two integer seeds
-**********************************************************************
-*/
-{
-#define numg 32L
-extern void gsrgs(long getset,long *qvalue);
-extern void gssst(long getset,long *qset);
-extern void gscgn(long getset,long *g);
-extern long Xm1,Xm2,Xa1vw,Xa2vw,Xig1[],Xig2[];
-static long T1;
-static long g,ocgn;
-static long qrgnin;
-    T1 = 1;
-/*
-     TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
-      HAS BEEN CALLED.
-*/
-    gssst(1,&T1);
-    gscgn(0L,&ocgn);
-/*
-     Initialize Common Block if Necessary
-*/
-    gsrgs(0L,&qrgnin);
-    if(!qrgnin) inrgcm();
-    *Xig1 = iseed1;
-    *Xig2 = iseed2;
-    initgn(-1L);
-    for(g=2; g<=numg; g++) {
-        *(Xig1+g-1) = mltmod(Xa1vw,*(Xig1+g-2),Xm1);
-        if (PyErr_Occurred ()) return;
-        *(Xig2+g-1) = mltmod(Xa2vw,*(Xig2+g-2),Xm2);
-        if (PyErr_Occurred ()) return;
-        gscgn(1L,&g);
-        initgn(-1L);
-    }
-    gscgn(1L,&ocgn);
-#undef numg
-}
-void setant(long qvalue)
-/*
-**********************************************************************
-     void setant(long qvalue)
-               SET ANTithetic
-     Sets whether the current generator produces antithetic values.  If
-     X   is  the value  normally returned  from  a uniform [0,1] random
-     number generator then 1  - X is the antithetic  value. If X is the
-     value  normally  returned  from a   uniform  [0,N]  random  number
-     generator then N - 1 - X is the antithetic value.
-     All generators are initialized to NOT generate antithetic values.
-     This is a transcription from Pascal to Fortran of routine
-     Set_Antithetic from the paper
-     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
-     with Splitting Facilities." ACM Transactions on Mathematical
-     Software, 17:98-111 (1991)
-                              Arguments
-     qvalue -> nonzero if generator G is to generating antithetic
-                    values, otherwise zero
-**********************************************************************
-*/
-{
-#define numg 32L
-extern void gsrgs(long getset,long *qvalue);
-extern void gscgn(long getset,long *g);
-extern long Xqanti[];
-static long g;
-static long qrgnin;
-/*
-     Abort unless random number generator initialized
-*/
-    gsrgs(0L,&qrgnin);
-    if(qrgnin) goto S10;
-    fprintf(stderr,"%s\n",
-      " SETANT called before random number generator  initialized -- abort!");
-    PyErr_SetString (PyExc_RuntimeError, "Described above.");
-    return;
-S10:
-    gscgn(0L,&g);
-    Xqanti[g-1] = qvalue;
-#undef numg
-}
-void setsd(long iseed1,long iseed2)
-/*
-**********************************************************************
-     void setsd(long iseed1,long iseed2)
-               SET S-ee-D of current generator
-     Resets the initial  seed of  the current  generator to  ISEED1 and
-     ISEED2. The seeds of the other generators remain unchanged.
-     This is a transcription from Pascal to Fortran of routine
-     Set_Seed from the paper
-     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
-     with Splitting Facilities." ACM Transactions on Mathematical
-     Software, 17:98-111 (1991)
-                              Arguments
-     iseed1 -> First integer seed
-     iseed2 -> Second integer seed
-**********************************************************************
-*/
-{
-#define numg 32L
-extern void gsrgs(long getset,long *qvalue);
-extern void gscgn(long getset,long *g);
-extern long Xig1[],Xig2[];
-static long g;
-static long qrgnin;
-/*
-     Abort unless random number generator initialized
-*/
-    gsrgs(0L,&qrgnin);
-    if(qrgnin) goto S10;
-    fprintf(stderr,"%s\n",
-      " SETSD called before random number generator  initialized -- abort!");
-    PyErr_SetString (PyExc_RuntimeError, "Described above.");
-    return;
-S10:
-    gscgn(0L,&g);
-    *(Xig1+g-1) = iseed1;
-    *(Xig2+g-1) = iseed2;
-    initgn(-1L);
-#undef numg
-}
-long Xm1,Xm2,Xa1,Xa2,Xcg1[32],Xcg2[32],Xa1w,Xa2w,Xig1[32],Xig2[32],Xlg1[32],
-    Xlg2[32],Xa1vw,Xa2vw;
-long Xqanti[32];
diff --git a/scipy/corelib/random_lite/linpack.c b/scipy/corelib/random_lite/linpack.c
deleted file mode 100644
index 9d5f3b44c..000000000
--- a/scipy/corelib/random_lite/linpack.c
+++ /dev/null
@@ -1,91 +0,0 @@
-#include <math.h>
-float sdot(long n,float *sx,long incx,float *sy,long incy)
-{
-static long i,ix,iy,m,mp1;
-static float sdot,stemp;
-    stemp = sdot = 0.0;
-    if(n <= 0) return sdot;
-    if(incx == 1 && incy == 1) goto S20;
-    ix = iy = 1;
-    if(incx < 0) ix = (-n+1)*incx+1;
-    if(incy < 0) iy = (-n+1)*incy+1;
-    for(i=1; i<=n; i++) {
-        stemp += (*(sx+ix-1)**(sy+iy-1));
-        ix += incx;
-        iy += incy;
-    }
-    sdot = stemp;
-    return sdot;
-S20:
-    m = n % 5L;
-    if(m == 0) goto S40;
-    for(i=0; i<m; i++) stemp += (*(sx+i)**(sy+i));
-    if(n < 5) goto S60;
-S40:
-    mp1 = m+1;
-    for(i=mp1; i<=n; i+=5) stemp += (*(sx+i-1)**(sy+i-1)+*(sx+i)**(sy+i)+*(sx+i
-      +1)**(sy+i+1)+*(sx+i+2)**(sy+i+2)+*(sx+i+3)**(sy+i+3));
-S60:
-    sdot = stemp;
-    return sdot;
-}
-void spofa(float *a,long lda,long n,long *info)
-/*
-     SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX.
-     SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED
-     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
-     (TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) .
-     ON ENTRY
-        A       REAL(LDA, N)
-                THE SYMMETRIC MATRIX TO BE FACTORED.  ONLY THE
-                DIAGONAL AND UPPER TRIANGLE ARE USED.
-        LDA     INTEGER
-                THE LEADING DIMENSION OF THE ARRAY  A .
-        N       INTEGER
-                THE ORDER OF THE MATRIX  A .
-     ON RETURN
-        A       AN UPPER TRIANGULAR MATRIX  R  SO THAT  A = TRANS(R)*R
-                WHERE  TRANS(R)  IS THE TRANSPOSE.
-                THE STRICT LOWER TRIANGLE IS UNALTERED.
-                IF  INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
-        INFO    INTEGER
-                = 0  FOR NORMAL RETURN.
-                = K  SIGNALS AN ERROR CONDITION.  THE LEADING MINOR
-                     OF ORDER  K  IS NOT POSITIVE DEFINITE.
-     LINPACK.  THIS VERSION DATED 08/14/78 .
-     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
-     SUBROUTINES AND FUNCTIONS
-     BLAS SDOT
-     FORTRAN SQRT
-     INTERNAL VARIABLES
-*/
-{
-extern float sdot(long n,float *sx,long incx,float *sy,long incy);
-static long j,jm1,k;
-static float t,s;
-/*
-     BEGIN BLOCK WITH ...EXITS TO 40
-*/
-    for(j=1; j<=n; j++) {
-        *info = j;
-        s = 0.0;
-        jm1 = j-1;
-        if(jm1 < 1) goto S20;
-        for(k=0; k<jm1; k++) {
-            t = *(a+k+(j-1)*lda)-sdot(k,(a+k*lda),1L,(a+(j-1)*lda),1L);
-            t /=  *(a+k+k*lda);
-            *(a+k+(j-1)*lda) = t;
-            s += (t*t);
-        }
-S20:
-        s = *(a+j-1+(j-1)*lda)-s;
-/*
-     ......EXIT
-*/
-        if(s <= 0.0) goto S40;
-        *(a+j-1+(j-1)*lda) = sqrt(s);
-    }
-    *info = 0;
-S40:
-    return;
-}
diff --git a/scipy/corelib/random_lite/pmath_rng.c b/scipy/corelib/random_lite/pmath_rng.c
deleted file mode 100644
index add310c50..000000000
--- a/scipy/corelib/random_lite/pmath_rng.c
+++ /dev/null
@@ -1,473 +0,0 @@
-/* pmath_rng.c
- *
- * This set of routines reproduces the Cray RANF family of random
- * number routines bit-for-bit.  Most of the routines here were
- * written by Jim Rathkopf, CP-Division, LLNL, and modified by
- * William C. Biester, B-Division. LLNL.  They were subsequently
- * modified for use in PMATH by Fred N. Fritsch, LC, LLNL.
- *
- * This version contains only the generator proper (equivalent to PM_RANF8
- * in the official PMATH version) and makes the seed and multiplier
- * directly accessible as arrays of doubles.  This was modified for use by
- * the Basis and Python module ranf.c, which contains its own seed and
- * multiplier management routines.  (FNF -- 9/4/96)
- */
-
-/* -----------------------------------------------------------------------
-
-Copyright (c) 1993,1994,1995,1996.
-The Regents of the University of California.  All rights reserved.
-
-This work was produced at the University of California, Lawrence
-Livermore National Laboratory (UC LLNL) under contract no. W-7405-ENG-48
-between the U.S. Department of Energy (DOE) and The Regents of the
-University of California (University) for the operation of UC LLNL.
-Copyright is reserved to the University for purposes of controlled
-dissemination, commercialization through formal licensing, or other
-disposition under terms of Contract 48; DOE policies, regulations and
-orders; and U.S. statutes.
-
-                         DISCLAIMER
- 
-This document was prepared as an account of work sponsored by an agency
-of the United States Government.  Neither the United States Government
-nor the University of California nor any of their employees, makes any
-warranty, express or implied, or assumes any liability or responsibility
-for the accuracy, completeness, or usefulness of any information,
-apparatus, product, or process disclosed, or represents that its use
-would not infringe privately-owned rights.  Reference herein to any
-specific commercial products, process, or service by trade name,
-trademark, manufacturer, or otherwise, does not necessarily constitute
-or imply its endorsement, recommendation, or favoring by the United
-States Government or the University of California.  The views and
-opinions of authors expressed herein do not necessarily state or reflect
-those of the United States Government or the University of California,
-and shall not be used for advertising or product endorsement purposes.
-*/
-
-/* -----------------------------------------------------------------------
- *
- * Service routines defined:
- *    PM_16to24 - Convert 3 16-bit shorts to 2 24-bit doubles
- *    PM_24to16 - Convert 2 24-bit doubles to 3 16-bit shorts
- *    PM_GSeed  - Get the seed as two 24-bit doubles
- *    PM_SSeed  - Set the seed from two 24-bit doubles (unsafe)
- *    PM_GMult  - Get the multiplier as two 24-bit doubles
- *    PM_SMult  - Set the multiplier from two 24-bit doubles (unsafe)
- *
- * User-callable routines defined:
- *    PM_RANF   - The generator itself
- *
- * Currently responsible person:
- *    Fred N. Fritsch
- *    Computer Applications Organization
- *    LCPD, ICF Group
- *    Lawrence Livermore National Laboratory
- *    fnf@llnl.gov
- *
- * Modifications:
- * (See individual SLATEC-format prologues for detailed modification records.)
- *    03-09-93  Integrated Jim's routines into my library system.  (WCB)
- *    10-21-93  Modified to provide RANF8 for PMATH.  (FNF)
- *    10-22-93  Added code for RNFCNT.  (FNF)
- *    11-09-93  Added code for RNMSET.  (FNF)
- *    12-15-93  Modified to process bytes of seeds and multipliers in
- *              in the reverse order, to agree with the way they come
- *              from the Cray.  (FNF)
- *    03-15-94  Changed general pm_binds.h to pm_rnfset.h.  (FNF)
- *    04-25-94  Added SLATEC-format prologues. (FNF)
- *    09-26-95  Simplified code and improved internal documentation.
- *              This included removing unnecessary internal procedures
- *              setseed and seed48. (FNF)
- *    09-27-95  Simplified further by eliminating internal procedure
- *              rnset and no longer used variable Multiplier.  Also,
- *              moved PM_RANF8 so that it appears first. (FNF)
- *    09-28-95  Modified behavior of RNMSET to agree with MATHLIB. (FNF)
- *    10-02-95  Eliminated unions via calls to CV16TO64 and CV64TO16.
- *              Most conditional compilations are now in the rand48_*
- *              procedures. (FNF)
- *    10-05-95  Corrected "overflow" problem to make code work on Cray.
- *              Removed some unnecessary code in the process. (FNF)
- *    10-09-95  Removed conditional compiles.  (No longer needed after
- *              change of 10-02-95.) (FNF)
- *    08-27-96  Removed code for PM_RNSSET, PM_RNMSET, PM_RNSGET, PM_RNFCNT.
- *              Changed name PM_RANF8 to PM_RANF.  Added simple interface
- *              routines to get/set seed and multiplier.  Removed dependency
- *              on CV16TO64 and CV64TO16.  Eliminated Fortran bindings.  (FNF)
- *    09-04-96  Improved descriptions of the get/set routines.  (FNF)
- *    09-05-96  Added definitions of internal procedures PM_16to24 and
- *              PM_24to16 (see ranf.h).  (FNF)
- * -----------------------------------------------------------------------
- */
-
-
-
-   /* Get includes, typedefs, and protypes for ranf.c and pmath_rng.c */
-#include "ranf.h"
-
-#define BASE24        16777216.0  /* 2^24 */
-
-#define TWO8               256.0  /* 2^8 */
-#define TWO16            65536.0  /* 2^16 */
-#define TWO48  281474976710656.0  /* 2^48 */
-
-/* ------------------------------------------------------
-   The following static global variables retain the state
-   of the random number generator between calls.
-   ------------------------------------------------------ */
-/* The default values to reproduce Cray values */
-static double dseed[2] = {16555217.0, 9732691.0} ;
-static double dmult[2] = {15184245.0, 2651554.0} ;
-static double dadder = 0.0;    /* This is a relic from drand48. */
-
-
-/* --------------------------------------
-      Define in-line mod function.
-   -------------------------------------- */
-#define MODF(r,m)  (r - (floor(r/m))*m)
-
-
-/* -----------------------------------------
-      Define PMATH RNG internal procedures.
-   -----------------------------------------
-*/
-
-/* PM_16to24 -- Take a number stored in 3 16-bit unsigned shorts and move it
- *              to 2 doubles each containing 24 bits of data.
- * 
- * This is a local function.
- *
- * Calling sequence:
- *   u16 x16[3];
- *   double x24[2];
- *   PM_16to24 (x16, x24);
- *
- *   x16 :   the 3 16-bit unsigned shorts.
- *   x24 :=  the 2 doubles of 24 bits returned.
- *
- * Return value:
- *   none.
- *
- * Note:
- *   This is the same as PMATH internal procedure rand48_16to24 except
- *   that the order of the entries in x16 have been reversed and its
- *   type declaration has been changed to u16.  It is intended for use
- *   by user-callable routines that set seeds or multipliers for users.
- */
-void PM_16to24(u16 x16[3], double x24[2])
-{
-   double t0, t1, t2, t1u, t1l;
-
-   t0 = (double) x16[0];
-   t1 = (double) x16[1];
-   t2 = (double) x16[2];
-
-   t1u = floor(t1/TWO8);
-   t1l = t1 -t1u*TWO8;
-   x24[0] = t0 + t1l*TWO16;
-   x24[1] = t1u + t2*TWO8;
-
-#ifdef RAN_DEBUG
-fprintf(stderr,"PM_16to24: x16 = %04x %04x %04x\n",x16[2],x16[1],x16[0]);
-fprintf(stderr,"PM_16to24: x24 = %.1f %.1f\n",x24[1],x24[0]);
-#endif
-   return;
-}
-
-
-/* PM_24to16 -- Take a number stored in 2 doubles each containing 24 bits
- *              of data and move it to 3 16-bit unsigned shorts. 
- *
- * Calling sequence:
- *   double x24[2];
- *   u16 x16[3];
- *   PM_24to16 (x24, x16);
- *
- *   x24 :   the 2 doubles of 24 bits.
- *   x16 :=  the 3 16-bit unsigned shorts returned.
- *
- * Return value:
- *   none.
- *
- * Note:
- *   This is the same as PMATH internal procedure rand48_24to16 except
- *   that the order of the entries in x16 have been reversed and its
- *   type declaration has been changed to u16.  It is intended for use
- *   by user-callable routines that get seeds or multipliers for users.
- */
-void PM_24to16(double x24[2], u16 x16[3])
-{
-   double t0u, t0l, t1u, t1l;
-
-   t0u = floor(x24[0]/TWO16);
-   t0l = x24[0] - t0u*TWO16;
-   x16[0] = (u16)t0l;
-
-   t1u = floor(x24[1]/TWO8);
-   x16[2] = (u16)t1u;
-   t1l = x24[1] - t1u*TWO8;
-   x16[1] = (u16)(t0u + t1l*TWO8);
-
-#ifdef RAN_DEBUG
-fprintf(stderr,"PM_24to16: x24 = %.1f %.1f\n",x24[1],x24[0]);
-fprintf(stderr,"PM_24to16: x16 = %04x %04x %04x\n",x16[2],x16[1],x16[0]);
-#endif
-   return;
-}
-
-/*   End of internal procedure definitions.
-   -----------------------------------------
-*/
-
-
-
-/* --------------------------------------
-      Define service routines.
-   --------------------------------------
-*/
-
-/* PM_GSeed -- Get the seed as 2 doubles each containing 24 bits of data.
- * 
- * Calling sequence:
- *   double myseed[2];
- *   PM_GSeed(myseed);
- *
- *   myseed :=  the 2 doubles of 24 bits returned.
- *
- * Return value:
- *   none.
- */
-void PM_GSeed(double seedout[2])
-{
-   seedout[0] = dseed[0];
-   seedout[1] = dseed[1];
-
-#ifdef RAN_DEBUG
-   fprintf(stderr,"PM_GSeed: seed = %.1f %.1f\n", dseed[1], dseed[0]);
-#endif
-   return;
-}
-
-/* PM_SSeed -- Set the seed from 2 doubles each containing 24 bits of data.
- * 
- * Calling sequence:
- *   double myseed[2];
- *   PM_SSeed(myseed);
- *
- *   myseed :  the 2 doubles of 24 bits to be used as the new seed.
- *
- * Return value:
- *   none.
- *
- * Caution:
- *   This routine does not check for valid seeds!  It is intended for
- *   use in a safe interface such as Setranf (in ranf.c).
- *   The elements of the myseed array are the coefficients in the base
- *   2**24 representation of the odd 48-bit integer
- *         seed =  myseed[1]*2**24 + myseed[0],
- *   so myseed must satisfy:
- *     (1) myseed[0] and myseed[1] are integers, stored in doubles;
- *     (2) 0 < myseed[0] < 2**24, 0 <= myseed[1] < 2**24;
- *     (3) myseed[0] is odd.
- */
-void PM_SSeed(double seedin[2])
-{
-   dseed[0] = seedin[0];
-   dseed[1] = seedin[1];
-
-#ifdef RAN_DEBUG
-   fprintf(stderr,"PM_SSeed: seed = %.1f %.1f\n", dseed[1], dseed[0]);
-#endif
-   return;
-}
-
-
-/* PM_GMult -- Get the multiplier as 2 doubles each containing 24 bits 
- *             of data.
- * 
- * Calling sequence:
- *   double mymult[2];
- *   PM_GMult(mymult);
- *
- *   mymult :=  the 2 doubles of 24 bits returned.
- *
- * Return value:
- *   none.
- */
-void PM_GMult(double multout[2])
-{
-   multout[0] = dmult[0];
-   multout[1] = dmult[1];
-
-#ifdef RAN_DEBUG
-   fprintf(stderr,"PM_GMult: mult = %.1f %.1f\n", dmult[1], dmult[0]);
-#endif
-   return;
-}
-
-/* PM_SMult -- Set the multiplier from 2 doubles each containing 24 bits 
- *             of data.
- * 
- * Calling sequence:
- *   double mymult[2];
- *   PM_SMult(mymult);
- *
- *   mymult :  the 2 doubles of 24 bits to be used as the new multiplier.
- *
- * Return value:
- *   none.
- *
- * Caution:
- *   This routine does not check for valid multipliers!  It is intended for
- *   use in a safe interface such as Setmult (in ranf.c).
- *   The elements of the mymult array are the coefficients in the base
- *   2**24 representation of the odd 48-bit integer
- *         mult =  mymult[1]*2**24 + mymult[0].
- *   Since we require  1 < mult < 2**46, mymult must satisfy:
- *     (1) mymult[0] and mymult[1] are integers, stored in doubles;
- *     (2) 1 < mymult[0] < 2**24, 0 <= mymult[1] < 2**22;
- *     (3) mymult[0] is odd.
- */
-void PM_SMult(double multin[2])
-{
-   dmult[0] = multin[0];
-   dmult[1] = multin[1];
-
-#ifdef RAN_DEBUG
-   fprintf(stderr,"PM_SMult: mult = %.1f %.1f\n", dmult[1], dmult[0]);
-#endif
-   return;
-}
-
-
-
-/* --------------------------------------
-      Define the generator itself.
-   --------------------------------------
-*/
-
-/*DECK PM_RANF
-C***BEGIN PROLOGUE  PM_RANF
-C***PURPOSE  Uniform random-number generator.
-C            The pseudorandom numbers generated by C function PM_RANF
-C            are uniformly distributed in the open interval (0,1).
-C***LIBRARY   PMATH
-C***CATEGORY  L6A21
-C***TYPE      REAL*8 (PM_RANF-8)
-C***KEYWORDS  RANDOM NUMBER GENERATION, UNIFORM DISTRIBUTION
-C***AUTHOR  Rathkopf, Jim, (LLNL/CP-Division)
-C           Fritsch, Fred N., (LLNL/LC/MSS)
-C***DESCRIPTION
-C     (Portable version of Cray MATHLIB routine RANF.)
-C     (Equivalent to Fortran-callable RANF8.)
-C *Usage:
-C        double r;
-C        r = PM_RANF();
-C
-C *Function Return Values:
-C     r        A random number between 0 and 1.
-C
-C *Description:
-C     PM_RANF generates pseudorandom numbers lying strictly between 0
-C     and 1. Each call to PM_RANF produces a different value, until the
-C     sequence cycles after 2**46 calls.
-C
-C     PM_RANF is a linear congruential pseudorandom-number generator.
-C     The default starting seed is
-C                SEED = 4510112377116321(oct) = 948253fc9cd1(hex).
-C     The multiplier is 1207264271730565(oct) = 2875a2e7b175(hex).
-C
-C *See Also:
-C     This is the same generator used by the Fortran generator family
-C     SRANF/DRANF/RANF8.
-C     The starting seed for PM_RANF may be set via PM_SSeed.
-C     The current PM_RANF seed may be obtained from PM_GSeed.
-C     The current PM_RANF multiplier may be obtained from PM_GMult.
-C     The PM_RANF multiplier may be set via PM_SMult (changing the
-C     multiplier is not recommended).
-C
-C***ROUTINES CALLED  (NONE)
-C***REVISION HISTORY  (YYMMDD)
-C   930308  DATE WRITTEN
-C           (Date from Biester's math_rnf.c.)
-C   931021  Changed name from B_RANF to PM_RANF8.  (FNF)
-C   940425  Added SLATEC-format prologue.  (FNF)
-C   950922  Re-ordered code to eliminate old EMULRAND routine. (FNF)
-C   951005  Eliminated unnecessary upper part computations.
-C           Added fmod invocations in t2_48 computation to avoid
-C           generating numbers too large for Cray floating point. (FNF)
-C   951020  Replaced fmod with macro MODF to avoid library calls. (FNF)
-C   951025  Removed MODF from first term of t2_48. (FNF)
-C   960827  Changed name from PM_RANF8 to PM_RANF, eliminated reference
-C           counter, and removed Fortran binding.  (FNF)
-C***END PROLOGUE  PM_RANF
-C
-C*Internal notes:
-C
-C    This routine generates pseudo-random numbers through a 48-bit
-C    linear congruential algorithm.  This emulates the drand48 library
-C    of random number generators available on many, but not all,
-C    UNIX machines.
-C
-C Algorithm:
-C    x(n+1) = (a*x(n) + c)  mod m
-C
-C   where the defaults for the standard UNIX rand48 are:
-C
-C                   double name    hexdecimal          decimal
-C    x: seed       -dseed[0],[1]  1234abcd330e      20017429951246
-C    a: multiplier -dmult[0],[1]     5deece66d         25214903917
-C    c: adder      -dadder                   b                  11
-C    m: 2**48                    1000000000000     281474976710656
-C
-C     24-bit defaults (decimal) (lower bits listed first)
-C    x: dseed[0],[1] = 13447950.0, 1193131.0
-C    a: dmult[0],[1] = 15525485.0, 1502.0
-C    c: dadder       = 11.0
-C
-C   The Cray defaults used in this code are:
-C
-C                   double name    hexdecimal           octal
-C    x: seed       -dseed[0],[1]  948253fc9cd1    4510112377116321
-C    a: multiplier -dmult[0],[1]  2875a2e7b175    1207264271730565
-C    c: adder      -dadder                   0                   0
-C    m: 2**48                    1000000000000   10000000000000000
-C
-C     24-bit defaults (decimal) (lower bits listed first)
-C    x: dseed[0],[1] = 16555217.0, 9732691.0
-C    a: dmult[0],[1] = 15184245.0, 2651554.0
-C    c: dadder       = 0.0
-C
-C Return value:
-C   double random number such that 0.0< d <1.0
-C
-C**End
- */
-double PM_RANF (void)
-{
-   double t1_48, t2_48, t1_24[2], t2_24;
-   double d;
-
-   /* perform 48-bit arithmetic using two part data */
-   t1_48 = dmult[0]*dseed[0] + dadder;
-   t2_48 = dmult[1]*dseed[0] + MODF(dmult[0]*dseed[1],BASE24);
-       /* First term safe, since dmult[1] < 2**22 for default mult */
-
-   t1_24[1] = floor(t1_48/BASE24);       /* upper part */
-   t1_24[0] = t1_48 - t1_24[1]*BASE24;   /* lower part */
-
-   t2_24 = MODF(t2_48, BASE24);      /* lower part */
-
-   t2_48 = t2_24 + t1_24[1];
-
-   t2_24 = MODF(t2_48, BASE24);      /* discard anything over 2**24 */
-
-   d = (dseed[0] + dseed[1]*BASE24)/TWO48;
-
-   dseed[0] = t1_24[0];
-   dseed[1] = t2_24;
-
-   return (d);
-}
-
-
diff --git a/scipy/corelib/random_lite/ranf.c b/scipy/corelib/random_lite/ranf.c
deleted file mode 100644
index 5ca7dc505..000000000
--- a/scipy/corelib/random_lite/ranf.c
+++ /dev/null
@@ -1,374 +0,0 @@
-/* ranf.c
-
-  The routines below implement an interface based on the drand48 family of
-  random number routines, but with seed handling and default starting seed
-  compatible with the CRI MATHLIB random number generator (RNG) family "ranf".
-*/
-
-/* -----------------------------------------------------------------------
- *
- * User-callable routines defined:
- *    Seedranf - Set seed from 32-bit integer
- *    Mixranf  - Set seed, with options; return seed
- *    Getranf  - Get 48-bit seed in integer array
- *    Setranf  - Set seed from 48-bit integer
- *    Getmult  - Get 48-bit multiplier in integer array
- *    Setmult  - Set multiplier from 48-bit integer
- *    Ranf     - The generator itself
- *
- * Currently responsible person:
- *    Fred N. Fritsch
- *    Computer Applications Organization
- *    LCPD, ICF Group
- *    Lawrence Livermore National Laboratory
- *    fnf@llnl.gov
- *
- * Revision history:
- *   (yymmdd)
- *    91????  DATE WRITTEN
- *            This started with the ranf.c that was checked out of the
- *            Basis repository in August 1996.  It was written by Lee
- *            Busby to provide an interface for Basis to the drand48
- *            family of generators that is compatible with the CRI
- *            MATHLIB ranf.  Following are the relevant cvs log entries.
- *    950511  Added back seedranf, setranf, getranf, mixranf from older
- *            Basis version. (LB)
- *    950530  Changed type of u32 from long to int, which is currently ok
- *            on all our supported platforms. (LB)
- *    960823  Revised to use the portable PMATH generator instead. (FNF)
- *    960903  Added new routines Getmult and Setmult. (FNF)
- *    960904  Changed declarations for 48-bit quantities from pointers to
- *            two-element arrays, to be consistent with actual usage, and
- *            brought internal documentation up to LDOC standard. (FNF)
- *    960905  Moved definitions of internal procedures PM_16to24 and
- *            PM_24to16 to pmath_rng.c (see ranf.h).  (FNF)
- *    960911  Corrected some problems with return values from Mixranf.
- *            Also improved calling sequence descriptions for Seedranf and
- *            Mixranf. (FNF)
- *    960916  Eliminated state variable ranf_started.  Since the PMATH
- *            family has the Cray default starting values built in, this
- *            is no longer needed. (FNF)
- *    961011  Modifed Setranf and Setmult to work OK on Cray's when given
- *            values saved on a 32-bit workstation. (FNF)
- *    961212  Corrected dimension error in Getmult.  (FNF per E. Brooks)
- *    
- * -----------------------------------------------------------------------
- */
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-   /* Get includes, typedefs, and protypes for ranf.c and pmath_rng.c */
-#include "ranf.h"
-
-
-/* Seedranf - Reset the RNG, given a new (32-bit) seed value.
- *
- * Calling sequence:
- *   u32 s;
- *   Seedranf(&s);
- *
- *   s :  (pointer to) the new (32-bit) seed value.
- *
- * Return value:
- *   none.
- *
- * Note:
- *   The upper 16 bits of the 48-bit seed are set to zero and then
- *   Setranf is called with the result.  This is provided primarily
- *   as a user convenience.
- */
-void Seedranf(u32 *s)
-{
-    u32 s48[2];
-
-    s48[0] = *s;
-    s48[1] = 0x0;
-#ifdef RAN_DEBUG
-fprintf(stderr,"Seedranf set 48-bit seed  %08x %08x\n",s48[1],s48[0]);
-#endif
-    (void)Setranf(s48);
-}
-
-
-/* Mixranf - Reset the RNG with a new 48-bit seed, with options.
- *
- * Calling sequence:
- *   int s;
- *   u32 s48[2];
- *   Mixranf(&s,s48);
- *
- *   s :     the value used to determine the new seed, as follows:
- *             s < 0  ==>  Use the default initial seed value.
- *             s = 0  ==>  Set a "random" value for the seed from the
- *                         system clock.
- *             s > 0  ==>  Set seed directly (32 bits only).
- *   s48 :=  the new 48-bit seed, returned in an array of u32's.
- *           (The upper 16 bits of s48[1] will be zero.)
- *
- * Return value:
- *   none.
- *
- * Note:
- *   In case of a positive s, all of s48[1] will be zero.  In order to set
- *   a seed with more than 32 bits, use Setranf directly.
- *
- * Portability:
- *   The structs timeval and timezone and routine gettimeofday are required
- *   when s=0.  These are generally available in <sys/time.h>.
- */
-void Mixranf(int *s,u32 s48[2])
-{
-    if(*s < 0){ /* Set default initial value */
-	s48[0] = s48[1] = 0;
-        Setranf(s48);
-        Getranf(s48);  /* Return default seed */
-
-    }else if(*s == 0){ /* Use the clock */
-	int i;
-#ifdef __MWERKS__
-/* use this for the Macintosh */
-	time_t theTime;
-	UnsignedWide	tv_usec;
-
-	(void)time(&theTime);
-	Microseconds(&tv_usec); /* the microsecs since start up */
-	s48[0] = (u32)theTime;
-	s48[1] = (u32)tv_usec.lo;
-#else
-#if defined(_WIN32)
-    // suggestions welcome for something better!
-    // one of these is from start of job, the other time of day
-	time_t long_time;
-	clock_t clock_time;
-	time(&long_time);
-	clock_time = clock();
-	s48[0] = (u32)long_time;
-	s48[1] = (u32)clock_time;
-#else
-	struct timeval tv;
-	struct timezone tz;
-#if !defined(__sgi)
-	int gettimeofday(struct timeval *, struct timezone *);
-#endif
-
-	(void)gettimeofday(&tv,&tz);
-	s48[0] = (u32)tv.tv_sec;
-	s48[1] = (u32)tv.tv_usec;
-#endif /* !_WIN32 */
-#endif /* !__MWERKS__ */
-	Setranf(s48);
-	for(i=0;i<10;i++) (void)Ranf();  /* Discard first 10 numbers */
-        Getranf(s48);  /* Return seed after these 10 calls */
-    }else{ /* s > 0, so set seed directly */
-        s48[0] = (u32)*s;
-        s48[1] = 0x0;
-        Setranf(s48);
-        Getranf(s48);  /* Return (possibly) corrected value of seed */
-    }
-#ifdef RAN_DEBUG
-fprintf(stderr,"Mixranf set 48-bit seed  %08x %08x\n",s48[1],s48[0]);
-#endif
-}
-
- 
-/* Getranf - Get the 48-bit seed value.
- *
- * Calling sequence:
- *   u32 s48[2];
- *   Getranf(s48);
- *
- *   s48 :=  the current 48-bit seed, stored in an array of u32's.
- *           (The upper 16 bits of s48[1] will be zero.)
- *
- * Return value:
- *   none.
- */
-void Getranf(u32 s48[2])
-{
-    u16 p[3];
-    double pm_seed[2];
-
-    /* Get the PMATH seed in an array of doubles. */
-    PM_GSeed(pm_seed);
-
-    /* Convert it to an array of u16's */
-    PM_24to16(pm_seed, p);
-
-/* Now put all the bits in the right place.  The picture below shows
-   s48[0], the rightmost 16 bits of s48[1], and all 16 bits of p[0],
-   p[1], and p[2] as aligned after assignment.
-
-  ---------------------------------------------------------------------
-  |        s48[1]        |                    s48[0]                  |
-  ---------------------------------------------------------------------
-  15                    0|31                  16|15                   0
-  ---------------------------------------------------------------------
-  |         p[2]         |         p[1]         |        p[0]         |
-  ---------------------------------------------------------------------
-
-*/
-
-    s48[0] = p[1];
-    s48[0] = (s48[0] << 16) + p[0];
-    s48[1] = p[2];
-#ifdef RAN_DEBUG
-fprintf(stderr,"Getranf returns %08x %08x\n",s48[1],s48[0]);
-#endif
-}
-
-
-/* Setranf - Reset the RNG seed from a 48-bit integer value.
- *
- * Calling sequence:
- *   u32 s48[2];
- *   Setranf(s48);
- *
- *   s48 :  the new 48-bit seed, contained in an array of u32's.
- *          If s48 is all zeros, set the default starting value,
- *          948253fc9cd1 (hex).
- *
- * Return value:
- *   none.
- *
- * Note:
- *   A 1 is masked into the lowest bit position to make sure the seed
- *   is odd.  (The upper 16 bits of s48[1] will be ignored if nonzero.)
- */
-void Setranf(u32 s48[2])
-{
-    u16 p[3];
-    double pm_seed[2];
-
-#ifdef RAN_DEBUG
-fprintf(stderr,"Setranf called with s48 = %08x %08x\n",s48[1],s48[0]);
-#endif
-
-    if(s48[0] == 0 && s48[1] == 0){ /* Set default starting value */
-	s48[0] = 0x53fc9cd1;
-	s48[1] = 0x9482;
-    }
-
-    /* Store seed as 3 u16's. */
-    p[0] = (s48[0] & 0xffff) | 0x1; /* Force an odd number */
-    p[1] = (s48[0] >> 16) & 0xffff; /* Protect against sign extension on Cray*/
-    p[2] = s48[1];
-
-    /* Convert to PMATH seed. */
-    PM_16to24(p, pm_seed);
-
-    /* Now reset the seed. */
-    PM_SSeed(pm_seed);
-
-#ifdef RAN_DEBUG
-fprintf(stderr,"Leaving Setranf\n");
-#endif
-}
-
- 
-/* Getmult - Get the current 48-bit multiplier.
- *
- * Calling sequence:
- *   u32 m48[2];
- *   Getmult(m48);
- *
- *   m48 :=  the current 48-bit multiplier, stored in an array of u32's.
- *           (The upper 16 bits of m48[1] will be zero.)
- *
- * Return value:
- *   none.
- */
-void Getmult(u32 m48[2])
-{
-    u16 p[3];
-    double pm_mult[2];
-
-    /* Get the PMATH multiplier in an array of doubles. */
-    PM_GMult(pm_mult);
-
-    /* Convert it to an array of u16's */
-    PM_24to16(pm_mult, p);
-
-/* Now put all the bits in the right place.  (See Getranf for picture.) */
-
-    m48[0] = p[1];
-    m48[0] = (m48[0] << 16) + p[0];
-    m48[1] = p[2];
-#ifdef RAN_DEBUG
-fprintf(stderr,"Getmult returns %08x %08x\n",m48[1],m48[0]);
-#endif
-}
-
-
-/* Setmult - Reset the RNG multiplier from a 48-bit integer value.
- *
- * Calling sequence:
- *   u32 m48[2];
- *   Setmult(m48);
- *
- *   m48 :  the new 48-bit multiplier, contained in an array of u32's.
- *          If m48 is all zeros, set the default starting value,
- *          2875a2e7b175 (hex).
- *
- * Return value:
- *   none.
- *
- * Note:
- *   A 1 is masked into the lowest bit position to make sure the value
- *   is odd.  (The upper 18 bits of m48[1] will be ignored if nonzero.)
- */
-void Setmult(u32 *m48)
-{
-    u16 p[3];
-    double pm_mult[2];
-
-#ifdef RAN_DEBUG
-fprintf(stderr,"Setmult called with m48 = %08x %08x\n",m48[1],m48[0]);
-#endif
-    if(m48[0] == 0 && m48[1] == 0){ /* Set default starting value */
-	m48[0] = 0xa2e7b175;
-	m48[1] = 0x2875;
-    }
-
-    /* Store multiplier as 3 u16's. */
-    p[0] = (m48[0] & 0xffff) | 0x1; /* Force an odd number */
-    p[1] = (m48[0] >> 16) & 0xffff; /* Protect against sign extension on Cray*/
-    p[2] = m48[1] & 0x3fff;       /* Force mult < 2**26 */
-
-    /* Convert to PMATH multiplier. */
-    PM_16to24(p, pm_mult);
-
-    /* Now reset the multiplier. */
-    PM_SMult(pm_mult);
-
-#ifdef RAN_DEBUG
-fprintf(stderr,"Leaving Setmult\n");
-#endif
-}
-
-
-/* Ranf - Uniform random number generator.
- *
- * Calling sequence:
- *   f64 value;
- *   value = Ranf();
- *
- * Return value:
- *   value :=  the next number in the ranf sequence.
- *
- * Note:
- *   The Ranf function is just a shell around PM_RANF.  Since the
- *   PMATH generator has the Cray default starting values built in,
- *   no initialization is needed.
- */
-f64 Ranf(void)
-{
-    return(PM_RANF());
-}
-
-
-#ifdef __cplusplus
-}
-#endif
-
diff --git a/scipy/corelib/random_lite/ranf.h b/scipy/corelib/random_lite/ranf.h
deleted file mode 100644
index 0c10d926d..000000000
--- a/scipy/corelib/random_lite/ranf.h
+++ /dev/null
@@ -1,41 +0,0 @@
-#include <stdlib.h>
-#include <math.h>
-#ifdef RAN_DEBUG
-#include <stdio.h>
-#endif
-
-#ifdef __MWERKS__
-/*#include <utime.h>*/
-#include <time.h>
-#include <Timer.h>
-#else
-#if defined(_WIN32)
-#include <time.h>
-#else
-#include <sys/time.h>
-#endif
-#endif
-
-typedef unsigned int u32;
-typedef unsigned short int u16;
-typedef double f64;
-
-/*  Prototypes for routines defined in ranf.c  */
-void Seedranf(u32 *s);            /* Set seed from 32-bit integer */
-void Mixranf(int *s, u32 s48[2]); /* Set seed, with options; return seed */
-void Getranf(u32 s48[2]);         /* Get 48-bit seed in integer array */
-void Setranf(u32 s48[2]);         /* Set seed from 48-bit integer */
-void Getmult(u32 m48[2]);         /* Get 48-bit multiplier in integer array */
-void Setmult(u32 m48[2]);         /* Set multiplier from 48-bit integer */
-f64  Ranf(void);                  /* The generator itself */
-
-/*  Prototypes for routines defined in pmath_rng.c  */
-void PM_16to24(u16 x16[3], double x24[2]);
-                              /* Convert 3 16-bit shorts to 2 24-bit doubles */
-void PM_24to16(double x24[2], u16 x16[3]);
-                              /* Convert 2 24-bit doubles to 3 16-bit shorts */
-void PM_GSeed(double seedout[2]);  /* Get the current seed */
-void PM_SSeed(double seedin[2]);   /* Reset the seed (unsafe) */
-void PM_GMult(double multout[2]);  /* Get the current multiplier */
-void PM_SMult(double multin[2]);   /* Reset the multiplier (unsafe) */
-f64  PM_RANF(void);                /* The generator itself */
diff --git a/scipy/corelib/random_lite/ranlib.c b/scipy/corelib/random_lite/ranlib.c
deleted file mode 100644
index 8e4f0ee00..000000000
--- a/scipy/corelib/random_lite/ranlib.c
+++ /dev/null
@@ -1,1941 +0,0 @@
-#include "Python.h"  /* for throwing exceptions */
-#include "ranlib.h"
-#include <stdio.h>
-#include <math.h>
-#include <stdlib.h>
-#define ABS(x) ((x) >= 0 ? (x) : -(x))
-#ifndef min
-#define min(a,b) ((a) <= (b) ? (a) : (b))
-#endif
-#ifndef max
-#define max(a,b) ((a) >= (b) ? (a) : (b))
-#endif
-void ftnstop(char*);
-float genbet(float aa,float bb)
-
-#ifdef _MSC_VER		/** disable extremly common MSVC warnings **/
-#pragma warning(disable:4305)
-#pragma warning(disable:4244)
-#endif
-
-
-/*
-**********************************************************************
-     float genbet(float aa,float bb)
-               GeNerate BETa random deviate
-                              Function
-     Returns a single random deviate from the beta distribution with
-     parameters A and B.  The density of the beta is
-               x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
-                              Arguments
-     aa --> First parameter of the beta distribution
-       
-     bb --> Second parameter of the beta distribution
-       
-                              Method
-     R. C. H. Cheng
-     Generating Beta Variatew with Nonintegral Shape Parameters
-     Communications of the ACM, 21:317-322  (1978)
-     (Algorithms BB and BC)
-**********************************************************************
-*/
-{
-#define expmax 89.0
-#define infnty 1.0E38
-static float olda = -1.0;
-static float oldb = -1.0;
-static float genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
-static long qsame;
-
-    qsame = olda == aa && oldb == bb;
-    if(qsame) goto S20;
-    if ((aa <= 0.0) || (bb <= 0.0)) {
-        char as[50], bs[50];
-        snprintf(as, 50, "%16.6E", aa);
-        snprintf(bs, 50, "%16.6E", bb);
-        PyErr_Format(PyExc_ValueError, "AA (%s) or BB (%s) <= 0 in GENBET",
-                     as, bs);
-        return 0.0;
-    }
-    olda = aa;
-    oldb = bb;
-S20:
-    if(!(min(aa,bb) > 1.0)) goto S100;
-/*
-     Alborithm BB
-     Initialize
-*/
-    if(qsame) goto S30;
-    a = min(aa,bb);
-    b = max(aa,bb);
-    alpha = a+b;
-    beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
-    gamma = a+1.0/beta;
-S30:
-S40:
-    u1 = ranf();
-/*
-     Step 1
-*/
-    u2 = ranf();
-    v = beta*log(u1/(1.0-u1));
-    if(!(v > expmax)) goto S50;
-    w = infnty;
-    goto S60;
-S50:
-    w = a*exp(v);
-S60:
-    z = pow(u1,2.0)*u2;
-    r = gamma*v-1.3862944;
-    s = a+r-w;
-/*
-     Step 2
-*/
-    if(s+2.609438 >= 5.0*z) goto S70;
-/*
-     Step 3
-*/
-    t = log(z);
-    if(s > t) goto S70;
-/*
-     Step 4
-*/
-    if(r+alpha*log(alpha/(b+w)) < t) goto S40;
-S70:
-/*
-     Step 5
-*/
-    if(!(aa == a)) goto S80;
-    genbet = w/(b+w);
-    goto S90;
-S80:
-    genbet = b/(b+w);
-S90:
-    goto S230;
-S100:
-/*
-     Algorithm BC
-     Initialize
-*/
-    if(qsame) goto S110;
-    a = max(aa,bb);
-    b = min(aa,bb);
-    alpha = a+b;
-    beta = 1.0/b;
-    delta = 1.0+a-b;
-    k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);
-    k2 = 0.25+(0.5+0.25/delta)*b;
-S110:
-S120:
-    u1 = ranf();
-/*
-     Step 1
-*/
-    u2 = ranf();
-    if(u1 >= 0.5) goto S130;
-/*
-     Step 2
-*/
-    y = u1*u2;
-    z = u1*y;
-    if(0.25*u2+z-y >= k1) goto S120;
-    goto S170;
-S130:
-/*
-     Step 3
-*/
-    z = pow(u1,2.0)*u2;
-    if(!(z <= 0.25)) goto S160;
-    v = beta*log(u1/(1.0-u1));
-    if(!(v > expmax)) goto S140;
-    w = infnty;
-    goto S150;
-S140:
-    w = a*exp(v);
-S150:
-    goto S200;
-S160:
-    if(z >= k2) goto S120;
-S170:
-/*
-     Step 4
-     Step 5
-*/
-    v = beta*log(u1/(1.0-u1));
-    if(!(v > expmax)) goto S180;
-    w = infnty;
-    goto S190;
-S180:
-    w = a*exp(v);
-S190:
-    if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120;
-S200:
-/*
-     Step 6
-*/
-    if(!(a == aa)) goto S210;
-    genbet = w/(b+w);
-    goto S220;
-S210:
-    genbet = b/(b+w);
-S230:
-S220:
-    return genbet;
-#undef expmax
-#undef infnty
-}
-float genchi(float df)
-/*
-**********************************************************************
-     float genchi(float df)
-                Generate random value of CHIsquare variable
-                              Function
-     Generates random deviate from the distribution of a chisquare
-     with DF degrees of freedom random variable.
-                              Arguments
-     df --> Degrees of freedom of the chisquare
-            (Must be positive)
-       
-                              Method
-     Uses relation between chisquare and gamma.
-**********************************************************************
-*/
-{
-static float genchi;
-
-    if (df <= 0.0) {
-        char dfs[50];
-        snprintf(dfs, 50, "%16.6E", df);
-        PyErr_Format(PyExc_ValueError, "DF (%s) <= 0 in GENCHI", dfs);
-        return 0.0;
-    }
-    genchi = 2.0*gengam(1.0,df/2.0);
-    return genchi;
-}
-float genexp(float av)
-/*
-**********************************************************************
-     float genexp(float av)
-                    GENerate EXPonential random deviate
-                              Function
-     Generates a single random deviate from an exponential
-     distribution with mean AV.
-                              Arguments
-     av --> The mean of the exponential distribution from which
-            a random deviate is to be generated.
-                              Method
-     Renames SEXPO from TOMS as slightly modified by BWB to use RANF
-     instead of SUNIF.
-     For details see:
-               Ahrens, J.H. and Dieter, U.
-               Computer Methods for Sampling From the
-               Exponential and Normal Distributions.
-               Comm. ACM, 15,10 (Oct. 1972), 873 - 882.
-**********************************************************************
-*/
-{
-static float genexp;
-
-    genexp = sexpo()*av;
-    return genexp;
-}
-float genf(float dfn,float dfd)
-/*
-**********************************************************************
-     float genf(float dfn,float dfd)
-                GENerate random deviate from the F distribution
-                              Function
-     Generates a random deviate from the F (variance ratio)
-     distribution with DFN degrees of freedom in the numerator
-     and DFD degrees of freedom in the denominator.
-                              Arguments
-     dfn --> Numerator degrees of freedom
-             (Must be positive)
-     dfd --> Denominator degrees of freedom
-             (Must be positive)
-                              Method
-     Directly generates ratio of chisquare variates
-**********************************************************************
-*/
-{
-static float genf,xden,xnum;
-
-    if((dfn <= 0.0 || dfd <= 0.0)) {
-        char dfns[50], dfds[50];
-        snprintf(dfns, 50, "%16.6E", dfn);
-        snprintf(dfds, 50, "%16.6E", dfd);
-        PyErr_Format(PyExc_ValueError,
-                     "Degrees of freedom nonpositive in GENF: DFN=%s DFD=%s",
-                     dfns,dfds);
-        return 0.0;
-    }
-    xnum = genchi(dfn)/dfn;
-/*
-      GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD )
-*/
-    xden = genchi(dfd)/dfd;
-    if(!(xden <= 9.999999999998E-39*xnum)) goto S20;
-    /*
-    fputs(" GENF - generated numbers would cause overflow",stderr);
-    fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden);
-    fputs(" GENF returning 1.0E38",stderr);
-    */
-    genf = 1.0E38;
-    goto S30;
-S20:
-    genf = xnum/xden;
-S30:
-    return genf;
-}
-float gengam(float a,float r)
-/*
-**********************************************************************
-     float gengam(float a,float r)
-           GENerates random deviates from GAMma distribution
-                              Function
-     Generates random deviates from the gamma distribution whose
-     density is
-          (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X)
-                              Arguments
-     a --> Location parameter of Gamma distribution
-     r --> Shape parameter of Gamma distribution
-                              Method
-     Renames SGAMMA from TOMS as slightly modified by BWB to use RANF
-     instead of SUNIF.
-     For details see:
-               (Case R >= 1.0)
-               Ahrens, J.H. and Dieter, U.
-               Generating Gamma Variates by a
-               Modified Rejection Technique.
-               Comm. ACM, 25,1 (Jan. 1982), 47 - 54.
-     Algorithm GD
-               (Case 0.0 <= R <= 1.0)
-               Ahrens, J.H. and Dieter, U.
-               Computer Methods for Sampling from Gamma,
-               Beta, Poisson and Binomial Distributions.
-               Computing, 12 (1974), 223-246/
-     Adapted algorithm GS.
-**********************************************************************
-*/
-{
-static float gengam;
-
-    gengam = sgamma(r);
-    gengam /= a;
-    return gengam;
-}
-void genmn(float *parm,float *x,float *work)
-/*
-**********************************************************************
-     void genmn(float *parm,float *x,float *work)
-              GENerate Multivariate Normal random deviate
-                              Arguments
-     parm --> Parameters needed to generate multivariate normal
-               deviates (MEANV and Cholesky decomposition of
-               COVM). Set by a previous call to SETGMN.
-               1 : 1                - size of deviate, P
-               2 : P + 1            - mean vector
-               P+2 : P*(P+3)/2 + 1  - upper half of cholesky
-                                       decomposition of cov matrix
-     x    <-- Vector deviate generated.
-     work <--> Scratch array
-                              Method
-     1) Generate P independent standard normal deviates - Ei ~ N(0,1)
-     2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM
-     3) trans(A)E + MEANV ~ N(MEANV,COVM)
-**********************************************************************
-*/
-{
-static long i,icount,j,p,D1,D2,D3,D4;
-static float ae;
-
-    p = (long) (*parm);
-/*
-     Generate P independent normal deviates - WORK ~ N(0,1)
-*/
-    for(i=1; i<=p; i++) *(work+i-1) = snorm();
-    for(i=1,D3=1,D4=(p-i+D3)/D3; D4>0; D4--,i+=D3) {
-/*
-     PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky
-      decomposition of the desired covariance matrix.
-          trans(A)(1,1) = PARM(P+2)
-          trans(A)(2,1) = PARM(P+3)
-          trans(A)(2,2) = PARM(P+2+P)
-          trans(A)(3,1) = PARM(P+4)
-          trans(A)(3,2) = PARM(P+3+P)
-          trans(A)(3,3) = PARM(P+2-1+2P)  ...
-     trans(A)*WORK + MEANV ~ N(MEANV,COVM)
-*/
-        icount = 0;
-        ae = 0.0;
-        for(j=1,D1=1,D2=(i-j+D1)/D1; D2>0; D2--,j+=D1) {
-            icount += (j-1);
-            ae += (*(parm+i+(j-1)*p-icount+p)**(work+j-1));
-        }
-        *(x+i-1) = ae+*(parm+i);
-    }
-}
-void genmul(long n,float *p,long ncat,long *ix)
-/*
-**********************************************************************
- 
-     void genmul(int n,float *p,int ncat,int *ix)
-     GENerate an observation from the MULtinomial distribution
-                              Arguments
-     N --> Number of events that will be classified into one of
-           the categories 1..NCAT
-     P --> Vector of probabilities.  P(i) is the probability that
-           an event will be classified into category i.  Thus, P(i)
-           must be [0,1]. Only the first NCAT-1 P(i) must be defined
-           since P(NCAT) is 1.0 minus the sum of the first
-           NCAT-1 P(i).
-     NCAT --> Number of categories.  Length of P and IX.
-     IX <-- Observation from multinomial distribution.  All IX(i)
-            will be nonnegative and their sum will be N.
-                              Method
-     Algorithm from page 559 of
- 
-     Devroye, Luc
- 
-     Non-Uniform Random Variate Generation.  Springer-Verlag,
-     New York, 1986.
- 
-**********************************************************************
-*/
-{
-static float prob,ptot,sum;
-static long i,icat,ntot;
-    if(n < 0) ftnstop("N < 0 in GENMUL");
-    if(ncat <= 1) ftnstop("NCAT <= 1 in GENMUL");
-    ptot = 0.0F;
-    for(i=0; i<ncat-1; i++) {
-        if(*(p+i) < 0.0F) ftnstop("Some P(i) < 0 in GENMUL");
-        if(*(p+i) > 1.0F) ftnstop("Some P(i) > 1 in GENMUL");
-        ptot += *(p+i);
-    }
-    if(ptot > 0.99999F) ftnstop("Sum of P(i) > 1 in GENMUL");
-/*
-     Initialize variables
-*/
-    ntot = n;
-    sum = 1.0F;
-    for(i=0; i<ncat; i++) ix[i] = 0;
-/*
-     Generate the observation
-*/
-    for(icat=0; icat<ncat-1; icat++) {
-        prob = *(p+icat)/sum;
-        *(ix+icat) = ignbin(ntot,prob);
-        ntot -= *(ix+icat);
-	if(ntot <= 0) return;
-        sum -= *(p+icat);
-    }
-    *(ix+ncat-1) = ntot;
-/*
-     Finished
-*/
-    return;
-}
-float gennch(float df,float xnonc)
-/*
-**********************************************************************
-     float gennch(float df,float xnonc)
-           Generate random value of Noncentral CHIsquare variable
-                              Function
-     Generates random deviate  from the  distribution  of a  noncentral
-     chisquare with DF degrees  of freedom and noncentrality  parameter
-     xnonc.
-                              Arguments
-     df --> Degrees of freedom of the chisquare
-            (Must be > 1.0)
-     xnonc --> Noncentrality parameter of the chisquare
-               (Must be >= 0.0)
-                              Method
-     Uses fact that  noncentral chisquare  is  the  sum of a  chisquare
-     deviate with DF-1  degrees of freedom plus the  square of a normal
-     deviate with mean XNONC and standard deviation 1.
-**********************************************************************
-*/
-{
-static float gennch;
-
-    if((df <= 1.0 || xnonc < 0.0)) {
-        char dfs[50], xnoncs[50];
-        snprintf(dfs, 50, "%16.6E", df);
-        snprintf(xnoncs, 50, "%16.6E", xnonc);
-        PyErr_Format(PyExc_ValueError,
-                     "DF (%s) <= 1 or XNONC (%s) < 0 in GENNCH", dfs, xnoncs);
-        return 0.0;
-    }
-    gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0);
-    return gennch;
-}
-float gennf(float dfn,float dfd,float xnonc)
-/*
-**********************************************************************
-     float gennf(float dfn,float dfd,float xnonc)
-           GENerate random deviate from the Noncentral F distribution
-                              Function
-     Generates a random deviate from the  noncentral F (variance ratio)
-     distribution with DFN degrees of freedom in the numerator, and DFD
-     degrees of freedom in the denominator, and noncentrality parameter
-     XNONC.
-                              Arguments
-     dfn --> Numerator degrees of freedom
-             (Must be >= 1.0)
-     dfd --> Denominator degrees of freedom
-             (Must be positive)
-     xnonc --> Noncentrality parameter
-               (Must be nonnegative)
-                              Method
-     Directly generates ratio of noncentral numerator chisquare variate
-     to central denominator chisquare variate.
-**********************************************************************
-*/
-{
-static float gennf,xden,xnum;
-static long qcond;
-
-    qcond = dfn <= 1.0 || dfd <= 0.0 || xnonc < 0.0;
-    if (qcond) {
-        char dfns[50], dfds[50], xnoncs[50];
-        snprintf(dfns, 50, "%16.6E", dfn);
-        snprintf(dfds, 50, "%16.6E", dfd);
-        snprintf(xnoncs, 50, "%16.16E", xnonc);
-        PyErr_Format(PyExc_ValueError,
-                     "either numerator (%s) <= 1.0 or denominator (%s) < 0.0 or noncentrality parameter (%s) < 0.0",
-                     dfns, dfds, xnoncs);
-        return 0.0;
-    }
-    xnum = gennch(dfn,xnonc)/dfn;
-/*
-      GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD )
-*/
-    xden = genchi(dfd)/dfd;
-    if(!(xden <= 9.999999999998E-39*xnum)) goto S20;
-    /*
-    fputs(" GENNF - generated numbers would cause overflow",stderr);
-    fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden);
-    fputs(" GENNF returning 1.0E38",stderr);
-    */
-    gennf = 1.0E38;
-    goto S30;
-S20:
-    gennf = xnum/xden;
-S30:
-    return gennf;
-}
-float gennor(float av,float sd)
-/*
-**********************************************************************
-     float gennor(float av,float sd)
-         GENerate random deviate from a NORmal distribution
-                              Function
-     Generates a single random deviate from a normal distribution
-     with mean, AV, and standard deviation, SD.
-                              Arguments
-     av --> Mean of the normal distribution.
-     sd --> Standard deviation of the normal distribution.
-                              Method
-     Renames SNORM from TOMS as slightly modified by BWB to use RANF
-     instead of SUNIF.
-     For details see:
-               Ahrens, J.H. and Dieter, U.
-               Extensions of Forsythe's Method for Random
-               Sampling from the Normal Distribution.
-               Math. Comput., 27,124 (Oct. 1973), 927 - 937.
-**********************************************************************
-*/
-{
-static float gennor;
-
-    gennor = sd*snorm()+av;
-    return gennor;
-}
-void genprm(long *iarray,int larray)
-/*
-**********************************************************************
-    void genprm(long *iarray,int larray)
-               GENerate random PeRMutation of iarray
-                              Arguments
-     iarray <--> On output IARRAY is a random permutation of its
-                 value on input
-     larray <--> Length of IARRAY
-**********************************************************************
-*/
-{
-static long i,itmp,iwhich,D1,D2;
-
-    for(i=1,D1=1,D2=(larray-i+D1)/D1; D2>0; D2--,i+=D1) {
-        iwhich = ignuin(i,larray);
-        itmp = *(iarray+iwhich-1);
-        *(iarray+iwhich-1) = *(iarray+i-1);
-        *(iarray+i-1) = itmp;
-    }
-}
-float genunf(float low,float high)
-/*
-**********************************************************************
-     float genunf(float low,float high)
-               GeNerate Uniform Real between LOW and HIGH
-                              Function
-     Generates a real uniformly distributed between LOW and HIGH.
-                              Arguments
-     low --> Low bound (exclusive) on real value to be generated
-     high --> High bound (exclusive) on real value to be generated
-**********************************************************************
-*/
-{
-static float genunf;
-
-    if((low > high)) {
-        char lows[50], highs[50];
-        snprintf(lows, 50, "%16.6E", low);
-        snprintf(highs, 50, "%16.6E", high);
-        PyErr_Format(PyExc_ValueError, "LOW (%s) > HIGH (%s) in GENUNF",
-                     lows, highs);
-        return 0.0;
-    }
-    genunf = low+(high-low)*ranf();
-    return genunf;
-}
-void gscgn(long getset,long *g)
-/*
-**********************************************************************
-     void gscgn(long getset,long *g)
-                         Get/Set GeNerator
-     Gets or returns in G the number of the current generator
-                              Arguments
-     getset --> 0 Get
-                1 Set
-     g <-- Number of the current random number generator (1..32)
-**********************************************************************
-*/
-{
-#define numg 32L
-static long curntg = 1;
-    if(getset == 0) *g = curntg;
-    else  {
-        if(*g < 0 || *g > numg) {
-            PyErr_SetString(PyExc_ValueError,
-                            "Generator number out of range in GSCGN");
-            return;
-        }
-        curntg = *g;
-    }
-#undef numg
-}
-void gsrgs(long getset,long *qvalue)
-/*
-**********************************************************************
-     void gsrgs(long getset,long *qvalue)
-               Get/Set Random Generators Set
-     Gets or sets whether random generators set (initialized).
-     Initially (data statement) state is not set
-     If getset is 1 state is set to qvalue
-     If getset is 0 state returned in qvalue
-**********************************************************************
-*/
-{
-static long qinit = 0;
-
-    if(getset == 0) *qvalue = qinit;
-    else qinit = *qvalue;
-}
-void gssst(long getset,long *qset)
-/*
-**********************************************************************
-     void gssst(long getset,long *qset)
-          Get or Set whether Seed is Set
-     Initialize to Seed not Set
-     If getset is 1 sets state to Seed Set
-     If getset is 0 returns T in qset if Seed Set
-     Else returns F in qset
-**********************************************************************
-*/
-{
-static long qstate = 0;
-    if(getset != 0) qstate = 1;
-    else  *qset = qstate;
-}
-long ignbin(long n,float pp)
-/*
-**********************************************************************
-     long ignbin(long n,float pp)
-                    GENerate BINomial random deviate
-                              Function
-     Generates a single random deviate from a binomial
-     distribution whose number of trials is N and whose
-     probability of an event in each trial is P.
-                              Arguments
-     n  --> The number of trials in the binomial distribution
-            from which a random deviate is to be generated.
-     p  --> The probability of an event in each trial of the
-            binomial distribution from which a random deviate
-            is to be generated.
-     ignbin <-- A random deviate yielding the number of events
-                from N independent trials, each of which has
-                a probability of event P.
-                              Method
-     This is algorithm BTPE from:
-         Kachitvichyanukul, V. and Schmeiser, B. W.
-         Binomial Random Variate Generation.
-         Communications of the ACM, 31, 2
-         (February, 1988) 216.
-**********************************************************************
-     SUBROUTINE BTPEC(N,PP,ISEED,JX)
-     BINOMIAL RANDOM VARIATE GENERATOR
-     MEAN .LT. 30 -- INVERSE CDF
-       MEAN .GE. 30 -- ALGORITHM BTPE:  ACCEPTANCE-REJECTION VIA
-       FOUR REGION COMPOSITION.  THE FOUR REGIONS ARE A TRIANGLE
-       (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
-       THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
-     BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
-     BTPEC REFERS TO BTPE AND "COMBINED."  THUS BTPE IS THE
-       RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
-       USABLE ALGORITHM.
-     REFERENCE:  VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
-       "BINOMIAL RANDOM VARIATE GENERATION,"
-       COMMUNICATIONS OF THE ACM, FORTHCOMING
-     WRITTEN:  SEPTEMBER 1980.
-       LAST REVISED:  MAY 1985, JULY 1987
-     REQUIRED SUBPROGRAM:  RAND() -- A UNIFORM (0,1) RANDOM NUMBER
-                           GENERATOR
-     ARGUMENTS
-       N : NUMBER OF BERNOULLI TRIALS            (INPUT)
-       PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
-       ISEED:  RANDOM NUMBER SEED                (INPUT AND OUTPUT)
-       JX:  RANDOMLY GENERATED OBSERVATION       (OUTPUT)
-     VARIABLES
-       PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
-       NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
-       XNP:  VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
-       P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
-       FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
-       M:  INTEGER VALUE OF THE CURRENT MODE
-       FM:  FLOATING POINT VALUE OF THE CURRENT MODE
-       XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
-       P1:  AREA OF THE TRIANGLE
-       C:  HEIGHT OF THE PARALLELOGRAMS
-       XM:  CENTER OF THE TRIANGLE
-       XL:  LEFT END OF THE TRIANGLE
-       XR:  RIGHT END OF THE TRIANGLE
-       AL:  TEMPORARY VARIABLE
-       XLL:  RATE FOR THE LEFT EXPONENTIAL TAIL
-       XLR:  RATE FOR THE RIGHT EXPONENTIAL TAIL
-       P2:  AREA OF THE PARALLELOGRAMS
-       P3:  AREA OF THE LEFT EXPONENTIAL TAIL
-       P4:  AREA OF THE RIGHT EXPONENTIAL TAIL
-       U:  A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
-           FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
-           FROM THE REGION
-       V:  A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
-           (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
-           REJECT THE CANDIDATE VALUE
-       IX:  INTEGER CANDIDATE VALUE
-       X:  PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
-           AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
-       K:  ABSOLUTE VALUE OF (IX-M)
-       F:  THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
-           ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
-           ALSO USED IN THE INVERSE TRANSFORMATION
-       R: THE RATIO P/Q
-       G: CONSTANT USED IN CALCULATION OF PROBABILITY
-       MP:  MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
-            OF F WHEN IX IS GREATER THAN M
-       IX1:  CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
-             CALCULATION OF F WHEN IX IS LESS THAN M
-       I:  INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
-       AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
-       YNORM: LOGARITHM OF NORMAL BOUND
-       ALV:  NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
-       X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
-       USED IN THE FINAL ACCEPT/REJECT TEST
-       QN: PROBABILITY OF NO SUCCESS IN N TRIALS
-     REMARK
-       IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
-       SAVE A MEMORY POSITION AND A LINE OF CODE.  HOWEVER, SOME
-       COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
-       ARE NOT INVOLVED.
-     ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
-     GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
-     TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
-**********************************************************************
-*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
-*/
-{
-static float psave = -1.0;
-static long nsave = -1;
-static long ignbin,i,ix,ix1,k,m,mp,T1;
-static float al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,q,qn,r,u,v,w,w2,x,x1,
-    x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2;
-
-    if(pp != psave) goto S10;
-    if(n != nsave) goto S20;
-    if(xnp < 30.0) goto S150;
-    goto S30;
-S10:
-/*
-*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
-*/
-    psave = pp;
-    p = min(psave,1.0-psave);
-    q = 1.0-p;
-S20:
-    xnp = n*p;
-    nsave = n;
-    if(xnp < 30.0) goto S140;
-    ffm = xnp+p;
-    m = ffm;
-    fm = m;
-    xnpq = xnp*q;
-    p1 = (long) (2.195*sqrt(xnpq)-4.6*q)+0.5;
-    xm = fm+0.5;
-    xl = xm-p1;
-    xr = xm+p1;
-    c = 0.134+20.5/(15.3+fm);
-    al = (ffm-xl)/(ffm-xl*p);
-    xll = al*(1.0+0.5*al);
-    al = (xr-ffm)/(xr*q);
-    xlr = al*(1.0+0.5*al);
-    p2 = p1*(1.0+c+c);
-    p3 = p2+c/xll;
-    p4 = p3+c/xlr;
-S30:
-/*
-*****GENERATE VARIATE
-*/
-    u = ranf()*p4;
-    v = ranf();
-/*
-     TRIANGULAR REGION
-*/
-    if(u > p1) goto S40;
-    ix = xm-p1*v+u;
-    goto S170;
-S40:
-/*
-     PARALLELOGRAM REGION
-*/
-    if(u > p2) goto S50;
-    x = xl+(u-p1)/c;
-    v = v*c+1.0-ABS(xm-x)/p1;
-    if(v > 1.0 || v <= 0.0) goto S30;
-    ix = x;
-    goto S70;
-S50:
-/*
-     LEFT TAIL
-*/
-    if(u > p3) goto S60;
-    ix = xl+log(v)/xll;
-    if(ix < 0) goto S30;
-    v *= ((u-p2)*xll);
-    goto S70;
-S60:
-/*
-     RIGHT TAIL
-*/
-    ix = xr-log(v)/xlr;
-    if(ix > n) goto S30;
-    v *= ((u-p3)*xlr);
-S70:
-/*
-*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
-*/
-    k = ABS(ix-m);
-    if(k > 20 && k < xnpq/2-1) goto S130;
-/*
-     EXPLICIT EVALUATION
-*/
-    f = 1.0;
-    r = p/q;
-    g = (n+1)*r;
-    T1 = m-ix;
-    if(T1 < 0) goto S80;
-    else if(T1 == 0) goto S120;
-    else  goto S100;
-S80:
-    mp = m+1;
-    for(i=mp; i<=ix; i++) f *= (g/i-r);
-    goto S120;
-S100:
-    ix1 = ix+1;
-    for(i=ix1; i<=m; i++) f /= (g/i-r);
-S120:
-    if(v <= f) goto S170;
-    goto S30;
-S130:
-/*
-     SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
-*/
-    amaxp = k/xnpq*((k*(k/3.0+0.625)+0.1666666666666)/xnpq+0.5);
-    ynorm = -(k*k/(2.0*xnpq));
-    alv = log(v);
-    if(alv < ynorm-amaxp) goto S170;
-    if(alv > ynorm+amaxp) goto S30;
-/*
-     STIRLING'S FORMULA TO MACHINE ACCURACY FOR
-     THE FINAL ACCEPTANCE/REJECTION TEST
-*/
-    x1 = ix+1.0;
-    f1 = fm+1.0;
-    z = n+1.0-fm;
-    w = n-ix+1.0;
-    z2 = z*z;
-    x2 = x1*x1;
-    f2 = f1*f1;
-    w2 = w*w;
-    if(alv <= xm*log(f1/x1)+(n-m+0.5)*log(z/w)+(ix-m)*log(w*p/(x1*q))+(13860.0-
-      (462.0-(132.0-(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0+(13860.0-(462.0-
-      (132.0-(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0+(13860.0-(462.0-(132.0-
-      (99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0+(13860.0-(462.0-(132.0-(99.0
-      -140.0/w2)/w2)/w2)/w2)/w/166320.0) goto S170;
-    goto S30;
-S140:
-/*
-     INVERSE CDF LOGIC FOR MEAN LESS THAN 30
-*/
-    qn = pow(q,(double)n);
-    r = p/q;
-    g = r*(n+1);
-S150:
-    ix = 0;
-    f = qn;
-    u = ranf();
-S160:
-    if(u < f) goto S170;
-    if(ix > 110) goto S150;
-    u -= f;
-    ix += 1;
-    f *= (g/ix-r);
-    goto S160;
-S170:
-    if(psave > 0.5) ix = n-ix;
-    ignbin = ix;
-    return ignbin;
-}
-long ignnbn(long n,float p)
-/*
-**********************************************************************
- 
-     long ignnbn(long n,float p)
-                GENerate Negative BiNomial random deviate
-                              Function
-     Generates a single random deviate from a negative binomial
-     distribution.
-                              Arguments
-     N  --> The number of trials in the negative binomial distribution
-            from which a random deviate is to be generated.
-     P  --> The probability of an event.
-                              Method
-     Algorithm from page 480 of
- 
-     Devroye, Luc
- 
-     Non-Uniform Random Variate Generation.  Springer-Verlag,
-     New York, 1986.
-**********************************************************************
-*/
-{
-static long ignnbn;
-static float y,a,r;
-/*
-     ..
-     .. Executable Statements ..
-*/
-/*
-     Check Arguments
-*/
-    if(n < 0) ftnstop("N < 0 in IGNNBN");
-    if(p <= 0.0F) ftnstop("P <= 0 in IGNNBN");
-    if(p >= 1.0F) ftnstop("P >= 1 in IGNNBN");
-/*
-     Generate Y, a random gamma (n,(1-p)/p) variable
-*/
-    r = (float)n;
-    a = p/(1.0F-p);
-    y = gengam(a,r);
-/*
-     Generate a random Poisson(y) variable
-*/
-    ignnbn = ignpoi(y);
-    return ignnbn;
-}
-long ignpoi(float mu)
-/*
-**********************************************************************
-     long ignpoi(float mu)
-                    GENerate POIsson random deviate
-                              Function
-     Generates a single random deviate from a Poisson
-     distribution with mean AV.
-                              Arguments
-     av --> The mean of the Poisson distribution from which
-            a random deviate is to be generated.
-     genexp <-- The random deviate.
-                              Method
-     Renames KPOIS from TOMS as slightly modified by BWB to use RANF
-     instead of SUNIF.
-     For details see:
-               Ahrens, J.H. and Dieter, U.
-               Computer Generation of Poisson Deviates
-               From Modified Normal Distributions.
-               ACM Trans. Math. Software, 8, 2
-               (June 1982),163-179
-**********************************************************************
-**********************************************************************
-                                                                      
-                                                                      
-     P O I S S O N  DISTRIBUTION                                      
-                                                                      
-                                                                      
-**********************************************************************
-**********************************************************************
-                                                                      
-     FOR DETAILS SEE:                                                 
-                                                                      
-               AHRENS, J.H. AND DIETER, U.                            
-               COMPUTER GENERATION OF POISSON DEVIATES                
-               FROM MODIFIED NORMAL DISTRIBUTIONS.                    
-               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. 
-                                                                      
-     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)  
-                                                                      
-**********************************************************************
-      INTEGER FUNCTION IGNPOI(IR,MU)
-     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
-             MU=MEAN MU OF THE POISSON DISTRIBUTION
-     OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
-     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
-     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
-     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
-     SEPARATION OF CASES A AND B
-*/
-{
-    extern float fsign( float num, float sign );
-    const float a0 = -0.5;
-    const float a1 = 0.3333333;
-    const float a2 = -0.2500068;
-    const float a3 = 0.2000118;
-    const float a4 = -0.1661269;
-    const float a5 = 0.1421878;
-    const float a6 = -0.1384794;
-    const float a7 = 0.125006;
-    static float muold = 0.0;
-    static float muprev = 0.0;
-    const float fact[10] = {
-        1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0
-    };
-    static long ignpoi=0,j=0,k=0,kflag=0,l=0,m=0;
-    static float b1,b2,c=0.0,c0=0.0,c1=0.0,c2=0.0,c3=0.0,d=0.0,del,difmuk=0.0,
-          e=0.0,fk=0.0,fx,fy,g,omega=0.0,p=0.0,p0=0.0,px,py,q=0.0,s=0.0,
-          t,u=0.0,v,x,xx,pp[35];
-
-    if (mu <= 0.0) return 0;
-    if(mu == muprev) goto S10;
-    if(mu < 10.0) goto S120;
-/*
-     C A S E  A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
-*/
-    muprev = mu;
-    s = sqrt(mu);
-    d = 6.0*mu*mu;
-/*
-             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
-             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
-             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
-*/
-    l = (long) (mu-1.1484);
-S10:
-/*
-     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
-*/
-    g = mu+s*snorm();
-    if(g < 0.0) goto S20;
-    ignpoi = (long) (g);
-/*
-     STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
-*/
-    if(ignpoi >= l) return ignpoi;
-/*
-     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
-*/
-    fk = (float)ignpoi;
-    difmuk = mu-fk;
-    u = ranf();
-    if(d*u >= difmuk*difmuk*difmuk) return ignpoi;
-S20:
-/*
-     STEP P. PREPARATIONS FOR STEPS Q AND H.
-             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
-             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
-             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
-             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
-             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
-*/
-    if(mu == muold) goto S30;
-    muold = mu;
-    omega = 0.3989423/s;
-    b1 = 4.166667E-2/mu;
-    b2 = 0.3*b1*b1;
-    c3 = 0.1428571*b1*b2;
-    c2 = b2-15.0*c3;
-    c1 = b1-6.0*b2+45.0*c3;
-    c0 = 1.0-b1+3.0*b2-15.0*c3;
-    c = 0.1069/mu;
-S30:
-    if(g < 0.0) goto S50;
-/*
-             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
-*/
-    kflag = 0;
-    goto S70;
-S40:
-/*
-     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
-*/
-    if(fy-u*fy <= py*exp(px-fx)) return ignpoi;
-S50:
-/*
-     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
-             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
-             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
-*/
-    e = sexpo();
-    u = ranf();
-    u += (u-1.0);
-    t = 1.8+fsign(e,u);
-    if(t <= -0.6744) goto S50;
-    ignpoi = (long) (mu+s*t);
-    fk = (float)ignpoi;
-    difmuk = mu-fk;
-/*
-             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
-*/
-    kflag = 1;
-    goto S70;
-S60:
-/*
-     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
-*/
-    if(c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50;
-    return ignpoi;
-S70:
-/*
-     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
-             CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
-*/
-    if(ignpoi >= 10) goto S80;
-    px = -mu;
-    py = pow(mu,(double)ignpoi)/ *(fact+ignpoi);
-    goto S110;
-S80:
-/*
-             CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
-             A0-A7 FOR ACCURACY WHEN ADVISABLE
-             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
-*/
-    del = 8.333333E-2/fk;
-    del -= (4.8*del*del*del);
-    v = difmuk/fk;
-    if(fabs(v) <= 0.25) goto S90;
-    px = fk*log(1.0+v)-difmuk-del;
-    goto S100;
-S90:
-    px = fk*v*v*(((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0)-del;
-S100:
-    py = 0.3989423/sqrt(fk);
-S110:
-    x = (0.5-difmuk)/s;
-    xx = x*x;
-    fx = -0.5*xx;
-    fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0);
-    if(kflag <= 0) goto S40;
-    goto S60;
-S120:
-/*
-     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
-*/
-    muprev = 0.0;
-    if(mu == muold) goto S130;
-    muold = mu;
-    m = max(1L,(long) (mu));
-    l = 0;
-    p = exp(-mu);
-    q = p0 = p;
-S130:
-/*
-     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
-*/
-    u = ranf();
-    ignpoi = 0;
-    if(u <= p0) return ignpoi;
-/*
-     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
-             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
-             (0.458=PP(9) FOR MU=10)
-*/
-    if(l == 0) goto S150;
-    j = 1;
-    if(u > 0.458) j = min(l,m);
-    for(k=j; k<=l; k++) {
-        if(u <= *(pp+k-1)) goto S180;
-    }
-    if(l == 35) goto S130;
-S150:
-/*
-     STEP C. CREATION OF NEW POISSON PROBABILITIES P
-             AND THEIR CUMULATIVES Q=PP(K)
-*/
-    l += 1;
-    for(k=l; k<=35; k++) {
-        p = p*mu/(float)k;
-        q += p;
-        *(pp+k-1) = q;
-        if(u <= q) goto S170;
-    }
-    l = 35;
-    goto S130;
-S170:
-    l = k;
-S180:
-    ignpoi = k;
-    return ignpoi;
-}
-long ignuin(long low,long high)
-/*
-**********************************************************************
-     long ignuin(long low,long high)
-               GeNerate Uniform INteger
-                              Function
-     Generates an integer uniformly distributed between LOW and HIGH.
-                              Arguments
-     low --> Low bound (inclusive) on integer value to be generated
-     high --> High bound (inclusive) on integer value to be generated
-                              Note
-     If (HIGH-LOW) > 2,147,483,561 prints error message on * unit and
-     stops the program.
-**********************************************************************
-     IGNLGI generates integers between 1 and 2147483562
-     MAXNUM is 1 less than maximum generable value
-*/
-{
-#define maxnum 2147483561L
-static long ignuin,ign,maxnow,range,ranp1;
-
-    if(!(low > high)) goto S10;
-    PyErr_SetString(PyExc_ValueError, "low > high in ignuin");
-    return 0;
-
-S10:
-    range = high-low;
-    if(!(range > maxnum)) goto S20;
-    PyErr_SetString(PyExc_ValueError, "high - low too large in ignuin");
-    return 0;
-
-S20:
-    if(!(low == high)) goto S30;
-    ignuin = low;
-    return ignuin;
-
-S30:
-/*
-     Number to be generated should be in range 0..RANGE
-     Set MAXNOW so that the number of integers in 0..MAXNOW is an
-     integral multiple of the number in 0..RANGE
-*/
-    ranp1 = range+1;
-    maxnow = maxnum/ranp1*ranp1;
-S40:
-    ign = ignlgi()-1;
-    if(!(ign <= maxnow)) goto S50;
-    ignuin = low+ign%ranp1;
-    return ignuin;
-S50:
-    goto S40;
-#undef maxnum
-#undef err1
-#undef err2
-}
-long lennob( char *str )
-/* 
-Returns the length of str ignoring trailing blanks but not 
-other white space.
-*/
-{
-long i, i_nb;
-
-for (i=0, i_nb= -1L; *(str+i); i++)
-    if ( *(str+i) != ' ' ) i_nb = i;
-return (i_nb+1);
-    }
-long mltmod(long a,long s,long m)
-/*
-**********************************************************************
-     long mltmod(long a,long s,long m)
-                    Returns (A*S) MOD M
-     This is a transcription from Pascal to Fortran of routine
-     MULtMod_Decompos from the paper
-     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
-     with Splitting Facilities." ACM Transactions on Mathematical
-     Software, 17:98-111 (1991)
-                              Arguments
-     a, s, m  -->
-**********************************************************************
-*/
-{
-#define h 32768L
-static long mltmod,a0,a1,k,p,q,qh,rh;
-/*
-     H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
-      machine. On a different machine recompute H
-*/
-    if ((a <= 0 || a >= m || s <= 0 || s >= m)) {
-        char as[50], ms[50], ss[50];
-        snprintf(as, 50, "%12ld", a);
-        snprintf(ms, 50, "%12ld", m);
-        snprintf(ss, 50, "%12ld", s);
-        PyErr_Format(PyExc_ValueError,
-                     "mltmod requires 0 < a (%s) < m (%s) and 0 < s (%s) < m",
-                     as, ms, ss);
-        return 0;
-    }
-    if(!(a < h)) goto S20;
-    a0 = a;
-    p = 0;
-    goto S120;
-S20:
-    a1 = a/h;
-    a0 = a-h*a1;
-    qh = m/h;
-    rh = m-h*qh;
-    if(!(a1 >= h)) goto S50;
-    a1 -= h;
-    k = s/qh;
-    p = h*(s-k*qh)-k*rh;
-S30:
-    if(!(p < 0)) goto S40;
-    p += m;
-    goto S30;
-S40:
-    goto S60;
-S50:
-    p = 0;
-S60:
-/*
-     P = (A2*S*H)MOD M
-*/
-    if(!(a1 != 0)) goto S90;
-    q = m/a1;
-    k = s/q;
-    p -= (k*(m-a1*q));
-    if(p > 0) p -= m;
-    p += (a1*(s-k*q));
-S70:
-    if(!(p < 0)) goto S80;
-    p += m;
-    goto S70;
-S90:
-S80:
-    k = p/qh;
-/*
-     P = ((A2*H + A1)*S)MOD M
-*/
-    p = h*(p-k*qh)-k*rh;
-S100:
-    if(!(p < 0)) goto S110;
-    p += m;
-    goto S100;
-S120:
-S110:
-    if(!(a0 != 0)) goto S150;
-/*
-     P = ((A2*H + A1)*H*S)MOD M
-*/
-    q = m/a0;
-    k = s/q;
-    p -= (k*(m-a0*q));
-    if(p > 0) p -= m;
-    p += (a0*(s-k*q));
-S130:
-    if(!(p < 0)) goto S140;
-    p += m;
-    goto S130;
-S150:
-S140:
-    mltmod = p;
-    return mltmod;
-#undef h
-}
-void phrtsd(char* phrase,long *seed1,long *seed2)
-/*
-**********************************************************************
-     void phrtsd(char* phrase,long *seed1,long *seed2)
-               PHRase To SeeDs
-
-                              Function
-
-     Uses a phrase (character string) to generate two seeds for the RGN
-     random number generator.
-                              Arguments
-     phrase --> Phrase to be used for random number generation
-      
-     seed1 <-- First seed for generator
-                        
-     seed2 <-- Second seed for generator
-                        
-                              Note
-
-     Trailing blanks are eliminated before the seeds are generated.
-     Generated seed values will fall in the range 1..2^30
-     (1..1,073,741,824)
-**********************************************************************
-*/
-{
-
-static char table[] =
-"abcdefghijklmnopqrstuvwxyz\
-ABCDEFGHIJKLMNOPQRSTUVWXYZ\
-0123456789\
-!@#$%^&*()_+[];:'\\\"<>?,./";
-
-long ix;
-
-static long twop30 = 1073741824L;
-static long shift[5] = {
-    1L,64L,4096L,262144L,16777216L
-};
-static long i,ichr,j,lphr,values[5];
-extern long lennob(char *str);
-
-    *seed1 = 1234567890L;
-    *seed2 = 123456789L;
-    lphr = lennob(phrase); 
-    if(lphr < 1) return;
-    for(i=0; i<=(lphr-1); i++) {
-	for (ix=0; table[ix]; ix++) if (*(phrase+i) == table[ix]) break; 
-        if (!table[ix]) ix = 0;
-        ichr = ix % 64;
-        if(ichr == 0) ichr = 63;
-        for(j=1; j<=5; j++) {
-            *(values+j-1) = ichr-j;
-            if(*(values+j-1) < 1) *(values+j-1) += 63;
-        }
-        for(j=1; j<=5; j++) {
-            *seed1 = ( *seed1+*(shift+j-1)**(values+j-1) ) % twop30;
-            *seed2 = ( *seed2+*(shift+j-1)**(values+6-j-1) )  % twop30;
-        }
-    }
-#undef twop30
-}
-double ranf(void)
-/*
-**********************************************************************
-     double ranf(void)
-                RANDom number generator as a Function
-     Returns a random floating point number from a uniform distribution
-     over 0 - 1 (endpoints of this interval are not returned) using the
-     current generator
-     This is a transcription from Pascal to Fortran of routine
-     Uniform_01 from the paper
-     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
-     with Splitting Facilities." ACM Transactions on Mathematical
-     Software, 17:98-111 (1991)
-**********************************************************************
-*/
-{
-    double ranf_;
-/*
-     4.656613057E-10 is 1/M1  M1 is set in a data statement in IGNLGI
-      and is currently 2147483563. If M1 changes, change this also.
-*/
-    ranf_ = ignlgi()*4.6566130573917691e-10;
-    return ranf_;
-}
-void setgmn(float *meanv,float *covm,long p,float *parm)
-/*
-**********************************************************************
-     void setgmn(float *meanv,float *covm,long p,float *parm)
-            SET Generate Multivariate Normal random deviate
-                              Function
-      Places P, MEANV, and the Cholesky factoriztion of COVM
-      in GENMN.
-                              Arguments
-     meanv --> Mean vector of multivariate normal distribution.
-     covm   <--> (Input) Covariance   matrix    of  the  multivariate
-                 normal distribution
-                 (Output) Destroyed on output
-     p     --> Dimension of the normal, or length of MEANV.
-     parm <-- Array of parameters needed to generate multivariate norma
-                deviates (P, MEANV and Cholesky decomposition of
-                COVM).
-                1 : 1                - P
-                2 : P + 1            - MEANV
-                P+2 : P*(P+3)/2 + 1  - Cholesky decomposition of COVM
-               Needed dimension is (p*(p+3)/2 + 1)
-**********************************************************************
-*/
-{
-extern void spofa(float *a,long lda,long n,long *info);
-static long T1;
-static long i,icount,info,j,D2,D3,D4,D5;
-    T1 = p*(p+3)/2+1;
-/*
-     TEST THE INPUT
-*/
-    if((p <= 0)) {
-        char ps[50];
-        snprintf(ps, 50, "%12ld", p);
-        PyErr_Format(PyExc_ValueError, "P=%s nonpositive in SETGMN", ps);
-        return;
-    }
-    *parm = p;
-/*
-     PUT P AND MEANV INTO PARM
-*/
-    for(i=2,D2=1,D3=(p+1-i+D2)/D2; D3>0; D3--,i+=D2) *(parm+i-1) = *(meanv+i-2);
-/*
-      Cholesky decomposition to find A s.t. trans(A)*(A) = COVM
-*/
-    spofa(covm,p,p,&info);
-    if (info != 0) {
-        PyErr_SetString(PyExc_ValueError, "COVM not positive definite in SETGMN");
-        return;
-    }
-    icount = p+1;
-/*
-     PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM
-          COVM(1,1) = PARM(P+2)
-          COVM(1,2) = PARM(P+3)
-                    :
-          COVM(1,P) = PARM(2P+1)
-          COVM(2,2) = PARM(2P+2)  ...
-*/
-    for(i=1,D4=1,D5=(p-i+D4)/D4; D5>0; D5--,i+=D4) {
-        for(j=i-1; j<p; j++) {
-            icount += 1;
-            *(parm+icount-1) = *(covm+i-1+j*p);
-        }
-    }
-}
-float sexpo(void)
-/*
-**********************************************************************
-                                                                      
-                                                                      
-     (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION                
-                                                                      
-                                                                      
-**********************************************************************
-**********************************************************************
-                                                                      
-     FOR DETAILS SEE:                                                 
-                                                                      
-               AHRENS, J.H. AND DIETER, U.                            
-               COMPUTER METHODS FOR SAMPLING FROM THE                 
-               EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  
-               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               
-                                                                      
-     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       
-     'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)       
-                                                                      
-     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
-     SUNIF.  The argument IR thus goes away.                          
-                                                                      
-**********************************************************************
-     Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
-     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
-*/
-{
-static float q[8] = {
-    0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0
-};
-static long i;
-static float sexpo,a,u,ustar,umin;
-static float *q1 = q;
-    a = 0.0;
-    u = ranf();
-    goto S30;
-S20:
-    a += *q1;
-S30:
-    u += u;
-    if(u <= 1.0) goto S20;
-    u -= 1.0;
-    if(u > *q1) goto S60;
-    sexpo = a+u;
-    return sexpo;
-S60:
-    i = 1;
-    ustar = ranf();
-    umin = ustar;
-S70:
-    ustar = ranf();
-    if(ustar < umin) umin = ustar;
-    i += 1;
-    if(u > *(q+i-1)) goto S70;
-    sexpo = a+umin**q1;
-    return sexpo;
-}
-float sgamma(float a)
-/*
-**********************************************************************
-                                                                      
-                                                                      
-     (STANDARD-)  G A M M A  DISTRIBUTION                             
-                                                                      
-                                                                      
-**********************************************************************
-**********************************************************************
-                                                                      
-               PARAMETER  A >= 1.0  !                                 
-                                                                      
-**********************************************************************
-                                                                      
-     FOR DETAILS SEE:                                                 
-                                                                      
-               AHRENS, J.H. AND DIETER, U.                            
-               GENERATING GAMMA VARIATES BY A                         
-               MODIFIED REJECTION TECHNIQUE.                          
-               COMM. ACM, 25,1 (JAN. 1982), 47 - 54.                  
-                                                                      
-     STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER     
-                                 (STRAIGHTFORWARD IMPLEMENTATION)     
-                                                                      
-     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
-     SUNIF.  The argument IR thus goes away.                          
-                                                                      
-**********************************************************************
-                                                                      
-               PARAMETER  0.0 < A < 1.0  !                            
-                                                                      
-**********************************************************************
-                                                                      
-     FOR DETAILS SEE:                                                 
-                                                                      
-               AHRENS, J.H. AND DIETER, U.                            
-               COMPUTER METHODS FOR SAMPLING FROM GAMMA,              
-               BETA, POISSON AND BINOMIAL DISTRIBUTIONS.              
-               COMPUTING, 12 (1974), 223 - 246.                       
-                                                                      
-     (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)    
-                                                                      
-**********************************************************************
-     INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
-     OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
-     COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
-     COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
-     COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
-     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
-     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
-*/
-{
-extern float fsign( float num, float sign );
-static float q1 = 4.166669E-2;
-static float q2 = 2.083148E-2;
-static float q3 = 8.01191E-3;
-static float q4 = 1.44121E-3;
-static float q5 = -7.388E-5;
-static float q6 = 2.4511E-4;
-static float q7 = 2.424E-4;
-static float a1 = 0.3333333;
-static float a2 = -0.250003;
-static float a3 = 0.2000062;
-static float a4 = -0.1662921;
-static float a5 = 0.1423657;
-static float a6 = -0.1367177;
-static float a7 = 0.1233795;
-static float e1 = 1.0;
-static float e2 = 0.4999897;
-static float e3 = 0.166829;
-static float e4 = 4.07753E-2;
-static float e5 = 1.0293E-2;
-static float aa = 0.0;
-static float aaa = 0.0;
-static float sqrt32 = 5.656854;
-static float sgamma,s2,s,d,t,x,u,r,q0,b,si,c,v,q,e,w,p;
-    if(a == aa) goto S10;
-    if(a < 1.0) goto S120;
-/*
-     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
-*/
-    aa = a;
-    s2 = a-0.5;
-    s = sqrt(s2);
-    d = sqrt32-12.0*s;
-S10:
-/*
-     STEP  2:  T=STANDARD NORMAL DEVIATE,
-               X=(S,1/2)-NORMAL DEVIATE.
-               IMMEDIATE ACCEPTANCE (I)
-*/
-    t = snorm();
-    x = s+0.5*t;
-    sgamma = x*x;
-    if(t >= 0.0) return sgamma;
-/*
-     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
-*/
-    u = ranf();
-    if(d*u <= t*t*t) return sgamma;
-/*
-     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
-*/
-    if(a == aaa) goto S40;
-    aaa = a;
-    r = 1.0/ a;
-    q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
-/*
-               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
-               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
-               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
-*/
-    if(a <= 3.686) goto S30;
-    if(a <= 13.022) goto S20;
-/*
-               CASE 3:  A .GT. 13.022
-*/
-    b = 1.77;
-    si = 0.75;
-    c = 0.1515/s;
-    goto S40;
-S20:
-/*
-               CASE 2:  3.686 .LT. A .LE. 13.022
-*/
-    b = 1.654+7.6E-3*s2;
-    si = 1.68/s+0.275;
-    c = 6.2E-2/s+2.4E-2;
-    goto S40;
-S30:
-/*
-               CASE 1:  A .LE. 3.686
-*/
-    b = 0.463+s+0.178*s2;
-    si = 1.235;
-    c = 0.195/s-7.9E-2+1.6E-1*s;
-S40:
-/*
-     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
-*/
-    if(x <= 0.0) goto S70;
-/*
-     STEP  6:  CALCULATION OF V AND QUOTIENT Q
-*/
-    v = t/(s+s);
-    if(fabs(v) <= 0.25) goto S50;
-    q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
-    goto S60;
-S50:
-    q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
-S60:
-/*
-     STEP  7:  QUOTIENT ACCEPTANCE (Q)
-*/
-    if(log(1.0-u) <= q) return sgamma;
-S70:
-/*
-     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
-               U= 0,1 -UNIFORM DEVIATE
-               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
-*/
-    e = sexpo();
-    u = ranf();
-    u += (u-1.0);
-    t = b+fsign(si*e,u);
-/*
-     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
-*/
-    if(t < -0.7187449) goto S70;
-/*
-     STEP 10:  CALCULATION OF V AND QUOTIENT Q
-*/
-    v = t/(s+s);
-    if(fabs(v) <= 0.25) goto S80;
-    q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
-    goto S90;
-S80:
-    q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
-S90:
-/*
-     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
-*/
-    if(q <= 0.0) goto S70;
-    if(q <= 0.5) goto S100;
-    w = exp(q)-1.0;
-    goto S110;
-S100:
-    w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
-S110:
-/*
-               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
-*/
-    if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
-    x = s+0.5*t;
-    sgamma = x*x;
-    return sgamma;
-S120:
-/*
-     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
-*/
-    aa = 0.0;
-    b = 1.0+0.3678794*a;
-S130:
-    p = b*ranf();
-    if(p >= 1.0) goto S140;
-    sgamma = exp(log(p)/ a);
-    if(sexpo() < sgamma) goto S130;
-    return sgamma;
-S140:
-    sgamma = -log((b-p)/ a);
-    if(sexpo() < (1.0-a)*log(sgamma)) goto S130;
-    return sgamma;
-}
-float snorm(void)
-/*
-**********************************************************************
-                                                                      
-                                                                      
-     (STANDARD-)  N O R M A L  DISTRIBUTION                           
-                                                                      
-                                                                      
-**********************************************************************
-**********************************************************************
-                                                                      
-     FOR DETAILS SEE:                                                 
-                                                                      
-               AHRENS, J.H. AND DIETER, U.                            
-               EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             
-               SAMPLING FROM THE NORMAL DISTRIBUTION.                 
-               MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          
-                                                                      
-     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  
-     (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  
-                                                                      
-     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
-     SUNIF.  The argument IR thus goes away.                          
-                                                                      
-**********************************************************************
-     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
-     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
-*/
-{
-static float a[32] = {
-    0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
-    0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
-    0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
-    1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
-    1.862732,2.153875
-};
-static float d[31] = {
-    0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
-    0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
-    0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
-    0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
-};
-static float t[31] = {
-    7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
-    1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
-    2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
-    4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
-    9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
-};
-static float h[31] = {
-    3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
-    4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
-    4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
-    5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
-    8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
-};
-static long i;
-static float snorm,u,s,ustar,aa,w,y,tt;
-    u = ranf();
-    s = 0.0;
-    if(u > 0.5) s = 1.0;
-    u += (u-s);
-    u = 32.0*u;
-    i = (long) (u);
-    if(i == 32) i = 31;
-    if(i == 0) goto S100;
-/*
-                                START CENTER
-*/
-    ustar = u-(float)i;
-    aa = *(a+i-1);
-S40:
-    if(ustar <= *(t+i-1)) goto S60;
-    w = (ustar-*(t+i-1))**(h+i-1);
-S50:
-/*
-                                EXIT   (BOTH CASES)
-*/
-    y = aa+w;
-    snorm = y;
-    if(s == 1.0) snorm = -y;
-    return snorm;
-S60:
-/*
-                                CENTER CONTINUED
-*/
-    u = ranf();
-    w = u*(*(a+i)-aa);
-    tt = (0.5*w+aa)*w;
-    goto S80;
-S70:
-    tt = u;
-    ustar = ranf();
-S80:
-    if(ustar > tt) goto S50;
-    u = ranf();
-    if(ustar >= u) goto S70;
-    ustar = ranf();
-    goto S40;
-S100:
-/*
-                                START TAIL
-*/
-    i = 6;
-    aa = *(a+31);
-    goto S120;
-S110:
-    aa += *(d+i-1);
-    i += 1;
-S120:
-    u += u;
-    if(u < 1.0) goto S110;
-    u -= 1.0;
-S140:
-    w = u**(d+i-1);
-    tt = (0.5*w+aa)*w;
-    goto S160;
-S150:
-    tt = u;
-S160:
-    ustar = ranf();
-    if(ustar > tt) goto S50;
-    u = ranf();
-    if(ustar >= u) goto S150;
-    u = ranf();
-    goto S140;
-}
-float fsign( float num, float sign )
-/* Transfers sign of argument sign to argument num */
-{
-if ( ( sign>0.0f && num<0.0f ) || ( sign<0.0f && num>0.0f ) )
-    return -num;
-else return num;
-}
-/************************************************************************
-FTNSTOP:
-Prints msg to standard error and then exits
-************************************************************************/
-void ftnstop(char* msg)
-/* msg - error message */
-{
-    if (msg == NULL) {
-        msg = "";
-    }
-    PyErr_SetString(PyExc_RuntimeError, msg);
-    return;
-}
diff --git a/scipy/corelib/random_lite/ranlib.h b/scipy/corelib/random_lite/ranlib.h
deleted file mode 100644
index 08af41329..000000000
--- a/scipy/corelib/random_lite/ranlib.h
+++ /dev/null
@@ -1,33 +0,0 @@
-/* Prototypes for all user accessible RANLIB routines */
-
-extern void advnst(long k);
-extern float genbet(float aa,float bb);
-extern float genchi(float df);
-extern float genexp(float av);
-extern float genf(float dfn, float dfd);
-extern float gengam(float a,float r);
-extern void genmn(float *parm,float *x,float *work);
-extern void genmul(long n,float *p,long ncat,long *ix);
-extern float gennch(float df,float xnonc);
-extern float gennf(float dfn, float dfd, float xnonc);
-extern float gennor(float av,float sd);
-extern void genprm(long *iarray,int larray);
-extern float genunf(float low,float high);
-extern void getsd(long *iseed1,long *iseed2);
-extern void gscgn(long getset,long *g);
-extern long ignbin(long n,float pp);
-extern long ignnbn(long n,float p);
-extern long ignlgi(void);
-extern long ignpoi(float mu);
-extern long ignuin(long low,long high);
-extern void initgn(long isdtyp);
-extern long mltmod(long a,long s,long m);
-extern void phrtsd(char* phrase,long* seed1,long* seed2);
-extern double ranf(void);
-extern void setall(long iseed1,long iseed2);
-extern void setant(long qvalue);
-extern void setgmn(float *meanv,float *covm,long p,float *parm);
-extern void setsd(long iseed1,long iseed2);
-extern float sexpo(void);
-extern float sgamma(float a);
-extern float snorm(void);
diff --git a/scipy/corelib/random_lite/ranlibmodule.c b/scipy/corelib/random_lite/ranlibmodule.c
deleted file mode 100644
index d8a913005..000000000
--- a/scipy/corelib/random_lite/ranlibmodule.c
+++ /dev/null
@@ -1,343 +0,0 @@
-#include "Python.h"
-#include "scipy/arrayobject.h"
-#include "ranlib.h"
-#include "stdio.h"
-
-/* ----------------------------------------------------- */
-
-static PyObject*
-get_continuous_random(int num_dist_params, PyObject* self, PyObject* args, void* fun) {
-  PyArrayObject *op;
-  double *out_ptr;
-  int i, n=-1;
-  float a, b, c;
-
-  switch(num_dist_params) {
-  case 0:
-    if( !PyArg_ParseTuple(args, "|i", &n) )
-      return NULL;
-    break;
-  case 1:
-    if( !PyArg_ParseTuple(args, "f|i", &a, &n) )
-      return NULL;
-    break;
-  case 2:
-    if( !PyArg_ParseTuple(args, "ff|i", &a, &b, &n) )
-      return NULL;
-    break;
-  case 3:
-    if( !PyArg_ParseTuple(args, "fff|i", &a, &b, &c, &n) )
-      return NULL;
-    break;
-  }
-  if( n == -1 )
-    n = 1;
-
-  /* Create a 1 dimensional array of length n of type double */
-  op = (PyArrayObject*) PyArray_FromDims(1, &n, PyArray_DOUBLE);
-  if( op == NULL )
-    return NULL;
-
-
-  out_ptr = (double *) op->data;
-  for(i=0; i<n; i++) {
-    switch(num_dist_params) {
-    case 0:
-      *out_ptr = (double) ((float (*)(void)) fun)();
-      break;
-    case 1:
-      *out_ptr = (double) ((float (*)(float)) fun)(a);
-      break;
-    case 2:
-      *out_ptr = (double) ((float (*)(float, float)) fun)(a,b);
-      break;
-    case 3:
-      *out_ptr = (double) ((float (*)(float, float, float)) fun)(a,b,c);
-      break;
-    }
-    out_ptr++;
-  }
-
-  return PyArray_Return(op);
-}
-
-
-static PyObject*
-get_discrete_scalar_random(int num_integer_args, PyObject* self, PyObject* args, void* fun) {
-  long int_arg;
-  int n=-1, i;
-  long* out_ptr;
-  PyArrayObject* op;
-  float float_arg;
-
-  switch( num_integer_args ) {
-  case 0:
-    if( !PyArg_ParseTuple(args, "f|i", &float_arg, &n) ) {
-      return NULL;
-    }
-    break;
-  case 1:
-    if( !PyArg_ParseTuple(args, "lf|i", &int_arg, &float_arg, &n) ) {
-      return NULL;
-    }
-    break;
-  }
-  if( n==-1 ) {
-    n = 1;
-  }
-
-  /* Create a 1 dimensional array of length n of type long */
-  op = (PyArrayObject*) PyArray_FromDims(1, &n, PyArray_LONG);
-  if( op == NULL ) {
-    return NULL;
-  }
-
-  out_ptr = (long*) op->data;
-  for(i=0; i<n; i++) {
-    switch( num_integer_args ) {
-    case 0:
-      *out_ptr = ((long (*)(float)) fun)(float_arg);
-      break;
-    case 1:
-      *out_ptr = ((long (*)(long, float)) fun)(int_arg, float_arg);
-      break;
-    }
-    out_ptr++;
-  }
-
-  return PyArray_Return(op);
-}
-
-
-static char random_sample__doc__[] ="";
-
-static PyObject *
-random_sample(PyObject *self, PyObject *args) {
-  PyArrayObject *op;
-  double *out_ptr;
-  int i, n=1;
-
-  if (!PyArg_ParseTuple(args, "|i", &n)) {
-      return NULL;
-  }
-
-  /* Create a 1 dimensional array of length n of type double */
-  op = (PyArrayObject*) PyArray_FromDims(1, &n, PyArray_DOUBLE);
-  if (op == NULL) {
-    return NULL;
-  }
-
-
-  out_ptr = (double *) op->data;
-  for (i=0; i<n; i++) {
-      *out_ptr = ranf();
-      out_ptr++;
-  }
-
-  return PyArray_Return(op);
-}
-
-
-static char standard_normal__doc__[] ="";
-
-static PyObject *
-standard_normal(PyObject *self, PyObject *args) {
-  return get_continuous_random(0, self, args, snorm);
-}
-
-
-static char beta__doc__[] ="";
-
-static PyObject *
-beta(PyObject *self, PyObject *args) {
-  return get_continuous_random(2, self, args, genbet);
-}
-
-
-static char gamma__doc__[] ="";
-
-static PyObject *
-/* there is a function named `gamma' in some libm's */
-_gamma(PyObject *self, PyObject *args) {
-  return get_continuous_random(2, self, args, gengam);
-}
-
-
-static char f__doc__[] ="";
-
-static PyObject *
-f(PyObject *self, PyObject *args) {
-  return get_continuous_random(2, self, args, genf);
-}
-
-
-static char noncentral_f__doc__[] ="";
-
-static PyObject *
-noncentral_f(PyObject *self, PyObject *args) {
-  return get_continuous_random(3, self, args, gennf);
-}
-
-
-static char noncentral_chisquare__doc__[] ="";
-
-static PyObject *
-noncentral_chisquare(PyObject *self, PyObject *args) {
-  return get_continuous_random(2, self, args, gennch);
-}
-
-
-static char chisquare__doc__[] ="";
-
-static PyObject *
-chisquare(PyObject *self, PyObject *args) {
-  return get_continuous_random(1, self, args, genchi);
-}
-
-
-static char binomial__doc__[] ="";
-
-static PyObject *
-binomial(PyObject *self, PyObject *args) {
-  return get_discrete_scalar_random(1, self, args, ignbin);
-}
-
-
-static char negative_binomial__doc__[] ="";
-
-static PyObject *
-negative_binomial(PyObject *self, PyObject *args) {
-  return get_discrete_scalar_random(1, self, args, ignnbn);
-}
-
-static char poisson__doc__[] ="";
-
-static PyObject *
-poisson(PyObject *self, PyObject *args) {
-  return get_discrete_scalar_random(0, self, args, ignpoi);
-}
-
-
-static char multinomial__doc__[] ="";
-
-static PyObject*
-multinomial(PyObject* self, PyObject* args) {
-  int n=-1, i;
-  long num_trials, num_categories;
-  char* out_ptr;
-  PyArrayObject* priors_array;
-  PyObject* priors_object;
-  PyArrayObject* op;
-  int out_dimensions[2];
-
-  if( !PyArg_ParseTuple(args, "lO|i", &num_trials, &priors_object, &n) ) {
-    return NULL;
-  }
-  priors_array = (PyArrayObject*) PyArray_ContiguousFromObject(priors_object, PyArray_FLOAT, 1, 1);
-  if( priors_array == NULL ) {
-    return NULL;
-  }
-  num_categories = priors_array->dimensions[0]+1;
-  if( n==-1 ) {
-    n = 1;
-  }
-
-  /* Create an n by num_categories array of long */
-  out_dimensions[0] = n;
-  out_dimensions[1] = num_categories;
-  op = (PyArrayObject*) PyArray_FromDims(2, out_dimensions, PyArray_LONG);
-  if( op == NULL ) {
-    return NULL;
-  }
-
-  out_ptr = op->data;
-  for(i=0; i<n; i++) {
-    genmul(num_trials, (float*)(priors_array->data), num_categories, (long*) out_ptr);
-    out_ptr += op->strides[0];
-  }
-
-  return PyArray_Return(op);
-}
-
-
-static PyObject *
-random_set_seeds(PyObject *self, PyObject *args)
-{
-  long seed1, seed2;
-
-  if (!PyArg_ParseTuple(args, "ll", &seed1, &seed2)) return NULL;
-
-
-  setall(seed1, seed2);
-  if (PyErr_Occurred ()) return NULL;
-  Py_INCREF(Py_None);
-  return (PyObject *)Py_None;
-}
-
-
-static PyObject *
-random_get_seeds(PyObject *self, PyObject *args)
-{
-  long seed1, seed2;
-
-  if (!PyArg_ParseTuple(args, "")) return NULL;
-
-  getsd(&seed1, &seed2);
-
-  return Py_BuildValue("ll", seed1, seed2);
-}
-
-
-/* Missing interfaces to */
-/* exponential (genexp), multivariate normal (genmn),
-   normal (gennor), permutation (genprm), uniform (genunf),
-   standard exponential (sexpo), standard gamma (sgamma) */
-
-/* List of methods defined in the module */
-
-static struct PyMethodDef random_methods[] = {
- {"sample",     random_sample,          1,      random_sample__doc__},
- {"standard_normal", standard_normal,   1,      standard_normal__doc__},
- {"beta",	beta,                   1,      beta__doc__},
- {"gamma",	_gamma,                  1,      gamma__doc__},
- {"f",	        f,                      1,      f__doc__},
- {"noncentral_f", noncentral_f,         1,      noncentral_f__doc__},
- {"chisquare",	chisquare,              1,      chisquare__doc__},
- {"noncentral_chisquare", noncentral_chisquare,
-                                        1,      noncentral_chisquare__doc__},
- {"binomial",	binomial,               1,      binomial__doc__},
- {"negative_binomial", negative_binomial,
-                                        1,      negative_binomial__doc__},
- {"multinomial", multinomial,           1,      multinomial__doc__},
- {"poisson",    poisson,                1,      poisson__doc__},
- {"set_seeds",  random_set_seeds,       1, },
- {"get_seeds",  random_get_seeds,       1, },
- {NULL,		NULL}		/* sentinel */
-};
-
-
-/* Initialization function for the module (*must* be called initranlib) */
-
-static char random_module_documentation[] =
-""
-;
-
-DL_EXPORT(void)
-initranlib(void)
-{
-	PyObject *m;
-
-	/* Create the module and add the functions */
-	m = Py_InitModule4("ranlib", random_methods,
-		random_module_documentation,
-		(PyObject*)NULL,PYTHON_API_VERSION);
-
-	/* Import the array object */
-	import_array();
-
-	/* XXXX Add constants here */
-
-	/* Check for errors */
-	if (PyErr_Occurred())
-		Py_FatalError("can't initialize module ranlib");
-}
diff --git a/scipy/corelib/random_lite/rngmodule.c b/scipy/corelib/random_lite/rngmodule.c
deleted file mode 100644
index d2673d7d8..000000000
--- a/scipy/corelib/random_lite/rngmodule.c
+++ /dev/null
@@ -1,641 +0,0 @@
-/*
- * Univariate random number generator. Based on (and upwards compatible
- * with) Paul Dubois' URNG module. This module provides other important
- * distributions in addition to the uniform [0., 1.) distribution.
- *
- * Written by Konrad Hinsen (based on the original URNG by Paul Dubois)
- * last revision: 1997-11-6
-
- * Modified 3/11/98 Source from P. Stoll to make it ANSI
- * , allow for C++, Windows. P. Dubois fix the &PyType_Type problems.
- */
-
-
-#include "Python.h"
-#include "scipy/arrayobject.h"
-
-#include "ranf.h"
-
-#ifndef M_PI
-#define M_PI    3.14159265358979323846
-#endif
-
-/* Define the Python interface */
-
-static PyObject *ErrorObject;
-
-static PyObject *ErrorReturn(char *message)
-{
-  PyErr_SetString(ErrorObject,message);
-  return NULL;
-}
-
-/* ---------------------------------------------------------------- */
-
-/* Declarations for objects of type distribution */
-
-typedef struct {
-	PyObject_HEAD
-	double (*density)(double x, double *param);
-	void (*sample)(double *buffer, int n, double *param);
-	PyArrayObject *parameters;
-} distributionobject;
-
-staticforward PyTypeObject distributiontype;
-
-
-static char dist_call__doc__[] = 
-"density(x) = The distribution function at x"
-;
-
-static PyObject *
-dist_call(distributionobject *self, PyObject *args, PyObject *kw)
-{
-	double x, p;
-
-	if (!PyArg_ParseTuple(args, "d", &x))
-		return NULL;
-
-	p = (*self->density)(x, (double *)self->parameters->data);
-
-	return PyFloat_FromDouble(p);
-}
-
-static void
-dist_sample(distributionobject *self, double *buffer, int n)
-{
-	(*self->sample)(buffer, n, (double *)self->parameters->data);
-}
-
-
-static struct PyMethodDef dist_methods[] = {
-	{"density",	(PyCFunction)dist_call,	1,	dist_call__doc__},
-	{NULL,		NULL}		/* sentinel */
-};
-
-
-static distributionobject *
-newdistributionobject(void)
-{
-	distributionobject *self;
-
-	self = PyObject_NEW(distributionobject, &distributiontype);
-	if (self == NULL)
-		return NULL;
-	self->density = NULL;
-	self->sample = NULL;
-	self->parameters = NULL;
-	return self;
-}
-
-static void
-dist_dealloc(distributionobject *self)
-{
-	Py_XDECREF(self->parameters);
-	PyMem_DEL(self);
-}
-
-
-static PyObject *
-dist_getattr(distributionobject *self, char *name)
-{
-	return Py_FindMethod(dist_methods, (PyObject *)self, name);
-}
-
-static char distributiontype__doc__[] = 
-"Random number distribution"
-;
-
-static PyTypeObject distributiontype = {
-  PyObject_HEAD_INIT(0)
-    0,				/*ob_size*/
-    "random_distribution",	/*tp_name*/
-    sizeof(distributionobject),	/*tp_basicsize*/
-    0,				/*tp_itemsize*/
-    /* methods */
-    (destructor) dist_dealloc,	/*tp_dealloc*/
-    (printfunc)0,		/*tp_print*/
-    (getattrfunc)dist_getattr,	/*tp_getattr*/
-    (setattrfunc)0,	/*tp_setattr*/
-    (cmpfunc)0,		/*tp_compare*/
-    (reprfunc)0,		/*tp_repr*/
-    0,			/*tp_as_number*/
-    0,		/*tp_as_sequence*/
-    0,		/*tp_as_mapping*/
-    (hashfunc)0,		/*tp_hash*/
-    (ternaryfunc)dist_call,	/*tp_call*/
-    (reprfunc)0,		/*tp_str*/
-    
-    /* Space for future expansion */
-    0L,0L,0L,0L,
-    distributiontype__doc__ /* Documentation string */
-};
-
-/* ---------------------------------------------------------------- */
-
-/* Specific distributions */
-
-static PyObject *default_distribution;
-
-/* ------------------------------------------------------ */
-
-/* Default: uniform in [0., 1.) */
-
-static double
-default_density(double x, double *param)
-{
-	return (x < 0. || x >= 1.) ? 0. : 1. ;
-}
-
-static void
-default_sample(double *buffer, int n, double *param)
-{
-	int i;
-	for(i = 0; i < n; i++)
-	  buffer[i] = Ranf();
-}
-
-static PyObject *
-create_default_distribution(void)
-{
-	distributionobject *self = newdistributionobject();
-
-	if (self != NULL) {
-	  int dims[1] = { 0 };
-	  self->density = default_density;
-	  self->sample = default_sample;
-	  self->parameters = 
-	    (PyArrayObject *)PyArray_FromDims(1, dims, PyArray_DOUBLE);
-	}
-
-	return (PyObject *)self;
-}
-
-
-/* ------------------------------------------------------ */
-
-/* Uniform in [a, b) */
-
-static double
-uniform_density(double x, double *param)
-{
-	return (x < param[0] || x >= param[1]) ? 0. : 1./(param[1]-param[0]);
-}
-
-static void
-uniform_sample(double *buffer, int n, double *param)
-{
-	double w = param[1]-param[0];
-	int i;
-	for (i = 0; i < n; i++)
-		buffer[i] = param[0]+w*Ranf();
-}
-
-static PyObject *
-RNG_UniformDistribution(PyObject *self, PyObject *args)
-{
-	distributionobject *dist;
-	double a, b;
-
-	if (!PyArg_ParseTuple(args, "dd", &a, &b))
-		return NULL;
-	if (a == b)
-		return ErrorReturn("width of uniform distribution must be > 0");
-
-	dist = newdistributionobject();
-	if (dist != NULL) {
-	  int dims[1] = { 2 };
-	  double *data;
-	  dist->density = uniform_density;
-	  dist->sample = uniform_sample;
-	  dist->parameters = 
-	    (PyArrayObject *)PyArray_FromDims(1, dims, PyArray_DOUBLE);
-	  data = (double *)dist->parameters->data;
-	  data[0] = a < b ? a : b;
-	  data[1] = a > b ? a : b;
-	}
-
-	return (PyObject *)dist;
-}
-
-static char RNG_UniformDistribution__doc__[] =
-"UniformDistribution(a, b) returns a uniform distribution between a and b.\n\
-";
-
-/* ------------------------------------------------------ */
-
-/* Normal (gaussian) with mean m and standard deviation s */
-
-static double
-normal_density(double x, double *param)
-{
-	double y = (x-param[0])/param[1];
-	double n = 1./sqrt(2*M_PI)/param[1];
-	return n*exp(-0.5*y*y);
-}
-
-static void
-normal_sample(double *buffer, int n, double *param)
-{
-	int i;
-	for (i = 0; i < n; i += 2) {
-		double v1, v2, s;
-		do {
-			v1 = 2.*Ranf()-1.;
-			v2 = 2.*Ranf()-1.;
-			s = v1*v1 + v2*v2;
-		} while (s >= 1. || s == 0.);
-		s = param[1]*sqrt(-2.*log(s)/s);
-		buffer[i] = param[0]+s*v1;
-		buffer[i+1] = param[0]+s*v2;
-	}
-}
-
-static PyObject *
-RNG_NormalDistribution(PyObject *self, PyObject *args)
-{
-	distributionobject *dist;
-	double m, s;
-
-	if (!PyArg_ParseTuple(args, "dd", &m, &s))
-		return NULL;
-	if (s <= 0.)
-		return ErrorReturn("standard deviation must be positive");
-
-	dist = newdistributionobject();
-	if (dist != NULL) {
-	  int dims[1] = { 2 };
-	  double *data;
-	  dist->density = normal_density;
-	  dist->sample = normal_sample;
-	  dist->parameters = 
-	    (PyArrayObject *)PyArray_FromDims(1, dims, PyArray_DOUBLE);
-	  data = (double *)dist->parameters->data;
-	  data[0] = m;
-	  data[1] = s;
-	}
-
-	return (PyObject *)dist;
-}
-
-static char RNG_NormalDistribution__doc__[] =
-"NormalDistribution(m, s) returns a normal distribution\n\
-with mean m and standard deviation s.\n\
-";
-
-/* ------------------------------------------------------ */
-
-/* Log Normal (gaussian) with mean m and standard deviation s */
-
-static double
-lognormal_density(double x, double *param)
-{
-	double y = (log(x)-param[2])/param[3];
-	double n = 1./sqrt(2*M_PI)/param[3];
-	return n*exp(-0.5*y*y)/x;
-}
-
-static void
-lognormal_sample(double *buffer, int n, double *param)
-{
-	int i;
-	for (i = 0; i < n; i += 2) {
-		double v1, v2, s;
-		do {
-			v1 = 2.*Ranf()-1.;
-			v2 = 2.*Ranf()-1.;
-			s = v1*v1 + v2*v2;
-		} while (s >= 1. || s == 0.);
-		s = param[3]*sqrt(-2.*log(s)/s);
-		buffer[i] = exp(param[2]+s*v1);
-		buffer[i+1] = exp(param[2]+s*v2);
-	}
-}
-
-static PyObject *
-RNG_LogNormalDistribution(PyObject *self, PyObject *args)
-{
-	distributionobject *dist;
-	double m, s, mn, sn;
-
-	if (!PyArg_ParseTuple(args, "dd", &m, &s))
-		return NULL;
-	if (s <= 0.)
-		return ErrorReturn("standard deviation must be positive");
-	sn = log(1. + s*s/(m*m));
-	mn = log(m) - 0.5*sn;
-	sn = sqrt(sn);
-
-	dist = newdistributionobject();
-	if (dist != NULL) {
-	  int dims[1] = { 4 };
-	  double *data;
-	  dist->density = lognormal_density;
-	  dist->sample = lognormal_sample;
-	  dist->parameters = 
-	    (PyArrayObject *)PyArray_FromDims(1, dims, PyArray_DOUBLE);
-	  data = (double *)dist->parameters->data;
-	  data[0] = m;
-	  data[1] = s;
-	  data[2] = mn;
-	  data[3] = sn;
-	}
-
-	return (PyObject *)dist;
-}
-
-static char RNG_LogNormalDistribution__doc__[] =
-"LogNormalDistribution(m, s) returns a log normal distribution\n\
-with mean m and standard deviation s.\n\
-";
-
-/* ------------------------------------------------------ */
-
-/* Exponential distribution */
-
-static double
-expo_density(double x, double *param)
-{
-	return (x < 0) ? 0. : param[0]*exp(-param[0]*x);
-}
-
-static void
-expo_sample(double *buffer, int n, double *param)
-{
-	int i;
-	double r;
-	for (i = 0; i < n; i++) {
-		do {
-			r = Ranf();
-		} while (r == 0.);
-		buffer[i] = -log(r)/param[0];
-	}
-}
-
-static PyObject *
-RNG_ExponentialDistribution(PyObject *self, PyObject *args)
-{
-	distributionobject *dist;
-	double l;
-
-	if (!PyArg_ParseTuple(args, "d", &l))
-		return NULL;
-	if (l <= 0.)
-		return ErrorReturn("parameter must be positive");
-
-	dist = newdistributionobject();
-	if (dist != NULL) {
-	  int dims[1] = { 1 };
-	  double *data;
-	  dist->density = expo_density;
-	  dist->sample = expo_sample;
-	  dist->parameters = 
-	    (PyArrayObject *)PyArray_FromDims(1, dims, PyArray_DOUBLE);
-	  data = (double *)dist->parameters->data;
-	  data[0] = l;
-	}
-
-	return (PyObject *)dist;
-}
-
-static char RNG_ExponentialDistribution__doc__[] =
-"ExponentialDistribution(l) returns an exponential distribution.\n\
-";
-
-/* ---------------------------------------------------------------- */
-
-/* Declarations for objects of type random generator */
-
-/* Make this bigger to reduce cost of making streams independent */
-#define SAMPLE_SIZE 128
-/* Make it smaller to save space */
-
-typedef struct {
-	PyObject_HEAD
-	distributionobject *distribution;
-	u32 seed[2];
-	int position;
-	double sample[SAMPLE_SIZE];
-} rngobject;
-
-staticforward PyTypeObject rngtype;
-
-
-static char rng_ranf__doc__[] = 
-"ranf() -- return next random number from this stream"
-;
-
-/* Get the next number in this stream */
-static double 
-rng_next(rngobject *self) 
-{
-  double d;
-
-  d = self->sample[self->position++];
-  if(self->position >= SAMPLE_SIZE) {
-    self->position = 0;
-    Setranf(self->seed);
-    dist_sample(self->distribution, self->sample, SAMPLE_SIZE);
-    Getranf(self->seed);
-  }
-  return d;
-}
-
-static PyObject *
-rng_ranf(rngobject *self, PyObject *args)
-{
-	if (!PyArg_ParseTuple(args, ""))
-		return NULL;
-	return Py_BuildValue("d", rng_next(self));
-}
-
-static char rng_sample__doc__[] = 
-"sample(n) = A vector of n random values"
-;
-
-static PyObject *
-rng_sample(rngobject *self, PyObject *args)
-{
-	 /* sample(n) : Return a vector of n random numbers */
-	int i;
-	double *x;
-	int dims[1];
-	PyArrayObject *result;
-
-	if (!PyArg_ParseTuple(args, "i", dims))
-		return NULL;
-
-	if (dims[0] <= 0) 
-	  return ErrorReturn("RNG sample length cannot be <= 0.");
-
-	result = (PyArrayObject *) PyArray_FromDims (1, dims, 'd');
-
-	if (result == NULL) 
-	  return ErrorReturn("RNG sample failed to create output array.");
-
-	x = (double *) result->data;
-
-	for(i=0; i < dims[0]; i ++) 
-	  x[i] = rng_next(self);
-
-	return PyArray_Return(result);
-}
-
-static struct PyMethodDef rng_methods[] = {
-	{"ranf",	(PyCFunction)rng_ranf,	1,	rng_ranf__doc__},
-	{"sample",	(PyCFunction)rng_sample,	1,	rng_sample__doc__},
- 
-	{NULL,		NULL}		/* sentinel */
-};
-
-
-static rngobject *
-newrngobject(int seed, distributionobject *distribution)
-{
-	rngobject *self;
-
-	self = PyObject_NEW(rngobject, &rngtype);
-	if (self == NULL)
-		return NULL;
-	self->distribution = distribution;
-	Py_INCREF(distribution);
-
-	Mixranf(&seed, self->seed);
-	self->position = 0;
-	dist_sample(self->distribution, self->sample, SAMPLE_SIZE);
-	Getranf(self->seed);
-
-
-#ifdef RAN_DEBUG
-	{
-		int i;
-/* Print first few elements of stored sample. */
-for(i = 0; i < 6; i++) {
-   fprintf(stderr,"sample[%i] = %f\n",i,self->sample[i]);
-}
-	}
-#endif
-	return self;
-}
-
-static void
-rng_dealloc(rngobject *self)
-{
-	Py_DECREF(self->distribution);
-	PyMem_DEL(self);
-}
-
-
-static PyObject *
-rng_getattr(rngobject *self, char *name)
-{
-	return Py_FindMethod(rng_methods, (PyObject *)self, name);
-}
-
-static char rngtype__doc__[] = 
-"Random number generator"
-;
-
-static PyTypeObject rngtype = {
-  PyObject_HEAD_INIT(0)
-    0,				/*ob_size*/
-    "random_number_generator",	/*tp_name*/
-    sizeof(rngobject),		/*tp_basicsize*/
-    0,				/*tp_itemsize*/
-    /* methods */
-    (destructor) rng_dealloc,	/*tp_dealloc*/
-    (printfunc)0,		/*tp_print*/
-    (getattrfunc)rng_getattr,	/*tp_getattr*/
-    (setattrfunc)0,	/*tp_setattr*/
-    (cmpfunc)0,		/*tp_compare*/
-    (reprfunc)0,		/*tp_repr*/
-    0,			/*tp_as_number*/
-    0,		/*tp_as_sequence*/
-    0,		/*tp_as_mapping*/
-    (hashfunc)0,		/*tp_hash*/
-    (ternaryfunc)0,		/*tp_call*/
-    (reprfunc)0,		/*tp_str*/
-    
-    /* Space for future expansion */
-    0L,0L,0L,0L,
-    rngtype__doc__ /* Documentation string */
-};
-
-/* End of code for randomgenerator objects */
-
-/* ---------------------------------------------------------------- */
-
-
-static char RNG_CreateGenerator__doc__[] =
-"CreateGenerator(s, d) returns an independent random number stream generator.\n\
-   s < 0  ==>  Use the default initial seed value.\n\
-   s = 0  ==>  Set a random value for the seed from the system clock.\n\
-   s > 0  ==>  Set seed directly (32 bits only).\n\
-   d (optional): distribution object.\n\
-";
-
-
-static PyObject *
-RNG_CreateGenerator(PyObject *self, PyObject *args)
-{
-	int seed;
-	PyObject *distribution = default_distribution;
-	PyObject *result;
-
-	if (!PyArg_ParseTuple(args, "i|O!", &seed,
-			      &distributiontype, &distribution))
-		return NULL;
-	result = (PyObject *)newrngobject(seed,(distributionobject *)distribution);
-
-	return result;
-}
-
-/* List of methods defined in the module */
-
-static struct PyMethodDef RNG_methods[] = {
-	{"CreateGenerator",	(PyCFunction) RNG_CreateGenerator,	1,	RNG_CreateGenerator__doc__},
-	{"UniformDistribution", (PyCFunction) RNG_UniformDistribution, 1,
-	 RNG_UniformDistribution__doc__},
-	{"NormalDistribution", (PyCFunction) RNG_NormalDistribution, 1,
-	 RNG_NormalDistribution__doc__},
-	{"LogNormalDistribution", (PyCFunction) RNG_LogNormalDistribution, 1,
-	 RNG_LogNormalDistribution__doc__},
-	{"ExponentialDistribution", (PyCFunction) RNG_ExponentialDistribution,
-	 1, RNG_ExponentialDistribution__doc__},
-
-	{NULL,		NULL}		/* sentinel */
-};
-
-
-/* Initialization function for the module (*must* be called initURNG) */
-
-static char RNG_module_documentation[] = 
-"Random number generator: independent random number streams."
-;
-
-DL_EXPORT(void)
-initrng(void)
-{
-	PyObject *m, *d;
-    distributiontype.ob_type = &PyType_Type;
-    rngtype.ob_type = &PyType_Type;
-
-	/* Create the module and add the functions */
-	m = Py_InitModule4("rng", RNG_methods,
-		RNG_module_documentation,
-		(PyObject*)NULL,PYTHON_API_VERSION);
-
-	import_array();
-
-	/* Add some symbolic constants to the module */
-	d = PyModule_GetDict(m);
-	ErrorObject = PyErr_NewException("random_lite.error", NULL, NULL);
-	PyDict_SetItemString(d, "error", ErrorObject);
-	default_distribution = create_default_distribution();
-	PyDict_SetItemString(d, "default_distribution", default_distribution);
-
-	/* Check for errors */
-	if (PyErr_Occurred())
-		Py_FatalError("can't initialize module random_lite");
-}
-
diff --git a/scipy/corelib/setup.py b/scipy/corelib/setup.py
index 7f60234fa..9db958af7 100644
--- a/scipy/corelib/setup.py
+++ b/scipy/corelib/setup.py
@@ -24,18 +24,6 @@ def configuration(parent_package='',top_path=None):
                                   ['fftpack_litemodule.c', 'fftpack.c']]
                          )
 
-    # Configure random_lite
-    config.add_extension('rng',
-                         sources=[join('random_lite',x) for x in \
-                                  ['rngmodule.c','ranf.c','pmath_rng.c']]
-                         )
-
-    config.add_extension('ranlib',
-                         sources=[join('random_lite',x) for x in \
-                                  ['ranlibmodule.c', 'ranlib.c', 'com.c',
-                                   'linpack.c']]
-                         )
-
     config.add_extension('mtrand',
                          sources=[join('mtrand', x) for x in 
                                   ['mtrand.c', 'randomkit.c', 'initarray.c',
diff --git a/scipy/stats/__init__.py b/scipy/stats/__init__.py
index 9ab09d4dc..2fa6484c4 100644
--- a/scipy/stats/__init__.py
+++ b/scipy/stats/__init__.py
@@ -17,6 +17,7 @@ import scipy.linalg as LinearAlgebra
 # some aliases
 ranf = random_sample
 random = random_sample
+sample = random_sample
 
 def rand(*args):
     """rand(d1,...,dn) returns a matrix of the given dimensions
@@ -80,4 +81,4 @@ def multivariate_normal(mean, cov, shape=[]):
        return x
 
 # XXX: should we also bring over mean_var_test() from random_lite.py? It seems 
-# out of place.
\ No newline at end of file
+# out of place.