Add RNG interface and clean up old-interfaces to be separate from newer ones.

9 files changed, 568 insertions, 120 deletions

diff --git a/numpy/core/include/numpy/arrayobject.h b/numpy/core/include/numpy/arrayobject.h
index 8c09823fe..ac8c391cb 100644
--- a/numpy/core/include/numpy/arrayobject.h
+++ b/numpy/core/include/numpy/arrayobject.h
@@ -164,7 +164,6 @@ enum PyArray_TYPES {    PyArray_BOOL=0,
 
 	/* basetype array priority */
 #define PyArray_PRIORITY 0.0
-#define PyArray_BIG_PRIORITY 0.1
 	/* default subtype priority */
 #define PyArray_SUBTYPE_PRIORITY 1.0
 
@@ -940,7 +939,6 @@ typedef struct _arr_descr {
 	PyObject *shape;       /* a tuple */
 } PyArray_ArrayDescr;
 
-
 /*
   The main array object structure. It is recommended to use the macros
   defined below (PyArray_DATA and friends) access fields here, instead
@@ -980,6 +978,8 @@ typedef struct {
         int flags;
 } PyArray_Chunk;
 
+typedef int (PyArray_FinalizeFunc)(PyArrayObject *, PyObject *);
+
 /* Array flags */
 /* For backward's compatibility only */
 #define OWN_DIMENSIONS 0
diff --git a/numpy/core/src/arraymethods.c b/numpy/core/src/arraymethods.c
index c83aea36e..79235e58c 100644
--- a/numpy/core/src/arraymethods.c
+++ b/numpy/core/src/arraymethods.c
@@ -548,14 +548,6 @@ array_wraparray(PyArrayObject *self, PyObject *args)
 	return ret;
 }
 
-/* NO-OP --- just so all subclasses will have one by default. */
-static PyObject *
-array_finalize(PyArrayObject *self, PyObject *args)
-{
-	Py_INCREF(Py_None);
-	return Py_None;
-}
-
 
 static char doc_array_getarray[] = "m.__array__(|dtype) just returns either a new reference to self if dtype is not given or a new array of provided data type if dtype is different from the current dtype of the array.";
 
@@ -1620,9 +1612,6 @@ static PyMethodDef array_methods[] = {
 	/* for subtypes */
 	{"__array__", (PyCFunction)array_getarray, 1, doc_array_getarray},
 	{"__array_wrap__", (PyCFunction)array_wraparray, 1, doc_wraparray},
-	/* default version so it is found... -- only used for subclasses */
-	{"__array_finalize__", (PyCFunction)array_finalize, 1, NULL},
-	
 	
 	/* for the copy module */
         {"__copy__", (PyCFunction)array_copy, 1, doc_copy},	 
diff --git a/numpy/core/src/arrayobject.c b/numpy/core/src/arrayobject.c
index 51f50a08b..ef25af71c 100644
--- a/numpy/core/src/arrayobject.c
+++ b/numpy/core/src/arrayobject.c
@@ -4587,26 +4587,34 @@ PyArray_NewFromDescr(PyTypeObject *subtype, PyArray_Descr *descr, int nd,
 		if (str == NULL) {
 			str = PyString_InternFromString("__array_finalize__");
 		}
-		if (strides != NULL) { /* did not allocate own data 
-					  or funny strides */
-			/* update flags before calling back into
-			   Python */
-			PyArray_UpdateFlags(self, UPDATE_ALL_FLAGS);
-		}
 		func = PyObject_GetAttr((PyObject *)self, str);
-		if (func) {
-			args = PyTuple_New(1);
-			if (obj == NULL) obj=Py_None;
-			Py_INCREF(obj);
-			PyTuple_SET_ITEM(args, 0, obj);
-			res = PyObject_Call(func, args, NULL);
-			Py_DECREF(args);
-			Py_DECREF(func);
-			if (res == NULL) goto fail;
-			else Py_DECREF(res);
-		}
-	}
-
+		if (func && func != Py_None) {
+                        if (strides != NULL) { /* did not allocate own data 
+                                                  or funny strides */
+                                /* update flags before finalize function */
+                                PyArray_UpdateFlags(self, UPDATE_ALL_FLAGS);
+                        }
+                        if PyCObject_Check(func) {
+                                PyArray_FinalizeFunc *cfunc;
+                                cfunc = PyCObject_AsVoidPtr(func);
+                                Py_DECREF(func);
+                                if (cfunc(self, obj) < 0) goto fail;
+                        }
+                        else {
+                                args = PyTuple_New(1);
+                                if (obj == NULL) obj=Py_None;
+                                Py_INCREF(obj);
+                                PyTuple_SET_ITEM(args, 0, obj);
+                                res = PyObject_Call(func, args, NULL);
+                                Py_DECREF(args);
+                                Py_DECREF(func);
+                                if (res == NULL) goto fail;
+                                else Py_DECREF(res);
+                        }
+                }
+                else Py_XDECREF(func);
+        }
+        
 	return (PyObject *)self;
 
  fail:
