diff --git a/test/testkiss.py b/test/testkiss.py
index 87edce2..89819ea 100755
--- a/test/testkiss.py
+++ b/test/testkiss.py
@@ -41,9 +41,9 @@ else:
 def dopack(x,cpx=1):
     x = Numeric.reshape( x, ( Numeric.size(x),) )
     
-    print 'packed=[',
-    print ' '.join([str(y) for y in x]),
-    print ']'
+    #print 'packed=[',
+    #print ' '.join([str(y) for y in x]),
+    #print ']'
 
     if cpx:
         s = ''.join( [ struct.pack(fmt*2,c.real,c.imag) for c in x ] )
@@ -80,9 +80,9 @@ def flatten(x):
 def randmat( ndims ):
     dims=[]
     for i in range( ndims ):
-        curdim = int( random.uniform(4,6) )
-        if doreal and i==0:
-            curdim = int(curdim/2)*2 # force even for first dimension of real
+        curdim = int( random.uniform(4,7) )
+        if doreal and i==(ndims-1):
+            curdim = int(curdim/2)*2 # force even last dimension if real
         dims.append( curdim )
     return make_random(dims )
 
@@ -94,9 +94,11 @@ def test_fft(ndims):
         xver = FFT.real_fftnd(x)
     else:
         xver = FFT.fftnd(x)
-    print 'x=',x
-    print 'xver=',xver
+    #print 'x=',x
+    #print 'xver=',xver
     
+    open('/tmp/fftexp.dat','w').write(dopack( flatten(xver) , True ) )
+
     x2=dofft(x)
     err = xver - x2
     errf = flatten(err)
@@ -132,22 +134,20 @@ def dofft(x):
 
     p = popen2.Popen3(cmd )
 
+    open('/tmp/fftin.dat','w').write(dopack( x , iscomp ) )
 
     p.tochild.write( dopack( x , iscomp ) )
     p.tochild.close()
 
     res = dounpack( p.fromchild.read() , 1 )
+    open('/tmp/fftout.dat','w').write(dopack( flatten(res) , True ) )
     if doreal:
         dims[-1] = int( dims[-1]/2 ) + 1
 
     res = scale * res
 
     p.wait()
-    try:
-        return Numeric.reshape(res,dims)
-    finally:
-        print 'len(res)=%d' % len(res)
-        print 'dims=%s' % str(dims)
+    return Numeric.reshape(res,dims)
 
 def main():
     opts,args = getopt.getopt(sys.argv[1:],'r')
@@ -156,11 +156,15 @@ def main():
     global doreal
     doreal = opts.has_key('-r')
 
-    print 'Testing multi-dimensional FFTs'
-    for dim in range(1,4):
+    if doreal:
+        print 'Testing multi-dimensional real FFTs'
+    else:
+        print 'Testing multi-dimensional FFTs'
+
+    for dim in range(1,5):
         test_fft( dim )
 
 if __name__ == "__main__":
-    random.seed(2);
+    #random.seed(2);
     main()
 
diff --git a/tools/fftutil.c b/tools/fftutil.c
index 73bb670..e80aea7 100644
--- a/tools/fftutil.c
+++ b/tools/fftutil.c
@@ -77,12 +77,12 @@ void fft_filend_real(FILE * fin,FILE * fout,int *dims,int ndims,int isinverse)
     for (i=0;i<ndims;++i) 
         dimprod *= dims[i];
     insize = outsize = dimprod;
+    int rdim = dims[ndims-1];
 
     if (isinverse)
-        insize = insize*2*(dims[0]/2+1)/dims[0];
+        insize = insize*2*(rdim/2+1)/rdim;
     else
-        outsize = outsize*2*(dims[0]/2+1)/dims[0];
-    fprintf(stderr,"insize=%d outsize=%d\n",insize,outsize);
+        outsize = outsize*2*(rdim/2+1)/rdim;
 
     ibuf = malloc(insize*sizeof(kiss_fft_scalar));
     obuf = malloc(outsize*sizeof(kiss_fft_scalar));
@@ -144,7 +144,7 @@ int get_dims(char * arg,int * dims)
         if (p0)
             *p0++ = '\0';
         dims[ndims++] = atoi(arg);
