Original Version

1 files changed, 1253 insertions, 0 deletions

diff --git a/src/c/signalProcessing/ifft/difftmx.c b/src/c/signalProcessing/ifft/difftmx.c
new file mode 100644
index 0000000..a2cce1b
--- /dev/null
+++ b/src/c/signalProcessing/ifft/difftmx.c
@@ -0,0 +1,1253 @@
+/*
+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ *  Copyright (C) 2008-2008 - INRIA - Allan SIMON
+ *
+ *  This file must be used under the terms of the CeCILL.
+ *  This source file is licensed as described in the file COPYING, which
+ *  you should have received as part of this distribution.  The terms
+ *  are also available at
+ *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+
+#include <stdlib.h>
+#include <math.h>
+#include "max.h"
+#include "min.h"
+#include "ifft_internal.h"
+
+
+/*c'est moche je sais */
+
+static double* a ;
+static double* b ;
+static int     ntot  ;
+static int     n     ;
+static int     nspan ;
+static int     isn   ;
+static int     m     ;
+static int     kt    ;
+static double*    wt    ;
+static double*    ck    ;
+static double*    bt    ;
+static double*    sk    ;
+static int*    np    ;
+static int*    nfac  ;
+static int     inc   ;
+static int     nt    ;
+static int     ks    ;
+
+static double  rad   ;
+static double  c72   ;
+static double  s72   ;
+static double  s120  ;
+
+static double  aa ;
+static double  ak ;
+static double  akm ;
+static double  akp ;
+static double  aj ;
+static double  ajp ;
+static double  ajm ;
+
+static double  bb ;
+static double  bk ;
+static double  bkm ;
+static double  bkp ;
+static double  bj ;
+static double  bjp ;
+static double  bjm ;
+static double  dr ;
+
+static double  cd ;
+static double  c1 ;
+static double  c2 ;
+static double  c3 ;
+
+static double  sd ;
+static double  s1 ;
+static double  s2 ;
+static double  s3 ;
+
+static int      kspan ;
+static int     nn ;
+static int     jc ;
+static int     klim ;
+static int     jf ;
+static int     maxf ;
+static int     kk ;
+static int     k ;
+static int     k1 ;
+static int     k2 ;
+static int     k3 ;
+static int     k4 ;
+static int     mm ;
+static int     kspnn ;
+static int      i ;
+static int      j ;
+static int      jj;
+
+/* Prototypes */
+
+static   void preliminaryWork (void);
+static    void permute_stage1 (void);
+static    void permute_stage2 (void);
+static    void f4t_150 (void);
+static    void factorOf3Transform (void) ;
+static    void factorOf5Transform (void) ;
+static    void preFOtherTransform (void);
+static    void factorOfOtherTransform (void);
+static    void pre_sqFactor2NormlOrder (void);
+static    void nonSqFactor2NormOrder (void) ;
+static    void detPermutCycles (void);
+static    void reorderMatrix (void ) ;
+
+static    int f4t_170 (void);
+static    int factorTransform (void);
+static    int pre_fOf2Trans (void);
+static    int factorOf2Transform (void);
+static    int factorOf4Transform (void);
+static    int  mulByRotationFactor (void );
+static    int  post_sqFactor2NormlOrder (void);
+static    void  single_sqFactor2NormlOrder (void);
+static    int  multi_sqFactor2NormlOrder (void);
+
+/* End Prototypes */
+
+
+
+int difftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan,
+             int _iIsn, int _iM, int _iKt, double* _pdblWt, double* _pdblCk,
+             double* _pdblBt, double* _pdblSk, int* _piNp, int* _piNfac)
+{
+
+   int retVal = 0 ;
+   int iii = 0 ;
+   a = _pdblA ;
+   b = _pdblB ;
+
+   ntot = _iNtot ;
+   n    = _iN ;
+   nspan= _iNspan ;
+   isn  = _iIsn;
+   m    = _iM ;
+   kt   = _iKt ;
+   wt   = _pdblWt ;
+   ck   = _pdblCk;
+   bt   = _pdblBt;
+   sk   = _pdblSk;
+   np   = _piNp;
+   nfac = _piNfac;
+
+   inc = abs ( isn ) ;
+   nt = inc*ntot ;
+   ks = inc*nspan;
+   rad = atan (1.0);
+
+   c72 = cos (rad/0.6250);
+   s72 = sin (rad/0.6250);
+   s120= sqrt(0.750);
+
+
+
+
+
+
+
+   preliminaryWork () ;
+
+
+   while ( retVal == 0 )
+      {
+
+        retVal = factorTransform ( ) ;
+      }
+
+
+   np[0] = ks ;
+
+   if ( kt != 0)
+    {
+
+     permute_stage1 ( ) ;
+    }
+
+   if ( 2*kt + 1 < m )
+    {
+
+      permute_stage2 ( ) ;
+    }
+
+
+/*  lines under are just for my own conveniance   */
+   for ( iii = 0 ; iii < 3 ; iii++)
+    {
+
+
+
+    }
+
+   _pdblA = a ;
+   _pdblB = b ;
+
+   return 0 ;
+}
+
+/** **************************************
+Sous-Fonctions
+******************************************/
+
+
+
+
+static void preliminaryWork (void)
+{
+
+    int lim ;
+
+
+
+       ak =(double) 1/n ;
+
+       /*scale by 1/n for isn > 0 */
+       for ( j = 1 ;  j <= nt ; j += inc )
+          {
+           a[j-1] *= ak ;
+           b[j-1] *= ak ;
+
+          }
+   kspan = ks ;
+   nn = nt -inc ;
+   jc = ks/n ;
+
+   /* sin , cos values are re-initialized each lim steps */
+
+   lim = 32 ;
+   klim = lim*jc ;
+   i = 0;
+   jf = 0 ;
+
+   maxf = m -kt ;
+   maxf = nfac[maxf-1] ;
+
+
+
+   if ( kt > 0 )
+      maxf = max ( nfac[kt-1] , maxf );
+
+
+}
+
+
+/*40*/
+/* this function is call as many time as dfftbi has determined factor for the size of the input vector
+   each time we call a transform function for each kind of factor , we begin by the smallest
+   factor are stored in nfac
+ */
+
+static int  factorTransform (void)
+{
+
+   int retVal = 42;
+
+   dr = 8 * (double)jc/(double)kspan ;
+   cd = 2 * sin(0.5*dr*rad)*sin(0.5*dr*rad);
+   sd = sin(dr*rad) ;
+   kk = 1 ;
+   i++ ;
+
+
+
+
+switch ( nfac[i-1] )
+   {
+      case 2 :
+         /*transform for factor of 2 (including rotation factor)*/
+
+         retVal = pre_fOf2Trans() ;
+         if ( retVal == 0 )  factorOf2Transform () ;
+
+         break ;
+
+      case 4 :
+         /*transform for factor of 4 */
+         kspnn = kspan ;
+         kspan = kspan >> 2 ; /*kspan /= 4 */
+
+         retVal = factorOf4Transform () ;
+         break ;
+
+      case 3 :
+         /*transform for factor of 3 */
+         k = nfac[i-1] ;
+         kspnn = kspan ;
+         kspan = kspan / k ;
+
+         factorOf3Transform ( ) ;
+         break ;
+
+      case 5 :
+   	/*transform for factor of 5 */
+         k = nfac[i-1] ;
+         kspnn = kspan ;
+         kspan = kspan / k ;
+
+         factorOf5Transform ( ) ;
+         break ;
+
+      default :
+
+         k = nfac[i-1] ;
+         kspnn = kspan ;
+         kspan = kspan / k ;
+
+         if ( nfac[i-1] != jf) preFOtherTransform ( ) ;
+
+         factorOfOtherTransform ( ) ;
+         break ;
+    }
+
+
+
+   if ( retVal == 42 )
+    {
+      if ( i !=  m)  retVal = mulByRotationFactor ( ) ;
+      else           retVal = 1 ;
+    }
+
+   if ( retVal == 1 ) return 1 ; /*goto permute */
+   else 	      return 0 ; /*goto factor_transform => once again*/
+
+}
+
+/* permutation for square factor of n */
+static void permute_stage1 (void)
+{
+
+ int retVal = 1 ;
+
+   pre_sqFactor2NormlOrder () ;
+
+   if ( n == ntot )
+      /*permutation for single-variate transform (optional code)*/
+      while ( retVal == 1)
+        {
+         single_sqFactor2NormlOrder () ;
+         retVal = post_sqFactor2NormlOrder () ;
+        }
+   else
+      /*permutation for multivariate transform*/
+      while ( retVal == 1) retVal = multi_sqFactor2NormlOrder ();
+
+}
+
+static void permute_stage2 (void)
+{
+	kspnn = np[kt] ;
+
+	/*permutation for square-free facotrs of n */
+	nonSqFactor2NormOrder () ;
+
+	/*determine the permutation cycles of length greater than 1*/
+	detPermutCycles ();
+
+	j = k3 + 1;
+	nt -= kspnn ;
+	i = nt - inc + 1 ;
+	while ( nt >= 0 )
+	{
+		reorderMatrix ( ) ;
+
+		j = k3 + 1 ;
+		nt -= kspnn ;
+		i = nt - inc + 1 ;
+         }
+}
+
+/*****************************************
+Sous-Sous-Fonctions
+******************************************/
+
+
+
+
+
+static int pre_fOf2Trans (void)
+{
+	kspan /= 2;
+	k1 = kspan + 2 ;
+	/*50*/
+	do{
+		do{
+			k2 = kk + kspan ;
+			ak = a[k2-1] ;
+			bk = b[k2-1] ;
+
+			a[k2-1] = a[kk-1] - ak;
+			b[k2-1] = b[kk-1] - bk;
+
+			a[kk-1] = a[kk-1] + ak;
+			b[kk-1] = b[kk-1] + bk;
+
+			kk = k2 + kspan ;
+		}while (kk <= nn);
+
+		kk -= nn ;
+	}while (kk <= jc);
+
+
+	if ( kk > kspan )  return 1 ; /*goto350*/
+	else               return 0 ; /*goto60*/
+
+
+}
+
+
+
+static int factorOf2Transform (void)
+{
+	do /*60*/ {/*while ( kk <= jc*2 )*/
+		c1 = 1 - cd ;
+		s1 = sd ;
+		mm = min( k1/2 , klim);
+
+		do/* do 80 */  {/*while ( kk <= mm || ( kk > mm && kk < k2 ))*/
+			do {/*while(kk > k2) */
+				do { /*while (  kk < nt )*/
+					k2 = kk + kspan;
+
+					ak = a[kk-1] - a[k2-1];
+					bk = b[kk-1] - b[k2-1];
+
+					a[kk-1] = a[kk-1] + a[k2-1];
+					b[kk-1] = b[kk-1] + b[k2-1];
+
+					a[k2-1] = c1*ak - s1*bk;
+					b[k2-1] = s1*ak + c1*bk;
+
+					kk = k2 + kspan;
+				}while (  kk < nt );
+
+				k2 = kk - nt;
+				c1 = -c1;
+				kk = k1 - k2;
+
+
+			}while (kk > k2);
+
+			kk += jc;
+
+			if ( kk <= mm ) /* 70 */
+			{
+				ak =  c1 - ( cd*c1+sd*s1) ;
+				s1 += (sd*c1-cd*s1) ;
+				/*c the following three statements compensate for truncation
+				c error.  if rounded arithmetic is used, substitute
+				c c1=ak*/
+				c1 = 0.5/(ak*ak+s1*s1) + 0.5 ;
+				s1 *= c1 ;
+				c1 *= ak ;
+			}
+			else {
+				if ( kk < k2 ) /*90*/ {
+					s1 = dr*rad*((double)(kk-1)/(double)jc);
+					c1 = cos(s1) ;
+					s1 = sin(s1) ;
+					mm = min(k1/2,mm+klim);
+				}
+			}
+
+		} while ( kk <= mm || ( kk > mm && kk < k2 ));
+
+		k1 += (inc+inc) ;
+		kk = (k1-kspan)/2 + jc;
+
+	} while ( kk <= jc*2 );
+
+
+	return 0 ; /*goto40*/
+}
+
+
+/* this one is just an optimisation of the factor of 2 transform , we compute more things each turn */
+
+static int factorOf4Transform (void)
+{
+   int return_value = 0 ;
+
+   /*120*/
+   do
+      {
+         c1 = 1 ;
+         s1 = 0 ;
+
+         mm = min ( kspan , klim ) ;
+
+         do
+            {
+            f4t_150 () ;
+            return_value = f4t_170 () ;
+            } while (  return_value == 0 );
+
+         kk += ( inc - kspan ) ;
+      } while ( kk <= jc ) ;
+
+
+
+   if ( kspan == jc )
+        return 1 ; /*goto350*/
+   else
+         return 0 ;/*goto40*/
+
+
+}
+
+/*this function and the following are just here for conveniance , they just do fourier transformation  for factor of 4
+  but as the code was a bit long in factorof4transform , we've created two sub-functions */
+
+static void f4t_150 (void)
+{
+
+   do{
+      k1 = kk + kspan ;
+      k2 = k1 + kspan ;
+      k3 = k2 + kspan ;
+
+      akp = a[kk-1] + a[k2-1] ;
+      akm = a[kk-1] - a[k2-1] ;
+
+      ajp = a[k1-1] + a[k3-1] ;
+      ajm = a[k1-1] - a[k3-1] ;
+
+      a[kk-1] = akp + ajp ;
+      ajp = akp - ajp ;
+
+      bkp = b[kk-1] + b[k2-1] ;
+      bkm = b[kk-1] - b[k2-1] ;
+
+      bjp = b[k1-1] + b[k3-1] ;
+      bjm = b[k1-1] - b[k3-1] ;
+
+      b[kk-1] = bkp + bjp ;
+      bjp = bkp - bjp ;
+
+      akp = akm - bjm ;
+      akm = akm + bjm ;
+
+      bkp = bkm + ajm ;
+      bkm = bkm - ajm ;
+
+      if ( s1 == 0 )/*190*/
+         {
+         a[k1-1] = akp ;
+         a[k2-1] = ajp ;
+         a[k3-1] = akm ;
+
+         b[k1-1] = bkp ;
+         b[k2-1] = bjp ;
+         b[k3-1] = bkm ;
+
+
+         }
+      else /*160*/
+         {
+
+         a[k1-1] = akp*c1 - bkp*s1 ;
+         a[k2-1] = ajp*c2 - bjp*s2 ;
+         a[k3-1] = akm*c3 - bkm*s3 ;
+
+         a[k1-1] = bkp*c1 + akp*s1 ;
+         a[k2-1] = bjp*c2 + ajp*s2 ;
+         a[k3-1] = bkm*c3 + akm*s3 ;
+         }
+       kk=k3+kspan;
+    }while  ( kk <= nt ) ;
+
+
+}
+
+static int  f4t_170 (void)
+{
+   kk += ( jc - nt ) ;
+
+   if ( kk <= mm )
+      {
+       c2 = c1 - (cd*c1 + sd*s1);
+       s1 = s1 + (sd*c1 - cd*s1);
+
+      /*
+      the following three statements compensate for truncation
+      error.  if rounded arithmetic is used, substitute
+      c1=c2
+      */
+
+      c1 = 0.5/(c2*c2+s1*s1) + 0.5 ;
+      s1 *= c1 ;
+      c1 *= c2 ;
+
+      /*140*/
+
+      c2 = c1*c1 - s1*s1 ;
+      s2 = c1*s1*2 ;
+      c3 = c2*c1 - s2*s1 ;
+      s3 = c2*s1 + s2*c1 ;
+
+
+      return 0 ;
+
+      }
+   else
+     {
+      if ( kk <= kspan )
+         {
+          s1 = dr*rad * (kk-1)/jc ;
+          c1 = cos (s1) ;
+          s1 = sin (s1) ;
+          mm = min ( kspan , mm  + klim );
+
+         /*140*/
+
+         c2 = c1*c1 - s1*s1 ;
+         s2 = c1*s1*2 ;
+      	 c3 = c2*c1 - s2*s1 ;
+         s3 = c2*s1 + s2*c1 ;
+
+         return 0 ;
+         }
+      }
+
+   return 1 ;
+}
+
+
+
+
+static void factorOf3Transform (void)
+{
+	do{
+		do{
+			k1 = kk + kspan ;
+			k2 = k1 + kspan ;
+
+			ak = a[kk-1] ;
+			bk = b[kk-1] ;
+
+			aj = a[k1-1] + a[k2-1] ;
+			bj = b[k1-1] + b[k2-1] ;
+
+			a[kk-1] = ak + aj ;
+			b[kk-1] = bk + bj ;
+
+			ak = -0.5*aj + ak ;
+			bk = -0.5*bj + bk ;
+
+			aj = (a[k1-1] - a[k2-1])*s120 ;
+			bj = (b[k1-1] - b[k2-1])*s120 ;
+
+			a[k1-1] = ak - bj ;
+			b[k1-1] = bk + aj ;
+			a[k2-1] = ak + bj ;
+			b[k2-1] = bk - aj ;
+
+			kk = k2 + kspan ;
+		} while (kk < nn);
+
+		kk -= nn ;
+	}while (kk <= kspan);
+
+}
+
+static void factorOf5Transform (void)
+{
+	c2 = c72*c72 - s72 *s72 ;
+	s2 = 2 * c72*s72;
+
+	do{
+		do{
+			k1 = kk + kspan ;
+			k2 = k1 + kspan ;
+			k3 = k2 + kspan ;
+			k4 = k3 + kspan ;
+
+
+
+			akp = a[k1-1] + a[k4-1] ;
+			akm = a[k1-1] - a[k4-1] ;
+
+			bkp = b[k1-1] + b[k4-1] ;
+			bkm = b[k1-1] - b[k4-1] ;
+
+			ajp = a[k2-1] + a[k3-1] ;
+			ajm = a[k2-1] - a[k3-1] ;
+
+			bjp = b[k2-1] + b[k3-1] ;
+			bjm = b[k2-1] - b[k3-1] ;
+
+			aa = a[kk-1] ;
+			bb = b[kk-1] ;
+
+			a[kk-1] = aa + akp + ajp;
+			b[kk-1] = bb + bkp + bjp;
+
+			ak = akp*c72 + ajp*c2 + aa ;
+			bk = bkp*c72 + bjp*c2 + bb ;
+
+			aj = akm*s72 + ajm*s2 ;
+			bj = bkm*s72 + bjm*s2 ;
+
+			a[k1-1] = ak - bj ;
+			a[k4-1] = ak + bj ;
+			b[k1-1] = bk + aj ;
+			b[k4-1] = bk - aj ;
+
+			ak = akp*c2 + ajp*c72 + aa ;
+			bk = bkp*c2 + bjp*c72 + bb ;
+
+			aj = akm*s2 - ajm*s72 ;
+
+			bj = bkm*s2 - bjm*s72 ;
+
+			a[k2-1] = ak - bj ;
+			a[k3-1] = ak + bj ;
+			b[k2-1] = bk + aj ;
+			b[k3-1] = bk - aj ;
+
+			kk = k4 + kspan;
+		}while (kk < nn);
+
+		kk -= nn ;
+	}while (kk <= kspan);
+}
+
+/* this function is the general case of non factor of 2 factor , the factorof3transform and factorof5trandform are just
+special case of this one */
+
+
+static void preFOtherTransform (void)
+{
+   jf = k ;
+   s1 = (rad*8)/k ;
+   c1 = cos (s1) ;
+   s1 = sin (s1) ;
+   ck[jf-1] = 1 ;
+   sk[jf-1] = 0 ;
+   j = 1 ;
+
+   do
+   {
+      ck[j-1] = ck[k-1] * c1 +  sk[k-1]*s1 ;
+      sk[j-1] = ck[k-1] * s1 -  sk[k-1]*c1 ;
+
+      k -- ;
+
+      ck[k-1] = ck[j-1] ;
+      sk[k-1] = - sk[j-1] ;
+
+      j++ ;
+
+   }while ( j < k );
+
+}
+
+static void factorOfOtherTransform (void)
+{
+int ktemp = 0 ;
+
+do
+   {
+   k1 = kk ;
+   k2 = kk + kspnn ;
+
+   aa = a[kk-1] ;
+   bb = b[kk-1] ;
+
+   ak = aa ;
+   bk = bb ;
+   j = 1 ;
+   k1 += kspan ;
+
+   do
+   {
+      k2 -= kspan ;
+      j++ ;
+
+      wt[j-1] = a[k1-1] + a[k2-1] ;
+      ak = wt[j-1] + ak ;
+
+      bt[j-1] = b[k1-1] + b[k2-1] ;
+      bk = bt[j-1] + bk ;
+      j++ ;
+
+      wt[j-1] = a[k1-1] - a[k2-1] ;
+      bt[j-1] = b[k1-1] - b[k2-1] ;
+
+      k1 += kspan;
+   }while ( k1 < k2 ) ;
+
+   a[kk-1] = ak ;
+   b[kk-1] = bk ;
+
+   k1 = kk ;
+   k2 = kk + kspnn ;
+   j = 1 ;
+
+   do
+      {
+      k1 += kspan ;
+      k2 -= kspan ;
+      jj = j  ;
+      ak = aa ;
+      bk  = bb ;
+      aj = 0 ;
+      bj = 0 ;
+      k  = 1 ;
+
+      do
+         {
+         k++ ;
+         ak += ( wt[k-1] * ck[jj-1] ) ;
+         bk += ( bt[k-1] * ck[jj-1] ) ;
+
+         k++ ;
+         aj += (wt[k-1] * sk[jj-1]) ;
+         bj += (bt[k-1] * sk[jj-1]) ;
+         jj += j ;
+
+         if ( jj > jf )
+            jj -= jf ;
+         } while ( k < jf ) ;
+
+      k = jf - j ;
+      a[k1-1] = ak - bj ;
+      b[k1-1] = bk + aj ;
+      a[k2-1] = ak + bj ;
+      b[k2-1] = bk - aj ;
+
+      j++ ;
+
+      }while ( j < k ) ;
+
+
+
+
+
+   kk += kspnn ;
+   ktemp = kk ;
+
+   if ( kk > nn )
+      {
+         kk -= nn;
+      }
+
+}while (  ktemp <= nn || (kk <= kspan &&  ( ktemp > nn))  );
+}
+
+
+
+
+static int mulByRotationFactor (void )
+{
+ int ktemp = 0 ;
+
+   if ( i != m )
+      {
+       kk = jc + 1 ;
+
+       /*300*/
+       do
+         {
+         c2 = 1 - cd ;
+         s1 = sd ;
+
+         mm = min ( kspan , klim ) ;
+
+         /*320 */
+         do
+           {
+
+            c1 = c2 ;
+            s2 = s1 ;
+            kk += kspan ;
+
+            do
+              {
+
+                ak = a[kk-1] ;
+                a[kk-1] = c2*ak - s2*b[kk-1] ;
+                b[kk-1] = s2*ak + c2*b[kk-1] ;
+
+                kk += kspnn ;
+                ktemp = kk ;
+
+
+                if ( kk > nt )
+                  {
+                   ak = s1*s2 ;
+                   s2 = s1*c2 + s2*c1 ;
+                   c2 = c1*c2 - ak ;
+                   kk += (kspan - nt ) ;
+
+
+
+                  }
+
+               }while (ktemp <= nt || ( kk <= kspnn && ktemp > nt  )) ;
+
+            kk += ( jc - kspnn );
+
+
+            if ( kk <= mm )
+               {
+
+               /* 310*/
+                c2 = c1 - ( cd*c1 + sd*s1 ) ;
+                s1 += (sd*c1 - cd*s1 ) ;
+
+               /*
+               the following three statements compensate for truncation
+               error.  if rounded arithmetic is used, substitute
+               c1=c2
+               */
+
+                c1 = 0.5/(c2*c2+s1*s1) + 0.5 ;
+                s1 *= c1 ;
+                c2 *= c1 ;
+               }
+            else
+               {
+                  if ( kk <= kspan )
+                     {
+
+                      s1 = dr*rad * (kk-1)/jc ;
+                      c2 = cos (s1) ;
+                      s1 = sin (s1) ;
+                      mm = min ( kspan , mm  + klim );
+                     }
+               }
+
+            }while ( kk <= mm  ||( kk <= kspan && kk > mm ) ) ;
+
+          kk += (jc + inc -kspan );
+
+
+         }while ( kk <= jc+jc);
+
+
+       return 0 ; /* goto40 */
+      }
+
+   return 1 ; /* goto350*/
+}
+
+
+
+
+static void pre_sqFactor2NormlOrder (void)
+{
+
+   k = kt + kt + 1 ;
+
+   if ( m < k )
+      k -- ;
+
+   j = 1 ;
+   np[k] = jc ;
+
+   do
+      {
+        np[j] = np[j-1]/nfac[j-1] ;
+        np[k-1] = np[k]*nfac[j-1] ;
+
+        j++ ;
+        k-- ;
+      }while ( j <  k ) ;
+
+   k3 = np[k] ;
+   kspan = np[1] ;
+   kk = jc + 1 ;
+   k2 = kspan + 1 ;
+   j = 1;
+
+
+}
+
+static int  post_sqFactor2NormlOrder (void)
+{
+
+    do
+     {
+      do
+        {
+          k2 -= np[j-1] ;
+          j++ ;
+          k2 += np[j]  ;
+
+         } while ( k2 >  np[j-1]);
+
+       j = 1 ;
+
+/* 390 */
+       do
+         {
+          if ( kk < k2 )
+            {
+
+
+             return 1 ;
+            }
+          else
+            {
+               kk += inc ;
+               k2 += kspan ;
+            }
+         }while( k2 < ks  );
+
+      }while ( kk < ks ) ;
+
+   jc = k3 ;
+
+   return 0;
+}
+
+
+/* appeler cetter fonction dans un do while valeur_retour != 1)*/
+static void  single_sqFactor2NormlOrder (void)
+{
+
+
+      do
+        {
+
+         ak = a[kk-1] ;
+         a[kk-1] = a[k2-1] ;
+         a[k2-1] = ak ;
+
+         bk = b[kk-1] ;
+         b[kk-1] = b[k2-1] ;
+         b[k2-1] = bk ;
+
+         kk += inc ;
+         k2 += kspan ;
+        } while ( k2 < ks );
+
+      /*380*/
+}
+
+/*idem que single_ */
+static int multi_sqFactor2NormlOrder (void)
+{
+
+
+
+      k = kk + jc ;
+
+      do /*410*/
+        {
+         ak = a[kk-1] ;
+         a[kk-1] = a[k2-1] ;
+         a[k2-1] = ak ;
+         bk = b[kk-1] ;
+         b[kk-1] = b[k2-1] ;
+         b[k2-1] = bk ;
+         kk += inc ;
+         k2 += kspan ;
+        } while ( kk < k );
+
+      kk += (ks - jc ) ;
+      k2 += (ks - jc ) ;
+
+      if ( kk < nt )
+         return 1 ;
+
+      k2 += (  kspan - nt );
+      kk += (  jc    - nt );
+
+      if ( k2 < ks )
+         {
+
+         return 1 ;
+         }
+      if( post_sqFactor2NormlOrder ( ) == 1 )
+         {
+
+         return 1 ;
+         }
+   jc = k3 ;
+
+   return 0;
+
+}
+
+
+
+static void nonSqFactor2NormOrder (void)
+{
+
+   j = m - kt ;
+   nfac[j] = 1 ;
+
+
+
+    do
+      {
+         nfac[j-1] *= nfac[j] ;
+
+         j-- ;
+
+
+      }while ( j != kt ) ;
+
+   kt ++ ;
+   nn = nfac[kt-1] - 1;
+
+   jj = 0 ;
+   j = 0;
+
+   /*480*/
+
+   k2 = nfac[kt-1] ;
+   k = kt + 1 ;
+   kk = nfac[k-1] ;
+   j ++ ;
+
+   while ( j <= nn )
+         {
+            jj += kk ;
+
+
+            while ( jj >=  k2 )
+               {
+                jj -= k2 ;
+                k2 = kk ;
+                k++ ;
+                kk = nfac[k-1] ;
+
+                jj += kk ;
+
+
+               }
+
+            np[j-1] = jj ;
+            k2 = nfac[kt-1] ;
+            k = kt + 1 ;
+            kk = nfac[k-1] ;
+            j ++ ;
+
+         }
+
+   j = 0 ;
+
+   return ;
+}
+
+/*  here we determine how many permutation cycles we need to do */
+static void detPermutCycles (void)
+{
+
+  do
+   {
+    do
+      {
+       j++ ;
+
+       kk = np[j-1] ;
+      }while ( kk < 0 ) ;
+
+
+    if ( kk != j )
+      {
+         do
+            {
+
+               k = kk ;
+               kk = np[k-1] ;
+               np[k-1] = -kk ;
+            }while ( kk != j ) ;
+         k3 = kk ;
+      }
+    else
+       np[j-1] = -j ;
+   }while ( j != nn );
+
+   maxf *= inc ;
+
+   return ;
+}
+
+static void  reorderMatrix (void)
+{
+do
+  {
+   do
+      {
+       j--  ;
+
+      }while (np[j-1] < 0 ) ;
+
+   jj = jc ;
+
+   /*520*/
+   do
+      {
+         kspan = jj ;
+
+         if ( jj > maxf )
+            kspan = maxf ;
+
+         jj -= kspan ;
+         k  = np [j-1];
+         kk = jc*k + i + jj ;
+         k1 = kk + kspan ;
+         k2 = 0 ;
+
+         do /*530*/
+            {
+               k2 ++ ;
+               wt[k2-1] = a[k1-1] ;
+               bt[k2-1] = b[k1-1] ;
+               k1 -= inc ;
+
+            }while ( k1 != kk );
+
+         do
+            {
+            k1 = kk + kspan ;
+            k2 = k1 - jc * (k + np[k-1]);
+            k = -np[k-1];
+
+
+            do
+               {
+                  a[k1-1] = a[k2-1] ;
+                  b[k1-1] = b[k2-1] ;
+
+                  k1 -= inc ;
+                  k2 -= inc ;
+
+                  }while ( k1 != kk ) ;
+
+            kk = k2 ;
+
+            }while ( k != j );
+
+            k1 = kk +kspan ;
+            k2 = 0 ;
+         /*560*/
+
+         do
+            {
+            k2 ++ ;
+            a[k1-1] = wt[k2-1] ;
+            b[k1-1] = bt[k2-1] ;
+            k1 -= inc ;
+
+
+            }while (  k1 != kk ) ;
+
+   } while ( jj != 0 ) ;
+}while ( j != 1 ) ;
+
+   return ;
+}
+
+