@@ -5598,6 +5606,14 @@ array_flat_set(PyArrayObject *self, PyObject *val)
 	return retval;
 }
 
+/* If this is None, no function call is made */
+static PyObject *
+array_finalize_get(PyArrayObject *self)
+{
+        Py_INCREF(Py_None);
+        return Py_None;
+}
+
 static PyGetSetDef array_getsetlist[] = {
         {"ndim",
 	 (getter)array_ndim_get,
@@ -5679,6 +5695,10 @@ static PyGetSetDef array_getsetlist[] = {
 	 (getter)array_priority_get,
 	 NULL,
 	 "Array priority"},
+        {"__array_finalize__",
+         (getter)array_finalize_get,
+         NULL,
+         "None"},
 	{NULL, NULL, NULL, NULL},  /* Sentinel */
 };
 
diff --git a/numpy/dft/fftpack.py b/numpy/dft/fftpack.py
index 6310ca071..70e76085c 100644
--- a/numpy/dft/fftpack.py
+++ b/numpy/dft/fftpack.py
@@ -5,26 +5,22 @@ The underlying code for these functions is an f2c translated and modified
 version of the FFTPACK routines.
 
 fft(a, n=None, axis=-1)
-inverse_fft(a, n=None, axis=-1)
-real_fft(a, n=None, axis=-1)
-inverse_real_fft(a, n=None, axis=-1)
-hermite_fft(a, n=None, axis=-1)
-inverse_hermite_fft(a, n=None, axis=-1)
-fftnd(a, s=None, axes=None)
-inverse_fftnd(a, s=None, axes=None)
-real_fftnd(a, s=None, axes=None)
-inverse_real_fftnd(a, s=None, axes=None)
-fft2d(a, s=None, axes=(-2,-1))
-inverse_fft2d(a, s=None, axes=(-2, -1))
-real_fft2d(a, s=None, axes=(-2,-1))
-inverse_real_fft2d(a, s=None, axes=(-2, -1))
+ifft(a, n=None, axis=-1)
+refft(a, n=None, axis=-1)
+irefft(a, n=None, axis=-1)
+hfft(a, n=None, axis=-1)
+ihfft(a, n=None, axis=-1)
+fftn(a, s=None, axes=None)
+ifftn(a, s=None, axes=None)
+refftn(a, s=None, axes=None)
+irefftn(a, s=None, axes=None)
+fft2(a, s=None, axes=(-2,-1))
+ifft2(a, s=None, axes=(-2, -1))
+refft2(a, s=None, axes=(-2,-1))
+irefft2(a, s=None, axes=(-2, -1))
 """
-__all__ = ['fft','inverse_fft', 'ifft', 'real_fft', 'refft',
-           'inverse_real_fft', 'irefft', 'hfft', 'ihfft', 'refftn',
-           'irefftn', 'refft2', 'irefft2', 'fft2', 'ifft2',
-           'hermite_fft','inverse_hermite_fft','fftnd','inverse_fftnd',
-           'fft2d', 'inverse_fft2d', 'real_fftnd', 'real_fft2d',
-           'inverse_real_fftnd', 'inverse_real_fft2d','fftn','ifftn']
+__all__ = ['fft','ifft', 'refft', 'irefft', 'hfft', 'ihfft', 'refftn',
+           'irefftn', 'refft2', 'irefft2', 'fft2', 'ifft2', 'fftn', 'ifftn']
 
 from numpy.core import asarray, zeros, swapaxes, shape, Complex, conjugate, \
      Float, take
@@ -90,8 +86,8 @@ def fft(a, n=None, axis=-1):
     return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf, _fft_cache)
 
 
-def inverse_fft(a, n=None, axis=-1):
-    """inverse_fft(a, n=None, axis=-1)
+def ifft(a, n=None, axis=-1):
+    """ifft(a, n=None, axis=-1)
 
     Will return the n point inverse discrete Fourier transform of a.  n
     defaults to the length of a. If n is larger than a, then a will be
@@ -101,7 +97,7 @@ def inverse_fft(a, n=None, axis=-1):
     The input array is expected to be packed the same way as the output of
     fft, as discussed in it's documentation.
 
-    This is the inverse of fft: inverse_fft(fft(a)) == a within numerical
+    This is the inverse of fft: ifft(fft(a)) == a within numerical
     accuracy.
 
     This is most efficient for n a power of two. This also stores a cache of
@@ -115,8 +111,8 @@ def inverse_fft(a, n=None, axis=-1):
     return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftb, _fft_cache) / n
 
 
-def real_fft(a, n=None, axis=-1):
-    """real_fft(a, n=None, axis=-1)
+def refft(a, n=None, axis=-1):
+    """refft(a, n=None, axis=-1)
 
     Will return the n point discrete Fourier transform of the real valued
     array a. n defaults to the length of a. n is the length of the input, not
@@ -136,8 +132,8 @@ def real_fft(a, n=None, axis=-1):
     return _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftf, _real_fft_cache)
 
 
-def inverse_real_fft(a, n=None, axis=-1):
-    """inverse_real_fft(a, n=None, axis=-1)
+def irefft(a, n=None, axis=-1):
+    """irefft(a, n=None, axis=-1)
 
     Will return the real valued n point inverse discrete Fourier transform of
     a, where a contains the nonnegative frequency terms of a Hermite-symmetric
@@ -148,10 +144,10 @@ def inverse_real_fft(a, n=None, axis=-1):
     If you specify an n such that a must be zero-padded or truncated, the
     extra/removed values will be added/removed at high frequencies. One can
     thus resample a series to m points via Fourier interpolation by: a_resamp
-    = inverse_real_fft(real_fft(a), m).
+    = irefft(refft(a), m).
 
-    This is the inverse of real_fft:
-    inverse_real_fft(real_fft(a), len(a)) == a
+    This is the inverse of refft:
+    irefft(refft(a), len(a)) == a
     within numerical accuracy."""
 
     a = asarray(a).astype(complex)
@@ -161,40 +157,40 @@ def inverse_real_fft(a, n=None, axis=-1):
                     _real_fft_cache) / n
 
 
