home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / permutation / permute_source.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-11-01  |  4.0 KB  |  164 lines
  1. /* permutation/permute_source.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* In-place Permutations 
  21.  
  22.    permute:    OUT[i]       = IN[perm[i]]     i = 0 .. N-1
  23.    invpermute: OUT[perm[i]] = IN[i]           i = 0 .. N-1
  24.  
  25.    PERM is an index map, i.e. a vector which contains a permutation of
  26.    the integers 0 .. N-1.
  27.  
  28.    From Knuth "Sorting and Searching", Volume 3 (3rd ed), Section 5.2
  29.    Exercise 10 (answers), p 617
  30.  
  31.    FIXME: these have not been fully tested.
  32. */
  33.  
  34. int
  35. TYPE (gsl_permute) (const size_t * p, ATOMIC * data, const size_t stride, const size_t n)
  36. {
  37.   size_t i, k, pk;
  38.  
  39.   for (i = 0; i < n; i++)
  40.     {
  41.       k = p[i];
  42.       
  43.       while (k > i) 
  44.         k = p[k];
  45.       
  46.       if (k < i)
  47.         continue ;
  48.       
  49.       /* Now have k == i, i.e the least in its cycle */
  50.       
  51.       pk = p[k];
  52.       
  53.       if (pk == i)
  54.         continue ;
  55.       
  56.       /* shuffle the elements of the cycle */
  57.       
  58.       {
  59.         unsigned int a;
  60.  
  61.         ATOMIC t[MULTIPLICITY];
  62.         
  63.         for (a = 0; a < MULTIPLICITY; a++)
  64.           t[a] = data[i*stride*MULTIPLICITY + a];
  65.       
  66.         while (pk != i)
  67.           {
  68.             for (a = 0; a < MULTIPLICITY; a++)
  69.               {
  70.                 ATOMIC r1 = data[pk*stride*MULTIPLICITY + a];
  71.                 data[k*stride*MULTIPLICITY + a] = r1;
  72.               }
  73.             k = pk;
  74.             pk = p[k];
  75.           };
  76.         
  77.         for (a = 0; a < MULTIPLICITY; a++)
  78.           data[k*stride*MULTIPLICITY + a] = t[a];
  79.       }
  80.     }
  81.  
  82.   return GSL_SUCCESS;
  83. }
  84.  
  85. int
  86. FUNCTION (gsl_permute,inverse) (const size_t * p, ATOMIC * data, const size_t stride, const size_t n)
  87. {
  88.   size_t i, k, pk;
  89.  
  90.   for (i = 0; i < n; i++)
  91.     {
  92.       k = p[i];
  93.           
  94.       while (k > i) 
  95.         k = p[k];
  96.  
  97.       if (k < i)
  98.         continue ;
  99.       
  100.       /* Now have k == i, i.e the least in its cycle */
  101.  
  102.       pk = p[k];
  103.  
  104.       if (pk == i)
  105.         continue ;
  106.       
  107.       /* shuffle the elements of the cycle in the inverse direction */
  108.       
  109.       {
  110.         unsigned int a;
  111.  
  112.         ATOMIC t[MULTIPLICITY];
  113.  
  114.         for (a = 0; a < MULTIPLICITY; a++)
  115.           t[a] = data[k*stride*MULTIPLICITY+a];
  116.  
  117.         while (pk != i)
  118.           {
  119.             for (a = 0; a < MULTIPLICITY; a++)
  120.               {
  121.                 ATOMIC r1 = data[pk*stride*MULTIPLICITY + a];
  122.                 data[pk*stride*MULTIPLICITY + a] = t[a];
  123.                 t[a] = r1;
  124.               }
  125.  
  126.             k = pk;
  127.             pk = p[k];
  128.           };
  129.  
  130.         for (a = 0; a < MULTIPLICITY; a++)
  131.           data[pk*stride*MULTIPLICITY+a] = t[a];
  132.       }
  133.     }
  134.  
  135.   return GSL_SUCCESS;
  136. }
  137.  
  138.  
  139. int
  140. TYPE (gsl_permute_vector) (const gsl_permutation * p, TYPE (gsl_vector) * v)
  141. {
  142.   if (v->size != p->size)
  143.     {
  144.       GSL_ERROR ("vector and permutation must be the same length", GSL_EBADLEN);
  145.     }
  146.  
  147.   TYPE (gsl_permute) (p->data, v->data, v->stride, v->size) ;
  148.  
  149.   return GSL_SUCCESS;
  150. }
  151.  
  152. int
  153. FUNCTION (gsl_permute_vector,inverse) (const gsl_permutation * p, TYPE (gsl_vector) * v)
  154. {
  155.   if (v->size != p->size)
  156.     {
  157.       GSL_ERROR ("vector and permutation must be the same length", GSL_EBADLEN);
  158.     }
  159.  
  160.   FUNCTION (gsl_permute,inverse) (p->data, v->data, v->stride, v->size) ;
  161.  
  162.   return GSL_SUCCESS;
  163. }
  164.