C++ Tensor Template Library (TTL)
If I just say, that here are templates, which implement full vector/matrix and partially tensor algebras over arbitrary class, I am afraid that I will hear...
Tensor?
- May be it has something to do with nuclear bomb?
- I have learned it ... years ago but since never used.
- Coordinates are tensors? Oh yes, DPtoLP, ScreenToClient... I hate it!
But sometimes mathematics makes complicated things easier. A lot of dull problems could be solved, just using four arithmetic operations +, -, *, /, applied to tensors.
At first - tensor is an array. Yes, just the array, which we use so often in programming. The array can be multi-dimensional. Amount of dimensions is equal to the number of tensor indexes. For example:
vi - one-dimensional array, or vector
aij - two-dimensional array, or matrix
tstuvwxyz - 8-dimensional array
The only difference to the array, as you have noticed, indexes can be upper or lower. Following Paul Dirac, we will call array with lower indexes Bra and with upper - Ket. Then we can represent in our program
vi - by class Ket<U> or Vector<U>, where U is class, placed into array
aij - by Ket<Bra<U>> or Matrix<U>
tstuvwxyz - by Ket<Ket<Bra<Ket<Ket<Bra<Bra<Bra<U>>>>>>>>
Forget now about tensors with more than two indexes (of course, if you are not solving equations of general relativity theory). We will find enough applications for vectors and matrices.
Do you know some programmer, who have got from the first attempt correct image on the screen after several successive rotations and mirrorings? It's like a nightmare - you always get one of 7 wrong images. I have started this libraryafter 11th attempt to fix a bug in such a program.
Any mirroring or rotation can be represented as a matrix. If you have such matrix and vector v, then vector x after transformation could be found from:
xi = mij * vj
If you apply subsequently another transformation with matrix n:
xi = nik * mkj * vj
Note that we always multiply new matrix from the left (the order of multiplication is important - for tensors a * b != b * a). You can also add any offset and scale your vector.
Now let us see how looks the program, which draws some image, stored in vector format after a set of such transformations:
int i;
int scale = 3;
Vector<int> image[100];
Matrix<int> mirrorX(0); mirrorX[1][1]=1; mirrorX[0][0]=-1;
Matrix<int> mirrorY(0); mirrorY[0][0]=1; mirrorY[1][1]=-1;
Matrix<int> rotateL(0); rotateL[1][0]=1; rotateL[0][1]=-1;
Matrix<int> rotateR(0); rotateR[0][1]=1; rotateR[1][0]=-1;
Vector<int> offset; offset[0]=20; offset[1]=-180;
Matrix<int> yAxeInWindowsDown=mirrorY;
Matrix<int> m = scale * yAxeInWindowsDown * rotateL * mirrorX * rotateR;
Vector<int> x = offset + m * image[0];
MoveToEx( hDC, x[0], x[1], NULL );
for ( i = 0; i < 100; i++ )
{
x = offset + m * image[i];
LineTo( hDC, x[0], x[1] );
}
Do I need any comments in the code above?
Following code allows to modify with the mouse your initial image. Look how we transform screen coordinates back to image coordinates:
case WM_LBUTTONUP:
{
Vector<int> mouse; mouse[0]=LOWORD(lParam); mouse[1]=HIWORD(lParam);
Matrix<int> inverted = Matrix<int>(1) / ( m / scale );
image[0] = inverted * ( mouse - offset ) / scale;
}
Stop. How it happened that inverted matrix was calculated for matrix of integers? In general, one should use Matrix<float> instead.
But look at our trick with ( m / scale ). This matrix contains only rotations and mirrorings, which means, that inverted matrix consists also only from rotations and mirrorings!
Any linear system of equations, like:
0.15*x0 + 3.7*x1 = 2.6 0.32*x0 + 2.4*x1 = 1.5
could be represented as vector/matrix equation:
/ 0.15 3.7 \ / x0 \ / 2.6 \ | | * | | = | | \ 0.32 2.4 / \ x1 / \ 1.5 /
or in tensor form:
mij * xj = ai
where mij - is a matrix and ai - resulting vector.
How do we solve this system of equations? Oh, if we can calculate an inverted matrix - it's very easy. We multiply inverted matrix by ai and find xj. We will denote calculation of inverted matrix as
m-1 jk = uik / mij
where uik is a unary matrix. But this divide operation can be applied not only to matrices. So we will get a solution even more simple and intuitive:
xj = ai / mij
And it will be calculated at least twice faster than using inverted matrix!
I will not continue the list of examples, because tensors represent rather deep level of abstraction and my aim was to transfer this deepness into the code. And it is known, that the deeper the abstraction is - the wider range where it can be applied. Will you ever list applications of numbers?
You can use any standard numeric types in tensor templates, as well as ANSI bcd and complex classes.
You can use also Fraction templates, which are supplied in header "fraction.h" together with "tensor.h". In this case (if you avoid overflow, selecting correct Fraction template argument) you get exact solutions (without rounding errors) for inverted matrices and systems of equations.
Tensors themselves could be also used as template arguments for tensors! If one really needs it - Vector<Matrix<double>> or Matrix<Matrix<Fraction<long>>> - will be also valid.
Designing own classes to be used as tensor arguments is also not prohibited. Minimum requirements (to create additive group):
Additional requirements (to create a ring or field):
The number of tensor indexes possible, tested with BC++ v 4.52 is not less than ten. May be more - I have not tried. Hope it's enough.
Strict. You will get an error message already at compile time, if you want to multiply Matrix<float,2,2> by Matrix<float,3,4> or subtract Vector<int,2> from Vector<int,3>. No need to wait for run-time "array bounds error" (even if it can be detected).
I have optimised programming in favour of speed, preserving deepest possible level of abstraction.