-def hermite_fft(a, n=None, axis=-1):
-    """hermite_fft(a, n=None, axis=-1)
-    inverse_hermite_fft(a, n=None, axis=-1)
+def hfft(a, n=None, axis=-1):
+    """hfft(a, n=None, axis=-1)
+    ihfft(a, n=None, axis=-1)
 
-    These are a pair analogous to real_fft/inverse_real_fft, but for the
+    These are a pair analogous to refft/irefft, but for the
     opposite case: here the signal is real in the frequency domain and has
     Hermite symmetry in the time domain. So here it's hermite_fft for which
     you must supply the length of the result if it is to be odd.
 
-    inverse_hermite_fft(hermite_fft(a), len(a)) == a
+    ihfft(hfft(a), len(a)) == a
     within numerical accuracy."""
 
     a = asarray(a).astype(Complex)
     if n == None:
         n = (shape(a)[axis] - 1) * 2
-    return inverse_real_fft(conjugate(a), n, axis) * n
+    return irefft(conjugate(a), n, axis) * n
 
 
-def inverse_hermite_fft(a, n=None, axis=-1):
-    """hermite_fft(a, n=None, axis=-1)
-    inverse_hermite_fft(a, n=None, axis=-1)
+def ihfft(a, n=None, axis=-1):
+    """hfft(a, n=None, axis=-1)
+    ihfft(a, n=None, axis=-1)
 
-    These are a pair analogous to real_fft/inverse_real_fft, but for the
+    These are a pair analogous to refft/irefft, but for the
     opposite case: here the signal is real in the frequency domain and has
-    Hermite symmetry in the time domain. So here it's hermite_fft for which
+    Hermite symmetry in the time domain. So here it's hfft for which
     you must supply the length of the result if it is to be odd.
 
-    inverse_hermite_fft(hermite_fft(a), len(a)) == a
+    ihfft(hfft(a), len(a)) == a
     within numerical accuracy."""
 
     a = asarray(a).astype(Float)
     if n == None:
         n = shape(a)[axis]
-    return conjugate(real_fft(a, n, axis))/n
+    return conjugate(refft(a, n, axis))/n
 
 
 def _cook_nd_args(a, s=None, axes=None, invreal=0):
@@ -226,8 +222,8 @@ def _raw_fftnd(a, s=None, axes=None, function=fft):
     return a
 
 
-def fftnd(a, s=None, axes=None):
-    """fftnd(a, s=None, axes=None)
+def fftn(a, s=None, axes=None):
+    """fftn(a, s=None, axes=None)
 
     The n-dimensional fft of a. s is a sequence giving the shape of the input
     an result along the transformed axes, as n for fft. Results are packed
@@ -243,16 +239,16 @@ def fftnd(a, s=None, axes=None):
 
     return _raw_fftnd(a,s,axes,fft)
 
-def inverse_fftnd(a, s=None, axes=None):
-    """inverse_fftnd(a, s=None, axes=None)
+def ifftn(a, s=None, axes=None):
+    """ifftn(a, s=None, axes=None)
 
-    The inverse of fftnd."""
+    The inverse of fftn."""
 
-    return _raw_fftnd(a, s, axes, inverse_fft)
+    return _raw_fftnd(a, s, axes, ifft)
 
 
-def fft2d(a, s=None, axes=(-2,-1)):
-    """fft2d(a, s=None, axes=(-2,-1))
+def fft2(a, s=None, axes=(-2,-1)):
+    """fft2(a, s=None, axes=(-2,-1))
 
     The 2d fft of a. This is really just fftnd with different default
     behavior."""
@@ -260,17 +256,17 @@ def fft2d(a, s=None, axes=(-2,-1)):
     return _raw_fftnd(a,s,axes,fft)
 
 
-def inverse_fft2d(a, s=None, axes=(-2,-1)):
-    """inverse_fft2d(a, s=None, axes=(-2, -1))
+def ifft2(a, s=None, axes=(-2,-1)):
+    """ifft2(a, s=None, axes=(-2, -1))
 
     The inverse of fft2d. This is really just inverse_fftnd with different
     default behavior."""
 
