htswanalysis/MACS/lib/gsl/gsl-1.11/permutation/canonical.c

   1 /* permutation/permutation.c
   2  *
   3  * Copyright (C) 2001, 2002 Nicolas Darnis
   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 3 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Modified for GSL by Brian Gough.
  21  * Use in-place algorithms, no need for workspace
  22  * Use conventions for canonical form given in Knuth (opposite of Sedgewick)
  23  */
  24
  25 #include <config.h>
  26 #include <gsl/gsl_errno.h>
  27 #include <gsl/gsl_permutation.h>
  28
  29 int
  30 gsl_permutation_linear_to_canonical (gsl_permutation * q,
  31                                      const gsl_permutation * p)
  32 {
  33   const size_t n = p->size;
  34   size_t i, k, s;
  35   size_t t = n;
  36
  37   const size_t *const pp = p->data;
  38   size_t *const qq = q->data;
  39
  40   if (q->size != p->size)
  41     {
  42       GSL_ERROR ("size of q does not match size of p", GSL_EINVAL);
  43     }
  44
  45   for (i = 0; i < n; i++)
  46     {
  47
  48       k = pp[i];
  49       s = 1;
  50
  51       while (k > i)
  52         {
  53           k = pp[k];
  54           s++;
  55         }
  56
  57       if (k < i)
  58         continue;
  59
  60       /* Now have k == i, i.e the least in its cycle, and s == cycle length */
  61
  62       t -= s;
  63
  64       qq[t] = i;
  65
  66       k = pp[i];
  67       s = 1;
  68
  69       while (k > i)
  70         {
  71           qq[t + s] = k;
  72           k = pp[k];
  73           s++;
  74         }
  75
  76       if (t == 0)
  77         break;
  78     }
  79
  80   return GSL_SUCCESS;
  81 }
  82
  83 int
  84 gsl_permutation_canonical_to_linear (gsl_permutation * p,
  85                                      const gsl_permutation * q)
  86 {
  87   size_t i, k, kk, first;
  88   const size_t n = p->size;
  89
  90   size_t *const pp = p->data;
  91   const size_t *const qq = q->data;
  92
  93   if (q->size != p->size)
  94     {
  95       GSL_ERROR ("size of q does not match size of p", GSL_EINVAL);
  96     }
  97
  98   for (i = 0; i < n; i++)
  99     {
 100       pp[i] = i;
 101     }
 102
 103   k = qq[0];
 104   first = pp[k];
 105
 106   for (i = 1; i < n; i++)
 107     {
 108       kk = qq[i];
 109
 110       if (kk > first)
 111         {
 112           pp[k] = pp[kk];
 113           k = kk;
 114         }
 115       else
 116         {
 117           pp[k] = first;
 118           k = kk;
 119           first = pp[kk];
 120         }
 121     }
 122
 123   pp[k] = first;
 124
 125   return GSL_SUCCESS;
 126 }
 127
 128
 129 size_t
 130 gsl_permutation_inversions (const gsl_permutation * p)
 131 {
 132   size_t count = 0;
 133   size_t i, j;
 134   const size_t size = p->size;
 135
 136   for (i = 0; i < size - 1; i++)
 137     {
 138       for (j = i + 1; j < size; j++)
 139         {
 140           if (p->data[i] > p->data[j])
 141             {
 142               count++;
 143             }
 144         }
 145     }
 146
 147   return count;
 148 }
 149
 150 size_t
 151 gsl_permutation_linear_cycles (const gsl_permutation * p)
 152 {
 153   size_t i, k;
 154   size_t count = 0;
 155   const size_t size = p->size;
 156
 157   for (i = 0; i < size; i++)
 158     {
 159
 160       k = p->data[i];
 161
 162       while (k > i)
 163         {
 164           k = p->data[k];
 165         }
 166
 167       if (k < i)
 168         continue;
 169
 170       count++;
 171     }
 172
 173   return count;
 174 }
 175
 176 size_t
 177 gsl_permutation_canonical_cycles (const gsl_permutation * p)
 178 {
 179   size_t i;
 180   size_t count = 1;
 181   size_t min = p->data[0];
 182
 183   for (i = 0; i < p->size; i++)
 184     {
 185       if (p->data[i] < min)
 186         {
 187           min = p->data[i];
 188           count++;
 189         }
 190     }
 191
 192   return count;
 193 }
 194