This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[PATCH] Adding multiset functionality atop GSL 1.13


Hi all,

Attached is a patch that adds multiset functionality atop GSL 1.13.
Multisets are like combinations but with replacement-- a single item
may appear multiple times.  Element order is insignificant and
enumeration uses lexicographical ordering.

The patch builds very, very, very heavily on the existing combinations
code.  Documentation, an example, and unit tests are included.
Please consider it for inclusion in the next feature release.

- Rhys

P.S. One weakness in the patch the amount of code duplication relative
to combinations.  The underlying data structures are the same, and
only the _init_first, _init_last, _valid, _next, and _prev functions
differ in a non-trivial way.  With some thought to naming, multisets
could be a combinations-with-replacement flag.  If you accept the
patch, please consider adding a TODO to both directories to merge the
functionality.
From 6effe93b67b84271ca21d4fde3f19c79eefbd49b Mon Sep 17 00:00:00 2001
From: Rhys Ulerich <rhys.ulerich@gmail.com>
Date: Mon, 21 Sep 2009 20:20:26 -0500
Subject: [PATCH] Added multiset functionality and documentation, based absurdly on combinations (Rhys Ulerich)

---
 Makefile.am               |    4 +-
 configure.ac              |    2 +-
 doc/Makefile.am           |    6 +-
 doc/examples/multiset.c   |   25 ++++
 doc/examples/multiset.out |   71 ++++++++++++
 doc/gsl-ref.texi          |    9 +-
 doc/multiset.texi         |  201 ++++++++++++++++++++++++++++++++
 multiset/.cvsignore       |    8 ++
 multiset/ChangeLog        |    3 +
 multiset/Makefile.am      |   21 ++++
 multiset/demo.c           |   72 ++++++++++++
 multiset/file.c           |  115 +++++++++++++++++++
 multiset/gsl_multiset.h   |   93 +++++++++++++++
 multiset/init.c           |  130 +++++++++++++++++++++
 multiset/inline.c         |   24 ++++
 multiset/multiset.c       |  181 +++++++++++++++++++++++++++++
 multiset/test.c           |  278 +++++++++++++++++++++++++++++++++++++++++++++
 17 files changed, 1235 insertions(+), 8 deletions(-)
 create mode 100644 doc/examples/multiset.c
 create mode 100644 doc/examples/multiset.out
 create mode 100644 doc/multiset.texi
 create mode 100644 multiset/.cvsignore
 create mode 100644 multiset/ChangeLog
 create mode 100644 multiset/Makefile.am
 create mode 100644 multiset/demo.c
 create mode 100644 multiset/file.c
 create mode 100644 multiset/gsl_multiset.h
 create mode 100644 multiset/init.c
 create mode 100644 multiset/inline.c
 create mode 100644 multiset/multiset.c
 create mode 100644 multiset/test.c

diff --git a/Makefile.am b/Makefile.am
index b303732..45b7c15 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -2,9 +2,9 @@
 
 # AUTOMAKE_OPTIONS = readme-alpha
 
-SUBDIRS = gsl utils sys test err const complex cheb block vector matrix permutation combination sort ieee-utils cblas blas linalg eigen specfunc dht qrng rng randist fft poly fit multifit statistics siman sum integration interpolation histogram ode-initval roots multiroots min multimin monte ntuple diff deriv cdf wavelet bspline doc
+SUBDIRS = gsl utils sys test err const complex cheb block vector matrix permutation combination multiset sort ieee-utils cblas blas linalg eigen specfunc dht qrng rng randist fft poly fit multifit statistics siman sum integration interpolation histogram ode-initval roots multiroots min multimin monte ntuple diff deriv cdf wavelet bspline doc
 
-SUBLIBS = block/libgslblock.la blas/libgslblas.la bspline/libgslbspline.la complex/libgslcomplex.la cheb/libgslcheb.la dht/libgsldht.la diff/libgsldiff.la deriv/libgslderiv.la eigen/libgsleigen.la err/libgslerr.la fft/libgslfft.la fit/libgslfit.la histogram/libgslhistogram.la ieee-utils/libgslieeeutils.la integration/libgslintegration.la interpolation/libgslinterpolation.la linalg/libgsllinalg.la matrix/libgslmatrix.la min/libgslmin.la monte/libgslmonte.la multifit/libgslmultifit.la multimin/libgslmultimin.la multiroots/libgslmultiroots.la ntuple/libgslntuple.la ode-initval/libgslodeiv.la permutation/libgslpermutation.la combination/libgslcombination.la poly/libgslpoly.la qrng/libgslqrng.la randist/libgslrandist.la rng/libgslrng.la roots/libgslroots.la siman/libgslsiman.la sort/libgslsort.la specfunc/libgslspecfunc.la statistics/libgslstatistics.la sum/libgslsum.la sys/libgslsys.la test/libgsltest.la utils/libutils.la vector/libgslvector.la cdf/libgslcdf.la wavelet/libgslwavelet.la
+SUBLIBS = block/libgslblock.la blas/libgslblas.la bspline/libgslbspline.la complex/libgslcomplex.la cheb/libgslcheb.la dht/libgsldht.la diff/libgsldiff.la deriv/libgslderiv.la eigen/libgsleigen.la err/libgslerr.la fft/libgslfft.la fit/libgslfit.la histogram/libgslhistogram.la ieee-utils/libgslieeeutils.la integration/libgslintegration.la interpolation/libgslinterpolation.la linalg/libgsllinalg.la matrix/libgslmatrix.la min/libgslmin.la monte/libgslmonte.la multifit/libgslmultifit.la multimin/libgslmultimin.la multiroots/libgslmultiroots.la ntuple/libgslntuple.la ode-initval/libgslodeiv.la permutation/libgslpermutation.la combination/libgslcombination.la multiset/libgslmultiset.la poly/libgslpoly.la qrng/libgslqrng.la randist/libgslrandist.la rng/libgslrng.la roots/libgslroots.la siman/libgslsiman.la sort/libgslsort.la specfunc/libgslspecfunc.la statistics/libgslstatistics.la sum/libgslsum.la sys/libgslsys.la test/libgsltest.la utils/libutils.la vector/libgslvector.la cdf/libgslcdf.la wavelet/libgslwavelet.la
 
 pkginclude_HEADERS = gsl_math.h gsl_pow_int.h gsl_nan.h gsl_machine.h gsl_mode.h gsl_precision.h gsl_types.h gsl_version.h gsl_minmax.h gsl_inline.h
 
