1
2 #include "comb/partition-asc-subset-lex.h"
3 #include "ds/bitarray.h"
4
5 //#define PART_DIST // for A111212
6 #ifdef PART_DIST
7 #include "comb/partition-dist-asc-subset-lex.h"
8 #endif
9
10 #include "comb/comb-print.h"
11
12 #include "fxtio.h"
13
14 #include "jjassert.h"
15
16 #include "fxttypes.h"
17 #include "nextarg.h"
18
19
20 //% OEIS sequence A069999:
21 //% Number of possible dimensions for commutators of n X n matrices.
22 //% Also sequence A111212.
23
24 // Cf. seq/A069999-2-demo.cc for computation using convolutions of sets
25 // Cf. comb/partition-asc-subset-lex-demo.cc
26 // Cf. comb/partition-dist-asc-subset-lex-demo.cc
27
28 //#define TIMING // uncomment to disable printing
29
30 class one_term
31 {
32 private:
33 #ifdef PART_DIST
34 partition_dist_asc_subset_lex P; // partitions
35 #else
36 partition_asc_subset_lex P; // partitions
37 #endif
38 bitarray B; // sums seen
39 ulong *S; // look up table for squares
40 ulong *T; // partial sums of squares
41
42 public:
one_term(ulong n)43 explicit one_term(ulong n)
44 : P(n), B(n*n+1), S(new ulong[n+1]), T(new ulong[n+1])
45 {
46 for (ulong j=0; j<=n; ++j) S[j] = j * j;
47 for (ulong j=0; j<=n; ++j) T[j] = 0;
48 }
49
~one_term()50 ~one_term()
51 {
52 delete [] S;
53 delete [] T;
54 }
55
56 #ifdef TIMING
sum_sqr()57 void sum_sqr()
58 // with subset-lex order, only the last two parts
59 // differ from the previous partition.
60 {
61 const ulong *D = P.data() - 1; // including D[0] = 0
62 const ulong m = P.num_parts();
63
64 jjassert( m >= 2 );
65 ulong ss = T[m-2];
66
67 ss += S[ D[m-1] ];
68 T[m-1] = ss;
69
70 ss += S[ D[m] ];
71 T[m] = ss;
72
73 B.set( ss );
74
75 // P.print(""); cout << " ==> " << ss << endl;
76 }
77 #else // TIMING
sum_sqr()78 void sum_sqr()
79 {
80 const ulong *D = P.data();
81 const ulong m = P.num_parts();
82 P.print("");
83 cout << " ==> ";
84 ulong ss = 0;
85 for (ulong j=0; j<m; ++j)
86 {
87 ulong s = S[ D[j] ];
88 ss += s;
89 cout << " " << s;
90 if ( j < m - 1 ) cout << " +";
91 }
92 B.set( ss );
93 cout << " = " << ss << endl;
94 }
95 #endif // TIMING
96
count()97 ulong count()
98 {
99 B.clear_all();
100 ulong n = P.n_;
101 B.set( n*n ); // subset-lex order, for first (skipped) partition
102 P.first();
103 // skip first partition (into one part)
104 while ( P.next() )
105 {
106 sum_sqr();
107 }
108
109 return B.count_ones();
110 }
111 };
112 // -------------------------
113
114
115
116 int
main(int argc,char ** argv)117 main(int argc, char **argv)
118 {
119 #ifdef TIMING
120 ulong n1 = 0;
121 ulong n2 = 13;
122 NXARG(n1, "partitions into n parts starting from n1");
123 NXARG(n2, "partitions into n parts up to n2");
124
125 for (ulong j=n1; j<=n2; ++j)
126 {
127 one_term X(j);
128 ulong ct = X.count();
129 cout << j << " " << ct << endl; // b-file format
130 cout.flush();
131 }
132
133 #else // TIMING
134 ulong n = 7;
135 NXARG(n, "partitions of n");
136
137 one_term X(n);
138 ulong ct = X.count();
139 cout << " ct = " << ct << endl;
140
141 #endif // TIMING
142
143 return 0;
144 }
145 // -------------------------
146
147 /*
148 for n in $(seq 0 10); do ./bin $n | grep 'ct' ; done
149
150 Timing: Intel(R) Core(TM) i7-8700K CPU @ 3.70GHz
151 GCC 8.3.0
152
153 time ./bin 0 100
154 ./bin 0 100 5.83s user 0.00s system 99% cpu 5.834 total
155
156 time ./bin 0 1000 > ~/b069999-in.txt
157 ## up to n = 153:
158 ./bin 1000 > ~/b069999-in.txt 2281.29s user 0.00s system 99% cpu 38:01.30 total
159 N=153;x='x+O('x^(N+1)); Vec( 1/prod(k=1,N,1-x^k)/(1-x) )[N+1] ==> 573,169,704,246
160 573169704246 / 2281.29 == 251,248,067 per second
161
162 time ./bin 0 1000 > ~/b111212-in.txt
163 ## up to n = 300:
164 ./bin 0 1000 > ~/b111212-in.txt 7239.64s user 0.28s system 99% cpu 2:00:39.98 total
165 N=300; x='x+O('x^(N+1)); Vec( prod(k=1,N,1+x^k)/(1-x) )[N+1] ==> 2,287,746,908,056
166 2287746908056 / 7239.64 == 316,002,854 per second
167
168 */
169
170
171 /// Emacs:
172 /// Local Variables:
173 /// MyRelDir: "demo/seq"
174 /// makefile-dir: "../../"
175 /// make-target: "1demo DSRC=demo/seq/A069999-demo.cc"
176 /// make-target2: "1demo DSRC=demo/seq/A069999-demo.cc DEMOFLAGS=-DTIMING"
177 /// End:
178
179