awdiagcub.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// A cubical version of the Alexander-Whitney diagonal.
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2009-2016 by Pawel Pilarczyk.
//
// This file is part of my research software package. This is free software:
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// This software is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this software; see the file "license.txt". If not,
// please, see <http://www.gnu.org/licenses/>.

// Started on March 24, 2009. Last revision: April 11, 2013.


#ifndef _CHAINCON_AWDIAGCUB_H_
#define _CHAINCON_AWDIAGCUB_H_


// include some standard C++ header files
#include <istream>
#include <ostream>

// include selected header files from the CHomP library
#include "chomp/system/config.h"

// include relevant local header files
#include "chaincon/cubcell.h"
#include "chaincon/combtensor.h"
#include "chaincon/tensor.h"


// --------------------------------------------------
// ------------------ AW diagonal -------------------
// --------------------------------------------------

/// Computes the Alexander-Whitney diagonal of a cubical cell
/// using the formula from Serre's paper. A purely combinatorial version.
template <class CoordT, class WrapT, class EmptyT>
inline void AWdiagonal (const tCubCell<CoordT,WrapT,EmptyT> &c,
        tCombTensor<tCubCell<CoordT,WrapT,EmptyT>,
        tCubCell<CoordT,WrapT,EmptyT> > &t)
{
        typedef tCubCell<CoordT,WrapT,EmptyT> CubeType;

        // determine the dimension of the cell and of the space
        int dim = c. dim ();
        if (dim > static_cast<int> (8 * sizeof (long) - 2))
                throw "Dimension too high for cubical A-W diagonal.";
        int spaceDim = c. spaceDim ();

        // prepare buffers for coordinates of pairs of new cubes
        CoordT *coordBuffer = new CoordT [6 * spaceDim + 1];
        CoordT *cLeft = coordBuffer;
        CoordT *cRight = cLeft + spaceDim;
        CoordT *leftLeft = cRight + spaceDim;
        CoordT *leftRight = leftLeft + spaceDim;
        CoordT *rightLeft = leftRight + spaceDim;
        CoordT *rightRight = rightLeft + spaceDim;

        // determine the indices of nonzero width of the cubical cell
        int *indices = new int [dim + 1];
        int indCur = 0;
        for (int i = 0; i < spaceDim; ++ i)
        {
                cLeft [i] = c. left (i);
                cRight [i] = c. right (i);
                if (cLeft [i] == cRight [i])
                        continue;
                indices [indCur ++] = i;
        }
        if (indCur != dim)
                throw "Wrong dimension of a cubical cell detected in AW.";

        // go through all the subsets and add up the tensor products
        const long subsetMax = 1 << dim;
        for (long subset = 0; subset < subsetMax; ++ subset)
        {
                // create the two cubical cells
                for (int i = 0; i < spaceDim; ++ i)
                {
                        leftLeft [i] = cLeft [i];
                        leftRight [i] = cRight [i];
                        rightLeft [i] = cLeft [i];
                        rightRight [i] = cRight [i];
                }
                for (int i = 0; i < dim; ++ i)
                {
                        int index = indices [i];
                        if (subset & (1 << i))
                        {
                                rightLeft [index] = rightRight [index];
                        }
                        else
                        {
                                leftRight [index] = leftLeft [index];
                        }
                }
                t. add (CubeType (spaceDim, leftLeft, leftRight),
                        CubeType (spaceDim, rightLeft, rightRight));
        }

        // release the allocated memory
        delete [] indices;
        delete [] coordBuffer;

        return;
} /* AWdiagonal */

/// Returns the sign rho of a term in the A-W diagonal
/// (see the paper for details).
inline int AW_rho (long subset, int dim)
{
        // TODO: more efficient
        int nu = 0;
        for (int i = 0; i < dim; ++ i)
        {
                for (int j = 0; j < i; ++ j)
                {
                        if (!(subset & (1 << i)))
                                continue;
                        if (subset & (1 << j))
                                continue;
                        ++ nu;
                }
        }
        return (nu & 1) ? -1 : 1;
} /* AW_rho */

/// Computes the Alexander-Whitney diagonal of a cubical cell
/// using the formula from Serre's paper. A general version.
template <class CoordT, class WrapT, class EmptyT, class CoefT>
inline void AWdiagonal (const tCubCell<CoordT,WrapT,EmptyT> &c,
        tTensor<tCubCell<CoordT,WrapT,EmptyT>,
        tCubCell<CoordT,WrapT,EmptyT>,CoefT> &t)
{
        typedef tCubCell<CoordT,WrapT,EmptyT> CubeType;

        // determine the dimension of the cell and of the space
        int dim = c. dim ();
        if (dim > static_cast<int> (8 * sizeof (long) - 2))
                throw "Dimension too high for cubical A-W diagonal.";
        int spaceDim = c. spaceDim ();

        // prepare buffers for coordinates of pairs of new cubes
        CoordT *coordBuffer = new CoordT [6 * spaceDim + 1];
        CoordT *cLeft = coordBuffer;
        CoordT *cRight = cLeft + spaceDim;
        CoordT *leftLeft = cRight + spaceDim;
        CoordT *leftRight = leftLeft + spaceDim;
        CoordT *rightLeft = leftRight + spaceDim;
        CoordT *rightRight = rightLeft + spaceDim;

        // determine the indices of nonzero width of the cubical cell
        int *indices = new int [dim + 1];
        int indCur = 0;
        for (int i = 0; i < spaceDim; ++ i)
        {
                cLeft [i] = c. left (i);
                cRight [i] = c. right (i);
                if (cLeft [i] == cRight [i])
                        continue;
                indices [indCur ++] = i;
        }
        if (indCur != dim)
                throw "Wrong dimension of a cubical cell detected in AW.";

        // go through all the subsets and add up the tensor products
        const long subsetMax = 1 << dim;
        for (long subset = 0; subset < subsetMax; ++ subset)
        {
                // create the two cubical cells
                for (int i = 0; i < spaceDim; ++ i)
                {
                        leftLeft [i] = cLeft [i];
                        leftRight [i] = cRight [i];
                        rightLeft [i] = cLeft [i];
                        rightRight [i] = cRight [i];
                }
                for (int i = 0; i < dim; ++ i)
                {
                        int index = indices [i];
                        if (subset & (1 << i))
                        {
                                rightLeft [index] = rightRight [index];
                        }
                        else
                        {
                                leftRight [index] = leftLeft [index];
                        }
                }
                CoefT coef (1);
                if (AW_rho (subset, dim) == -1)
                        coef. negate ();
                t. add (CubeType (spaceDim, leftLeft, leftRight),
                        CubeType (spaceDim, rightLeft, rightRight), coef);
        }

        // release the allocated memory
        delete [] indices;
        delete [] coordBuffer;

        return;
} /* AWdiagonal */


#endif // _CHAINCON_AWDIAGCUB_H_