Moving source code

1 files changed, 0 insertions, 1211 deletions

diff --git a/src/signalProcessing/fft/dfftmx.c b/src/signalProcessing/fft/dfftmx.c
deleted file mode 100644
index f7d6ce30..00000000
--- a/src/signalProcessing/fft/dfftmx.c
+++ /dev/null
@@ -1,1211 +0,0 @@
-/*
- *  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 "fft_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 */
-
-/*note on this code all numbers alone in comment is
- a reference to the corresponding goto in the original fotran code */
-
-int dfftmx ( 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 ;
-   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((double)1);
-   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 ( ) ;
-
-   _pdblA = a ;
-   _pdblB = b ;
-
-   return 0 ;
-}
-
-/** **************************************
-Sous-Fonctions
-******************************************/
-
-
-
-/* this function only set the value of variable */
-static void preliminaryWork (void)
-{
-   s72 = -s72 ;
-   s120= -s120;
-   rad = -rad ;
-
-   kspan = ks ;
-   nn = nt -inc ;
-   jc = ks/n ;
-
-   /* sin , cos values are re-initialized each lim steps */
-
-   klim = 32*jc ;
-   i = 0;
-   jf = 0 ;
-
-   if ( kt > 0 )  maxf = max ( nfac[kt-1] , nfac[m-kt-1] );
-   else 	  maxf = nfac[m-kt-1] ;
-
-}
-
-
-
-/*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 :
-
-         kspnn = kspan ;
-         kspan = kspan >> 2 ; /*kspan /= 4 */
-
-         retVal = factorOf4Transform () ;
-         break ;
-
-      case 3 :
-
-         k = nfac[i-1] ;
-         kspnn = kspan ;
-         kspan = kspan / k ;
-
-         factorOf3Transform ( ) ;
-         break ;
-
-      case 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 ;
-}
-
-