gsl_permutation__canonical.c 3.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195
  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. /* Modified for GSL by Brian Gough.
  20. * Use in-place algorithms, no need for workspace
  21. * Use conventions for canonical form given in Knuth (opposite of Sedgewick)
  22. */
  23. #include "gsl__config.h"
  24. #include "gsl_errno.h"
  25. #include "gsl_permutation.h"
  26. int
  27. gsl_permutation_linear_to_canonical (gsl_permutation * q,
  28. const gsl_permutation * p)
  29. {
  30. const size_t n = p->size;
  31. size_t i, k, s;
  32. size_t t = n;
  33. const size_t *const pp = p->data;
  34. size_t *const qq = q->data;
  35. if (q->size != p->size)
  36. {
  37. GSL_ERROR ("size of q does not match size of p", GSL_EINVAL);
  38. }
  39. for (i = 0; i < n; i++)
  40. {
  41. k = pp[i];
  42. s = 1;
  43. while (k > i)
  44. {
  45. k = pp[k];
  46. s++;
  47. }
  48. if (k < i)
  49. continue;
  50. /* Now have k == i, i.e the least in its cycle, and s == cycle length */
  51. t -= s;
  52. qq[t] = i;
  53. k = pp[i];
  54. s = 1;
  55. while (k > i)
  56. {
  57. qq[t + s] = k;
  58. k = pp[k];
  59. s++;
  60. }
  61. if (t == 0)
  62. break;
  63. }
  64. return GSL_SUCCESS;
  65. }
  66. int
  67. gsl_permutation_canonical_to_linear (gsl_permutation * p,
  68. const gsl_permutation * q)
  69. {
  70. size_t i, k, kk, first;
  71. const size_t n = p->size;
  72. size_t *const pp = p->data;
  73. const size_t *const qq = q->data;
  74. if (q->size != p->size)
  75. {
  76. GSL_ERROR ("size of q does not match size of p", GSL_EINVAL);
  77. }
  78. for (i = 0; i < n; i++)
  79. {
  80. pp[i] = i;
  81. }
  82. k = qq[0];
  83. first = pp[k];
  84. for (i = 1; i < n; i++)
  85. {
  86. kk = qq[i];
  87. if (kk > first)
  88. {
  89. pp[k] = pp[kk];
  90. k = kk;
  91. }
  92. else
  93. {
  94. pp[k] = first;
  95. k = kk;
  96. first = pp[kk];
  97. }
  98. }
  99. pp[k] = first;
  100. return GSL_SUCCESS;
  101. }
  102. size_t
  103. gsl_permutation_inversions (const gsl_permutation * p)
  104. {
  105. size_t count = 0;
  106. size_t i, j;
  107. const size_t size = p->size;
  108. for (i = 0; i < size - 1; i++)
  109. {
  110. for (j = i + 1; j < size; j++)
  111. {
  112. if (p->data[i] > p->data[j])
  113. {
  114. count++;
  115. }
  116. }
  117. }
  118. return count;
  119. }
  120. size_t
  121. gsl_permutation_linear_cycles (const gsl_permutation * p)
  122. {
  123. size_t i, k;
  124. size_t count = 0;
  125. const size_t size = p->size;
  126. for (i = 0; i < size; i++)
  127. {
  128. k = p->data[i];
  129. while (k > i)
  130. {
  131. k = p->data[k];
  132. }
  133. if (k < i)
  134. continue;
  135. count++;
  136. }
  137. return count;
  138. }
  139. size_t
  140. gsl_permutation_canonical_cycles (const gsl_permutation * p)
  141. {
  142. size_t i;
  143. size_t count = 1;
  144. size_t min = p->data[0];
  145. for (i = 0; i < p->size; i++)
  146. {
  147. if (p->data[i] < min)
  148. {
  149. min = p->data[i];
  150. count++;
  151. }
  152. }
  153. return count;
  154. }