Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F19251233
Tensor.cc
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Award Token
Flag For Later
Size
24 KB
Referenced Files
None
Subscribers
None
Tensor.cc
View Options
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include
"HEJ/Tensor.hh"
#include
<array>
#include
<iostream>
#include
<valarray>
#include
<vector>
#include
<CLHEP/Vector/LorentzVector.h>
namespace
HEJ
{
namespace
{
// Tensor Template definitions
short
int
sigma_index5
[
1024
];
short
int
sigma_index3
[
64
];
std
::
valarray
<
COM
>
permfactor5
;
std
::
valarray
<
COM
>
permfactor3
;
short
int
helfactor5
[
2
][
1024
];
short
int
helfactor3
[
2
][
64
];
// map from a list of tensor lorentz indices onto a single index
// in d dimensional spacetime
// TODO: basically the same as detail::accumulate_idx
template
<
std
::
size_t
N
>
std
::
size_t
tensor_to_list_idx
(
std
::
array
<
int
,
N
>
const
&
indices
)
{
std
::
size_t
list_idx
=
0
;
for
(
std
::
size_t
i
=
1
,
factor
=
1
;
i
<=
N
;
++
i
)
{
list_idx
+=
factor
*
indices
[
N
-
i
];
factor
*=
detail
::
d
;
}
#ifndef NDEBUG
constexpr
std
::
size_t
max_possible_idx
=
detail
::
power
(
detail
::
d
,
N
);
assert
(
list_idx
<
max_possible_idx
);
#endif
return
list_idx
;
}
// generate all unique perms of vectors {a,a,a,a,b}, return in perms
// set_permfactor is a bool which encodes the anticommutation relations of the
// sigma matrices namely if we have sigma0, set_permfactor= false because it
// commutes with all others otherwise we need to assign a minus sign to odd
// permutations, set in permfactor
// note, inital perm is always even
void
perms41
(
int
same4
,
int
diff
,
std
::
vector
<
std
::
array
<
int
,
5
>>
&
perms
){
bool
set_permfactor
(
true
);
if
(
same4
==
0
||
diff
==
0
)
set_permfactor
=
false
;
for
(
int
diffpos
=
0
;
diffpos
<
5
;
diffpos
++
){
std
::
array
<
int
,
5
>
tempvec
=
{
same4
,
same4
,
same4
,
same4
,
same4
};
tempvec
[
diffpos
]
=
diff
;
perms
.
push_back
(
tempvec
);
if
(
set_permfactor
){
if
(
diffpos
%
2
==
1
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
;
}
}
}
// generate all unique perms of vectors {a,a,a,b,b}, return in perms
// note, inital perm is always even
void
perms32
(
int
same3
,
int
diff
,
std
::
vector
<
std
::
array
<
int
,
5
>>
&
perms
){
bool
set_permfactor
(
true
);
if
(
same3
==
0
||
diff
==
0
)
set_permfactor
=
false
;
for
(
int
diffpos1
=
0
;
diffpos1
<
5
;
diffpos1
++
){
for
(
int
diffpos2
=
diffpos1
+
1
;
diffpos2
<
5
;
diffpos2
++
){
std
::
array
<
int
,
5
>
tempvec
=
{
same3
,
same3
,
same3
,
same3
,
same3
};
tempvec
[
diffpos1
]
=
diff
;
tempvec
[
diffpos2
]
=
diff
;
perms
.
push_back
(
tempvec
);
if
(
set_permfactor
){
if
((
diffpos2
-
diffpos1
)
%
2
==
0
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
;
}
}
}
}
// generate all unique perms of vectors {a,a,a,b,c}, return in perms
// we have two bools since there are three distinct type of sigma matrices to
// permute, so need to test if the 3xrepeated = sigma0 or if one of the
// singles is sigma0
// if diff1/diff2 are distinct, flipping them results in a change of perm,
// otherwise it'll be symmetric under flip -> encode this in set_permfactor2
// as long as diff2!=0 can ensure inital perm is even
// if diff2=0 then it is not possible to ensure inital perm even -> encode in
// bool signflip
void
perms311
(
int
same3
,
int
diff1
,
int
diff2
,
std
::
vector
<
std
::
array
<
int
,
5
>>
&
perms
){
bool
set_permfactor2
(
true
);
bool
same0
(
false
);
bool
diff0
(
false
);
bool
signflip
(
false
);
// if true, inital perm is odd
if
(
same3
==
0
)
// is the repeated matrix sigma0?
same0
=
true
;
else
if
(
diff2
==
0
){
// is one of the single matrices sigma0
diff0
=
true
;
if
((
diff1
%
3
-
same3
)
!=-
1
)
signflip
=
true
;
}
else
if
(
diff1
==
0
){
std
::
cerr
<<
"Note, only first and last argument may be zero"
<<
std
::
endl
;
return
;
}
// possible outcomes: tt, ft, tf
for
(
int
diffpos1
=
0
;
diffpos1
<
5
;
diffpos1
++
){
for
(
int
diffpos2
=
diffpos1
+
1
;
diffpos2
<
5
;
diffpos2
++
){
std
::
array
<
int
,
5
>
tempvec
=
{
same3
,
same3
,
same3
,
same3
,
same3
};
tempvec
[
diffpos1
]
=
diff1
;
tempvec
[
diffpos2
]
=
diff2
;
perms
.
push_back
(
tempvec
);
if
(
!
same0
&&
!
diff0
){
// full antisymmetric under exchange of same3,diff1,diff2
if
((
diffpos2
-
diffpos1
)
%
2
==
0
){
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd perm
// if this perm is odd, swapping diff1<->diff2 automatically even
set_permfactor2
=
false
;
}
else
{
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=
COM
(
0
,
1
);
// even perm
// if this perm is even, swapping diff1<->diff2 automatically odd
set_permfactor2
=
true
;
}
}
else
if
(
diff0
){
// one of the single matrices is sigma0
if
(
signflip
){
// initial config is odd!
if
(
diffpos1
%
2
==
1
){
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=
COM
(
0
,
1
);
// even perm
// initally symmetric under diff1<->diff2 =>
// if this perm is even, automatically even for first diffpos2
set_permfactor2
=
false
;
}
else
{
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd perm
// initally symmetric under diff1<->diff2 =>
// if this perm is odd, automatically odd for first diffpos2
set_permfactor2
=
true
;
}
}
else
{
// diff0=true, initial config is even
if
(
diffpos1
%
2
==
1
){
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd perm
// initally symmetric under diff1<->diff2 =>
// if this perm is odd, automatically odd for first diffpos2
set_permfactor2
=
true
;
}
else
{
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=
COM
(
0
,
1
);
// even perm
// initally symmetric under diff1<->diff2 =>
// if this perm is even, automatically even for first diffpos2
set_permfactor2
=
false
;
}
}
if
((
diffpos2
-
diffpos1
-
1
)
%
2
==
1
)
set_permfactor2
=!
set_permfactor2
;
// change to account for diffpos2
}
else
if
(
same0
){
// the repeated matrix is sigma0
// => only relative positions of diff1, diff2 matter.
// always even before flip because diff1pos<diff2pos
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=
COM
(
0
,
1
);
// if this perm is odd, swapping diff1<->diff2 automatically odd
set_permfactor2
=
true
;
}
tempvec
[
diffpos1
]
=
diff2
;
tempvec
[
diffpos2
]
=
diff1
;
perms
.
push_back
(
tempvec
);
if
(
set_permfactor2
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
else
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=
COM
(
0
,
1
);
}
}
}
// end perms311
// generate all unique perms of vectors {a,a,b,b,c}, return in perms
void
perms221
(
int
same2a
,
int
same2b
,
int
diff
,
std
::
vector
<
std
::
array
<
int
,
5
>>
&
perms
){
bool
set_permfactor1
(
true
);
bool
set_permfactor2
(
true
);
if
(
same2b
==
0
){
std
::
cerr
<<
"second entry in perms221() shouldn't be zero"
<<
std
::
endl
;
return
;
}
else
if
(
same2a
==
0
)
set_permfactor1
=
false
;
else
if
(
diff
==
0
)
set_permfactor2
=
false
;
for
(
int
samepos
=
0
;
samepos
<
5
;
samepos
++
){
int
permcount
=
0
;
for
(
int
samepos2
=
samepos
+
1
;
samepos2
<
5
;
samepos2
++
){
for
(
int
diffpos
=
0
;
diffpos
<
5
;
diffpos
++
){
if
(
diffpos
==
samepos
||
diffpos
==
samepos2
)
continue
;
std
::
array
<
int
,
5
>
tempvec
=
{
same2a
,
same2a
,
same2a
,
same2a
,
same2a
};
tempvec
[
samepos
]
=
same2b
;
tempvec
[
samepos2
]
=
same2b
;
tempvec
[
diffpos
]
=
diff
;
perms
.
push_back
(
tempvec
);
if
(
set_permfactor1
){
if
(
set_permfactor2
){
// full anti-symmetry
if
(
permcount
%
2
==
1
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
;
}
else
{
// diff is sigma0
if
(
((
samepos2
-
samepos
-
1
)
%
3
>
0
)
&&
(
abs
(
abs
(
samepos2
-
diffpos
)
-
abs
(
samepos
-
diffpos
))
%
3
>
0
)
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
;
}
}
else
{
// same2a is sigma0
if
((
diffpos
>
samepos
)
&&
(
diffpos
<
samepos2
))
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
;
}
permcount
++
;
}
}
}
}
// generate all unique perms of vectors {a,a,b,b,c}, return in perms
// there must be a sigma zero if we have 4 different matrices in string
// bool is true if sigma0 is the repeated matrix
void
perms2111
(
int
same2
,
int
diff1
,
int
diff2
,
int
diff3
,
std
::
vector
<
std
::
array
<
int
,
5
>>
&
perms
){
bool
twozero
(
false
);
if
(
same2
==
0
)
twozero
=
true
;
else
if
(
diff1
!=
0
){
std
::
cerr
<<
"One of first or second argurments must be a zero"
<<
std
::
endl
;
return
;
}
else
if
(
diff2
==
0
||
diff3
==
0
){
std
::
cerr
<<
"Only first and second arguments may be a zero."
<<
std
::
endl
;
return
;
}
int
permcount
=
0
;
for
(
int
diffpos1
=
0
;
diffpos1
<
5
;
diffpos1
++
){
for
(
int
diffpos2
=
0
;
diffpos2
<
5
;
diffpos2
++
){
if
(
diffpos2
==
diffpos1
)
continue
;
for
(
int
diffpos3
=
0
;
diffpos3
<
5
;
diffpos3
++
){
if
(
diffpos3
==
diffpos2
||
diffpos3
==
diffpos1
)
continue
;
std
::
array
<
int
,
5
>
tempvec
=
{
same2
,
same2
,
same2
,
same2
,
same2
};
tempvec
[
diffpos1
]
=
diff1
;
tempvec
[
diffpos2
]
=
diff2
;
tempvec
[
diffpos3
]
=
diff3
;
perms
.
push_back
(
tempvec
);
if
(
twozero
){
// don't care about exact positions of singles, just order
if
(
diffpos2
>
diffpos3
&&
diffpos3
>
diffpos1
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd
else
if
(
diffpos1
>
diffpos2
&&
diffpos2
>
diffpos3
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd
else
if
(
diffpos3
>
diffpos1
&&
diffpos1
>
diffpos2
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd
else
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=
COM
(
0
,
1
);
// evwn
}
else
{
if
(
permcount
%
2
==
1
)
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
else
permfactor5
[
tensor_to_list_idx
(
tempvec
)]
=
COM
(
0
,
1
);
}
permcount
++
;
}
}
}
}
void
perms21
(
int
same
,
int
diff
,
std
::
vector
<
std
::
array
<
int
,
3
>>
&
perms
){
bool
set_permfactor
(
true
);
if
(
same
==
0
||
diff
==
0
)
set_permfactor
=
false
;
for
(
int
diffpos
=
0
;
diffpos
<
3
;
diffpos
++
){
std
::
array
<
int
,
3
>
tempvec
=
{
same
,
same
,
same
};
tempvec
[
diffpos
]
=
diff
;
perms
.
push_back
(
tempvec
);
if
(
set_permfactor
&&
diffpos
==
1
)
permfactor3
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
;
}
}
void
perms111
(
int
diff1
,
int
diff2
,
int
diff3
,
std
::
vector
<
std
::
array
<
int
,
3
>>
&
perms
){
bool
sig_zero
(
false
);
if
(
diff1
==
0
)
sig_zero
=
true
;
else
if
(
diff2
==
0
||
diff3
==
0
){
std
::
cerr
<<
"Only first argument may be a zero."
<<
std
::
endl
;
return
;
}
int
permcount
=
0
;
for
(
int
pos1
=
0
;
pos1
<
3
;
pos1
++
){
for
(
int
pos2
=
pos1
+
1
;
pos2
<
3
;
pos2
++
){
std
::
array
<
int
,
3
>
tempvec
=
{
diff1
,
diff1
,
diff1
};
tempvec
[
pos1
]
=
diff2
;
tempvec
[
pos2
]
=
diff3
;
perms
.
push_back
(
tempvec
);
if
(
sig_zero
){
permfactor3
[
tensor_to_list_idx
(
tempvec
)]
=
1.
*
COM
(
0
,
1
);
// even
}
else
if
(
permcount
%
2
==
1
){
permfactor3
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd
}
else
{
permfactor3
[
tensor_to_list_idx
(
tempvec
)]
=
1.
*
COM
(
0
,
1
);
// even
}
tempvec
[
pos1
]
=
diff3
;
tempvec
[
pos2
]
=
diff2
;
perms
.
push_back
(
tempvec
);
if
(
sig_zero
){
permfactor3
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd
}
else
if
(
permcount
%
2
==
1
){
permfactor3
[
tensor_to_list_idx
(
tempvec
)]
=
1.
*
COM
(
0
,
1
);
// even
}
else
{
permfactor3
[
tensor_to_list_idx
(
tempvec
)]
=-
1.
*
COM
(
0
,
1
);
// odd
}
permcount
++
;
}
}
}
COM
perp
(
CLHEP
::
HepLorentzVector
const
&
p
)
{
return
COM
{
p
.
x
(),
p
.
y
()};
}
COM
perp_hat
(
CLHEP
::
HepLorentzVector
const
&
p
)
{
const
COM
p_perp
=
perp
(
p
);
if
(
p_perp
==
COM
{
0.
,
0.
})
return
-
1.
;
return
p_perp
/
std
::
abs
(
p_perp
);
}
}
// anonymous namespace
// "angle" product angle(pi, pj) = <i j>
// see eq. (53) (\eqref{eq:angle_product}) in developer manual
COM
angle
(
CLHEP
::
HepLorentzVector
const
&
pi
,
CLHEP
::
HepLorentzVector
const
&
pj
)
{
return
+
std
::
sqrt
(
COM
{
pi
.
minus
()
*
pj
.
plus
()})
*
perp_hat
(
pi
)
-
std
::
sqrt
(
COM
{
pi
.
plus
()
*
pj
.
minus
()})
*
perp_hat
(
pj
);
}
// "square" product square(pi, pj) = [i j]
// see eq. (54) (\eqref{eq:square_product}) in developer manual
COM
square
(
CLHEP
::
HepLorentzVector
const
&
pi
,
CLHEP
::
HepLorentzVector
const
&
pj
)
{
return
-
std
::
conj
(
angle
(
pi
,
pj
));
}
Tensor
<
2
>
metric
(){
static
const
Tensor
<
2
>
g
=
[](){
Tensor
<
2
>
g
(
0.
);
g
(
0
,
0
)
=
1.
;
g
(
1
,
1
)
=
-
1.
;
g
(
2
,
2
)
=
-
1.
;
g
(
3
,
3
)
=
-
1.
;
return
g
;
}();
return
g
;
}
// <1|mu|2>
// see eqs. (48), (49) in developer manual
Tensor
<
1
>
current
(
CLHEP
::
HepLorentzVector
const
&
p
,
CLHEP
::
HepLorentzVector
const
&
q
,
Helicity
h
){
using
namespace
helicity
;
const
COM
p_perp_hat
=
perp_hat
(
p
);
const
COM
q_perp_hat
=
perp_hat
(
q
);
std
::
array
<
std
::
array
<
COM
,
2
>
,
2
>
sqrt_pq
;
sqrt_pq
[
minus
][
minus
]
=
std
::
sqrt
(
COM
{
p
.
minus
()
*
q
.
minus
()})
*
p_perp_hat
*
conj
(
q_perp_hat
);
sqrt_pq
[
minus
][
plus
]
=
std
::
sqrt
(
COM
{
p
.
minus
()
*
q
.
plus
()})
*
p_perp_hat
;
sqrt_pq
[
plus
][
minus
]
=
std
::
sqrt
(
COM
{
p
.
plus
()
*
q
.
minus
()})
*
conj
(
q_perp_hat
);
sqrt_pq
[
plus
][
plus
]
=
std
::
sqrt
(
COM
{
p
.
plus
()
*
q
.
plus
()});
Tensor
<
1
>
j
;
j
(
0
)
=
sqrt_pq
[
plus
][
plus
]
+
sqrt_pq
[
minus
][
minus
];
j
(
1
)
=
sqrt_pq
[
plus
][
minus
]
+
sqrt_pq
[
minus
][
plus
];
j
(
2
)
=
COM
{
0.
,
1.
}
*
(
sqrt_pq
[
plus
][
minus
]
-
sqrt_pq
[
minus
][
plus
]);
j
(
3
)
=
sqrt_pq
[
plus
][
plus
]
-
sqrt_pq
[
minus
][
minus
];
return
(
h
==
minus
)
?
j
:
j
.
complex_conj
();
}
// <1|mu nu sigma|2>
// TODO: how does this work?
Tensor
<
3
>
rank3_current
(
CLHEP
::
HepLorentzVector
const
&
p1
,
CLHEP
::
HepLorentzVector
const
&
p2
,
Helicity
h
){
const
CLHEP
::
HepLorentzVector
p1new
{
p1
.
e
()
<
0
?-
p1
:
p1
};
const
CLHEP
::
HepLorentzVector
p2new
{
p2
.
e
()
<
0
?-
p2
:
p2
};
auto
j
=
HEJ
::
current
(
p1new
,
p2new
,
h
);
if
(
h
==
helicity
::
minus
)
{
j
*=
-
1
;
j
(
0
)
*=
-
1
;
}
Tensor
<
3
>
j3
;
for
(
std
::
size_t
tensor_idx
=
0
;
tensor_idx
<
j3
.
size
();
++
tensor_idx
){
const
double
hfact
=
helfactor3
[
h
][
tensor_idx
];
j3
.
components
[
tensor_idx
]
=
j
(
sigma_index3
[
tensor_idx
])
*
hfact
*
permfactor3
[
tensor_idx
];
}
return
j3
;
}
// <1|mu nu sigma tau rho|2>
Tensor
<
5
>
rank5_current
(
CLHEP
::
HepLorentzVector
const
&
p1
,
CLHEP
::
HepLorentzVector
const
&
p2
,
Helicity
h
){
const
CLHEP
::
HepLorentzVector
p1new
{
p1
.
e
()
<
0
?-
p1
:
p1
};
const
CLHEP
::
HepLorentzVector
p2new
{
p2
.
e
()
<
0
?-
p2
:
p2
};
auto
j
=
HEJ
::
current
(
p1new
,
p2new
,
h
);
if
(
h
==
helicity
::
minus
)
{
j
*=
-
1
;
j
(
0
)
*=
-
1
;
}
Tensor
<
5
>
j5
;
for
(
std
::
size_t
tensor_idx
=
0
;
tensor_idx
<
j5
.
size
();
++
tensor_idx
){
const
double
hfact
=
helfactor5
[
h
][
tensor_idx
];
j5
.
components
[
tensor_idx
]
=
j
(
sigma_index5
[
tensor_idx
])
*
hfact
*
permfactor5
[
tensor_idx
];
}
return
j5
;
}
Tensor
<
1
>
to_tensor
(
CLHEP
::
HepLorentzVector
const
&
p
){
Tensor
<
1
>
newT
;
newT
.
components
=
{
p
.
e
(),
p
.
x
(),
p
.
y
(),
p
.
z
()};
return
newT
;
}
// polarisation vector according to eq. (55) in developer manual
Tensor
<
1
>
eps
(
CLHEP
::
HepLorentzVector
const
&
pg
,
CLHEP
::
HepLorentzVector
const
&
pr
,
Helicity
pol
){
if
(
pol
==
helicity
::
plus
)
{
return
current
(
pr
,
pg
,
pol
)
/
(
std
::
sqrt
(
2
)
*
square
(
pg
,
pr
));
}
return
current
(
pr
,
pg
,
pol
)
/
(
std
::
sqrt
(
2
)
*
angle
(
pg
,
pr
));
}
namespace
{
// slow function! - but only need to evaluate once.
bool
init_sigma_index_helper
(){
// initialize permfactor(3) to list of ones (change to minus one for each odd
// permutation and multiply by i for all permutations in perms2111, perms311,
// perms111)
permfactor5
.
resize
(
1024
,
1.
);
permfactor3
.
resize
(
64
,
1.
);
// first set sigma_index (5)
std
::
vector
<
std
::
array
<
int
,
5
>>
sigma0indices
;
std
::
vector
<
std
::
array
<
int
,
5
>>
sigma1indices
;
std
::
vector
<
std
::
array
<
int
,
5
>>
sigma2indices
;
std
::
vector
<
std
::
array
<
int
,
5
>>
sigma3indices
;
// need to generate all possible permuations of {i,j,k,l,m}
// where each index can be {0,1,2,3,4}
// 1024 possibilities
// perms with 5 same
sigma0indices
.
push_back
({
0
,
0
,
0
,
0
,
0
});
sigma1indices
.
push_back
({
1
,
1
,
1
,
1
,
1
});
sigma2indices
.
push_back
({
2
,
2
,
2
,
2
,
2
});
sigma3indices
.
push_back
({
3
,
3
,
3
,
3
,
3
});
// perms with 4 same
perms41
(
1
,
0
,
sigma0indices
);
perms41
(
2
,
0
,
sigma0indices
);
perms41
(
3
,
0
,
sigma0indices
);
perms41
(
0
,
1
,
sigma1indices
);
perms41
(
2
,
1
,
sigma1indices
);
perms41
(
3
,
1
,
sigma1indices
);
perms41
(
0
,
2
,
sigma2indices
);
perms41
(
1
,
2
,
sigma2indices
);
perms41
(
3
,
2
,
sigma2indices
);
perms41
(
0
,
3
,
sigma3indices
);
perms41
(
1
,
3
,
sigma3indices
);
perms41
(
2
,
3
,
sigma3indices
);
// perms with 3 same, 2 same
perms32
(
0
,
1
,
sigma0indices
);
perms32
(
0
,
2
,
sigma0indices
);
perms32
(
0
,
3
,
sigma0indices
);
perms32
(
1
,
0
,
sigma1indices
);
perms32
(
1
,
2
,
sigma1indices
);
perms32
(
1
,
3
,
sigma1indices
);
perms32
(
2
,
0
,
sigma2indices
);
perms32
(
2
,
1
,
sigma2indices
);
perms32
(
2
,
3
,
sigma2indices
);
perms32
(
3
,
0
,
sigma3indices
);
perms32
(
3
,
1
,
sigma3indices
);
perms32
(
3
,
2
,
sigma3indices
);
// perms with 3 same, 2 different
perms311
(
1
,
2
,
3
,
sigma0indices
);
perms311
(
2
,
3
,
1
,
sigma0indices
);
perms311
(
3
,
1
,
2
,
sigma0indices
);
perms311
(
0
,
2
,
3
,
sigma1indices
);
perms311
(
2
,
3
,
0
,
sigma1indices
);
perms311
(
3
,
2
,
0
,
sigma1indices
);
perms311
(
0
,
3
,
1
,
sigma2indices
);
perms311
(
1
,
3
,
0
,
sigma2indices
);
perms311
(
3
,
1
,
0
,
sigma2indices
);
perms311
(
0
,
1
,
2
,
sigma3indices
);
perms311
(
1
,
2
,
0
,
sigma3indices
);
perms311
(
2
,
1
,
0
,
sigma3indices
);
perms221
(
1
,
2
,
0
,
sigma0indices
);
perms221
(
1
,
3
,
0
,
sigma0indices
);
perms221
(
2
,
3
,
0
,
sigma0indices
);
perms221
(
0
,
2
,
1
,
sigma1indices
);
perms221
(
0
,
3
,
1
,
sigma1indices
);
perms221
(
2
,
3
,
1
,
sigma1indices
);
perms221
(
0
,
1
,
2
,
sigma2indices
);
perms221
(
0
,
3
,
2
,
sigma2indices
);
perms221
(
1
,
3
,
2
,
sigma2indices
);
perms221
(
0
,
1
,
3
,
sigma3indices
);
perms221
(
0
,
2
,
3
,
sigma3indices
);
perms221
(
1
,
2
,
3
,
sigma3indices
);
perms2111
(
0
,
1
,
2
,
3
,
sigma0indices
);
perms2111
(
1
,
0
,
2
,
3
,
sigma1indices
);
perms2111
(
2
,
0
,
3
,
1
,
sigma2indices
);
perms2111
(
3
,
0
,
1
,
2
,
sigma3indices
);
if
(
sigma0indices
.
size
()
!=
256
){
std
::
cerr
<<
"sigma_index not set: "
;
std
::
cerr
<<
"sigma0indices has "
<<
sigma0indices
.
size
()
<<
" components"
<<
std
::
endl
;
return
false
;
}
else
if
(
sigma1indices
.
size
()
!=
256
){
std
::
cerr
<<
"sigma_index not set: "
;
std
::
cerr
<<
"sigma1indices has "
<<
sigma0indices
.
size
()
<<
" components"
<<
std
::
endl
;
return
false
;
}
else
if
(
sigma2indices
.
size
()
!=
256
){
std
::
cerr
<<
"sigma_index not set: "
;
std
::
cerr
<<
"sigma2indices has "
<<
sigma0indices
.
size
()
<<
" components"
<<
std
::
endl
;
return
false
;
}
else
if
(
sigma3indices
.
size
()
!=
256
){
std
::
cerr
<<
"sigma_index not set: "
;
std
::
cerr
<<
"sigma3indices has "
<<
sigma0indices
.
size
()
<<
" components"
<<
std
::
endl
;
return
false
;
}
for
(
int
i
=
0
;
i
<
256
;
i
++
){
// map each unique set of tensor indices to its position in a list
int
index0
=
tensor_to_list_idx
(
sigma0indices
.
at
(
i
));
int
index1
=
tensor_to_list_idx
(
sigma1indices
.
at
(
i
));
int
index2
=
tensor_to_list_idx
(
sigma2indices
.
at
(
i
));
int
index3
=
tensor_to_list_idx
(
sigma3indices
.
at
(
i
));
sigma_index5
[
index0
]
=
0
;
sigma_index5
[
index1
]
=
1
;
sigma_index5
[
index2
]
=
2
;
sigma_index5
[
index3
]
=
3
;
short
int
sign
[
4
]
=
{
1
,
-
1
,
-
1
,
-
1
};
// plus->true->1
helfactor5
[
1
][
index0
]
=
sign
[
sigma0indices
.
at
(
i
)[
1
]]
*
sign
[
sigma0indices
.
at
(
i
)[
3
]];
helfactor5
[
1
][
index1
]
=
sign
[
sigma1indices
.
at
(
i
)[
1
]]
*
sign
[
sigma1indices
.
at
(
i
)[
3
]];
helfactor5
[
1
][
index2
]
=
sign
[
sigma2indices
.
at
(
i
)[
1
]]
*
sign
[
sigma2indices
.
at
(
i
)[
3
]];
helfactor5
[
1
][
index3
]
=
sign
[
sigma3indices
.
at
(
i
)[
1
]]
*
sign
[
sigma3indices
.
at
(
i
)[
3
]];
// minus->false->0
helfactor5
[
0
][
index0
]
=
sign
[
sigma0indices
.
at
(
i
)[
0
]]
*
sign
[
sigma0indices
.
at
(
i
)[
2
]]
*
sign
[
sigma0indices
.
at
(
i
)[
4
]];
helfactor5
[
0
][
index1
]
=
sign
[
sigma1indices
.
at
(
i
)[
0
]]
*
sign
[
sigma1indices
.
at
(
i
)[
2
]]
*
sign
[
sigma1indices
.
at
(
i
)[
4
]];
helfactor5
[
0
][
index2
]
=
sign
[
sigma2indices
.
at
(
i
)[
0
]]
*
sign
[
sigma2indices
.
at
(
i
)[
2
]]
*
sign
[
sigma2indices
.
at
(
i
)[
4
]];
helfactor5
[
0
][
index3
]
=
sign
[
sigma3indices
.
at
(
i
)[
0
]]
*
sign
[
sigma3indices
.
at
(
i
)[
2
]]
*
sign
[
sigma3indices
.
at
(
i
)[
4
]];
}
// now set sigma_index3
std
::
vector
<
std
::
array
<
int
,
3
>>
sigma0indices3
;
std
::
vector
<
std
::
array
<
int
,
3
>>
sigma1indices3
;
std
::
vector
<
std
::
array
<
int
,
3
>>
sigma2indices3
;
std
::
vector
<
std
::
array
<
int
,
3
>>
sigma3indices3
;
// perms with 3 same
sigma0indices3
.
push_back
({
0
,
0
,
0
});
sigma1indices3
.
push_back
({
1
,
1
,
1
});
sigma2indices3
.
push_back
({
2
,
2
,
2
});
sigma3indices3
.
push_back
({
3
,
3
,
3
});
// 2 same
perms21
(
1
,
0
,
sigma0indices3
);
perms21
(
2
,
0
,
sigma0indices3
);
perms21
(
3
,
0
,
sigma0indices3
);
perms21
(
0
,
1
,
sigma1indices3
);
perms21
(
2
,
1
,
sigma1indices3
);
perms21
(
3
,
1
,
sigma1indices3
);
perms21
(
0
,
2
,
sigma2indices3
);
perms21
(
1
,
2
,
sigma2indices3
);
perms21
(
3
,
2
,
sigma2indices3
);
perms21
(
0
,
3
,
sigma3indices3
);
perms21
(
1
,
3
,
sigma3indices3
);
perms21
(
2
,
3
,
sigma3indices3
);
// none same
perms111
(
1
,
2
,
3
,
sigma0indices3
);
perms111
(
0
,
2
,
3
,
sigma1indices3
);
perms111
(
0
,
3
,
1
,
sigma2indices3
);
perms111
(
0
,
1
,
2
,
sigma3indices3
);
if
(
sigma0indices3
.
size
()
!=
16
){
std
::
cerr
<<
"sigma_index3 not set: "
;
std
::
cerr
<<
"sigma0indices3 has "
<<
sigma0indices3
.
size
()
<<
" components"
<<
std
::
endl
;
return
false
;
}
else
if
(
sigma1indices3
.
size
()
!=
16
){
std
::
cerr
<<
"sigma_index3 not set: "
;
std
::
cerr
<<
"sigma1indices3 has "
<<
sigma0indices3
.
size
()
<<
" components"
<<
std
::
endl
;
return
false
;
}
else
if
(
sigma2indices3
.
size
()
!=
16
){
std
::
cerr
<<
"sigma_index3 not set: "
;
std
::
cerr
<<
"sigma2indices3 has "
<<
sigma0indices3
.
size
()
<<
" components"
<<
std
::
endl
;
return
false
;
}
else
if
(
sigma3indices3
.
size
()
!=
16
){
std
::
cerr
<<
"sigma_index3 not set: "
;
std
::
cerr
<<
"sigma3indices3 has "
<<
sigma0indices3
.
size
()
<<
" components"
<<
std
::
endl
;
return
false
;
}
for
(
int
i
=
0
;
i
<
16
;
i
++
){
int
index0
=
tensor_to_list_idx
(
sigma0indices3
.
at
(
i
));
int
index1
=
tensor_to_list_idx
(
sigma1indices3
.
at
(
i
));
int
index2
=
tensor_to_list_idx
(
sigma2indices3
.
at
(
i
));
int
index3
=
tensor_to_list_idx
(
sigma3indices3
.
at
(
i
));
sigma_index3
[
index0
]
=
0
;
sigma_index3
[
index1
]
=
1
;
sigma_index3
[
index2
]
=
2
;
sigma_index3
[
index3
]
=
3
;
short
int
sign
[
4
]
=
{
1
,
-
1
,
-
1
,
-
1
};
// plus->true->1
helfactor3
[
1
][
index0
]
=
sign
[
sigma0indices3
.
at
(
i
)[
1
]];
helfactor3
[
1
][
index1
]
=
sign
[
sigma1indices3
.
at
(
i
)[
1
]];
helfactor3
[
1
][
index2
]
=
sign
[
sigma2indices3
.
at
(
i
)[
1
]];
helfactor3
[
1
][
index3
]
=
sign
[
sigma3indices3
.
at
(
i
)[
1
]];
// minus->false->0
helfactor3
[
0
][
index0
]
=
sign
[
sigma0indices3
.
at
(
i
)[
0
]]
*
sign
[
sigma0indices3
.
at
(
i
)[
2
]];
helfactor3
[
0
][
index1
]
=
sign
[
sigma1indices3
.
at
(
i
)[
0
]]
*
sign
[
sigma1indices3
.
at
(
i
)[
2
]];
helfactor3
[
0
][
index2
]
=
sign
[
sigma2indices3
.
at
(
i
)[
0
]]
*
sign
[
sigma2indices3
.
at
(
i
)[
2
]];
helfactor3
[
0
][
index3
]
=
sign
[
sigma3indices3
.
at
(
i
)[
0
]]
*
sign
[
sigma3indices3
.
at
(
i
)[
2
]];
}
return
true
;
}
// end init_sigma_index
}
void
init_sigma_index
(){
static
const
bool
initialised
=
init_sigma_index_helper
();
(
void
)
initialised
;
// shut up compiler warnings
}
}
File Metadata
Details
Attached
Mime Type
text/x-c
Expires
Tue, Sep 30, 5:48 AM (1 d, 8 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6503374
Default Alt Text
Tensor.cc (24 KB)
Attached To
Mode
rHEJ HEJ
Attached
Detach File
Event Timeline
Log In to Comment