-    return _raw_fftnd(a, s, axes, inverse_fft)
+    return _raw_fftnd(a, s, axes, ifft)
 
 
-def real_fftnd(a, s=None, axes=None):
-    """real_fftnd(a, s=None, axes=None)
+def refftn(a, s=None, axes=None):
+    """refftn(a, s=None, axes=None)
 
     The n-dimensional discrete Fourier transform of a real array a. A real
     transform as real_fft is performed along the axis specified by the last
@@ -279,24 +275,24 @@ def real_fftnd(a, s=None, axes=None):
 
     a = asarray(a).astype(Float)
     s, axes = _cook_nd_args(a, s, axes)
-    a = real_fft(a, s[-1], axes[-1])
+    a = refft(a, s[-1], axes[-1])
     for ii in range(len(axes)-1):
         a = fft(a, s[ii], axes[ii])
     return a
 
-def real_fft2d(a, s=None, axes=(-2,-1)):
-    """real_fft2d(a, s=None, axes=(-2,-1))
+def refft2(a, s=None, axes=(-2,-1)):
+    """refft2(a, s=None, axes=(-2,-1))
 
-    The 2d fft of the real valued array a. This is really just real_fftnd with
+    The 2d fft of the real valued array a. This is really just refftn with
     different default behavior."""
 
     return real_fftnd(a, s, axes)
 
 
-def inverse_real_fftnd(a, s=None, axes=None):
-    """inverse_real_fftnd(a, s=None, axes=None)
+def irefftn(a, s=None, axes=None):
+    """irefftn(a, s=None, axes=None)
 
-    The inverse of real_fftnd. The transform implemented in inverse_fft is
+    The inverse of refftn. The transform implemented in inverse_fft is
     applied along all axes but the last, then the transform implemented in
     inverse_real_fft is performed along the last axis. As with
     inverse_real_fft, the length of the result along that axis must be
@@ -305,31 +301,15 @@ def inverse_real_fftnd(a, s=None, axes=None):
     a = asarray(a).astype(Complex)
     s, axes = _cook_nd_args(a, s, axes, invreal=1)
     for ii in range(len(axes)-1):
-        a = inverse_fft(a, s[ii], axes[ii])
-    a = inverse_real_fft(a, s[-1], axes[-1])
+        a = ifft(a, s[ii], axes[ii])
+    a = irefft(a, s[-1], axes[-1])
     return a
 
+def irefft2(a, s=None, axes=(-2,-1)):
+    """irefft2(a, s=None, axes=(-2, -1))
 
-def inverse_real_fft2d(a, s=None, axes=(-2,-1)):
-    """inverse_real_fft2d(a, s=None, axes=(-2, -1))
-
-    The inverse of real_fft2d. This is really just inverse_real_fftnd with
+    The inverse of refft2. This is really just irefftn with
     different default behavior."""
 
-    return inverse_real_fftnd(a, s, axes)
-
-ifft = inverse_fft
-refft = real_fft
-irefft = inverse_real_fft
-hfft = hermite_fft
-ihfft = inverse_hermite_fft
-
-fftn = fftnd
-ifftn = inverse_fftnd
-refftn = real_fftnd
-irefftn = inverse_real_fftnd
+    return irefftn(a, s, axes)
 
-fft2 = fft2d
-ifft2 = inverse_fft2d
-refft2 = real_fft2d
-irefft2 = inverse_real_fft2d
diff --git a/numpy/dft/old.py b/numpy/dft/old.py
new file mode 100644
index 000000000..2ed5321d9
--- /dev/null
+++ b/numpy/dft/old.py
@@ -0,0 +1,19 @@
+
+__all__ = ['fft', 'fft2d', 'fftnd', 'hermite_fft', 'inverse_fft', 'inverse_fft2d',
+           'inverse_fftnd', 'inverse_hermite_fft', 'inverse_real_fft', 'inverse_real_fft2d',
+           'inverse_real_fftnd', 'real_fft', 'real_fft2d', 'real_fftnd']
+
+from fftpack import fft
+from fftpack import fft2 as fft2d
+from fftpack import fftn as fftnd
+from fftpack import hfft as hermite_fft
+from fftpack import ifft as inverse_fft
+from fftpack import ifft2 as inverse_fft2d
+from fftpack import ifftn as inverse_fftnd
+from fftpack import ihfft as inverse_hermite_fft
+from fftpack import irefft as inverse_real_fft
+from fftpack import irefft2 as inverse_real_fft2d
+from fftpack import irefftn as inverse_real_fftnd
+from fftpack import refft as real_fft
+from fftpack import refft2 as real_fft2d
+from fftpack import refftn as real_fftnd
diff --git a/numpy/lib/convertcode.py b/numpy/lib/convertcode.py
index 67a98cb1d..d45bd2fe3 100644
--- a/numpy/lib/convertcode.py
+++ b/numpy/lib/convertcode.py
@@ -98,13 +98,14 @@ def fromstr(filestr):
                                           'numpy.core.umath')
     filestr, fromall1 = changeimports(filestr, 'Precision', 'numpy')
     filestr, fromall2 = changeimports(filestr, 'numerix', 'numpy')