diff --git a/configure.ac b/configure.ac
index 5de7d8f..b1eb03f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -543,5 +543,5 @@ AH_VERBATIM([GSL_DISABLE_DEPRECATED],
 #define GSL_DISABLE_DEPRECATED 1])
 
 dnl
-AC_CONFIG_FILES([gsl_version.h gsl.spec gsl/Makefile test/Makefile err/Makefile sys/Makefile utils/Makefile const/Makefile min/Makefile multimin/Makefile ieee-utils/Makefile fft/Makefile specfunc/Makefile dht/Makefile fit/Makefile multifit/Makefile bspline/Makefile statistics/Makefile sum/Makefile roots/Makefile multiroots/Makefile ntuple/Makefile poly/Makefile qrng/Makefile rng/Makefile randist/Makefile siman/Makefile integration/Makefile interpolation/Makefile doc/Makefile block/Makefile vector/Makefile matrix/Makefile histogram/Makefile monte/Makefile ode-initval/Makefile cblas/Makefile blas/Makefile linalg/Makefile eigen/Makefile permutation/Makefile combination/Makefile sort/Makefile complex/Makefile diff/Makefile deriv/Makefile cheb/Makefile cdf/Makefile wavelet/Makefile Makefile])
+AC_CONFIG_FILES([gsl_version.h gsl.spec gsl/Makefile test/Makefile err/Makefile sys/Makefile utils/Makefile const/Makefile min/Makefile multimin/Makefile ieee-utils/Makefile fft/Makefile specfunc/Makefile dht/Makefile fit/Makefile multifit/Makefile bspline/Makefile statistics/Makefile sum/Makefile roots/Makefile multiroots/Makefile ntuple/Makefile poly/Makefile qrng/Makefile rng/Makefile randist/Makefile siman/Makefile integration/Makefile interpolation/Makefile doc/Makefile block/Makefile vector/Makefile matrix/Makefile histogram/Makefile monte/Makefile ode-initval/Makefile cblas/Makefile blas/Makefile linalg/Makefile eigen/Makefile permutation/Makefile combination/Makefile multiset/Makefile sort/Makefile complex/Makefile diff/Makefile deriv/Makefile cheb/Makefile cdf/Makefile wavelet/Makefile Makefile])
 AC_OUTPUT
diff --git a/doc/Makefile.am b/doc/Makefile.am
index 41578e6..c172250 100644
--- a/doc/Makefile.am
+++ b/doc/Makefile.am
@@ -2,15 +2,15 @@
 
 info_TEXINFOS = gsl-ref.texi
 noinst_TEXINFOS = gsl-design.texi
-gsl_ref_TEXINFOS = err.texi cblas.texi blas.texi min.texi fft.texi rng.texi randist.texi roots.texi statistics.texi specfunc.texi specfunc-airy.texi specfunc-bessel.texi specfunc-clausen.texi specfunc-coulomb.texi specfunc-coupling.texi specfunc-dawson.texi specfunc-debye.texi specfunc-dilog.texi specfunc-elementary.texi specfunc-ellint.texi specfunc-elljac.texi specfunc-erf.texi specfunc-exp.texi specfunc-expint.texi specfunc-fermi-dirac.texi specfunc-gamma.texi specfunc-gegenbauer.texi specfunc-hyperg.texi specfunc-lambert.texi specfunc-laguerre.texi specfunc-legendre.texi specfunc-log.texi specfunc-mathieu.texi specfunc-pow-int.texi specfunc-psi.texi specfunc-synchrotron.texi specfunc-transport.texi specfunc-trig.texi specfunc-zeta.texi siman.texi vectors.texi debug.texi histogram.texi ode-initval.texi integration.texi ieee754.texi montecarlo.texi sum.texi intro.texi usage.texi dwt.texi dht.texi interp.texi poly.texi linalg.texi eigen.texi multiroots.texi sort.texi permutation.texi combination.texi complex.texi math.texi fitting.texi multifit.texi const.texi ntuple.texi diff.texi qrng.texi cheb.texi multimin.texi gpl.texi fdl.texi autoconf.texi freemanuals.texi bspline.texi
+gsl_ref_TEXINFOS = err.texi cblas.texi blas.texi min.texi fft.texi rng.texi randist.texi roots.texi statistics.texi specfunc.texi specfunc-airy.texi specfunc-bessel.texi specfunc-clausen.texi specfunc-coulomb.texi specfunc-coupling.texi specfunc-dawson.texi specfunc-debye.texi specfunc-dilog.texi specfunc-elementary.texi specfunc-ellint.texi specfunc-elljac.texi specfunc-erf.texi specfunc-exp.texi specfunc-expint.texi specfunc-fermi-dirac.texi specfunc-gamma.texi specfunc-gegenbauer.texi specfunc-hyperg.texi specfunc-lambert.texi specfunc-laguerre.texi specfunc-legendre.texi specfunc-log.texi specfunc-mathieu.texi specfunc-pow-int.texi specfunc-psi.texi specfunc-synchrotron.texi specfunc-transport.texi specfunc-trig.texi specfunc-zeta.texi siman.texi vectors.texi debug.texi histogram.texi ode-initval.texi integration.texi ieee754.texi montecarlo.texi sum.texi intro.texi usage.texi dwt.texi dht.texi interp.texi poly.texi linalg.texi eigen.texi multiroots.texi sort.texi permutation.texi combination.texi multiset.texi complex.texi math.texi fitting.texi multifit.texi const.texi ntuple.texi diff.texi qrng.texi cheb.texi multimin.texi gpl.texi fdl.texi autoconf.texi freemanuals.texi bspline.texi
 
 man_MANS = gsl.3 gsl-config.1 gsl-randist.1 gsl-histogram.1
 
 figures = multimin.eps siman-test.eps siman-energy.eps 12-cities.eps initial-route.eps final-route.eps fft-complex-radix2-f.eps fft-complex-radix2-t.eps fft-complex-radix2.eps fft-real-mixedradix.eps roots-bisection.eps roots-false-position.eps roots-newtons-method.eps roots-secant-method.eps histogram.eps histogram2d.eps min-interval.eps fit-wlinear.eps fit-wlinear2.eps fit-exp.eps ntuple.eps qrng.eps cheb.eps vdp.eps interp2.eps rand-beta.tex rand-binomial.tex rand-cauchy.tex rand-chisq.tex rand-erlang.tex rand-exponential.tex rand-fdist.tex rand-flat.tex rand-gamma.tex rand-gaussian.tex rand-geometric.tex rand-laplace.tex rand-logarithmic.tex rand-logistic.tex rand-lognormal.tex rand-pareto.tex rand-poisson.tex rand-hypergeometric.tex rand-nbinomial.tex rand-pascal.tex rand-bivariate-gaussian.tex rand-rayleigh.tex rand-rayleigh-tail.tex rand-tdist.tex rand-weibull.tex random-walk.tex randplots.gnp rand-exppow.tex rand-levy.tex rand-levyskew.tex rand-gumbel.tex rand-bernoulli.tex rand-gaussian-tail.tex rand-gumbel1.tex rand-gumbel2.tex landau.dat rand-landau.tex dwt-orig.eps dwt-samp.eps interpp2.eps bspline.eps
 
-examples_src = examples/blas.c examples/block.c examples/cblas.c examples/cdf.c examples/cheb.c examples/combination.c examples/const.c examples/demo_fn.c examples/diff.c examples/eigen.c examples/fft.c examples/fftmr.c examples/fftreal.c examples/fitting.c examples/fitting2.c examples/fitting3.c examples/histogram.c examples/histogram2d.c examples/ieee.c examples/ieeeround.c examples/integration.c examples/interp.c examples/intro.c examples/linalglu.c examples/matrix.c examples/matrixw.c examples/min.c examples/monte.c examples/ntupler.c examples/ntuplew.c examples/ode-initval.c examples/odefixed.c examples/permseq.c examples/permshuffle.c examples/polyroots.c examples/qrng.c examples/randpoisson.c examples/randwalk.c examples/rng.c examples/rngunif.c examples/rootnewt.c examples/roots.c examples/siman.c examples/sortsmall.c examples/specfun.c examples/specfun_e.c examples/stat.c examples/statsort.c examples/sum.c examples/vector.c examples/vectorr.c examples/vectorview.c examples/vectorw.c examples/demo_fn.h examples/dwt.c examples/expfit.c examples/nlfit.c examples/interpp.c examples/eigen_nonsymm.c examples/bspline.c examples/multimin.c examples/multiminfn.c examples/nmsimplex.c
+examples_src = examples/blas.c examples/block.c examples/cblas.c examples/cdf.c examples/cheb.c examples/combination.c examples/multiset.c examples/const.c examples/demo_fn.c examples/diff.c examples/eigen.c examples/fft.c examples/fftmr.c examples/fftreal.c examples/fitting.c examples/fitting2.c examples/fitting3.c examples/histogram.c examples/histogram2d.c examples/ieee.c examples/ieeeround.c examples/integration.c examples/interp.c examples/intro.c examples/linalglu.c examples/matrix.c examples/matrixw.c examples/min.c examples/monte.c examples/ntupler.c examples/ntuplew.c examples/ode-initval.c examples/odefixed.c examples/permseq.c examples/permshuffle.c examples/polyroots.c examples/qrng.c examples/randpoisson.c examples/randwalk.c examples/rng.c examples/rngunif.c examples/rootnewt.c examples/roots.c examples/siman.c examples/sortsmall.c examples/specfun.c examples/specfun_e.c examples/stat.c examples/statsort.c examples/sum.c examples/vector.c examples/vectorr.c examples/vectorview.c examples/vectorw.c examples/demo_fn.h examples/dwt.c examples/expfit.c examples/nlfit.c examples/interpp.c examples/eigen_nonsymm.c examples/bspline.c examples/multimin.c examples/multiminfn.c examples/nmsimplex.c
 
-examples_out = examples/blas.out examples/block.out examples/cblas.out examples/cdf.out examples/combination.out examples/const.out examples/diff.out examples/integration.out examples/intro.out examples/linalglu.out examples/min.out examples/polyroots.out examples/randpoisson.2.out examples/randpoisson.out examples/rng.out examples/rngunif.2.out examples/rngunif.out examples/sortsmall.out examples/specfun.out examples/specfun_e.out examples/stat.out examples/statsort.out examples/sum.out examples/vectorview.out examples/ecg.dat examples/dwt.dat examples/multimin.out examples/nmsimplex.out
+examples_out = examples/blas.out examples/block.out examples/cblas.out examples/cdf.out examples/combination.out examples/multiset.out examples/const.out examples/diff.out examples/integration.out examples/intro.out examples/linalglu.out examples/min.out examples/polyroots.out examples/randpoisson.2.out examples/randpoisson.out examples/rng.out examples/rngunif.2.out examples/rngunif.out examples/sortsmall.out examples/specfun.out examples/specfun_e.out examples/stat.out examples/statsort.out examples/sum.out examples/vectorview.out examples/ecg.dat examples/dwt.dat examples/multimin.out examples/nmsimplex.out
 
 noinst_DATA = $(examples_src) $(examples_out) $(figures)
 
diff --git a/doc/examples/multiset.c b/doc/examples/multiset.c
new file mode 100644
index 0000000..6ac1e7c
--- /dev/null
+++ b/doc/examples/multiset.c
@@ -0,0 +1,25 @@
+#include <stdio.h>
+#include <gsl/gsl_multiset.h>
+
+int
+main (void)
+{
+  gsl_multiset * c;
+  size_t i;
+
+  printf ("All multisets of {0,1,2,3} by size:\n") ;
+  for (i = 0; i <= 4; i++)
+    {
+      c = gsl_multiset_calloc (4, i);
+      do
+        {
+          printf ("{");
+          gsl_multiset_fprintf (stdout, c, " %u");
+          printf (" }\n");
+        }
+      while (gsl_multiset_next (c) == GSL_SUCCESS);
+      gsl_multiset_free (c);
+    }
+
+  return 0;
+}
diff --git a/doc/examples/multiset.out b/doc/examples/multiset.out
new file mode 100644
index 0000000..b071e1b
--- /dev/null
+++ b/doc/examples/multiset.out
@@ -0,0 +1,71 @@
+all multisets of {0,1,2,3} by size:
+{ }
+{ 0 }
+{ 1 }
+{ 2 }
+{ 3 }
+{ 0 0 }
+{ 0 1 }
+{ 0 2 }
+{ 0 3 }
+{ 1 1 }
+{ 1 2 }
+{ 1 3 }
+{ 2 2 }
+{ 2 3 }
+{ 3 3 }
+{ 0 0 0 }
+{ 0 0 1 }
+{ 0 0 2 }
+{ 0 0 3 }
+{ 0 1 1 }
+{ 0 1 2 }
+{ 0 1 3 }
+{ 0 2 2 }
+{ 0 2 3 }
+{ 0 3 3 }
+{ 1 1 1 }
+{ 1 1 2 }
+{ 1 1 3 }
+{ 1 2 2 }
+{ 1 2 3 }
+{ 1 3 3 }
+{ 2 2 2 }
+{ 2 2 3 }
+{ 2 3 3 }
+{ 3 3 3 }
+{ 0 0 0 0 }
+{ 0 0 0 1 }
+{ 0 0 0 2 }
+{ 0 0 0 3 }
+{ 0 0 1 1 }
+{ 0 0 1 2 }
+{ 0 0 1 3 }
+{ 0 0 2 2 }
+{ 0 0 2 3 }
+{ 0 0 3 3 }
+{ 0 1 1 1 }
+{ 0 1 1 2 }
+{ 0 1 1 3 }
+{ 0 1 2 2 }
+{ 0 1 2 3 }
+{ 0 1 3 3 }
+{ 0 2 2 2 }
+{ 0 2 2 3 }
+{ 0 2 3 3 }
+{ 0 3 3 3 }
+{ 1 1 1 1 }
+{ 1 1 1 2 }
+{ 1 1 1 3 }
+{ 1 1 2 2 }
+{ 1 1 2 3 }
+{ 1 1 3 3 }
+{ 1 2 2 2 }
+{ 1 2 2 3 }
+{ 1 2 3 3 }
+{ 1 3 3 3 }
+{ 2 2 2 2 }
+{ 2 2 2 3 }
+{ 2 2 3 3 }
+{ 2 3 3 3 }
+{ 3 3 3 3 }
diff --git a/doc/gsl-ref.texi b/doc/gsl-ref.texi
index 6bba7c9..13b57d2 100644
--- a/doc/gsl-ref.texi
+++ b/doc/gsl-ref.texi
@@ -248,6 +248,7 @@ project homepage thanks to Daisuke Tominaga.
 * Vectors and Matrices::        
 * Permutations::                
 * Combinations::                
+* Multisets::
 * Sorting::                     
 * BLAS Support::                
 * Linear Algebra::              
@@ -329,11 +330,15 @@ project homepage thanks to Daisuke Tominaga.
 @chapter Permutations
 @include permutation.texi
 
-@node   Combinations, Sorting, Permutations, Top
+@node   Combinations, Multisets, Permutations, Top
 @chapter Combinations
 @include combination.texi
 
-@node  Sorting, BLAS Support, Combinations, Top
+@node   Multisets, Sorting, Combinations, Top
+@chapter Multisets
+@include multiset.texi
+
+@node  Sorting, BLAS Support, Multisets, Top
 @chapter Sorting
 @include sort.texi
 
diff --git a/doc/multiset.texi b/doc/multiset.texi
new file mode 100644
index 0000000..000f532
--- /dev/null
+++ b/doc/multiset.texi
@@ -0,0 +1,201 @@
+@cindex multisets
+
+This chapter describes functions for creating and manipulating multisets. A
+multiset @math{c} is represented by an array of @math{k} integers in the range
+0 to @math{n-1}, where each value @math{c_i} may occur more than once.  The
+multiset @math{c} corresponds to indices of @math{k} elements chosen from an
+@math{n} element vector with replacement.  In mathematical terms, @math{n} is
+the cardinality of the multiset while @math{k} is the maximum multiplicity of
+any value.  Multisets are useful, for example, when iterating over the indices
+of a @math{k}-th order symmetric tensor in @math{n}-space.
+
+The functions described in this chapter are defined in the header file
+@file{gsl_multiset.h}.
+
+@menu
+* The Multiset struct::
+* Multiset allocation::
+* Accessing multiset elements::
+* Multiset properties::
+* Multiset functions::
+* Reading and writing multisets::
+* Multiset Examples::
+@end menu
+
+@node The Multiset struct
+@section The Multiset struct
+@tpindex gsl_multiset
+A multiset is defined by a structure containing three components, the
+values of @math{n} and @math{k}, and a pointer to the multiset array.
+The elements of the multiset array are all of type @code{size_t}, and
+are stored in increasing order.  The @code{gsl_multiset} structure
+looks like this,
+
+@example
+typedef struct
+@{
+  size_t n;
+  size_t k;
+  size_t *data;
+@} gsl_multiset;
+@end example
+@comment
+
+@noindent
+
+@node Multiset allocation
+@section Multiset allocation
+
+@deftypefun {gsl_multiset *} gsl_multiset_alloc (size_t @var{n}, size_t @var{k})
+This function allocates memory for a new multiset with parameters @var{n},
+@var{k}.  The multiset is not initialized and its elements are undefined.  Use
+the function @code{gsl_multiset_calloc} if you want to create a multiset which
+is initialized to the lexicographically first multiset element. A null pointer
+is returned if insufficient memory is available to create the multiset.
+@end deftypefun
+
+@deftypefun {gsl_multiset *} gsl_multiset_calloc (size_t @var{n}, size_t @var{k})
+This function allocates memory for a new multiset with parameters @var{n},
+@var{k} and initializes it to the lexicographically first multiset element. A
+null pointer is returned if insufficient memory is available to create the
+multiset.
+@end deftypefun
+
+@deftypefun void gsl_multiset_init_first (gsl_multiset * @var{c})
+This function initializes the multiset @var{c} to the lexicographically first
+multiset element, i.e. @math{0} repeated @math{k} times.
+@end deftypefun
+
+@deftypefun void gsl_multiset_init_last (gsl_multiset * @var{c})
+This function initializes the multiset @var{c} to the lexicographically last
+multiset element, i.e. @math{n-1} repeated @math{k} times.
+@end deftypefun
+
+@deftypefun void gsl_multiset_free (gsl_multiset * @var{c})
+This function frees all the memory used by the multiset @var{c}.
+@end deftypefun
+
+@deftypefun int gsl_multiset_memcpy (gsl_multiset * @var{dest}, const gsl_multiset * @var{src})
+This function copies the elements of the multiset @var{src} into the
+multiset @var{dest}.  The two multisets must have the same size.
+@end deftypefun
+
+
+@node Accessing multiset elements
+@section Accessing multiset elements
+
+The following function can be used to access the elements of a multiset.
+
+@deftypefun size_t gsl_multiset_get (const gsl_multiset * @var{c}, const size_t @var{i})
+This function returns the value of the @var{i}-th element of the
+multiset @var{c}.  If @var{i} lies outside the allowed range of 0 to
+@math{@var{k}-1} then the error handler is invoked and 0 is returned.  @inlinefn{}
+@end deftypefun
+
+@node Multiset properties
+@section Multiset properties
+
+@deftypefun size_t gsl_multiset_n (const gsl_multiset * @var{c})
+This function returns the range (@math{n}) of the multiset @var{c}.
+@end deftypefun
+
+@deftypefun size_t gsl_multiset_k (const gsl_multiset * @var{c})
+This function returns the number of elements (@math{k}) in the multiset @var{c}.
+@end deftypefun
+
+@deftypefun {size_t *} gsl_multiset_data (const gsl_multiset * @var{c})
+This function returns a pointer to the array of elements in the
+multiset @var{c}.
+@end deftypefun
+
+@deftypefun int gsl_multiset_valid (gsl_multiset * @var{c})
+@cindex checking multiset for validity
+@cindex testing multiset for validity
+This function checks that the multiset @var{c} is valid.  The @var{k}
+elements should lie in the range 0 to @math{@var{n}-1}, with each
+value occurring in nondecreasing order.
+@end deftypefun
+
+@node Multiset functions
+@section Multiset functions
+
+@deftypefun int gsl_multiset_next (gsl_multiset * @var{c})
+@cindex iterating through multisets
+This function advances the multiset @var{c} to the next multiset element in
+lexicographic order and returns @code{GSL_SUCCESS}.  If no further multisets
+elements are available it returns @code{GSL_FAILURE} and leaves @var{c}
+unmodified.  Starting with the first multiset and repeatedly applying this
+function will iterate through all possible multisets of a given order.
+@end deftypefun
+
+@deftypefun int gsl_multiset_prev (gsl_multiset * @var{c})
+This function steps backwards from the multiset @var{c} to the previous
+multiset element in lexicographic order, returning @code{GSL_SUCCESS}.  If no
+previous multiset is available it returns @code{GSL_FAILURE} and leaves @var{c}
+unmodified.
+@end deftypefun
+
+@node Reading and writing multisets
+@section Reading and writing multisets
+
+The library provides functions for reading and writing multisets to a
+file as binary data or formatted text.
+
+@deftypefun int gsl_multiset_fwrite (FILE * @var{stream}, const gsl_multiset * @var{c})
+This function writes the elements of the multiset @var{c} to the
+stream @var{stream} in binary format.  The function returns
+@code{GSL_EFAILED} if there was a problem writing to the file.  Since the
+data is written in the native binary format it may not be portable
+between different architectures.
+@end deftypefun
+
+@deftypefun int gsl_multiset_fread (FILE * @var{stream}, gsl_multiset * @var{c})
+This function reads elements from the open stream @var{stream} into the
+multiset @var{c} in binary format.  The multiset @var{c} must be
+preallocated with correct values of @math{n} and @math{k} since the
+function uses the size of @var{c} to determine how many bytes to read.
+The function returns @code{GSL_EFAILED} if there was a problem reading
+from the file.  The data is assumed to have been written in the native
+binary format on the same architecture.
+@end deftypefun
+
+@deftypefun int gsl_multiset_fprintf (FILE * @var{stream}, const gsl_multiset * @var{c}, const char * @var{format})
+This function writes the elements of the multiset @var{c}
+line-by-line to the stream @var{stream} using the format specifier
+@var{format}, which should be suitable for a type of @var{size_t}.
+In ISO C99 the type modifier @code{z} represents @code{size_t}, so
+@code{"%zu\n"} is a suitable format.@footnote{In versions of the
+GNU C library prior to the ISO C99 standard,
+the type modifier @code{Z} was used instead.}  The function returns
+@code{GSL_EFAILED} if there was a problem writing to the file.
+@end deftypefun
+
+@deftypefun int gsl_multiset_fscanf (FILE * @var{stream}, gsl_multiset * @var{c})
+This function reads formatted data from the stream @var{stream} into the
+multiset @var{c}.  The multiset @var{c} must be preallocated with
+correct values of @math{n} and @math{k} since the function uses the size of @var{c} to
+determine how many numbers to read.  The function returns
+@code{GSL_EFAILED} if there was a problem reading from the file.
+@end deftypefun
+
+
+@node Multiset Examples
+@section Examples
+The example program below prints all multisets elements containing the values
+@math{@{0,1,2,3@}} ordered by size.  Multiset elements of the same size are
+ordered lexicographically.
+
+@example
+@verbatiminclude examples/multiset.c
+@end example
+
+@noindent
+Here is the output from the program,
+
+@example
+$ ./a.out
+@verbatiminclude examples/multiset.out
+@end example
+
+@noindent
+All 70 multisets are generated and sorted lexicographically.
diff --git a/multiset/.cvsignore b/multiset/.cvsignore
new file mode 100644
index 0000000..52fdab8
--- /dev/null
+++ b/multiset/.cvsignore
@@ -0,0 +1,8 @@
+*.la
+*.lo
+.deps
+.libs
+demo
+Makefile
+Makefile.in
+test
diff --git a/multiset/ChangeLog b/multiset/ChangeLog
new file mode 100644
index 0000000..3289c0e
--- /dev/null
+++ b/multiset/ChangeLog
@@ -0,0 +1,3 @@
+2009-09-21  Rhys Ulerich  <rhys.ulerich@gmail.com>
+
+	* added multiset support to GSL based on combinations
diff --git a/multiset/Makefile.am b/multiset/Makefile.am
new file mode 100644
index 0000000..2408d3f
--- /dev/null
+++ b/multiset/Makefile.am
@@ -0,0 +1,21 @@
+noinst_LTLIBRARIES = libgslmultiset.la
+
+pkginclude_HEADERS = gsl_multiset.h
+
+INCLUDES = -I$(top_srcdir)
+
+libgslmultiset_la_SOURCES = init.c file.c multiset.c inline.c
+
+noinst_HEADERS =
+
+TESTS = $(check_PROGRAMS)
+
+check_PROGRAMS = test
+
+test_SOURCES = test.c
+
+test_LDADD = libgslmultiset.la ../vector/libgslvector.la  ../block/libgslblock.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la
+
+noinst_PROGRAMS = demo
+demo_SOURCES = demo.c
+demo_LDADD = libgslmultiset.la ../vector/libgslvector.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la
diff --git a/multiset/demo.c b/multiset/demo.c
new file mode 100644
index 0000000..3c4749a
--- /dev/null
+++ b/multiset/demo.c
@@ -0,0 +1,72 @@
+#include <stdio.h>
+#include <gsl/gsl_multiset.h>
+
+void
+print_all_multisets ( size_t n, size_t k );
+
+int
+main (void)
+{
+  gsl_multiset * c;
+  size_t i;
+
+  printf("all multisets of {0,1,2,3} by size (lex. order)\n") ;
+  for(i = 0; i <= 4; i++)
+    {
+      c = gsl_multiset_calloc (4, i);
+      do
+        {
+          printf("{");
+          gsl_multiset_fprintf (stdout, c, " %u");
+          printf(" }\n");
+        }
+      while (gsl_multiset_next(c) == GSL_SUCCESS);
+      gsl_multiset_free(c);
+    }
+  printf("all multisets of {1,2,3,4} by size (reverse lex. order)\n") ;
+  for(i = 0; i <= 4; i++)
+    {
+      c = gsl_multiset_alloc (4, i) ;
+      gsl_multiset_init_last(c) ;
+      do
+        {
+          printf("{");
+          gsl_multiset_fprintf (stdout, c, " %u");
+          printf(" }\n");
+        }
+      while (gsl_multiset_prev(c) == GSL_SUCCESS);
+      gsl_multiset_free(c);
+    }
+  printf("\n");
+
+  print_all_multisets(5, 3);
+  print_all_multisets(5, 0);
+  print_all_multisets(5, 5);
+  print_all_multisets(1, 1);
+  print_all_multisets(3, 1);
+
+  return 0;
+}
+
+void
+print_all_multisets (size_t n, size_t k)
+{
+  gsl_multiset * c = gsl_multiset_calloc (n, k);
+
+  printf("multisets %u choose %u (with replacement)\n", n, k);
+  do
+    {
+      gsl_multiset_fprintf (stdout, c, " %u");
+      printf("\n");
+    }
+  while (gsl_multiset_next(c) == GSL_SUCCESS);
+  while (gsl_multiset_next(c) == GSL_SUCCESS);
+  do
+    {
+      gsl_multiset_fprintf (stdout, c, " %u");
+      printf("\n");
+    }
+  while (gsl_multiset_prev(c) == GSL_SUCCESS);
+  printf("\n");
+  gsl_multiset_free(c);
+}
diff --git a/multiset/file.c b/multiset/file.c
new file mode 100644
index 0000000..97be1f1
--- /dev/null
+++ b/multiset/file.c
@@ -0,0 +1,115 @@
+/* multiset/file.c
+ * based on combination/file.c by Szymon Jaroszewicz
+ * based on permutation/file.c by Brian Gough
+ *
+ * Copyright (C) 2001 Szymon Jaroszewicz
+ * Copyright (C) 2009 Rhys Ulerich
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <config.h>
+#include <stdio.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_multiset.h>
+
+#define IN_FORMAT "%lu"
+
+int
+gsl_multiset_fread (FILE * stream, gsl_multiset * c)
+{
+  size_t k = c->k ;
+
+  size_t * data = c->data ;
+
+  size_t items = fread (data, sizeof (size_t), k, stream);
+
+  if (items != k)
+    {
+      GSL_ERROR ("fread failed", GSL_EFAILED);
+    }
+
+  return GSL_SUCCESS;
+}
+
+int
+gsl_multiset_fwrite (FILE * stream, const gsl_multiset * c)
+{
+  size_t k = c->k ;
+
+  size_t * data = c->data ;
+
+  size_t items = fwrite (data, sizeof (size_t), k, stream);
+
+  if (items != k)
+    {
+      GSL_ERROR ("fwrite failed", GSL_EFAILED);
+    }
+
+  return GSL_SUCCESS;
+}
+
+int
+gsl_multiset_fprintf (FILE * stream, const gsl_multiset * c, const char *format)
+{
+  size_t k = c->k ;
+
+  size_t * data = c->data ;
+
+  size_t i;
+
+  for (i = 0; i < k; i++)
+    {
+      int status = fprintf (stream, format, data[i]);
+
+      if (status < 0)
+        {
+          GSL_ERROR ("fprintf failed", GSL_EFAILED);
+        }
+    }
+
+  return GSL_SUCCESS;
+}
+
+int
+gsl_multiset_fscanf (FILE * stream, gsl_multiset * c)
+{
+  size_t k = c->k ;
+
+  size_t * data = c->data ;
+
+  size_t i;
+
+  for (i = 0; i < k; i++)
+    {
+      unsigned long j ;
+
+      /* FIXME: what if size_t != unsigned long ???
+
+         want read in size_t but have to read in unsigned long to avoid
+         error from compiler */
+
+      int status = fscanf (stream, IN_FORMAT, &j);
+
+      if (status != 1)
+        {
+          GSL_ERROR ("fscanf failed", GSL_EFAILED);
+        }
+
+      data[i] = j;
+    }
+
+  return GSL_SUCCESS;
+}
diff --git a/multiset/gsl_multiset.h b/multiset/gsl_multiset.h
new file mode 100644
index 0000000..1ba5d6c
--- /dev/null
+++ b/multiset/gsl_multiset.h
@@ -0,0 +1,93 @@
+/* multiset/gsl_multiset.h
+ * based on combination/gsl_combination.h by Szymon Jaroszewicz
+ * based on permutation/gsl_permutation.h by Brian Gough
+ *
+ * Copyright (C) 2009 Rhys Ulerich
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __GSL_MULTISET_H__
+#define __GSL_MULTISET_H__
+
+#include <stdlib.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_types.h>
+#include <gsl/gsl_inline.h>
+#include <gsl/gsl_check_range.h>
+
+#undef __BEGIN_DECLS
+#undef __END_DECLS
+#ifdef __cplusplus
+# define __BEGIN_DECLS extern "C" {
+# define __END_DECLS }
+#else
+# define __BEGIN_DECLS /* empty */
+# define __END_DECLS /* empty */
+#endif
+
+__BEGIN_DECLS
+
+struct gsl_multiset_struct
+{
+  size_t n;
+  size_t k;
+  size_t *data;
+};
+
+typedef struct gsl_multiset_struct gsl_multiset;
+
+gsl_multiset *gsl_multiset_alloc (const size_t n, const size_t k);
+gsl_multiset *gsl_multiset_calloc (const size_t n, const size_t k);
+void gsl_multiset_init_first (gsl_multiset * c);
+void gsl_multiset_init_last (gsl_multiset * c);
+void gsl_multiset_free (gsl_multiset * c);
+int gsl_multiset_memcpy (gsl_multiset * dest, const gsl_multiset * src);
+
+int gsl_multiset_fread (FILE * stream, gsl_multiset * c);
+int gsl_multiset_fwrite (FILE * stream, const gsl_multiset * c);
+int gsl_multiset_fscanf (FILE * stream, gsl_multiset * c);
+int gsl_multiset_fprintf (FILE * stream, const gsl_multiset * c, const char *format);
+
+size_t gsl_multiset_n (const gsl_multiset * c);
+size_t gsl_multiset_k (const gsl_multiset * c);
+size_t * gsl_multiset_data (const gsl_multiset * c);
+
+int gsl_multiset_valid (gsl_multiset * c);
+int gsl_multiset_next (gsl_multiset * c);
+int gsl_multiset_prev (gsl_multiset * c);
+
+INLINE_DECL size_t gsl_multiset_get (const gsl_multiset * c, const size_t i);
+
+#ifdef HAVE_INLINE
+
+INLINE_FUN
+size_t
+gsl_multiset_get (const gsl_multiset * c, const size_t i)
+{
+#if GSL_RANGE_CHECK
+  if (GSL_RANGE_COND(i >= c->k)) /* size_t is unsigned, can't be negative */
+    {
+      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
+    }
+#endif
+  return c->data[i];
+}
+
+#endif /* HAVE_INLINE */
+
+__END_DECLS
+
+#endif /* __GSL_MULTISET_H__ */
diff --git a/multiset/init.c b/multiset/init.c
new file mode 100644
index 0000000..b71b736
--- /dev/null
+++ b/multiset/init.c
@@ -0,0 +1,130 @@
+/* multiset/init.c
+ * based on combination/init.c by Szymon Jaroszewicz
+ * based on permutation/init.c by Brian Gough
+ *
+ * Copyright (C) 2001 Szymon Jaroszewicz
+ * Copyright (C) 2009 Brian Gough
+ * Copyright (C) 2009 Rhys Ulerich
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <config.h>
+#include <stdlib.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_multiset.h>
+
+gsl_multiset *
+gsl_multiset_alloc (const size_t n, const size_t k)
+{
+  gsl_multiset * c;
+
+  if (n == 0)
+    {
+      GSL_ERROR_VAL ("multiset parameter n must be positive integer",
+                        GSL_EDOM, 0);
+    }
+  if (k > n)
+    {
+      GSL_ERROR_VAL ("multiset length k must be an integer less than or equal to n",
+                        GSL_EDOM, 0);
+    }
+  c = (gsl_multiset *) malloc (sizeof (gsl_multiset));
+
+  if (c == 0)
+    {
+      GSL_ERROR_VAL ("failed to allocate space for multiset struct",
+                        GSL_ENOMEM, 0);
+    }
+
+  if (k > 0)
+    {
+      c->data = (size_t *) malloc (k * sizeof (size_t));
+
+      if (c->data == 0)
+        {
+          free (c);             /* exception in constructor, avoid memory leak */
+
+          GSL_ERROR_VAL ("failed to allocate space for multiset data",
+                         GSL_ENOMEM, 0);
+        }
+    }
+  else
+    {
+      c->data = 0;
+    }
+
+  c->n = n;
+  c->k = k;
+
+  return c;
+}
+
+gsl_multiset *
+gsl_multiset_calloc (const size_t n, const size_t k)
+{
+  size_t i;
+
+  gsl_multiset * c =  gsl_multiset_alloc (n, k);
+
+  if (c == 0)
+    return 0;
+
+  /* initialize multiset to repeated first element */
+
+  for (i = 0; i < k; i++)
+    {
+      c->data[i] = 0;
+    }
+
+  return c;
+}
+
+void
+gsl_multiset_init_first (gsl_multiset * c)
+{
+  const size_t k = c->k ;
+  size_t i;
+
+  /* initialize multiset to repeated first element */
+
+  for (i = 0; i < k; i++)
+    {
+      c->data[i] = 0;
+    }
+}
+
+void
+gsl_multiset_init_last (gsl_multiset * c)
+{
+  const size_t k = c->k ;
+  size_t i;
+  size_t n = c->n;
+
+  /* initialize multiset to repeated last element */
+
+  for (i = 0; i < k; i++)
+    {
+      c->data[i] = n - 1;
+    }
+}
+
+void
+gsl_multiset_free (gsl_multiset * c)
+{
+  RETURN_IF_NULL (c);
+  if (c->k > 0) free (c->data);
+  free (c);
+}
diff --git a/multiset/inline.c b/multiset/inline.c
new file mode 100644
index 0000000..f59e9f3
--- /dev/null
+++ b/multiset/inline.c
@@ -0,0 +1,24 @@
+/* multiset/inline.c
+ *
+ * Copyright (C) 2008 Brian Gough
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* Compile all the inline functions */
+
+#define COMPILE_INLINE_STATIC
+#include "build.h"
+#include <gsl/gsl_multiset.h>
diff --git a/multiset/multiset.c b/multiset/multiset.c
new file mode 100644
index 0000000..9b1e145
--- /dev/null
+++ b/multiset/multiset.c
@@ -0,0 +1,181 @@
+/* multiset/multiset.c
+ * based on combination/combination.c by Szymon Jaroszewicz
+ * based on permutation/permutation.c by Brian Gough
+ *
+ * Copyright (C) 2001 Szymon Jaroszewicz
+ * Copyright (C) 2009 Rhys Ulerich
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <config.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_multiset.h>
+
+size_t
+gsl_multiset_n (const gsl_multiset * c)
+{
+  return c->n ;
+}
+
+size_t
+gsl_multiset_k (const gsl_multiset * c)
+{
+  return c->k ;
+}
+
+size_t *
+gsl_multiset_data (const gsl_multiset * c)
+{
+  return c->data ;
+}
+
+int
+gsl_multiset_valid (gsl_multiset * c)
+{
+  const size_t n = c->n ;
+  const size_t k = c->k ;
+
+  size_t i, j ;
+
+  if( k > n )
+    {
+      GSL_ERROR("multiset has k greater than n", GSL_FAILURE) ;
+    }
+  for (i = 0; i < k; i++)
+    {
+      const size_t ci = c->data[i];
+
+      if (ci >= n)
+        {
+          GSL_ERROR("multiset index outside range", GSL_FAILURE) ;
+        }
+
+      for (j = 0; j < i; j++)
+        {
+          if (c->data[j] > ci)
+            {
+              GSL_ERROR("multiset indices not in increasing order",
+                        GSL_FAILURE) ;
+            }
+        }
+    }
+
+  return GSL_SUCCESS;
+}
+
+
+int
+gsl_multiset_next (gsl_multiset * c)
+{
+  /* Replaces c with the next multiset (in the standard lexicographical
+   * ordering).  Returns GSL_FAILURE if there is no next multiset.
+   */
+  const size_t n = c->n;
+  const size_t k = c->k;
+  size_t *data = c->data;
+  size_t i;
+
+  if(k == 0)
+    {
+      return GSL_FAILURE;
+    }
+  i = k - 1;
+
+  while(i > 0 && data[i] == n-1)
+    {
+      --i;
+    }
+
+  if (i == 0 && data[0] == n-1)
+    {
+      return GSL_FAILURE;
+    }
+
+  ++data[i];
+
+  while(i < k-1)
+    {
+      data[i+1] = data[i];
+      ++i;
+    }
+
+  return GSL_SUCCESS;
+}
+
+int
+gsl_multiset_prev (gsl_multiset * c)
+{
+  /* Replaces c with the previous multiset (in the standard
+   * lexicographical ordering).  Returns GSL_FAILURE if there is no
+   * previous multiset.
+   */
+  const size_t n = c->n;
+  const size_t k = c->k;
+  size_t *data = c->data;
+  size_t i;
+
+  if(k == 0)
+    {
+      return GSL_FAILURE;
+    }
+  i = k - 1;
+
+  while(i > 0 && data[i-1] == data[i])
+    {
+      --i;
+    }
+
+  if(i == 0 && data[i] == 0)
+    {
+      return GSL_FAILURE;
+    }
+
+  data[i]--;
+
+  if (data[i] < n-1)
+    {
+      while (i < k-1) {
+        data[++i] = n - 1;
+      }
+    }
+
+  return GSL_SUCCESS;
+}
+
+int
+gsl_multiset_memcpy (gsl_multiset * dest, const gsl_multiset * src)
+{
+   const size_t src_n = src->n;
+   const size_t src_k = src->k;
+   const size_t dest_n = dest->n;
+   const size_t dest_k = dest->k;
+
+   if (src_n != dest_n || src_k != dest_k)
+     {
+       GSL_ERROR ("multiset lengths are not equal", GSL_EBADLEN);
+     }
+
+   {
+     size_t j;
+
+     for (j = 0; j < src_k; j++)
+       {
+         dest->data[j] = src->data[j];
+       }
+   }
+
+   return GSL_SUCCESS;
+}
diff --git a/multiset/test.c b/multiset/test.c
new file mode 100644
index 0000000..4548554
--- /dev/null
+++ b/multiset/test.c
@@ -0,0 +1,278 @@
+/* multiset/test.c
+ * based on combination/test.c by Szymon Jaroszewicz
+ * based on permutation/test.c by Brian Gough
+ *
+ * Copyright (C) 2001 Szymon Jaroszewicz
+ * Copyright (C) 2009 Rhys Ulerich
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <config.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <gsl/gsl_multiset.h>
+#include <gsl/gsl_test.h>
+#include <gsl/gsl_ieee_utils.h>
+
+size_t c63[56][3] = {
+  { 0,0,0 },  { 0,0,1 },  { 0,0,2 },  { 0,0,3 }, { 0,0,4 },  { 0,0,5 },
+              { 0,1,1 },  { 0,1,2 },  { 0,1,3 }, { 0,1,4 },  { 0,1,5 },
+                          { 0,2,2 },  { 0,2,3 }, { 0,2,4 },  { 0,2,5 },
+                                      { 0,3,3 }, { 0,3,4 },  { 0,3,5 },
+                                                 { 0,4,4 },  { 0,4,5 },
+                                                             { 0,5,5 },
+              { 1,1,1 },  { 1,1,2 },  { 1,1,3 }, { 1,1,4 },  { 1,1,5 },
+                          { 1,2,2 },  { 1,2,3 }, { 1,2,4 },  { 1,2,5 },
+                                      { 1,3,3 }, { 1,3,4 },  { 1,3,5 },
+                                                 { 1,4,4 },  { 1,4,5 },
+                                                             { 1,5,5 },
+                          { 2,2,2 },  { 2,2,3 }, { 2,2,4 },  { 2,2,5 },
+                                      { 2,3,3 }, { 2,3,4 },  { 2,3,5 },
+                                                 { 2,4,4 },  { 2,4,5 },
+                                                             { 2,5,5 },
+                                      { 3,3,3 }, { 3,3,4 },  { 3,3,5 },
+                                                 { 3,4,4 },  { 3,4,5 },
+                                                             { 3,5,5 },
+                                                 { 4,4,4 },  { 4,4,5 },
+                                                             { 4,5,5 },
+                                                             { 5,5,5 }
+} ;
+
+void my_error_handler (const char *reason, const char *file, int line, int err);
+
+int
+main (void)
+{
+  size_t i, j;
+  int status = 0, s;
+  gsl_multiset * c ;
+
+  gsl_ieee_env_setup ();
+
+  c = gsl_multiset_alloc (6,3);
+
+  /* Test multisets in forward order */
+
+  gsl_multiset_init_first (c);
+
+  i = 0;
+
+  do
+    {
+      if ( i >= 56 )
+        {
+          status = 1;
+          break;
+        }
+      for (j = 0; j < 3; j++)
+        {
+          status |= (c->data[j] != c63[i][j]);
+        }
+
+      {
+        int s1 = gsl_multiset_valid (c);
+        gsl_test (s1, "gsl_multiset_valid (%u)", i);
+      }
+
+      i++;
+    }
+  while (gsl_multiset_next(c) == GSL_SUCCESS);
+
+  gsl_test(status, "gsl_multiset_next, 6 choose 3 multiset, 56 steps");
+
+  gsl_multiset_next(c);
+  gsl_multiset_next(c);
+  gsl_multiset_next(c);
+  for (j = 0; j < 3; j++)
+    {
+      status |= (c->data[j] != c63[55][j]);
+    }
+  gsl_test(status, "gsl_multiset_next on the last multiset");
+
+  {
+    int s1 = gsl_multiset_valid (c);
+    gsl_test (s1, "gsl_multiset_valid on the last multiset");
+  }
+
+  {
+    gsl_multiset * d = gsl_multiset_alloc (6,3);
+    gsl_multiset_memcpy (d, c);
+
+    status = 0;
+
+    for (j = 0; j < 3; j++)
+      {
+        status |= (d->data[j] != c->data[j]);
+      }
+
+    gsl_test (status, "gsl_multiset_memcpy, 6 choose 3 multiset");
+    gsl_multiset_free(d);
+  }
+
+
+  /* Now test multisets in reverse order */
+
+  gsl_multiset_init_last (c);
+
+  i = 56;
+  do
+    {
+      if ( i == 0 )
+        {
+          status = 1;
+          break;
+        }
+
+      i--;
+
+      for (j = 0; j < 3; j++)
+        {
+          status |= (c->data[j] != c63[i][j]);
+        }
+
+      {
+        int s1 = gsl_multiset_valid (c);
+        gsl_test (s1, "gsl_multiset_valid (%u)", i);
+      }
+    }
+  while (gsl_multiset_prev(c) == GSL_SUCCESS);
+
+  gsl_test(status, "gsl_multiset_prev, 6 choose 3 multiset, 20 steps");
+
+  gsl_multiset_prev(c);
+  gsl_multiset_prev(c);
+  gsl_multiset_prev(c);
+  for (j = 0; j < 3; j++)
+    {
+      status |= (c->data[j] != c63[0][j]);
+    }
+  gsl_test(status, "gsl_multiset_prev on the first multiset");
+
+  {
+    int s1 = gsl_multiset_valid (c);
+    gsl_test (s1, "gsl_multiset_valid on the first multiset");
+  }
+
+  {
+    gsl_multiset * d = gsl_multiset_alloc (6,3);
+    gsl_multiset_memcpy (d, c);
+
+    status = 0;
+
+    for (j = 0; j < 3; j++)
+      {
+        status |= (d->data[j] != c->data[j]);
+      }
+
+    gsl_test (status, "gsl_multiset_memcpy, 6 choose 3 multiset");
+    gsl_multiset_free(d);
+  }
+
+  gsl_multiset_free (c);
+
+  c = gsl_multiset_calloc(7, 0);
+  /* should return GSL_FAILURE every time */
+  status |= (gsl_multiset_next(c) != GSL_FAILURE);
+  status |= (gsl_multiset_next(c) != GSL_FAILURE);
+  status |= (gsl_multiset_prev(c) != GSL_FAILURE);
+  status |= (gsl_multiset_prev(c) != GSL_FAILURE);
+  gsl_test(status, "gsl_multiset 7 choose 0");
+  gsl_multiset_free (c);
+
+  c = gsl_multiset_calloc(1, 1);
+  /* should return GSL_FAILURE every time */
+  for(j = 0; j < 1; j++)
+  {
+    status |= (gsl_multiset_get(c, j) != j);
+  }
+  status |= (gsl_multiset_next(c) != GSL_FAILURE);
+  for(j = 0; j < 1; j++)
+  {
+    status |= (gsl_multiset_get(c, j) != j);
+  }
+  status |= (gsl_multiset_next(c) != GSL_FAILURE);
+  for(j = 0; j < 1; j++)
+  {
+    status |= (gsl_multiset_get(c, j) != j);
+  }
+  status |= (gsl_multiset_prev(c) != GSL_FAILURE);
+  for(j = 0; j < 1; j++)
+  {
+    status |= (gsl_multiset_get(c, j) != j);
+  }
+  status |= (gsl_multiset_prev(c) != GSL_FAILURE);
+  for(j = 0; j < 1; j++)
+  {
+    status |= (gsl_multiset_get(c, j) != j);
+  }
+  gsl_test(status, "gsl_multiset 7 choose 7");
+  gsl_multiset_free (c);
+
+  c = gsl_multiset_calloc(6, 3);
+
+  gsl_set_error_handler (&my_error_handler);
+
+  c->data[0] = 1;
+  c->data[1] = 2;
+  c->data[2] = 1;
+  s = gsl_multiset_valid (c);
+  gsl_test (!s, "gsl_multiset_valid on an invalid multiset (1,1,2)");
+
+  c->data[0] = 2;
+  c->data[1] = 1;
+  c->data[2] = 0;
+  s = gsl_multiset_valid (c);
+  gsl_test (!s, "gsl_multiset_valid on an invalid multiset (2,1,0)");
+
+  c->data[0] = 1;
+  c->data[1] = 2;
+  c->data[2] = 0;
+  s = gsl_multiset_valid (c);
+  gsl_test (!s, "gsl_multiset_valid on an invalid multiset (1,2,0)");
+
+  {
+    gsl_multiset * d = gsl_multiset_alloc (6,4);
+    int s = gsl_multiset_memcpy (d, c);
+    gsl_test (!s, "gsl_multiset_memcpy, (6,4) vs (6,3)");
+    gsl_multiset_free(d);
+  }
+
+  {
+    gsl_multiset * d = gsl_multiset_alloc (7,3);
+    int s = gsl_multiset_memcpy (d, c);
+    gsl_test (!s, "gsl_multiset_memcpy, (7,3) vs (6,3)");
+    gsl_multiset_free(d);
+  }
+
+  {
+    gsl_multiset * d = gsl_multiset_alloc (7,2);
+    int s = gsl_multiset_memcpy (d, c);
+    gsl_test (!s, "gsl_multiset_memcpy, (7,2) vs (6,3)");
+    gsl_multiset_free(d);
+  }
+
+
+  gsl_multiset_free (c);
+
+  exit (gsl_test_summary());
+}
+
+void
+my_error_handler (const char *reason, const char *file, int line, int err)
+{
+  if (0) printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err) ;
+}
-- 
1.5.4.3


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]