check/exam_clifford.cpp

   1 /** @file exam_clifford.cpp
   2  *
   3  *  Here we test GiNaC's Clifford algebra objects. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  21  */
  22
  23 #include "exams.h"
  24
  25 const numeric half(1, 2);
  26
  27 static unsigned check_equal(const ex &e1, const ex &e2)
  28 {
  29         ex e = normal(e1 - e2);
  30         if (!e.is_zero()) {
  31                 clog << "(" << e1 << ") - (" << e2 << ") erroneously returned "
  32                      << e << " instead of 0" << endl;
  33                 return 1;
  34         }
  35         return 0;
  36 }
  37
  38 static unsigned check_equal_simplify(const ex &e1, const ex &e2)
  39 {
  40         ex e = normal(simplify_indexed(e1) - e2);
  41         if (!e.is_zero()) {
  42                 clog << "simplify_indexed(" << e1 << ") - (" << e2 << ") erroneously returned "
  43                          << e << " instead of 0" << endl;
  44                 return 1;
  45         }
  46         return 0;
  47 }
  48
  49 static unsigned check_equal_lst(const ex & e1, const ex & e2)
  50 {
  51         for (unsigned int i = 0; i < e1.nops(); i++) {
  52                 ex e = e1.op(i) - e2.op(i);
  53                 if (!e.normal().is_zero()) {
  54                         clog << "(" << e1 << ") - (" << e2 << ") erroneously returned "
  55                              << e << " instead of 0 (in the entry " << i  << ")" << endl;
  56                         return 1;
  57                 }
  58         }
  59         return 0;
  60 }
  61
  62 static unsigned check_equal_simplify_term(const ex & e1, const ex & e2, idx & mu)
  63 {
  64         ex e = expand_dummy_sum(normal(simplify_indexed(e1) - e2), true);
  65
  66         for (int j=0; j<4; j++) {
  67                 ex esub = e.subs(lst(is_a<varidx>(mu) ?
  68                                                          mu == idx(j, mu.get_dim()), ex_to<varidx>(mu).toggle_variance() == idx(j, mu.get_dim())
  69                                                          : mu == idx(j, mu.get_dim())));
  70                 if (!(canonicalize_clifford(esub).is_zero())) {
  71                         clog << "simplify_indexed(" << e1 << ") - (" << e2 << ") erroneously returned "
  72                                  << canonicalize_clifford(esub) << " instead of 0 for mu=" << j << endl;
  73                         return 1;
  74                 }
  75         }
  76         return 0;
  77 }
  78
  79 static unsigned check_equal_simplify_term2(const ex & e1, const ex & e2)
  80 {
  81         ex e = expand_dummy_sum(normal(simplify_indexed(e1) - e2), true);
  82         if (!(canonicalize_clifford(e).is_zero())) {
  83                 clog << "simplify_indexed(" << e1 << ") - (" << e2 << ") erroneously returned "
  84                          << canonicalize_clifford(e) << " instead of 0" << endl;
  85                 return 1;
  86         }
  87         return 0;
  88 }
  89
  90
  91 static unsigned clifford_check1()
  92 {
  93         // checks general identities and contractions
  94
  95         unsigned result = 0;
  96
  97         symbol dim("D");
  98         varidx mu(symbol("mu"), dim), nu(symbol("nu"), dim), rho(symbol("rho"), dim);
  99         ex e;
 100
 101         e = dirac_ONE() * dirac_ONE();
 102         result += check_equal(e, dirac_ONE());
 103
 104         e = dirac_ONE() * dirac_gamma(mu) * dirac_ONE();
 105         result += check_equal(e, dirac_gamma(mu));
 106
 107         e = dirac_gamma(varidx(2, dim)) * dirac_gamma(varidx(1, dim)) *
 108             dirac_gamma(varidx(1, dim)) * dirac_gamma(varidx(2, dim));
 109         result += check_equal(e, dirac_ONE());
 110
 111         e = dirac_gamma(mu) * dirac_gamma(nu) *
 112             dirac_gamma(nu.toggle_variance()) * dirac_gamma(mu.toggle_variance());
 113         result += check_equal_simplify(e, pow(dim, 2) * dirac_ONE());
 114
 115         e = dirac_gamma(mu) * dirac_gamma(nu) *
 116             dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance());
 117         result += check_equal_simplify(e, 2*dim*dirac_ONE()-pow(dim, 2)*dirac_ONE());
 118
 119         e = dirac_gamma(nu.toggle_variance()) * dirac_gamma(rho.toggle_variance()) *
 120             dirac_gamma(mu) * dirac_gamma(rho) * dirac_gamma(nu);
 121         e = e.simplify_indexed().collect(dirac_gamma(mu));
 122         result += check_equal(e, pow(2 - dim, 2).expand() * dirac_gamma(mu));
 123
 124         return result;
 125 }
 126
 127 static unsigned clifford_check2()
 128 {
 129         // checks identities relating to gamma5
 130
 131         unsigned result = 0;
 132
 133         symbol dim("D");
 134         varidx mu(symbol("mu"), dim), nu(symbol("nu"), dim);
 135         ex e;
 136
 137         e = dirac_gamma(mu) * dirac_gamma5() + dirac_gamma5() * dirac_gamma(mu);
 138         result += check_equal(e, 0);
 139
 140         e = dirac_gamma5() * dirac_gamma(mu) * dirac_gamma5() + dirac_gamma(mu);
 141         result += check_equal(e, 0);
 142
 143         return result;
 144 }
 145
 146 static unsigned clifford_check3()
 147 {
 148         // checks traces
 149
 150         unsigned result = 0;
 151
 152         symbol dim("D"), m("m"), q("q"), l("l"), ldotq("ldotq");
 153         varidx mu(symbol("mu"), dim), nu(symbol("nu"), dim), rho(symbol("rho"), dim),
 154                sig(symbol("sig"), dim), kap(symbol("kap"), dim), lam(symbol("lam"), dim);
 155         ex e;
 156
 157         e = dirac_gamma(mu);
 158         result += check_equal(dirac_trace(e), 0);
 159
 160         e = dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho);
 161         result += check_equal(dirac_trace(e), 0);
 162
 163         e = dirac_gamma5() * dirac_gamma(mu);
 164         result += check_equal(dirac_trace(e), 0);
 165
 166         e = dirac_gamma5() * dirac_gamma(mu) * dirac_gamma(nu);
 167         result += check_equal(dirac_trace(e), 0);
 168
 169         e = dirac_gamma5() * dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho);
 170         result += check_equal(dirac_trace(e), 0);
 171
 172         scalar_products sp;
 173         sp.add(q, q, pow(q, 2));
 174         sp.add(l, l, pow(l, 2));
 175         sp.add(l, q, ldotq);
 176
 177         e = pow(m, 2) * dirac_slash(q, dim) * dirac_slash(q, dim);
 178         e = dirac_trace(e).simplify_indexed(sp);
 179         result += check_equal(e, 4*pow(m, 2)*pow(q, 2));
 180
 181         // cyclicity without gamma5
 182         e = dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(sig)
 183           - dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(sig) * dirac_gamma(mu);
 184         e = dirac_trace(e);
 185         result += check_equal(e, 0);
 186
 187         e = dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(sig) * dirac_gamma(kap) * dirac_gamma(lam)
 188           - dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(sig) * dirac_gamma(kap) * dirac_gamma(lam) * dirac_gamma(mu);
 189         e = dirac_trace(e).expand();
 190         result += check_equal(e, 0);
 191
 192         // cyclicity of gamma5 * S_4
 193         e = dirac_gamma5() * dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(sig)
 194           - dirac_gamma(sig) * dirac_gamma5() * dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho);
 195         e = dirac_trace(e);
 196         result += check_equal(e, 0);
 197
 198         // non-cyclicity of order D-4 of gamma5 * S_6
 199         e = dirac_gamma5() * dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(sig) * dirac_gamma(kap) * dirac_gamma(mu.toggle_variance())
 200           + dim * dirac_gamma5() * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(sig) * dirac_gamma(kap);
 201         e = dirac_trace(e).simplify_indexed();
 202         e = (e / (dim - 4)).normal();
 203         result += check_equal(e, 8 * I * lorentz_eps(nu.replace_dim(4), rho.replace_dim(4), sig.replace_dim(4), kap.replace_dim(4)));
 204
 205         // one-loop vacuum polarization in QED
 206         e = dirac_gamma(mu) *
 207             (dirac_slash(l, dim) + dirac_slash(q, 4) + m * dirac_ONE()) *
 208             dirac_gamma(mu.toggle_variance()) *
 209             (dirac_slash(l, dim) + m * dirac_ONE());
 210         e = dirac_trace(e).simplify_indexed(sp);
 211         result += check_equal(e, 4*((2-dim)*l*l + (2-dim)*ldotq + dim*m*m).expand());
 212
 213         e = dirac_slash(q, 4) *
 214             (dirac_slash(l, dim) + dirac_slash(q, 4) + m * dirac_ONE()) *
 215             dirac_slash(q, 4) *
 216             (dirac_slash(l, dim) + m * dirac_ONE());
 217         e = dirac_trace(e).simplify_indexed(sp);
 218         result += check_equal(e, 4*(2*ldotq*ldotq + q*q*ldotq - q*q*l*l + q*q*m*m).expand());
 219
 220         // stuff that had problems in the past
 221         ex prop = dirac_slash(q, dim) - m * dirac_ONE();
 222         e = dirac_slash(l, dim) * dirac_gamma5() * dirac_slash(l, dim) * prop;
 223         e = dirac_trace(dirac_slash(q, dim) * e) - dirac_trace(m * e)
 224           - dirac_trace(prop * e);
 225         result += check_equal(e, 0);
 226
 227         e = (dirac_gamma5() + dirac_ONE()) * dirac_gamma5();
 228         e = dirac_trace(e);
 229         result += check_equal(e, 4);
 230
 231         // traces with multiple representation labels
 232         e = dirac_ONE(0) * dirac_ONE(1) / 16;
 233         result += check_equal(dirac_trace(e, 0), dirac_ONE(1) / 4);
 234         result += check_equal(dirac_trace(e, 1), dirac_ONE(0) / 4);
 235         result += check_equal(dirac_trace(e, 2), e);
 236         result += check_equal(dirac_trace(e, lst(0, 1)), 1);
 237
 238         e = dirac_gamma(mu, 0) * dirac_gamma(mu.toggle_variance(), 1) * dirac_gamma(nu, 0) * dirac_gamma(nu.toggle_variance(), 1);
 239         result += check_equal_simplify(dirac_trace(e, 0), 4 * dim * dirac_ONE(1));
 240         result += check_equal_simplify(dirac_trace(e, 1), 4 * dim * dirac_ONE(0));
 241         // Fails with new tinfo mechanism because the order of gamme matrices with different rl depends on luck.
 242         // TODO: better check.
 243         //result += check_equal_simplify(dirac_trace(e, 2), canonicalize_clifford(e)); // e will be canonicalized by the calculation of the trace
 244         result += check_equal_simplify(dirac_trace(e, lst(0, 1)), 16 * dim);
 245
 246         return result;
 247 }
 248
 249 static unsigned clifford_check4()
 250 {
 251         // simplify_indexed()/dirac_trace() cross-checks
 252
 253         unsigned result = 0;
 254
 255         symbol dim("D");
 256         varidx mu(symbol("mu"), dim), nu(symbol("nu"), dim), rho(symbol("rho"), dim),
 257                sig(symbol("sig"), dim), lam(symbol("lam"), dim);
 258         ex e, t1, t2;
 259
 260         e = dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(mu.toggle_variance());
 261         t1 = dirac_trace(e).simplify_indexed();
 262         t2 = dirac_trace(e.simplify_indexed());
 263         result += check_equal((t1 - t2).expand(), 0);
 264
 265         e = dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(sig) * dirac_gamma(mu.toggle_variance()) * dirac_gamma(lam);
 266         t1 = dirac_trace(e).simplify_indexed();
 267         t2 = dirac_trace(e.simplify_indexed());
 268         result += check_equal((t1 - t2).expand(), 0);
 269
 270         e = dirac_gamma(sig) * dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(nu.toggle_variance()) * dirac_gamma(mu.toggle_variance());
 271         t1 = dirac_trace(e).simplify_indexed();
 272         t2 = dirac_trace(e.simplify_indexed());
 273         result += check_equal((t1 - t2).expand(), 0);
 274
 275         e = dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(rho) * dirac_gamma(mu.toggle_variance()) * dirac_gamma(sig) * dirac_gamma(nu.toggle_variance());
 276         t1 = dirac_trace(e).simplify_indexed();
 277         t2 = dirac_trace(e.simplify_indexed());
 278         result += check_equal((t1 - t2).expand(), 0);
 279
 280         return result;
 281 }
 282
 283 static unsigned clifford_check5()
 284 {
 285         // canonicalize_clifford() checks
 286
 287         unsigned result = 0;
 288
 289         symbol dim("D");
 290         varidx mu(symbol("mu"), dim), nu(symbol("nu"), dim), lam(symbol("lam"), dim);
 291         ex e;
 292
 293         e = dirac_gamma(mu) * dirac_gamma(nu) + dirac_gamma(nu) * dirac_gamma(mu);
 294         result += check_equal(canonicalize_clifford(e), 2*dirac_ONE()*lorentz_g(mu, nu));
 295
 296         e = (dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(lam)
 297            + dirac_gamma(nu) * dirac_gamma(lam) * dirac_gamma(mu)
 298            + dirac_gamma(lam) * dirac_gamma(mu) * dirac_gamma(nu)
 299            - dirac_gamma(nu) * dirac_gamma(mu) * dirac_gamma(lam)
 300            - dirac_gamma(lam) * dirac_gamma(nu) * dirac_gamma(mu)
 301            - dirac_gamma(mu) * dirac_gamma(lam) * dirac_gamma(nu)) / 6
 302           + lorentz_g(mu, nu) * dirac_gamma(lam)
 303           - lorentz_g(mu, lam) * dirac_gamma(nu)
 304           + lorentz_g(nu, lam) * dirac_gamma(mu)
 305           - dirac_gamma(mu) * dirac_gamma(nu) * dirac_gamma(lam);
 306         result += check_equal(canonicalize_clifford(e), 0);
 307
 308         return result;
 309 }
 310
 311 /* We make two identical checks with metrics defined through a matrix in
 312  * the cases when used indexes have or have not variance.
 313  * To this end we recycle the code through the following macros */
 314
 315 #define CHECK6(IDX,TOGGLE) {IDX v(symbol("v"), 4), nu(symbol("nu"), 4), mu(symbol("mu"), 4), \
 316                psi(symbol("psi"),4), lam(symbol("lambda"), 4),\
 317                xi(symbol("xi"), 4),  rho(symbol("rho"),4);\
 318 \
 319 /* checks general identities and contractions for clifford_unit*/\
 320         e = dirac_ONE(2) * clifford_unit(mu, A, 2) * dirac_ONE(2);\
 321         result += check_equal(e, clifford_unit(mu, A, 2));\
 322 \
 323         e = clifford_unit(idx(2, 4), A) * clifford_unit(idx(1, 4), A)\
 324           * clifford_unit(idx(1, 4), A) * clifford_unit(idx(2, 4), A);\
 325         result += check_equal(e, A(1, 1) * A(2, 2) * dirac_ONE());\
 326 \
 327         e = clifford_unit(IDX(2, 4), A) * clifford_unit(IDX(1, 4), A)\
 328           * clifford_unit(IDX(1, 4), A) * clifford_unit(IDX(2, 4), A);\
 329         result += check_equal(e, A(1, 1) * A(2, 2) * dirac_ONE());\
 330 \
 331         e = clifford_unit(nu, A) * clifford_unit(nu TOGGLE, A);\
 332         result += check_equal_simplify(e, A.trace() * dirac_ONE());\
 333 \
 334         e = clifford_unit(nu, A) * clifford_unit(nu, A);\
 335         result += check_equal_simplify(e, indexed(A_symm, sy_symm(), nu, nu) * dirac_ONE());\
 336 \
 337         e = clifford_unit(nu, A) * clifford_unit(nu TOGGLE, A) * clifford_unit(mu, A);\
 338         result += check_equal_simplify(e, A.trace() * clifford_unit(mu, A));\
 339 \
 340         e = clifford_unit(nu, A) * clifford_unit(mu, A) * clifford_unit(nu TOGGLE, A);\
 341         \
 342         result += check_equal_simplify_term(e,  2 * indexed(A_symm, sy_symm(), nu TOGGLE, mu) *clifford_unit(nu, A)-A.trace()*clifford_unit(mu, A), mu);\
 343 \
 344         e = clifford_unit(nu, A) * clifford_unit(nu TOGGLE, A)\
 345           * clifford_unit(mu, A) * clifford_unit(mu TOGGLE, A);\
 346         result += check_equal_simplify(e, pow(A.trace(), 2) * dirac_ONE());\
 347 \
 348         e = clifford_unit(mu, A) * clifford_unit(nu, A)\
 349           * clifford_unit(nu TOGGLE, A) * clifford_unit(mu TOGGLE, A);\
 350         result += check_equal_simplify(e, pow(A.trace(), 2)  * dirac_ONE());\
 351 \
 352         e = clifford_unit(mu, A) * clifford_unit(nu, A)\
 353           * clifford_unit(mu TOGGLE, A) * clifford_unit(nu TOGGLE, A);\
 354 \
 355         result += check_equal_simplify_term2(e, 2*indexed(A_symm, sy_symm(), nu TOGGLE, mu TOGGLE) * clifford_unit(mu, A) * clifford_unit(nu, A) - pow(A.trace(), 2)*dirac_ONE());\
 356 \
 357         e = clifford_unit(mu TOGGLE, A) * clifford_unit(nu, A)\
 358           * clifford_unit(mu, A) * clifford_unit(nu TOGGLE, A);\
 359 \
 360         result += check_equal_simplify_term2(e, 2*indexed(A_symm, nu, mu) * clifford_unit(mu TOGGLE, A) * clifford_unit(nu TOGGLE, A) - pow(A.trace(), 2)*dirac_ONE());\
 361 \
 362         e = clifford_unit(nu TOGGLE, A) * clifford_unit(rho TOGGLE, A)\
 363           * clifford_unit(mu, A) * clifford_unit(rho, A) * clifford_unit(nu, A);\
 364         e = e.simplify_indexed().collect(clifford_unit(mu, A));\
 365         \
 366         result += check_equal_simplify_term(e, 4* indexed(A_symm, sy_symm(), nu TOGGLE,  rho)*indexed(A_symm, sy_symm(), rho TOGGLE, mu) *clifford_unit(nu, A) \
 367                                             - 2*A.trace() * (clifford_unit(rho, A) * indexed(A_symm, sy_symm(), rho TOGGLE, mu) \
 368                                                              + clifford_unit(nu, A) * indexed(A_symm, sy_symm(), nu TOGGLE, mu)) + pow(A.trace(),2)* clifford_unit(mu, A), mu);\
 369 \
 370         e = clifford_unit(nu TOGGLE, A) * clifford_unit(rho, A)\
 371           * clifford_unit(mu, A) * clifford_unit(rho TOGGLE, A) * clifford_unit(nu, A);\
 372         e = e.simplify_indexed().collect(clifford_unit(mu, A));\
 373         \
 374         result += check_equal_simplify_term(e, 4* indexed(A_symm, sy_symm(), nu TOGGLE,  rho)*indexed(A_symm, sy_symm(), rho TOGGLE, mu) *clifford_unit(nu, A) \
 375                                             - 2*A.trace() * (clifford_unit(rho, A) * indexed(A_symm, sy_symm(), rho TOGGLE, mu) \
 376                                                              + clifford_unit(nu, A) * indexed(A_symm, sy_symm(), nu TOGGLE, mu)) + pow(A.trace(),2)* clifford_unit(mu, A), mu);\
 377 \
 378         e = clifford_unit(mu, A) * clifford_unit(nu, A) + clifford_unit(nu, A) * clifford_unit(mu, A);\
 379         result += check_equal(canonicalize_clifford(e), 2*dirac_ONE()*indexed(A_symm, sy_symm(), mu, nu));\
 380 \
 381         e = (clifford_unit(mu, A) * clifford_unit(nu, A) * clifford_unit(lam, A)\
 382                  + clifford_unit(nu, A) * clifford_unit(lam, A) * clifford_unit(mu, A)\
 383                  + clifford_unit(lam, A) * clifford_unit(mu, A) * clifford_unit(nu, A)\
 384                  - clifford_unit(nu, A) * clifford_unit(mu, A) * clifford_unit(lam, A)\
 385                  - clifford_unit(lam, A) * clifford_unit(nu, A) * clifford_unit(mu, A)\
 386                  - clifford_unit(mu, A) * clifford_unit(lam, A) * clifford_unit(nu, A)) / 6\
 387                 + indexed(A_symm, sy_symm(), mu, nu) * clifford_unit(lam, A)\
 388                 - indexed(A_symm, sy_symm(), mu, lam) * clifford_unit(nu, A)\
 389                 + indexed(A_symm, sy_symm(), nu, lam) * clifford_unit(mu, A)\
 390                 - clifford_unit(mu, A) * clifford_unit(nu, A) * clifford_unit(lam, A);\
 391         result += check_equal(canonicalize_clifford(e), 0);\
 392 \
 393 /* lst_to_clifford() and clifford_inverse()  check*/\
 394         realsymbol x("x"), y("y"), t("t"), z("z");\
 395         \
 396         ex c = clifford_unit(nu, A, 1);\
 397         e = lst_to_clifford(lst(t, x, y, z), mu, A, 1) * lst_to_clifford(lst(1, 2, 3, 4), c);\
 398         e1 = clifford_inverse(e);\
 399         result += check_equal_simplify_term2((e*e1).simplify_indexed(), dirac_ONE(1));\
 400 \
 401 /* Moebius map (both forms) checks for symmetric metrics only */\
 402         matrix M1(2, 2),  M2(2, 2);\
 403         c = clifford_unit(nu, A);\
 404         \
 405         e = clifford_moebius_map(0, dirac_ONE(), \
 406                                                          dirac_ONE(), 0, lst(t, x, y, z), A); \
 407 /* this is just the inversion*/\
 408         M1 = 0, dirac_ONE(),\
 409                 dirac_ONE(), 0;\
 410         e1 = clifford_moebius_map(M1, lst(t, x, y, z), A); \
 411 /* the inversion again*/\
 412         result += check_equal_lst(e, e1);\
 413         \
 414         e1 = clifford_to_lst(clifford_inverse(lst_to_clifford(lst(t, x, y, z), mu, A)), c);\
 415         result += check_equal_lst(e, e1);\
 416         \
 417         e = clifford_moebius_map(dirac_ONE(), lst_to_clifford(lst(1, 2, 3, 4), nu, A), \
 418                                                          0, dirac_ONE(), lst(t, x, y, z), A); \
 419 /*this is just a shift*/\
 420         M2 = dirac_ONE(), lst_to_clifford(lst(1, 2, 3, 4), c),\
 421                 0, dirac_ONE();\
 422         e1 = clifford_moebius_map(M2, lst(t, x, y, z), c); \
 423 /* the same shift*/\
 424         result += check_equal_lst(e, e1);\
 425                 \
 426         result += check_equal(e, lst(t+1, x+2, y+3, z+4));\
 427         \
 428 /* Check the group law for Moebius maps */\
 429         e = clifford_moebius_map(M1, ex_to<lst>(e1), c); \
 430 /*composition of M1 and M2*/\
 431         e1 = clifford_moebius_map(M1.mul(M2), lst(t, x, y, z), c); \
 432 /* the product M1*M2*/\
 433         result += check_equal_lst(e, e1);}
 434
 435 static unsigned clifford_check6(const matrix & A)
 436 {
 437         matrix A_symm(4,4), A2(4, 4);
 438         A_symm = A.add(A.transpose()).mul(half);
 439         A2 = A_symm.mul(A_symm);
 440
 441         ex e, e1;
 442         int result = 0;
 443
 444         CHECK6(varidx,.toggle_variance())
 445
 446         return result;
 447 }
 448
 449 static unsigned clifford_check6a(const matrix & A)
 450 {
 451         matrix A_symm(4,4), A2(4, 4);
 452         A_symm = A.add(A.transpose()).mul(half);
 453         A2 = A_symm.mul(A_symm);
 454
 455         ex e, e1;
 456         int result = 0;
 457
 458         CHECK6(idx,)
 459
 460         return result;
 461 }
 462
 463 static unsigned clifford_check7(const ex & G, const symbol & dim)
 464 {
 465         // checks general identities and contractions
 466
 467         unsigned result = 0;
 468
 469         varidx mu(symbol("mu"), dim), nu(symbol("nu"), dim), rho(symbol("rho"), dim),
 470                psi(symbol("psi"),dim), lam(symbol("lambda"), dim), xi(symbol("xi"), dim);
 471
 472         ex e;
 473         clifford unit = ex_to<clifford>(clifford_unit(mu, G));
 474         ex scalar = unit.get_metric(varidx(0, dim), varidx(0, dim));
 475
 476         e = dirac_ONE() * dirac_ONE();
 477         result += check_equal(e, dirac_ONE());
 478
 479         e = dirac_ONE() * clifford_unit(mu, G) * dirac_ONE();
 480         result += check_equal(e, clifford_unit(mu, G));
 481
 482         e = clifford_unit(varidx(2, dim), G) * clifford_unit(varidx(1, dim), G)
 483           * clifford_unit(varidx(1, dim), G) * clifford_unit(varidx(2, dim), G);
 484         result += check_equal(e, dirac_ONE()*pow(scalar, 2));
 485
 486         e = clifford_unit(mu, G) * clifford_unit(nu, G)
 487           * clifford_unit(nu.toggle_variance(), G) * clifford_unit(mu.toggle_variance(), G);
 488         result += check_equal_simplify(e, pow(dim*scalar, 2) * dirac_ONE());
 489
 490         e = clifford_unit(mu, G) * clifford_unit(nu, G)
 491           * clifford_unit(mu.toggle_variance(), G) * clifford_unit(nu.toggle_variance(), G);
 492         result += check_equal_simplify(e, (2*dim - pow(dim, 2))*pow(scalar,2)*dirac_ONE());
 493
 494         e = clifford_unit(nu.toggle_variance(), G) * clifford_unit(rho.toggle_variance(), G)
 495           * clifford_unit(mu, G) * clifford_unit(rho, G) * clifford_unit(nu, G);
 496         e = e.simplify_indexed().collect(clifford_unit(mu, G));
 497         result += check_equal(e, pow(scalar*(dim-2), 2).expand() * clifford_unit(mu, G));
 498
 499         // canonicalize_clifford() checks, only for symmetric metrics
 500         if (ex_to<symmetry>(ex_to<indexed>(ex_to<clifford>(clifford_unit(mu, G)).get_metric()).get_symmetry()).has_symmetry()) {
 501                 e = clifford_unit(mu, G) * clifford_unit(nu, G) + clifford_unit(nu, G) * clifford_unit(mu, G);
 502                 result += check_equal(canonicalize_clifford(e), 2*dirac_ONE()*unit.get_metric(nu, mu));
 503
 504                 e = (clifford_unit(mu, G) * clifford_unit(nu, G) * clifford_unit(lam, G)
 505                          + clifford_unit(nu, G) * clifford_unit(lam, G) * clifford_unit(mu, G)
 506                          + clifford_unit(lam, G) * clifford_unit(mu, G) * clifford_unit(nu, G)
 507                          - clifford_unit(nu, G) * clifford_unit(mu, G) * clifford_unit(lam, G)
 508                          - clifford_unit(lam, G) * clifford_unit(nu, G) * clifford_unit(mu, G)
 509                          - clifford_unit(mu, G) * clifford_unit(lam, G) * clifford_unit(nu, G)) / 6
 510                         + unit.get_metric(mu, nu) * clifford_unit(lam, G)
 511                         - unit.get_metric(mu, lam) * clifford_unit(nu, G)
 512                         + unit.get_metric(nu, lam) * clifford_unit(mu, G)
 513                         - clifford_unit(mu, G) * clifford_unit(nu, G) * clifford_unit(lam, G);
 514                 result += check_equal(canonicalize_clifford(e), 0);
 515         } else {
 516                 e = clifford_unit(mu, G) * clifford_unit(nu, G) + clifford_unit(nu, G) * clifford_unit(mu, G);
 517                 result += check_equal(canonicalize_clifford(e), dirac_ONE()*(unit.get_metric(mu, nu) + unit.get_metric(nu, mu)));
 518
 519                 e = (clifford_unit(mu, G) * clifford_unit(nu, G) * clifford_unit(lam, G)
 520                          + clifford_unit(nu, G) * clifford_unit(lam, G) * clifford_unit(mu, G)
 521                          + clifford_unit(lam, G) * clifford_unit(mu, G) * clifford_unit(nu, G)
 522                          - clifford_unit(nu, G) * clifford_unit(mu, G) * clifford_unit(lam, G)
 523                          - clifford_unit(lam, G) * clifford_unit(nu, G) * clifford_unit(mu, G)
 524                          - clifford_unit(mu, G) * clifford_unit(lam, G) * clifford_unit(nu, G)) / 6
 525                         + half * (unit.get_metric(mu, nu) + unit.get_metric(nu, mu)) * clifford_unit(lam, G)
 526                         - half * (unit.get_metric(mu, lam) + unit.get_metric(lam, mu)) * clifford_unit(nu, G)
 527                         + half * (unit.get_metric(nu, lam) + unit.get_metric(lam, nu)) * clifford_unit(mu, G)
 528                         - clifford_unit(mu, G) * clifford_unit(nu, G) * clifford_unit(lam, G);
 529                 result += check_equal(canonicalize_clifford(e), 0);
 530         }
 531         return result;
 532 }
 533
 534 unsigned exam_clifford()
 535 {
 536         unsigned result = 0;
 537
 538         cout << "examining clifford objects" << flush;
 539         clog << "----------clifford objects:" << endl;
 540
 541         result += clifford_check1(); cout << '.' << flush;
 542         result += clifford_check2(); cout << '.' << flush;
 543         result += clifford_check3(); cout << '.' << flush;
 544         result += clifford_check4(); cout << '.' << flush;
 545         result += clifford_check5(); cout << '.' << flush;
 546
 547         // anticommuting, symmetric examples
 548         result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, 1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, 1))));; cout << '.' << flush;
 549         result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, -1, -1, -1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, -1, -1, -1))));; cout << '.' << flush;
 550         result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, -1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, -1))));; cout << '.' << flush;
 551         result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, 0, 1, -1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, 0, 1, -1))));; cout << '.' << flush;
 552         result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-3, 0, 2, -1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-3, 0, 2, -1))));; cout << '.' << flush;
 553
 554         realsymbol s("s"), t("t"); // symbolic entries in matric
 555         result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, 1, s, t))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, 1, s, t))));; cout << '.' << flush;
 556
 557         matrix A(4, 4);
 558         A = 1, 0, 0, 0, // anticommuting, not symmetric, Tr=0
 559                 0, -1, 0, 0,
 560                 0, 0, 0, -1,
 561                 0, 0, 1, 0;
 562         result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
 563
 564         A = 1, 0, 0, 0, // anticommuting, not symmetric, Tr=2
 565                 0, 1, 0, 0,
 566                 0, 0, 0, -1,
 567                 0, 0, 1, 0;
 568         result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
 569
 570         A = 1, 0, 0, 0, // not anticommuting, symmetric, Tr=0
 571                 0, -1, 0, 0,
 572                 0, 0, 0, -1,
 573                 0, 0, -1, 0;
 574         result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
 575
 576         A = 1, 0, 0, 0, // not anticommuting, symmetric, Tr=2
 577                 0, 1, 0, 0,
 578                 0, 0, 0, -1,
 579                 0, 0, -1, 0;
 580         result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
 581
 582         A = 1, 1, 0, 0, // not anticommuting, not symmetric, Tr=4
 583                 0, 1, 1, 0,
 584                 0, 0, 1, 1,
 585                 0, 0, 0, 1;
 586         result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
 587
 588         symbol dim("D");
 589         result += clifford_check7(minkmetric(), dim); cout << '.' << flush;
 590
 591         varidx chi(symbol("chi"), dim), xi(symbol("xi"), dim);
 592         result += clifford_check7(delta_tensor(xi, chi), dim); cout << '.' << flush;
 593
 594         result += clifford_check7(lorentz_g(xi, chi), dim); cout << '.' << flush;
 595
 596         result += clifford_check7(indexed(-2*minkmetric(), sy_symm(), xi, chi), dim); cout << '.' << flush;
 597         result += clifford_check7(-2*delta_tensor(xi, chi), dim); cout << '.' << flush;
 598
 599         if (!result) {
 600                 cout << " passed " << endl;
 601                 clog << "(no output)" << endl;
 602         } else {
 603                 cout << " failed " << endl;
 604         }
 605
 606         return result;
 607 }