-         fprintf(stderr,"dims[%d] = %d\n",ndims-1,dims[ndims-1]); 
+//         fprintf(stderr,"dims[%d] = %d\n",ndims-1,dims[ndims-1]); 
         arg = p0;
     }while (p0);
     return ndims;
diff --git a/tools/kiss_fftnd.c b/tools/kiss_fftnd.c
index 00ec783..d6c9124 100644
--- a/tools/kiss_fftnd.c
+++ b/tools/kiss_fftnd.c
@@ -177,7 +177,6 @@ void kiss_fftnd(kiss_fftnd_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
     for ( k=0; k < st->ndims; ++k) {
         int curdim = st->dims[k];
         int stride = st->dimprod / curdim;
-        fprintf(stderr,"curdim = %d, stride = %d\n",curdim,stride);
 
         for ( i=0 ; i<stride ; ++i ) 
             kiss_fft_stride( st->states[k], bufin+i , bufout+i*curdim, stride );
diff --git a/tools/kiss_fftndr.c b/tools/kiss_fftndr.c
index 7a638e5..10fc643 100644
--- a/tools/kiss_fftndr.c
+++ b/tools/kiss_fftndr.c
@@ -14,18 +14,7 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 
 #include "kiss_fftndr.h"
 #include "_kiss_fft_guts.h"
-
-
-// copied from kiss_fftnd.c
-struct kiss_fftnd_state{
-    int dimprod; 
-    int ndims;
-    int *dims;
-    kiss_fft_cfg *states; /* cfg states for each dimension */
-    kiss_fft_cpx * tmpbuf; /*buffer capable of hold the entire input */
-};
-
-
+#define MAX(x,y) ( ( (x)<(y) )?(y):(x) )
 
 struct kiss_fftndr_state
 {
@@ -44,22 +33,19 @@ static int prod(const int *dims, int ndims)
     return x;
 }
 
-#define CHK fprintf(stderr,"line=%d\t",__LINE__)
 kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem)
 {
     kiss_fftndr_cfg st = NULL;
     size_t nr=0 , nd=0,ntmp=0;
-    int dimReal = dims[0];
-    int dimOther = prod(dims+1,ndims-1);
+    int dimReal = dims[ndims-1];
+    int dimOther = prod(dims,ndims-1);
     
-    (void)kiss_fftr_alloc(dims[0],inverse_fft,NULL,&nr);
-    (void)kiss_fftnd_alloc(dims+1,ndims-1,inverse_fft,NULL,&nd);
+    (void)kiss_fftr_alloc(dimReal,inverse_fft,NULL,&nr);
+    (void)kiss_fftnd_alloc(dims,ndims-1,inverse_fft,NULL,&nd);
     ntmp =
-        dimReal * sizeof(kiss_fft_scalar)   // time buffer for one pass
-        + (dimReal+2) * sizeof(kiss_fft_scalar)  // freq buffer for one pass
+        MAX( 2*dimOther , dimReal+2) * sizeof(kiss_fft_scalar)  // freq buffer for one pass
         + dimOther*(dimReal+2) * sizeof(kiss_fft_scalar);  // large enough to hold entire input in case of in-place
 
-    
     size_t memneeded = sizeof( struct kiss_fftndr_state ) + nr + nd + ntmp;
 
     if (lenmem==NULL) {
@@ -75,12 +61,8 @@ kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void
     
     st->dimReal = dimReal;
     st->dimOther = dimOther;
-    st->cfg_r = kiss_fftr_alloc( dims[0],inverse_fft,st+1,&nr);
-    st->cfg_nd = kiss_fftnd_alloc(dims+1,ndims-1,inverse_fft, ((char*) st->cfg_r)+nr,&nd);
-
-    // tell the N-D complex FFT to stride the input
-    st->cfg_nd->dimprod *= (dims[0]/2+1);
-    
+    st->cfg_r = kiss_fftr_alloc( dimReal,inverse_fft,st+1,&nr);
+    st->cfg_nd = kiss_fftnd_alloc(dims,ndims-1,inverse_fft, ((char*) st->cfg_r)+nr,&nd);
     st->tmpbuf = (char*)st->cfg_nd + nd;
 
     return st;
@@ -92,44 +74,50 @@ void kiss_fftndr(kiss_fftndr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx
     int dimReal = st->dimReal;
     int dimOther = st->dimOther;
 
-    kiss_fft_scalar * t1 = (kiss_fft_scalar*)st->tmpbuf; // enough to hold dimReal scalars
-    kiss_fft_cpx * t2 = (kiss_fft_cpx*)(t1+dimReal); // enough to hold the entire input (only used if in==out)
+    kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf; 
+    kiss_fft_cpx * tmp2;
+    if (dimReal/2+1 > dimOther)
+        tmp2 = tmp1+dimReal/2+1;
+    else
+        tmp2 = tmp1+dimOther;
 
-    fprintf(stderr,"dimReal=%d\n",dimReal);
-    fprintf(stderr,"dimOther=%d\n",dimOther);
-    fprintf(stderr,"t1=%p\n",t1);
-    fprintf(stderr,"t2=%p\n",t2);
+    // timedata is N0 x N1 x ... x Nk real
 
-    
-    // take the real data, fft it and transpose
+    // take a real chunk of data, fft it and place the output at correct intervals
     for (k1=0;k1<dimOther;++k1) {
-        for (k2=0;k2<dimReal;++k2)
-            t1[k2] = timedata[k1 + k2*dimOther]; // select the appropriate samples into a contiguous buffer
-
-        //then fft the N real samples to create N/2+1 complex points
-        kiss_fftr( st->cfg_r, t1, t2 + k1*(dimReal/2+1) );
-
-        fprintf(stderr,"t1=fft([ " );
-        for (k2=0;k2<dimReal;++k2) 
-            fprintf(stderr,"%f ",t1[k2] );
-        fprintf(stderr,"])\n" );
+        kiss_fftr( st->cfg_r, timedata + k1*dimReal , tmp1 ); // tmp1 now holds dimReal/2+1 complex points
+        for (k2=0;k2<dimReal/2+1;++k2)
+           tmp2[ k2*dimOther+k1 ] = tmp1[k2];
     }
 
-        // fft the remaining dimensions just like the N-D complex case
-        kiss_fftnd(st->cfg_nd, t2,freqdata);
-
-    fprintf(stderr,"out=[ " );
-    for (k1=0;k1<dimOther;++k1) {
-        for (k2=0;k2<dimReal/2+1;++k2) {
-            kiss_fft_cpx c = freqdata[k1*(dimReal/2+1) +k2];
-            fprintf(stderr,"%f%+fi ",c.r , c.i );
-        }
-        fprintf(stderr,"\n" );
+    for (k2=0;k2<dimReal/2+1;++k2) {
+        kiss_fftnd(st->cfg_nd, tmp2+k2*dimOther, tmp1);  // tmp1 now holds dimOther complex points
+        for (k1=0;k1<dimOther;++k1) 
+            freqdata[ k1*(dimReal/2+1) + k2] = tmp1[k1];
     }
-    fprintf(stderr,"]\n" );
 }
 
 void kiss_fftndri(kiss_fftndr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
 {
-    memset(timedata,0,sizeof(kiss_fft_scalar) * st->dimReal * st->dimOther);
+    int k1,k2;
+    int dimReal = st->dimReal;
+    int dimOther = st->dimOther;
+    kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf; 
+    kiss_fft_cpx * tmp2;
+    if (dimReal/2+1 > dimOther)
+        tmp2 = tmp1+dimReal/2+1;
+    else
+        tmp2 = tmp1+dimOther;
+
+    for (k2=0;k2<dimReal/2+1;++k2) {
+        for (k1=0;k1<dimOther;++k1) 
+            tmp1[k1] = freqdata[ k1*(dimReal/2+1) + k2 ];
+        kiss_fftnd(st->cfg_nd, tmp1, tmp2+k2*dimOther);
+    }
+
+    for (k1=0;k1<dimOther;++k1) {
+        for (k2=0;k2<dimReal/2+1;++k2)
+            tmp1[k2] = tmp2[ k2*dimOther+k1 ];
+        kiss_fftri( st->cfg_r,tmp1,timedata + k1*dimReal);
+    }
 }