-    filestr, fromall3 = changeimports(filestr, 'numpy_base', 'numpy')
+    filestr, fromall3 = changeimports(filestr, 'scipy_base', 'numpy')
     filestr, fromall3 = changeimports(filestr, 'MLab', 'numpy.lib.mlab')
     filestr, fromall3 = changeimports(filestr, 'LinearAlgebra',
                                       'numpy.linalg.old')
-    filestr, fromall3 = changeimports(filestr, 'RNG', 'numpy.random')
-    filestr, fromall3 = changeimports(filestr, 'RandomArray', 'numpy.random')
-    filestr, fromall3 = changeimports(filestr, 'FFT', 'numpy.dft')
+    filestr, fromall3 = changeimports(filestr, 'RNG', 'numpy.random.oldrng')
+    filestr, fromall3 = changeimports(filestr, 'RNG.Statistics', 'numpy.random.oldrngstats')
+    filestr, fromall3 = changeimports(filestr, 'RandomArray', 'numpy.random.oldrandomarray')
+    filestr, fromall3 = changeimports(filestr, 'FFT', 'numpy.dft.oldfft')
     filestr, fromall3 = changeimports(filestr, 'MA', 'numpy.core.ma')
     fromall = fromall1 or fromall2 or fromall3
     filestr = replaceattr(filestr)
