← Back to branch summary

~ubuntu-branches/ubuntu/wily/octave-miscellaneous/wily

« back to all changes in this revision

Viewing changes to src/partint.cc

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2012-10-17 13:40:55 UTC
  • mfrom: (1.1.6)
  • mto: (5.1.2 experimental)
  • mto: This revision was merged to the branch mainline in revision 13.
  • Revision ID: package-import@ubuntu.com-20121017134055-e8lrxjd3qgcd3kmt
Tags: upstream-1.2.0
Import upstream version 1.2.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/partint.cc
 
1
// Copyright (C) 2006 Torsten Finke <fi@igh-essen.com>
 
2
//
 
3
// This program is free software; you can redistribute it and/or modify it under
 
4
// the terms of the GNU General Public License as published by the Free Software
 
5
// Foundation; either version 3 of the License, or (at your option) any later
 
6
// version.
 
7
//
 
8
// This program is distributed in the hope that it will be useful, but WITHOUT
 
9
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 
10
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
 
11
// details.
 
12
//
 
13
// You should have received a copy of the GNU General Public License along with
 
14
// this program; if not, see <http://www.gnu.org/licenses/>.
 
15
 
 
16
#include <octave/oct.h>
 
17
#include <octave/lo-ieee.h>
 
18
 
 
19
#include "partint.h"
 
20
 
 
21
unsigned long * pcnt(unsigned long n) 
 
22
{
 
23
    unsigned long *s = new unsigned long[n];
 
24
    unsigned long *x = new unsigned long[n*n];
 
25
    unsigned long **p = new unsigned long*[n];
 
26
    for (unsigned long k=0; k<n; ++k)  {
 
27
        p[k] = x + (k*n);
 
28
        s[k] = 0;
 
29
    }
 
30
    for (unsigned long k=0; k<n*n; ++k)  x[k] = 0;
 
31
    p[0][0] = 1;
 
32
    for (unsigned long k=1; k<n; ++k)
 
33
    {
 
34
        // p[N][j] == numpart of N with max summand j
 
35
        for (unsigned long j=1; j<=k; ++j) {
 
36
            p[k][j] = p[k-1][j-1] + p[k-j][j];
 
37
        }
 
38
        for (unsigned long j=1; j<=k; ++j) {
 
39
            s[k] += p[k][j]; // S(k) = numpart(n)
 
40
        }
 
41
    }
 
42
    return s;
 
43
}
 
