Skip to content

Taylor Indices

Jump to bottom Edit New page
Parsa Amini edited this page Apr 27, 2017 · 2 revisions

The indicies defined in the arrays taylor_consts::map2, taylor_consts::map3, and taylor_consts::map4 are set up in a way, such that using them in loops which are structured as

taylor<5, T> x;
x() = ...
for (int a = 0; a != NDIM; ++a) {
    x(a) = ...;
    for (int b = a; b != NDIM; ++b) {
        x(a, b) = ...;
        for (int c = b; c != NDIM; ++c) {
            x(a, b, c) = ...;
            for (int d = c; d != NDIM; ++d) {
                x(a, b, c, d) = ...;
        }
    }
}

the expressions x(a), x(a, b), x(a, b, c), and x(a, b, c, d) evaluate to monotonically increasing indicies into the one dimensional array representing the storage inside x.data:

Zero-dimensional index:

x() x.data[0]

One-dimensional indicies:

a x(a)
0 x.data[1]
1 x.data[2]
2 x.data[3]

Two-dimensional indicies:

a b x(a, b) a b x(a, b) a b x(a, b)
0 0 x.data[4] * 1 0 * 2 0
0 1 x.data[5] 1 1 x.data[7] * 2 1
0 2 x.data[6] 1 2 x.data[8] 2 2 x.data[9]

Three-dimensional indicies:

a b c x(a, b, c) a b c x(a, b, c) a b c x(a, b, c)
0 0 0 x.data[10] * 1 0 0 * 2 0 0
0 0 1 x.data[11] * 1 0 1 * 2 0 1
0 0 2 x.data[12] * 1 0 2 * 2 0 2
* 0 1 0 * 1 1 0 * 2 1 0
0 1 1 x.data[13] 1 1 1 x.data[16] * 2 1 1
0 1 2 x.data[14] 1 1 2 x.data[17] * 2 1 2
* 0 2 0 * 1 2 0 * 2 2 0
* 0 2 1 * 1 2 1 * 2 2 1
0 2 2 x.data[15] 1 2 2 x.data[18] 2 2 2 x.data[19]

Four-dimensional indicies:

a b c d x(a, b, c, d) a b c d x(a, b, c, d) a b c d x(a, b, c, d)
0 0 0 0 x.data[20] * 1 0 0 0 * 2 0 0 0
0 0 0 1 x.data[21] * 1 0 0 1 * 2 0 0 1
0 0 0 2 x.data[22] * 1 0 0 2 * 2 0 0 2
* 0 0 1 0 * 1 0 1 0 * 2 0 1 0
0 0 1 1 x.data[23] * 1 0 1 1 * 2 0 1 1
0 0 1 2 x.data[24] * 1 0 1 2 * 2 0 1 2
* 0 0 2 0 * 1 0 2 0 * 2 0 2 0
* 0 0 2 1 * 1 0 2 1 * 2 0 2 1
0 0 2 2 x.data[25] * 1 0 2 2 * 2 0 2 2
* 0 1 0 0 * 1 1 0 0 * 2 1 0 0
* 0 1 0 1 * 1 1 0 1 * 2 1 0 1
* 0 1 0 2 * 1 1 0 2 * 2 1 0 2
* 0 1 1 0 * 1 1 1 0 * 2 1 1 0
0 1 1 1 x.data[26] 1 1 1 1 x.data[30] * 2 1 1 1
0 1 1 2 x.data[27] 1 1 1 2 x.data[31] * 2 1 1 2
* 0 1 2 0 * 1 1 2 0 * 2 1 2 0
* 0 1 2 1 * 1 1 2 1 * 2 1 2 1
0 1 2 2 x.data[28] 1 1 2 2 x.data[32] * 2 1 2 2
* 0 2 0 0 * 1 2 0 0 * 2 2 0 0
* 0 2 0 1 * 1 2 0 1 * 2 2 0 1
* 0 2 0 2 * 1 2 0 2 * 2 2 0 2
* 0 2 1 0 * 1 2 1 0 * 2 2 1 0
* 0 2 1 1 * 1 2 1 1 * 2 2 1 1
* 0 2 1 2 * 1 2 1 2 * 2 2 1 2
* 0 2 2 0 * 1 2 2 0 * 2 2 2 0
* 0 2 2 1 * 1 2 2 1 * 2 2 2 1
0 2 2 2 x.data[29] 1 2 2 2 x.data[33] 2 2 2 2 x.data[34]

The flipside of all of this is that now we can convert any loops which are structured as shown above:

taylor<5, T> x;
x() = ...
for (int a = 0; a != NDIM; ++a) {
    x(a) = ...;
    for (int b = a; b != NDIM; ++b) {
        x(a, b) = ...;
        for (int c = b; c != NDIM; ++c) {
            x(a, b, c) = ...;
            for (int d = c; d != NDIM; ++d) {
                x(a, b, c, d) = ...;
        }
    }
}

into a equivalent set of loops structured as:

for (integer i = taylor_sizes[0]; i != taylor_sizes[1]; ++i) {
    x[i] = ...;
}
for (integer i = taylor_sizes[1]; i != taylor_sizes[2]; ++i) {
    x[i] = ...;
}
for (integer i = taylor_sizes[2]; i != taylor_sizes[3]; ++i) {
    x[i] = ...;
}
for (integer i = taylor_sizes[3]; i != taylor_sizes[4]; ++i) {
    x[i] = ...;
}

Note that the taylor<>::operator[]() directly accesses the underlying one-dimensional data store encapsulated by this type.

Additionally we now can convert any more complex index-based referral into a simpler one, for instance:

taylor<4, T> XX;
std::array<T, NDIM> X;
for (integer a = 0; a != NDIM; a++) {
    XX(a) = X[a];
    for (integer b = a; b != NDIM; b++) {
        XX(a, b) = X[a] * X[b];
        for (integer c = b; c != NDIM; c++) {
            XX(a, b, c) = X[a] * X[b] * X[c];
        }
    }
}

into:

taylor<4, T> XX;
std::array<T, NDIM> X;
for (integer i = taylor_sizes[0]; i != taylor_sizes[1]; ++i) {
    XX[to_a(i)] = X[to_a(i)-1];
}
for (integer i = taylor_sizes[1]; i != taylor_sizes[2]; ++i) {
    XX[to_ab(i)] = X[to_a(i) - 1] * X[to_b(i) - 1];
}
for (integer i = taylor_sizes[2]; i != taylor_sizes[3]; ++i) {
    XX[to_abc(i)] = X[to_a(i) - 1] * X[to_b(i) - 1] * X[to_c(i) - 1];
}

where the index manipulation functions to_a, to_a_base, etc. are accessing pre-calculated index-mapping tables to retrieve the proper index into the underlying one-dimensional data store.

The array taylor_sizes is defined as follows:

constexpr integer taylor_sizes[] = { 1, 4, 10, 20, 35 };
One dimensional iteration

This:

for (integer a = 0; a != NDIM; a++) {
    XX(a) = X[a];
}

should be equivalent to:

for (integer i = taylor_sizes[0]; i != taylor_sizes[1]; ++i) {
    XX[to_a(i)] = X[to_a(i) - 1];
}
i to_a
1 1
2 2
3 3
Two dimensional iteration

This:

for (integer a = 0; a != NDIM; a++) {
    XX(a) = X[a];
    for (integer b = a; b != NDIM; b++) {
        XX(a, b) = X[a] * X[b];
    }
}

should be equivalent to:

for (integer i = taylor_sizes[0]; i != taylor_sizes[1]; ++i) {
    XX[to_a(i)] = X[to_a(i) - 1];
}
for (integer i = taylor_sizes[1]; i != taylor_sizes[2]; ++i) {
    XX[to_ab(i)] = X[to_a(i) - 1] * X[to_b(i) - 1];
}
i to_a to_b to_ab
1 1 * *
2 2 * *
3 3 * *
4 1 1 4
5 1 2 5
6 1 3 6
7 2 2 7
8 2 3 8
9 3 3 9
Three dimensional iteration

This:

for (integer a = 0; a != NDIM; a++) {
    XX(a) = X[a];
    for (integer b = a; b != NDIM; b++) {
        XX(a, b) = X[a] * X[b];
        for (integer c = b; c != NDIM; c++) {
            XX(a, b, c) = X[a] * X[b] * X[c];
        }
    }
}

should be equivalent to:

for (integer i = taylor_sizes[0]; i != taylor_sizes[1]; ++i) {
    XX[to_a(i)] = X[to_a(i) - 1];
}
for (integer i = taylor_sizes[1]; i != taylor_sizes[2]; ++i) {
    XX[to_ab(i)] = X[to_a(i) - 1] * X[to_b(i) - 1];
}
for (integer i = taylor_sizes[2]; i != taylor_sizes[3]; ++i) {
    XX[to_abc(i)] = X[to_a(i) - 1] * X[to_b(i) - 1] * X[to_c(i) - 1];
}
i to_a to_b to_ab to_c to_abc
1 1 * * * *
2 2 * * * *
3 3 * * * *
4 1 1 4 * *
5 1 2 5 * *
6 1 3 6 * *
7 2 2 7 * *
8 2 3 8 * *
9 3 3 9 * *
10 1 1 4 1 10
11 1 1 4 2 11
12 1 1 4 3 12
13 1 2 5 2 13
14 1 2 5 3 14
15 1 3 6 3 15
16 2 2 7 2 16
17 2 2 7 3 17
18 2 3 8 3 18
19 3 3 9 3 19

General note: with the index mapping as described by the tables above the following is true:

A(a, b, c) === A(to_a(i), to_b(i), to_c(i)) === A[to_abc(i)]
Add a custom footer
Add a custom sidebar
Clone this wiki locally