diff --git a/numpy/random/old.py b/numpy/random/old.py
new file mode 100644
index 000000000..338779ac8
--- /dev/null
+++ b/numpy/random/old.py
@@ -0,0 +1,268 @@
+
+__all__ = ['ArgumentError','F','beta','binomial','chi_square', 'exponential', 'gamma', 'get_seed',
+           'mean_var_test', 'multinomial', 'multivariate_normal', 'negative_binomial',
+           'noncentral_F', 'noncentral_chi_square', 'normal', 'permutation', 'poisson', 'randint',
+           'random', 'random_integers', 'seed', 'standard_normal', 'uniform']
+
+ArgumentError = ValueError
+
+import numpy.random.mtrand as mt
+import numpy as Numeric
+
+from types import IntType
+
+def seed(x=0, y=0):
+    if (x == 0 or y == 0):
+        mt.seed()
+    else:
+        mt.seed((x,y))        
+
+def get_seed():
+    raise NotImplementedError, \
+          "If you want to save the state of the random number generator.\n"\
+          "Then you should use obj = numpy.random.get_state() followed by.\n"\
+          "numpy.random.set_state(obj)."
+
+def random(shape=[]):
+    "random(n) or random([n, m, ...]) returns array of random numbers"
+    if shape == []:
+        shape = None
+    return mt.random_sample(shape)
+
+def uniform(minimum, maximum, shape=[]):
+    """uniform(minimum, maximum, shape=[]) returns array of given shape of random reals
+    in given range"""
+    if shape == []:
+        shape = None
+    return mt.uniform(minimum, maximum, shape)
+
+def randint(minimum, maximum=None, shape=[]):
+    """randint(min, max, shape=[]) = random integers >=min, < max
+    If max not given, random integers >= 0, <min"""
+    if not isinstance(minimum, IntType):
+        raise ArgumentError, "randint requires first argument integer"
+    if maximum is None:
+        maximum = minimum
+        minimum = 0
+    if not isinstance(maximum, IntType):
+        raise ArgumentError, "randint requires second argument integer"
+    a = ((maximum-minimum)* random(shape))
+    if isinstance(a, Numeric.ArrayType):
+        return minimum + a.astype(Numeric.Int)
+    else:
+        return minimum + int(a)
+
+def random_integers(maximum, minimum=1, shape=[]):
+    """random_integers(max, min=1, shape=[]) = random integers in range min-max inclusive"""
+    return randint(minimum, maximum+1, shape)
+
+def permutation(n):
+    "permutation(n) = a permutation of indices range(n)"
+    return mt.permutation(n)
+
+def standard_normal(shape=[]):
+    """standard_normal(n) or standard_normal([n, m, ...]) returns array of
+           random numbers normally distributed with mean 0 and standard
+           deviation 1"""
+    if shape == []:
+        shape = None
+    return mt.standard_normal(shape)
+
+def normal(mean, std, shape=[]):
+        """normal(mean, std, n) or normal(mean, std, [n, m, ...]) returns
+           array of random numbers randomly distributed with specified mean and
+           standard deviation"""
+        if shape == []:
+            shape = None
+        return mt.normal(mean, std, shape)
+
+def multivariate_normal(mean, cov, shape=[]):
+       """multivariate_normal(mean, cov) or multivariate_normal(mean, cov, [m, n, ...])
+          returns an array containing multivariate normally distributed random numbers
+          with specified mean and covariance.
+
+          mean must be a 1 dimensional array. cov must be a square two dimensional
+          array with the same number of rows and columns as mean has elements.
+
+          The first form returns a single 1-D array containing a multivariate
+          normal.
+
+          The second form returns an array of shape (m, n, ..., cov.shape[0]).
+          In this case, output[i,j,...,:] is a 1-D array containing a multivariate
+          normal."""
+       if shape == []:
+           shape = None
+       return mt.multivariate_normal(mean, cov, shape)
+
+def exponential(mean, shape=[]):
+    """exponential(mean, n) or exponential(mean, [n, m, ...]) returns array
+      of random numbers exponentially distributed with specified mean"""
+    if shape == []:
+        shape = None
+    return mt.exponential(mean, shape)
+
+def beta(a, b, shape=[]):
+    """beta(a, b) or beta(a, b, [n, m, ...]) returns array of beta distributed random numbers."""
+    if shape == []:
+        shape = None
+    return mt.beta(a, b, shape)
+
+def gamma(a, r, shape=[]):
+    """gamma(a, r) or gamma(a, r, [n, m, ...]) returns array of gamma distributed random numbers."""
+    if shape == []:
+        shape = None
+    return mt.gamma(a, r, shape)
+
+def F(dfn, dfd, shape=[]):
+    """F(dfn, dfd) or F(dfn, dfd, [n, m, ...]) returns array of F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator."""
+    if shape == []:
+        shape == None
+    return mt.f(dfn, dfd, shape)
+
+def noncentral_F(dfn, dfd, nconc, shape=[]):
+    """noncentral_F(dfn, dfd, nonc) or noncentral_F(dfn, dfd, nonc, [n, m, ...]) returns array of noncentral F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator, and noncentrality parameter nconc."""
+    if shape == []:
+        shape = None
+    return mt.noncentral_f(dfn, dfd, nconc, shape)
+
+def chi_square(df, shape=[]):
+    """chi_square(df) or chi_square(df, [n, m, ...]) returns array of chi squared distributed random numbers with df degrees of freedom."""
+    if shape == []:
+        shape = None
+    return mt.chisquare(df, shape)
+
+def noncentral_chi_square(df, nconc, shape=[]):
+    """noncentral_chi_square(df, nconc) or chi_square(df, nconc, [n, m, ...]) returns array of noncentral chi squared distributed random numbers with df degrees of freedom and noncentrality parameter."""
+    if shape == []:
+        shape = None
+    return mt.noncentral_chisquare(df, nconc, shape)
+
+def binomial(trials, p, shape=[]):
+    """binomial(trials, p) or binomial(trials, p, [n, m, ...]) returns array of binomially distributed random integers.
+
+           trials is the number of trials in the binomial distribution.
+           p is the probability of an event in each trial of the binomial distribution."""
+    if shape == []:
+        shape = None
+    return mt.binomial(trials, p, shape)
+
+def negative_binomial(trials, p, shape=[]):
+    """negative_binomial(trials, p) or negative_binomial(trials, p, [n, m, ...]) returns
+           array of negative binomially distributed random integers.
+
+           trials is the number of trials in the negative binomial distribution.
+           p is the probability of an event in each trial of the negative binomial distribution."""
+    if shape == []:
+        shape = None
+    return mt.negative_binomial(trials, p, shape)
+
+def multinomial(trials, probs, shape=[]):
+    """multinomial(trials, probs) or multinomial(trials, probs, [n, m, ...]) returns
+           array of multinomial distributed integer vectors.
+
+           trials is the number of trials in each multinomial distribution.
+           probs is a one dimensional array. There are len(prob)+1 events.
+           prob[i] is the probability of the i-th event, 0<=i<len(prob).
+           The probability of event len(prob) is 1.-Numeric.sum(prob).
+
+       The first form returns a single 1-D array containing one multinomially
+           distributed vector.
+
+           The second form returns an array of shape (m, n, ..., len(probs)).
+           In this case, output[i,j,...,:] is a 1-D array containing a multinomially
+           distributed integer 1-D array."""
+    if shape == []:
+        shape = None
+    return mt.multinomial(trials, probs, shape)
+
+def poisson(mean, shape=[]):
+    """poisson(mean) or poisson(mean, [n, m, ...]) returns array of poisson
+           distributed random integers with specified mean."""
+    if shape == []:
+        shape = None
+    return mt.poisson(mean, shape)
+
+
+def mean_var_test(x, type, mean, var, skew=[]):
+    n = len(x) * 1.0
+    x_mean = Numeric.sum(x)/n
+    x_minus_mean = x - x_mean
+    x_var = Numeric.sum(x_minus_mean*x_minus_mean)/(n-1.0)
+    print "\nAverage of ", len(x), type
+    print "(should be about ", mean, "):", x_mean
+    print "Variance of those random numbers (should be about ", var, "):", x_var
+    if skew != []:
+       x_skew = (Numeric.sum(x_minus_mean*x_minus_mean*x_minus_mean)/9998.)/x_var**(3./2.)
+       print "Skewness of those random numbers (should be about ", skew, "):", x_skew
+
+def test():
+    from types import *
+    obj = mt.get_state()
+    mt.set_state(obj)
+    obj2 = mt.get_state()
+    if (obj2[1] - obj[1]).any():
+        raise SystemExit, "Failed seed test."
+    print "First random number is", random()
+    print "Average of 10000 random numbers is", Numeric.sum(random(10000))/10000.
+    x = random([10,1000])
+    if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000:
+        raise SystemExit, "random returned wrong shape"
+    x.shape = (10000,)
+    print "Average of 100 by 100 random numbers is", Numeric.sum(x)/10000.
+    y = uniform(0.5,0.6, (1000,10))
+    if len(y.shape) !=2 or y.shape[0] != 1000 or y.shape[1] != 10:
+        raise SystemExit, "uniform returned wrong shape"
+    y.shape = (10000,)
+    if Numeric.minimum.reduce(y) <= 0.5 or Numeric.maximum.reduce(y) >= 0.6:
+        raise SystemExit, "uniform returned out of desired range"
+    print "randint(1, 10, shape=[50])"
+    print randint(1, 10, shape=[50])
+    print "permutation(10)", permutation(10)
+    print "randint(3,9)", randint(3,9)
+    print "random_integers(10, shape=[20])"
+    print random_integers(10, shape=[20])
+    s = 3.0
+    x = normal(2.0, s, [10, 1000])
+    if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000:
+        raise SystemExit, "standard_normal returned wrong shape"
+    x.shape = (10000,)
+    mean_var_test(x, "normally distributed numbers with mean 2 and variance %f"%(s**2,), 2, s**2, 0)
+    x = exponential(3, 10000)
+    mean_var_test(x, "random numbers exponentially distributed with mean %f"%(s,), s, s**2, 2)
+    x = multivariate_normal(Numeric.array([10,20]), Numeric.array(([1,2],[2,4])))
+    print "\nA multivariate normal", x
+    if x.shape != (2,): raise SystemExit, "multivariate_normal returned wrong shape"
+    x = multivariate_normal(Numeric.array([10,20]), Numeric.array([[1,2],[2,4]]), [4,3])
+    print "A 4x3x2 array containing multivariate normals"
+    print x
+    if x.shape != (4,3,2): raise SystemExit, "multivariate_normal returned wrong shape"
+    x = multivariate_normal(Numeric.array([-100,0,100]), Numeric.array([[3,2,1],[2,2,1],[1,1,1]]), 10000)
+    x_mean = Numeric.sum(x)/10000.
+    print "Average of 10000 multivariate normals with mean [-100,0,100]"
+    print x_mean
+    x_minus_mean = x - x_mean
+    print "Estimated covariance of 10000 multivariate normals with covariance [[3,2,1],[2,2,1],[1,1,1]]"
+    print Numeric.dot(Numeric.transpose(x_minus_mean),x_minus_mean)/9999.
+    x = beta(5.0, 10.0, 10000)
+    mean_var_test(x, "beta(5.,10.) random numbers", 0.333, 0.014)
+    x = gamma(.01, 2., 10000)
+    mean_var_test(x, "gamma(.01,2.) random numbers", 2*100, 2*100*100)
+    x = chi_square(11., 10000)
+    mean_var_test(x, "chi squared random numbers with 11 degrees of freedom", 11, 22, 2*Numeric.sqrt(2./11.))
+    x = F(5., 10., 10000)
+    mean_var_test(x, "F random numbers with 5 and 10 degrees of freedom", 1.25, 1.35)
+    x = poisson(50., 10000)
+    mean_var_test(x, "poisson random numbers with mean 50", 50, 50, 0.14)
+    print "\nEach element is the result of 16 binomial trials with probability 0.5:"
+    print binomial(16, 0.5, 16)
+    print "\nEach element is the result of 16 negative binomial trials with probability 0.5:"
+    print negative_binomial(16, 0.5, [16,])
+    print "\nEach row is the result of 16 multinomial trials with probabilities [0.1, 0.5, 0.1 0.3]:"
+    x = multinomial(16, [0.1, 0.5, 0.1], 8)
+    print x
+    print "Mean = ", Numeric.sum(x)/8.
+
+if __name__ == '__main__':
+    test()
+
+
diff --git a/numpy/random/oldrng.py b/numpy/random/oldrng.py
new file mode 100644
index 000000000..1e6b2fbd6
--- /dev/null
+++ b/numpy/random/oldrng.py
@@ -0,0 +1,136 @@
+# This module re-creates the RNG interface from Numeric
+# Replace import RNG with import numpy.random.oldrng as RNG
+#
+# It is for backwards compatibility only. 
+
+
+__all__ = ['CreateGenerator','ExponentialDistribution','LogNormalDistribution','NormalDistribution',
+           'UniformDistribution', 'error', 'default_distribution', 'random_sample', 'ranf',
+           'standard_generator']
+
+import numpy.random.mtrand as mt
+import math
+
+class error(Exception):
+    pass
+
+class Distribution(object):
+    def __init__(self, meth, *args):
+        self._meth = meth
+        self._args = args
+
+    def density(self,x):
+        raise NotImplementedError
+    
+    def __call__(self, x):
+        return self.density(x)
+
+    def _onesample(self, rng):
+        return getattr(rng, self._meth)(*self._args)
+
+    def _sample(self, rng, n):
+        kwds = {'size' : n}
+        return getattr(rng, self._meth)(*self._args, **kwds)
+
+
+class ExponentialDistribution(Distribution):
+    def __init__(self, lambda_):
+        if (lambda_ <= 0):
+            raise error, "parameter must be positive"  
+        Distribution.__init__(self, 'exponential', lambda_)
+
+    def density(x):
+        if x < 0:
+            return 0.0
+        else:
+            lambda_ = self._args[0]
+            return lambda_*exp(-lambda_*x)
+
+class LogNormalDistribution(Distribution):
+    def __init__(self, m, s):
+        m = float(m)
+        s = float(s)
+        if (s <= 0):
+            raise error, "standard deviation must be positive"
+        Distribution.__init__(self, 'lognormal', m, s)
+        sn = math.log(1.0+s*s/(m*m));
+        self._mn = math.log(m)-0.5*sn
+        self._sn = math.sqrt(sn)
+        self._fac = 1.0/math.sqrt(2*math.pi)/self._sn
+
+    def density(x):
+        m,s = self._args
+        y = (math.log(x)-self._mn)/self._sn
+        return self._fac*exp(-0.5*y*y)/x
+
+
+class NormalDistribution(Distribution):
+    def __init__(self, m, s):
+        m = float(m)
+        s = float(s)
+        if (s <= 0):
+            raise error, "standard deviation must be positive"
+        Distribution.__init__(self, 'normal', m, s)
+        self._fac = 1.0/math.sqrt(2*math.pi)/s
+
+    def density(x):
+        m,s = self._args
+        y = (x-m)/s
+        return self._fac*exp(-0.5*y*y)
+
+class UniformDistribution(Distribution):
+    def __init__(self, a, b):
+        a = float(a)
+        b = float(b)
+        width = b-a
+        if (width <=0):
+            raise error, "width of uniform distribution must be > 0"
+        Distribution.__init__(self, 'uniform', a, b)
+        self._fac = 1.0/width
+
+    def density(x):
+        a, b = self._args
+        if (x < a) or (x >= b):
+            return 0.0
+        else:
+            return self._fac
+
+default_distribution = UniformDistribution(0.0,1.0)
+
+class CreateGenerator(object):
+    def __init__(self, seed, dist=None):
+        if seed <= 0:
+            self._rng = mt.RandomState()
+        elif seed > 0:
+            self._rng = mt.RandomState(seed)
+        if dist is None:
+            dist = default_distribution
+        if not isinstance(dist, Distribution):
+            raise error, "Not a distribution object"
+        self._dist = dist
+
+    def ranf(self):
+        return self._dist._onesample(self._rng)
+
+    def sample(self, n):
+        return self._dist._sample(self._rng, n)
+
+    
+standard_generator = CreateGenerator(-1)
+
+def ranf():
+    "ranf() = a random number from the standard generator."
+    return standard_generator.ranf()
+
+def random_sample(*n):
+    """random_sample(n) = array of n random numbers;
+
+    random_sample(n1, n2, ...)= random array of shape (n1, n2, ..)"""
+    
+    if not n:
+        return standard_generator.ranf()
+    m = 1
+    for i in n:
+        m = m * i
+    return standard_generator.sample(m).reshape(*n)
+    
diff --git a/numpy/random/oldrngstats.py b/numpy/random/oldrngstats.py
new file mode 100644
index 000000000..b795eee09
--- /dev/null
+++ b/numpy/random/oldrngstats.py
@@ -0,0 +1,35 @@
+
+__all__ = ['average', 'histogram', 'standardDeviation', 'variance']
+
+import numpy as Numeric
+
+def average(data):
+    data = Numeric.array(data)
+    return Numeric.add.reduce(data)/len(data)
+
+def variance(data):
+    data = Numeric.array(data)
+    return Numeric.add.reduce((data-average(data))**2)/(len(data)-1)
+
+def standardDeviation(data):
+    data = Numeric.array(data)
+    return Numeric.sqrt(variance(data))
+
+def histogram(data, nbins, range = None):
+    data = Numeric.array(data, Numeric.Float)
+    if range is None:
+        min = Numeric.minimum.reduce(data)
+        max = Numeric.maximum.reduce(data)
+    else:
+        min, max = range
+        data = Numeric.repeat(data,
+                              Numeric.logical_and(Numeric.less_equal(data, max),
+                                                  Numeric.greater_equal(data,
+                                                                        min)))
+    bin_width = (max-min)/nbins
+    data = Numeric.floor((data - min)/bin_width).astype(Numeric.Int)
+    histo = Numeric.add.reduce(Numeric.equal(
+        Numeric.arange(nbins)[:,Numeric.NewAxis], data), -1)
+    histo[-1] = histo[-1] + Numeric.add.reduce(Numeric.equal(nbins, data))
+    bins = min + bin_width*(Numeric.arange(nbins)+0.5)
+    return Numeric.transpose(Numeric.array([bins, histo]))