44
 
 
45
DEFUN_DLD (partcnt, args, ,
 
46
"-*- texinfo -*-\n\
 
47
@deftypefn{Loadable Function} {@var{p} =} partcnt(@var{n})\n\
 
48
\n\
 
49
@cindex partition count\n\
 
50
\n\
 
51
Calculate integer partition count. Argument @var{n} be scalar, vector\n\
 
52
or matrix of non-negative numbers. A partition is the sum decomposition \n\
 
53
of integers. \n\
 
54
\n\
 
55
Example:\n\
 
56
@example \n\
 
57
partcnt([1, 5; 17 -3])\n\
 
58
@end example\n\
 
59
@noindent\n\
 
60
returns\n\
 
61
@example\n\
 
62
ans =\n\
 
63
    1     7\n\
 
64
  297   NaN\n\
 
65
@end example\n\
 
66
\n\
 
67
Reference:\n\
 
68
Joerg Arndt: Algorithms for programmers (http://www.jjj.de), 2006.\n\n\
 
69
@end deftypefn\n\
 
70
@seealso{partint}\n\
 
71
")
 
72
{
 
73
    octave_value r;
 
74
    
 
75
    int nargin = args.length ();
 
76
    if (nargin != 1) {
 
77
        error("partcnt accepts exactly one argument");
 
78
        return r; 
 
79
    }
 
80
    if ( ! args(0).is_numeric_type()) {
 
81
        error("partcnt only accepts a numeric argument");
 
82
        return r;
 
83
    }
 
84
 
 
85
    NDArray f(args(0).matrix_value());
 
86
    RowVector m(f.max());
 
87
    double mmax = m.max();
 
88
    if ( mmax < 1 ) {
 
89
        error("partcnt is only defined for non-negative arguments");
 
90
        return r;
 
91
    }
 
92
 
 
93
    unsigned long n = (unsigned long) mmax + 1;
 
94
    unsigned long *s = pcnt(n);
 
95
    unsigned long fr = (unsigned long) f.rows();
 
96
    unsigned long fc = (unsigned long) f.columns();
 
97
    for (unsigned long i=0; i<fr; i++) {
 
98
        for (unsigned long k=0; k<fc; k++) {
 
99
            unsigned long idx = (unsigned long) f(i, k);
 
100
            if (0 < idx && idx < n) {
 
101
                f(i, k) = s[idx];
 
102
            } else {
 
103
                f(i, k) = lo_ieee_nan_value();
 
104
            }
 
105
        }
 
106
    }
 
107
    r = f;
 
108
    return r;
 
109
}
 
110
 
 
111
/*
 
112
 
 
113
%!assert(partcnt(1), 1);
 
114
%!assert(partcnt(17), 297);
 
115
%!fail("partcnt()", "partcnt");
 
116
%!fail("partcnt(1,2)", "partcnt");
 
117
%!fail("partcnt('xyz')", "partcnt");
 
118
%!demo
 
119
%! p = partcnt([1, 5; 17 -5])
 
120
 
 
121
*/
 
122
 
 
123
DEFUN_DLD (partint, args, ,
 
124
"-*- texinfo -*-\n\
 
125
@deftypefn{Loadable Function} {@var{p} =} partint(@var{n})\n\
 
126
\n\
 
127
@cindex partition\n\
 
128
\n\
 
129
Calculate all integer partitions. Argument @var{n} be \n\
 
130
a positive number. A partition is the sum decomposition \n\
 
131
of integers. \n\
 
132
\n\
 
133
Example:\n\
 
134
@example \n\
 
135
partint(4)\n\
 
136
@end example\n\
 
137
@noindent\n\
 
138
returns\n\
 
139
@example\n\
 
140
ans =\n\
 
141
  4  0  0  0\n\
 
142
  2  1  0  0\n\
 
143
  0  2  0  0\n\
 
144
  1  0  1  0\n\
 
145
  0  0  0  1\n\
 
146
@end example\n\
 
147
\n\
 
148
This means\n\n\
 
149
@iftex\n\
 
150
@tex\n\
 
151
$$\eqalign{4 & = 4 \\cdot 1 \\cr\n\
 
152
 & = 2 \\cdot 1 + 1 \\cdot 2 \\cr\n\
 
153
 & = 2 \\cdot 2 \\cr\n\
 
154
 & = 1 \\cdot 1 + 1 \\cdot 3 \\cr\n\
 
155
 & = 1 \\cdot 1 \\cr\n\
 
156
\\cr}$$\n\
 
157
@end tex\n\
 
158
@end iftex\n\
 
159
@ifinfo\n\
 
160
@example\n\
 
161
4 = 4 * 1\n\
 
162
  = 2 * 1 + 1 * 2\n\
 
163
  =         2 * 2\n\
 
164
  = 1 * 1         + 1 * 3\n\
 
165
  =                 1 * 4\n\
 
166
@end example\n\
 
167
@end ifinfo\n\
 
168
\n\
 
169
Note:\n\
 
170
\n\
 
171
partint(n) * [1:n]' == n\n\
 
172
\n\
 
173
Reference:\n\
 
174
Joerg Arndt: Algorithms for programmers (http://www.jjj.de), 2006.\n\n\
 
175
@end deftypefn\n\
 
176
@seealso{partcnt}\n\
 
177
")
 
178
{
 
179
    octave_value r;
 
180
    
 
181
    int nargin = args.length ();
 
182
    if (nargin != 1 || 
 
183
        ! args(0).is_scalar_type() ||
 
184
        ! args(0).is_numeric_type()
 
185
        ) {
 
186
        error("partint only accepts one scalar positive integer argument");
 
187
        return r;
 
188
    }
 
189
    double arg0 = args(0).double_value();
 
190
    if ( arg0 < 1 ) {
 
191
        error("partint is only defined for positive integer arguments");
 
192
        return r;
 
193
    }
 
194
 
 
195
    unsigned long n = (unsigned long) arg0;
 
196
    unsigned long *s = pcnt(n+1);
 
197
    unsigned long k = s[n];
 
198
    Matrix pa(k, n, 0);
 
199
    int_partition p(n);
 
200
    unsigned long i = 0;
 
201
    do {
 
202
        const unsigned long *d = p.data();
 
203
        for (unsigned long j=0; j<n; j++) {
 
204
            pa(i, j) = (unsigned long)d[j+1];
 
205
        }
 
206
        i ++;
 
207
    } while (p.next());
 
208
    r = pa;
 
209
    return r;
 
210
}
 
211
 
 
212
/*
 
213
 
 
214
%!assert(partint(1), 1);
 
215
%!assert(all(partint(n=17) * [1:n]' == n) - 1, 0); 
 
216
%!test
 
217
%! expected = [4,0,0,0; 2,1,0,0; 0,2,0,0; 1,0,1,0; 0,0,0,1];
 
218
%! assert(partint(4), expected);
 
219
%!fail("partint()", "partint");
 
220
%!fail("partint(1,2)", "partint");
 
221
%!fail("partint('xyz')", "partint");
 
222
%!demo
 
223
%! p = partint(4)
 
224
 
 
225
*/