OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
FTpsData.h
Go to the documentation of this file.
1#ifndef CLASSIC_FTpsData_H
2#define CLASSIC_FTpsData_H
3
4// ------------------------------------------------------------------------
5// $RCSfile: FTpsData.h,v $
6// ------------------------------------------------------------------------
7// $Revision: 1.2.2.10 $
8// ------------------------------------------------------------------------
9// Copyright: see Copyright.readme
10// ------------------------------------------------------------------------
11//
12// Template class: FTpsData
13//
14// ------------------------------------------------------------------------
15// Class category: FixedAlgebra
16// ------------------------------------------------------------------------
17//
18// $Date: 2004/11/12 18:57:54 $
19// $Author: adelmann $
20//
21// ------------------------------------------------------------------------
22
23#include "Algebra/Array1D.h"
27
28// Template class FTpsData<N>
29// ------------------------------------------------------------------------
31
32template <int N>
33class FTpsData {
34
35public:
36
37 FTpsData();
38 ~FTpsData();
39
40 // Return exponents for the monomial with given Giorgilli index.
41 static inline const FMonomial<N> &getExponents(int index);
42
43 // Return the Giorgilli index for a given monomial.
44 static int getIndex(const FMonomial<N> &);
45
46 // Order of monomial which has Giorgilli index "index".
47 static inline int getOrder(int index);
48
49 // Number of terms in an FTps<T,N> object of order "order".
50 static inline int getSize(int order);
51
52 // Index at which order "order" begins; this equals the number
53 // of monomials of degree <= (order - 1) in an FTps<T,N> object.
54 static inline int orderStart(int order);
55
56 // Index at which order "order" begins for polynomials in nv (<= N) variables.
57 static inline int orderStart(int order, int nv);
58
59 // One plus the index at which order "order" ends.
60 static inline int orderEnd(int order);
61
62 // One plus index at which order "order" ends for polynomials in nv (<= N) variables.
63 static inline int orderEnd(int order, int nv);
64
65 // Number of monomials of degree "order".
66 static inline int orderLength(int order);
67
68 // Number of monomials such that "orderL" <= degree <= "orderH".
69 static inline int orderLength(int orderL, int orderH);
70
71 // Return the product index array for monomial "index".
72 static inline const Array1D<int> &getProductArray(int index);
73
74 // Return the variable list for monomial "index".
75 static inline const Array1D<int> &getVariableList(int index);
76
77 // Return the substitution table.
78 static inline const Array1D<TpsSubstitution> &getSubTable();
79
80 // To ensure proper table setup, this method must be
81 // called by FTps whenever a new object is created.
82 static void setup(int order);
83
84private:
85
86 // Not implemented.
88 void operator=(const FTpsData &);
89
90 // Build tables for given order.
91 void build(int order);
92
93 // Initialise Substitution table.
94 void fillSubst(int var, int order, FMonomial<N> &pow, int &next);
95
96 // Binomial coefficients: binom[i-1][k] = (i+k)!/(i! k!); used in the indexing
97 // methods getIndex(), getSize(), orderStart(), orderEnd, and orderLength().
99
100 // Boundary indices for each order.
102
103 // Exponents: expon[index] = set of exponents in monomial "index".
105
106 // Indexing information, used for Tps multiplication.
107 // prod[i][j] = index of monomial that results from
108 // multiplying monomials i and j
110
111 // Indexing information, used for Tps multiplication
112 // using the look-back method. These two arrays are
113 // used to identify which monomial pairs contribute
114 // to a given monomial in the result.
117
118 // Indexing information, used for Tps multiplication
119 // using the compact index method.
121
122 // Array of variable lists; used for Tps substitution.
124
125 // The substitution table.
127
128 // The current maximum order for which the tables exist.
129 static int topOrder;
130
131 // The current maximum size of any FTps object.
132 static int topSize;
133
134 // The singleton object for local book-keeping.
136};
137
138
139// Class FTpsData, public methods.
140// ------------------------------------------------------------------------
141
142// The singleton object for local book-keeping.
143template <int N>
145theBook = new FTpsData<N>();
146
147// The current maximum order for which the tables exist.
148template <int N>
150
151// The current maximum size of any FTps object.
152template <int N>
154
155
156template <int N> inline
158getExponents(int index) {
159 return theBook->expon[index];
160}
161
162
163template <int N>
165getIndex(const FMonomial<N> &pows) {
166 int order = 0;
167 int index = 0;
168
169 for(int var = N; var-- > 0;) {
170 order += pows[var];
171 index += theBook->binom[order][var];
172 }
173
174 return index;
175}
176
177
178template <int N> inline
180getOrder(int index) {
181 return theBook->expon[index].getOrder();
182}
183
184
185template <int N> inline
187getSize(int order) {
188 return theBook->bin[order+1];
189}
190
191
192template <int N> inline
194orderStart(int order) {
195 return theBook->bin[order];
196}
197
198
199template <int N> inline
201orderStart(int order, int nv) {
202 return theBook->binom[order][N-nv];
203}
204
205
206template <int N> inline
208orderEnd(int order) {
209 return theBook->bin[order+1];
210}
211
212
213template <int N> inline
215orderEnd(int order, int nv) {
216 return theBook->binom[order+1][N-nv];
217}
218
219
220template <int N> inline
222orderLength(int order) {
223 // return theBook->bin[order+1] - theBook->bin[order];
224 return theBook->binom[order+1][1];
225}
226
227
228template <int N> inline
230orderLength(int orderL, int orderH) {
231 return theBook->bin[orderH+1] - theBook->bin[orderL];
232}
233
234
235template <int N> inline
237getProductArray(int index) {
238 return theBook->prod[index];
239}
240
241
242template <int N> inline
244getVariableList(int index) {
245 return theBook->vrblList[index];
246}
247
248
249//template <int N> inline
250//const Array1D<int> &FTpsData<N>::
251//getVariableList(int index)
252//{
253// FMonomial<N> jlist = theBook->expon[index];
254// Array1D<int> *result = new Array1D<int> (jlist.getOrder(),0);
255//
256// int k = 0;
257// for (int v = 0; v < N; ++v) {
258// int jv = jlist[v];
259// while (jv-- > 0) (*result)[k++] = v;
260// }
262// return *result;
263//}
264
265
266template <int N> inline
268getSubTable() {
269 return theBook->subTable;
270}
271
272
273template <int N>
275setup(int order) {
276 theBook->build(order);
277}
278
279
280// Class FTpsData, protected methods.
281// ------------------------------------------------------------------------
282
283template <int N>
285FTpsData() {
286 // Initialize with order 4. May expand later.
287 topOrder = 0;
288 build(4);
289}
290
291
292template <int N>
294~FTpsData()
295{}
296
297
298// Class FTpsData, private methods.
299// ------------------------------------------------------------------------
300
301template <int N>
303build(int order) {
304 if(topOrder < order) {
305 topOrder = order;
306 int prodSize = 0;
307
308 // Build array containing binomial coefficients:
309 // binom[0 ... topOrder+1][0 ... N].
310 binom = Array1D < FArray1D < int, N + 1 > > (topOrder + 2);
311 // Here var runs from N to 0.
312 for(int var = N + 1; var-- > 0;) binom[0][var] = 0;
313 for(int ord = 1; ord <= topOrder + 1; ord++) {
314 binom[ord][N] = 1;
315 // Here var runs from N-1 to 0.
316 for(int var = N; var-- > 0;)
317 binom[ord][var] = binom[ord-1][var] + binom[ord][var+1];
318 }
319
320 // Build array of indices that mark order boundaries.
321 bin = Array1D<int>(topOrder + 2);
322 for(int i = 0; i < topOrder + 2; i++) bin[i] = binom[i][0];
323
324 // Build table of monomial exponents for orders zero through topOrder.
325 // (algorithm due to Liam Healey (Marylie 3.0)).
326 {
327 topSize = bin[topOrder+1];
328 expon = Array1D < FMonomial<N> >(topSize);
329 FMonomial<N> power;
330
331 if(N == 1)
332 for(int index = 1; index < topSize; index++) expon[index][0] = index;
333 else
334 for(int index = 1; index < topSize; index++) {
335 int carry = power[N-1];
336 power[N-1] = 0;
337 int lastnz = N - 2;
338 while(power[lastnz] == 0 && lastnz-- > 0) {}
339
340 if(lastnz == -1) power[0] = 1 + carry;
341 else {
342 power[lastnz]--;
343 power[lastnz+1] += 1 + carry;
344 }
345
346 expon[index] = power;
347 }
348 }
349
350 // Build table of products for orders zero through topOrder.
351 {
352 FMonomial<N> power;
353 prod = Array1D < Array1D<int> >(topSize);
354 for(int xord = 0; xord <= topOrder; xord++) {
355 int yord = topOrder - xord;
356 int ysize = bin[yord+1];
357 for(int i = bin[xord]; i < bin[xord+1]; i++) {
358 prod[i] = Array1D<int>(ysize);
359 // use symmetry for LL half of prod array
360 for(int j = 0; j < std::min(i, ysize); j++) {
361 prod[i][j] = prod[j][i];
362 ++prodSize;
363 }
364 for(int j = i; j < ysize; j++) {
365 power = expon[i] * expon[j];
366 int ord = 0;
367 int ind = 0;
368 for(int vv = N; vv-- > 0;) {
369 ord += power[vv];
370 ind += binom[ord][vv];
371 }
372 prod[i][j] = ind;
373 ++prodSize;
374 }
375 }
376 }
377 }
378
379 // Build table of variable lists.
380 {
381 // Allocate variable list pointers and working array.
382 vrblList = Array1D< Array1D<int> >(topSize);
383 int *vars = new int[topOrder];
384 // Initialise counter.
385 int j = 1, N1 = N - 1;
386 // Loop over the orders.
387 for(int ord = 1; ord <= topOrder; ord++) {
388 // Allocate first variable list; initialise both it and vars.
389 vrblList[j] = Array1D<int>(ord, 0);
390 std::fill(vars, vars + ord, int(0));
391 // Define end points.
392 int last_j = bin[ord+1];
393 int *vlast = vars + ord;
394 // Build remaining variable lists at this order.
395 while(++j < last_j) {
396 // From last variable list (in vars), construct next one.
397 int *vi = vlast;
398 while(*--vi == N1) {};
399 int k = *vi + 1;
400 std::fill(vi, vlast, k);
401 // Allocate j-th variable list and copy from vars.
402 vrblList[j] = Array1D<int>(ord);
403 std::copy(vars, vars + ord, vrblList[j].begin());
404 }
405 }
406 delete [] vars;
407 }
408
409 // Build the substitution table.
410 {
411 subTable = Array1D<TpsSubstitution>(topSize);
412 FMonomial<N> power;
413 int next = 1;
414 fillSubst(0, 1, power, next);
415 }
416 }
417}
418
419
420template <int N>
422fillSubst(int var, int order, FMonomial<N> &power, int &next) {
423 for(int v = var; v < N; v++) {
424 TpsSubstitution &s = subTable[next];
425 power[v]++;
426 int ord = 0;
427 int ind = 0;
428 for(int vv = N; vv-- > 0;) {
429 ord += power[vv];
430 ind += binom[ord][vv];
431 }
432 s.index = ind;
433 s.order = order;
434 s.variable = v;
435 next++;
436 if(order < topOrder) fillSubst(v, order + 1, power, next);
437 s.skip = next;
438 power[v]--;
439 }
440}
441
442#endif // CLASSIC_FTpsData_H
PartBunchBase< T, Dim >::ConstIterator begin(PartBunchBase< T, Dim > const &bunch)
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Definition: TpsMath.h:76
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:76
T::PETE_Expr_t::PETE_Return_t prod(const PETE_Expr< T > &expr)
Definition: PETE.h:1121
Substitution for Tps<T>.
Representation of the exponents for a monomial with fixed dimension.
Definition: FMonomial.h:32
Internal utility class for FTps<T,N> class.
Definition: FTpsData.h:33
static int orderEnd(int order)
Definition: FTpsData.h:208
static const FMonomial< N > & getExponents(int index)
Definition: FTpsData.h:158
~FTpsData()
Definition: FTpsData.h:294
void fillSubst(int var, int order, FMonomial< N > &pow, int &next)
Definition: FTpsData.h:422
FTpsData()
Definition: FTpsData.h:285
static int getIndex(const FMonomial< N > &)
Definition: FTpsData.h:165
static const Array1D< TpsSubstitution > & getSubTable()
Definition: FTpsData.h:268
Array1D< Array1D< int > > prod
Definition: FTpsData.h:109
static FTpsData< N > * theBook
Definition: FTpsData.h:135
Array1D< int > lookBack
Definition: FTpsData.h:115
FTpsData(const FTpsData &)
Array1D< Array1D< int > > vrblList
Definition: FTpsData.h:123
void build(int order)
Definition: FTpsData.h:303
static int getOrder(int index)
Definition: FTpsData.h:180
Array1D< int > lbBound
Definition: FTpsData.h:116
static int getSize(int order)
Definition: FTpsData.h:187
Array1D< TpsSubstitution > subTable
Definition: FTpsData.h:126
static int topSize
Definition: FTpsData.h:132
static int topOrder
Definition: FTpsData.h:129
static const Array1D< int > & getVariableList(int index)
Definition: FTpsData.h:244
static void setup(int order)
Definition: FTpsData.h:275
static int orderLength(int order)
Definition: FTpsData.h:222
Array1D< FArray1D< int, N+1 > > binom
Definition: FTpsData.h:98
void operator=(const FTpsData &)
Array1D< int > Giorgilli2ExponIndex
Definition: FTpsData.h:120
static int orderStart(int order)
Definition: FTpsData.h:194
Array1D< int > bin
Definition: FTpsData.h:101
static const Array1D< int > & getProductArray(int index)
Definition: FTpsData.h:237
Array1D< FMonomial< N > > expon
Definition: FTpsData.h:104