YoungTab.hh

Go to the documentation of this file.

/* 

    Cadabra: a field-theory motivated computer algebra system.
    Copyright (C) 2001-2011  Kasper Peeters <kasper.peeters@aei.mpg.de>

   This program 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 program 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 program.  If not, see <http://www.gnu.org/licenses/>.
 
*/

/*  
    - TODO: has_nullifying trace is wrong, but needs to be merged with the
            input_asym code in order to be more useful.
    
*/

#pragma once

#include <cstddef>
#include <iostream>
#include <iterator>
#include <vector>
#include <list>
#include <gmpxx.h>
#include "Combinatorics.hh"
#include <cstddef>

typedef mpz_class yngint_t;
typedef mpq_class yngrat_t;

namespace yngtab {

// The tableau_base is the abstract interface; does not depend on the
// actual storage format.

class tableau_base {
    public:
        tableau_base();
        virtual ~tableau_base();
        virtual unsigned int number_of_rows() const=0;
        virtual unsigned int row_size(unsigned int row) const=0;
        virtual unsigned int column_size(unsigned int col) const; // FIXME: maybe make pure virt too
        virtual void         add_box(unsigned int row)=0;
        virtual void         remove_box(unsigned int row)=0;
        virtual void         add_row(unsigned int row_size);
        virtual void         clear()=0;
        
        yngrat_t             multiplicity;    // also keeps track of signs
        int                  selfdual_column; // -n, 0, n  for antiselfdual, no, selfdual (count from 1)
       yngint_t             dimension(unsigned int) const;
        unsigned long        hook_length(unsigned int row, unsigned int col) const;
        yngint_t             hook_length_prod() const;
};

class tableau : public tableau_base {
    public:
        virtual ~tableau();
        virtual unsigned int number_of_rows() const;
        virtual unsigned int row_size(unsigned int row) const;
        virtual void         add_box(unsigned int row);
        virtual void         remove_box(unsigned int row);
        virtual void         clear();

        tableau& operator=(const tableau&); 
    private:
        std::vector<int> rows;
};

template<class T>
class tableaux;

template<class T>
class filled_tableau : public tableau {
    public:
        typedef T value_type;

        virtual ~filled_tableau();
        virtual unsigned int number_of_rows() const;
        virtual unsigned int row_size(unsigned int row) const;
        virtual void         add_box(unsigned int row);
        virtual void         remove_box(unsigned int row);
        std::pair<int, int>  find(const T&) const;
        virtual void         clear();
        
        void                 copy_shape(const tableau&);
        
        T&                   operator()(unsigned int row, unsigned int col);
        const T&             operator()(unsigned int row, unsigned int col) const;
        const T&             operator[](unsigned int boxnum) const;
        void                 add_box(unsigned int rownum, T val);
        void                 swap_columns(unsigned int c1, unsigned int c2);

        bool                 compare_without_multiplicity(const filled_tableau<T>& other) const;
        bool                 has_nullifying_trace() const;
        void                 sort_within_columns();
        void                 sort_columns();
        void                 canonicalise();
        std::pair<int, int>  nonstandard_loc() const;
        template<class StrictWeakOrdering> void sort_within_columns(StrictWeakOrdering comp);
        template<class StrictWeakOrdering> void sort_columns(StrictWeakOrdering comp);
        template<class StrictWeakOrdering> void canonicalise(StrictWeakOrdering comp, bool only_col_ex=false);
        void                 projector(combin::symmetriser<T>&, bool modulo_monoterm=false) const;
        void                 projector(combin::symmetriser<T>&, combin::range_vector_t&) const;
        yngrat_t             projector_normalisation() const;

        filled_tableau<T>& operator=(const filled_tableau<T>&); 

        class iterator_base {
            public:
                typedef T                               value_type;
                typedef T*                              pointer;
                typedef T&                              reference;
                typedef size_t                          size_type;
                typedef ptrdiff_t                       difference_type;
                typedef std::random_access_iterator_tag iterator_category;
        };

        class in_column_iterator : public iterator_base {
            public:
                in_column_iterator(unsigned int r, unsigned int c, filled_tableau<T> *);
                T&                  operator*() const;
                T*                  operator->() const;
                in_column_iterator& operator++();
                in_column_iterator  operator++(int);
                in_column_iterator& operator--();
                in_column_iterator  operator--(int);
                in_column_iterator  operator+(unsigned int);
                in_column_iterator  operator-(unsigned int);
                in_column_iterator& operator+=(unsigned int);
                in_column_iterator& operator-=(unsigned int);
                bool                operator<(const in_column_iterator& other) const;
                bool                operator>(const in_column_iterator& other) const;
                bool                operator<=(const in_column_iterator& other) const;
                bool                operator>=(const in_column_iterator& other) const;
                ptrdiff_t           operator-(const in_column_iterator&) const;
                bool                operator==(const in_column_iterator&) const;
                bool                operator!=(const in_column_iterator&) const;
                
                friend class filled_tableau<T>;
            private:
                filled_tableau<T> *tab;
                unsigned int       column_number, row_number;
        };

        class iterator : public iterator_base {
            public:
                iterator(unsigned int r, unsigned int c, filled_tableau<T> *);
                T&                  operator*() const;
                T*                  operator->() const;
                iterator&           operator++();
                iterator            operator++(int);
                iterator&           operator--();
                iterator            operator--(int);
                iterator            operator+(unsigned int);
                iterator            operator-(unsigned int);
                iterator&           operator+=(unsigned int);
                iterator&           operator-=(unsigned int);
                bool                operator<(const iterator& other) const;
                bool                operator>(const iterator& other) const;
                ptrdiff_t           operator-(const iterator&) const;
                bool                operator==(const iterator&) const;
                bool                operator!=(const iterator&) const;
                
                friend class filled_tableau<T>;
            private:
                filled_tableau<T> *tab;
                unsigned int       column_number, row_number;
        };


        in_column_iterator   begin_column(unsigned int column_number);
        in_column_iterator   end_column(unsigned int column_number);
        iterator             begin() const; 
        iterator             end() const;
        
        template<class OutputIterator>
        OutputIterator       Garnir_set(OutputIterator, unsigned int, unsigned int) const;
    private:
        typedef std::vector<T>       box_row;
        typedef std::vector<box_row> row_stack;
        row_stack rows;
};

template<class T>
class tableaux {
    public:
        yngint_t         total_dimension(unsigned int dim);
        void             remove_nullifying_traces();
        bool             standard_form();
        void             add_tableau(const T&);
        void             symmetrise(const T& tabsym);

        typedef std::list<T> tableau_container_t;
        tableau_container_t  storage;

        typedef std::back_insert_iterator<tableau_container_t> back_insert_iterator;

        back_insert_iterator get_back_insert_iterator();
};

bool legal_box(const std::vector<std::pair<int,int> >& prev, 
                    const std::vector<std::pair<int,int> >& ths, 
                    int colpos, int trypos);

// --------------------------------------


template<class T>
typename tableaux<T>::back_insert_iterator tableaux<T>::get_back_insert_iterator()
    {
    return back_insert_iterator(storage);
    }

template<class T>
void tableaux<T>::remove_nullifying_traces()
    {
    typename tableau_container_t::iterator it=storage.begin();
    while(it!=storage.end()) {
        if(it->has_nullifying_trace()) 
            it=storage.erase(it);
        else ++it;
        }
    }

template<class T>
void tableaux<T>::symmetrise(const T& symtab)
    {
//
//  typename tableau_container_t::iterator thetab=storage.begin();
//  while(thetab!=storage.end()) {
//      (*thetab).sort_columns();
//      std::pair<int,int> where=(*thetab).nonstandard_loc();
//      if(where.first!=-1) {
//          combinations<typename T::value_type> com;
//  
    
/*
    FIXME: we should have two LR_tensor routines, because if you do 'alltabs', you should
    keep track of which boxes came from tableau 2. So do a LR_tensor with numbered boxes,
   and then after the LR_tensor apply the symmetries of the original tableaux, put back
   the original index names, sort columns and determine whether the tableau is identically
    non-zero. Then add to the product.

    Another issue: adding to tableaux should have an option to not insert doubles.

    There was something third, forgotten...
*/  
    }

template<class T>
void filled_tableau<T>::copy_shape(const tableau& other)
    {
    rows.clear();
    for(unsigned int r=0; r<other.number_of_rows(); ++r) {
        rows.push_back(box_row(other.row_size(r)));
        }
    tableau::operator=(other);
    }

template<class T>
bool filled_tableau<T>::compare_without_multiplicity(const filled_tableau<T>& other) const
    {
    return (rows==other.rows);
    }

template<class T>
bool filled_tableau<T>::has_nullifying_trace() const
    {
    return false;

// Old, probably incorrect code:
//
//  for(unsigned int r1=0; r1<number_of_rows(); ++r1) {
//      for(unsigned c1=0; c1<row_size(r1); ++c1) {
//          for(unsigned int c2=c1+1; c2<row_size(r1); ++c2) {
//              // (r1,c1) and (r1,c2)
//              for(unsigned int c3=0; c3<row_size(0); ++c3) {
//                  unsigned int r3=0;
//                  while(r3<number_of_rows()-1 && c3<row_size(r3)) {
//                      unsigned int r4=r3+1;
//                      while(r4<number_of_rows() && c3<row_size(r4)) {
//                          if((rows[r1][c1]==rows[r3][c3] && rows[r1][c2]==rows[r4][c3]) ||
//                              (rows[r1][c1]==rows[r4][c3] && rows[r1][c2]==rows[r3][c3]) )
//                              return true;
//                          ++r4;
//                          }
//                      ++r3;
//                      }
//                  }
//              }
//          }
//      }
//  return false;
    }

template<class T>
std::pair<int, int> filled_tableau<T>::find(const T& obj) const
    {
    for(unsigned int ir=0; ir<rows.size(); ++ir) {
        for(unsigned int ic=0; ic<rows[ir].size(); ++ic) {
            if(rows[ir][ic]==obj)
                return std::pair<int,int>(ir, ic);
            }
        }
    return std::pair<int,int>(-1,-1);
    }

template<class T>
void filled_tableau<T>::sort_within_columns()
    {
    std::less<T> comp;
    sort_within_columns(comp);
    }

template<class T>
void filled_tableau<T>::sort_columns()
    {
    std::less<T> comp;
    sort_columns(comp);
    }

template<class T>
void filled_tableau<T>::canonicalise()
    {
    std::less<T> comp;
    canonicalise(comp);
    }

template<class T>
template<class StrictWeakOrdering>
void filled_tableau<T>::sort_within_columns(StrictWeakOrdering comp)
    {
    filled_tableau<T> tmp(*this);
    if(number_of_rows()==0) return;
    for(unsigned int c=0; c<row_size(0); ++c) {
        std::sort(begin_column(c), end_column(c), comp);
        multiplicity*=combin::ordersign(begin_column(c), end_column(c), tmp.begin_column(c), tmp.end_column(c));
        }
    }

template<class T>
template<class StrictWeakOrdering>
void filled_tableau<T>::sort_columns(StrictWeakOrdering comp)
    {
    for(unsigned int c1=0; c1<row_size(0); ++c1) {
        for(unsigned int c2=c1; c2<row_size(0); ++c2) {
            if(column_size(c1)==column_size(c2)) {
                if(comp((*this)(0,c2), (*this)(0,c1)))
                    swap_columns(c1,c2);
                }
            }
        }
    }

template<class T>
template<class StrictWeakOrdering>
void filled_tableau<T>::canonicalise(StrictWeakOrdering comp, bool only_col_ex)
    {
    if(!only_col_ex)
        sort_within_columns(comp);
    sort_columns(comp);
    }

//---------------------------------------------------------------------------
// in_column_iterator

template<class T>
filled_tableau<T>::in_column_iterator::in_column_iterator(unsigned int r, unsigned int c, filled_tableau<T> *t)
    : tab(t), column_number(c), row_number(r)
    {
    }

template<class T>
typename filled_tableau<T>::in_column_iterator filled_tableau<T>::in_column_iterator::operator+(unsigned int n)
    {
    typename filled_tableau<T>::in_column_iterator it2(*this);
    it2+=n;
    return it2;
    }

template<class T>
typename filled_tableau<T>::in_column_iterator filled_tableau<T>::in_column_iterator::operator-(unsigned int n)
    {
    typename filled_tableau<T>::in_column_iterator it2(*this);
    it2-=n;
    return it2;
    }

template<class T>
ptrdiff_t filled_tableau<T>::in_column_iterator::operator-(const in_column_iterator& other) const
    {
    return row_number-other.row_number;
    }

template<class T>
T& filled_tableau<T>::in_column_iterator::operator*() const
    {
    return (*tab)(row_number,column_number);
    }

template<class T>
T* filled_tableau<T>::in_column_iterator::operator->() const
    {
    return &((*tab)(row_number,column_number));
    }

template<class T>
typename filled_tableau<T>::in_column_iterator& filled_tableau<T>::in_column_iterator::operator++()
    {
    ++row_number;
    return (*this);
    }

template<class T>
typename filled_tableau<T>::in_column_iterator& filled_tableau<T>::in_column_iterator::operator+=(unsigned int n)
    {
    row_number+=n;
    return (*this);
    }

template<class T>
typename filled_tableau<T>::in_column_iterator& filled_tableau<T>::in_column_iterator::operator--()
    {
    --row_number;
    return (*this);
    }

template<class T>
typename filled_tableau<T>::in_column_iterator filled_tableau<T>::in_column_iterator::operator--(int)
    {
    in_column_iterator tmp(*this);
    --row_number;
    return tmp;
    }

template<class T>
typename filled_tableau<T>::in_column_iterator filled_tableau<T>::in_column_iterator::operator++(int)
    {
    in_column_iterator tmp(*this);
    ++row_number;
    return tmp;
    }

template<class T>
typename filled_tableau<T>::in_column_iterator& filled_tableau<T>::in_column_iterator::operator-=(unsigned int n)
    {
    row_number-=n;
    return (*this);
    }

template<class T>
bool filled_tableau<T>::in_column_iterator::operator==(const in_column_iterator& other) const
    {
    if(tab==other.tab && row_number==other.row_number && column_number==other.column_number)
        return true;
    return false;
    }

template<class T>
bool filled_tableau<T>::in_column_iterator::operator<=(const in_column_iterator& other) const
    {
    if(row_number<=other.row_number) return true;
    return false;
    }

template<class T>
bool filled_tableau<T>::in_column_iterator::operator>=(const in_column_iterator& other) const
    {
    if(row_number>=other.row_number) return true;
    return false;
    }

template<class T>
bool filled_tableau<T>::in_column_iterator::operator<(const in_column_iterator& other) const
    {
    if(row_number<other.row_number) return true;
    return false;
    }

template<class T>
bool filled_tableau<T>::in_column_iterator::operator>(const in_column_iterator& other) const
    {
    if(row_number>other.row_number) return true;
    return false;
    }

template<class T>
bool filled_tableau<T>::in_column_iterator::operator!=(const in_column_iterator& other) const
    {
    return !((*this)==other);
    }

//---------------------------------------------------------------------------
// iterator

template<class T>
filled_tableau<T>::iterator::iterator(unsigned int r, unsigned int c, filled_tableau<T> *t)
    : tab(t), column_number(c), row_number(r)
    {
    }

template<class T>
typename filled_tableau<T>::iterator filled_tableau<T>::iterator::operator+(unsigned int n)
    {
    typename filled_tableau<T>::iterator it2(*this);
    it2+=n;
    return it2;
    }

template<class T>
typename filled_tableau<T>::iterator filled_tableau<T>::iterator::operator-(unsigned int n)
    {
    typename filled_tableau<T>::iterator it2(*this);
    it2-=n;
    return it2;
    }

template<class T>
ptrdiff_t filled_tableau<T>::iterator::operator-(const iterator& other) const
    {
    return row_number-other.row_number;
    }

template<class T>
T& filled_tableau<T>::iterator::operator*() const
    {
    return (*tab)(row_number,column_number);
    }

template<class T>
T* filled_tableau<T>::iterator::operator->() const
    {
    return &((*tab)(row_number,column_number));
    }

template<class T>
typename filled_tableau<T>::iterator& filled_tableau<T>::iterator::operator++()
    {
    if(++column_number==tab->rows[row_number].size()) {
        column_number=0;
        ++row_number;
        }
    return (*this);
    }

template<class T>
typename filled_tableau<T>::iterator& filled_tableau<T>::iterator::operator+=(unsigned int n)
    {
    while(n>0) {
        if(++column_number==tab->rows[row_number]) {
            column_number=0;
            ++row_number;
            }
        --n;
        }
    return (*this);
    }

template<class T>
typename filled_tableau<T>::iterator& filled_tableau<T>::iterator::operator--()
    {
    if(column_number==0) {
        --row_number;
        column_number=tab->rows[row_number].size()-1;
        }
    else --column_number;
    return (*this);
    }

template<class T>
typename filled_tableau<T>::iterator filled_tableau<T>::iterator::operator--(int)
    {
    iterator tmp(*this);
    if(column_number==0) {
        --row_number;
        column_number=tab->rows[row_number].size()-1;
        }
    else --column_number;
    return tmp;
    }

template<class T>
typename filled_tableau<T>::iterator filled_tableau<T>::iterator::operator++(int)
    {
    iterator tmp(*this);
    while(this->n>0) {
        if(++column_number==tab->rows[row_number]) {
            column_number=0;
            ++row_number;
            }
        --this->n;
        }
    return tmp;
    }

template<class T>
typename filled_tableau<T>::iterator& filled_tableau<T>::iterator::operator-=(unsigned int n)
    {
    while(n>0) {
        if(column_number==0) {
            --row_number;
            column_number=tab->rows[row_number].size()-1;
            }
        else --column_number;
        --n;
        }
    return (*this);
    }

template<class T>
bool filled_tableau<T>::iterator::operator==(const iterator& other) const
    {
    if(tab==other.tab && row_number==other.row_number && column_number==other.column_number)
        return true;
    return false;
    }

template<class T>
bool filled_tableau<T>::iterator::operator<(const iterator& other) const
    {
    if(row_number<other.row_number) return true;
    return false;
    }

template<class T>
bool filled_tableau<T>::iterator::operator>(const iterator& other) const
    {
    if(row_number>other.row_number) return true;
    return false;
    }

template<class T>
bool filled_tableau<T>::iterator::operator!=(const iterator& other) const
    {
    return !((*this)==other);
    }

//---
// other

template<class T>
typename filled_tableau<T>::iterator filled_tableau<T>::begin() const
    {
    return iterator(0,0,const_cast<filled_tableau<T> *>(this));
    }

template<class T>
typename filled_tableau<T>::iterator filled_tableau<T>::end() const
    {
    return iterator(rows.size(), 0, const_cast<filled_tableau<T> *>(this));
    }

template<class T>
typename filled_tableau<T>::in_column_iterator filled_tableau<T>::begin_column(unsigned int column)
    {
    typename filled_tableau<T>::in_column_iterator it(0,column,this);
    assert(number_of_rows()>0);
    assert(column<row_size(0));
    return it;
    }

template<class T>
typename filled_tableau<T>::in_column_iterator filled_tableau<T>::end_column(unsigned int column)
    {
    unsigned int r=0;
    while(r<number_of_rows()) {
        if(row_size(r)<=column)
            break;
        ++r;
        }
    typename filled_tableau<T>::in_column_iterator it(r,column,this);
    return it;
    }

template<class T>
template<class OutputIterator>
OutputIterator filled_tableau<T>::Garnir_set(OutputIterator it, unsigned int row, unsigned int col) const
    {
    assert(col>0);
    unsigned int r=row, c=col;
    *it=(*this)(r,c);
    ++it;
    while(r>0) {
        --r;
        *it=(*this)(r,c);
        ++it;
        }
    r=row;
    --c;
    *it=(*this)(r,c);
    ++it;
    while(r+1<column_size(c)) {
        ++r;
        *it=(*this)(r,c);
        ++it;
        }
    return it;
    }

template<class T>
std::pair<int, int> filled_tableau<T>::nonstandard_loc() const
    {
    unsigned int r=number_of_rows();
    assert(r>0);
    do {
        --r;
        for(unsigned int c=0; c<row_size(r)-1; ++c) {
            if((*this)(r,c) > (*this)(r,c+1) )
                return std::pair<int, int>(r,c);
            }
        } while(r>0);
    return std::pair<int,int>(-1,-1);
    }
    
template<class T>
bool tableaux<T>::standard_form()
    {
    bool already_standard=true;

    typename tableau_container_t::iterator thetab=storage.begin();
    while(thetab!=storage.end()) {
        (*thetab).sort_within_columns();
        std::pair<int,int> where=(*thetab).nonstandard_loc();
        if(where.first!=-1) {
            already_standard=false;
            combin::combinations<typename T::value_type> com;
            for(unsigned int i1=where.first; i1<(*thetab).column_size(where.second); ++i1)
                com.original.push_back((*thetab)(i1,where.second));
            for(unsigned int i1=0; i1<=(unsigned int)(where.first); ++i1)
                com.original.push_back((*thetab)(i1,where.second+1));
            com.sublengths.push_back((*thetab).column_size(where.second)-where.first);
            com.sublengths.push_back(where.first+1);
            com.permute();
            for(unsigned int tabi=1; tabi<com.size(); ++tabi) {
                T ntab((*thetab));
                unsigned int offset=0;
                for(unsigned int i1=where.first; i1<(*thetab).column_size(where.second); ++i1, ++offset)
                    ntab(i1,where.second)=com[tabi][offset];
                for(unsigned int i1=0; i1<=(unsigned int)(where.first); ++i1, ++offset)
                    ntab(i1,where.second+1)=com[tabi][offset];
                ntab.multiplicity*=-1*com.ordersign(tabi);
                add_tableau(ntab);
                }
            thetab=storage.erase(thetab);
            }
        else ++thetab;
        }
   return already_standard;
    }

template<class T>
void tableaux<T>::add_tableau(const T& ntab)
    {
    typename tableau_container_t::iterator it=storage.begin();
    while(it!=storage.end()) {
        if((*it).compare_without_multiplicity(ntab)) {
            (*it).multiplicity+=ntab.multiplicity;
            if((*it).multiplicity==0)
                storage.erase(it);
            return;
            }
        ++it;
        }
    storage.push_back(ntab);
    }


template<class T>
yngrat_t filled_tableau<T>::projector_normalisation() const
    {
    yngrat_t norm=1;
    norm/=hook_length_prod();
    return norm;
    }

template<class T>
void filled_tableau<T>::projector(combin::symmetriser<T>& sym, bool modulo_monoterm) const
    {
    for(unsigned int r=0; r<number_of_rows(); ++r) 
        for(unsigned int c=0; c<row_size(r); ++c)
            sym.original.push_back(rows[r][c]);

    unsigned int offset=0;
    // symmetrise over boxes in rows
    for(unsigned int r=0; r<number_of_rows(); ++r) {
        sym.permutation_sign=1;
        sym.permute_blocks.clear();
        sym.block_length=1;
        sym.input_asym.clear();
        for(unsigned int c=0; c<row_size(r); ++c) 
            sym.permute_blocks.push_back(offset++);
        sym.apply_symmetry();
        }
//  sym.collect();
    // anti-symmetrise over boxes in columns
    if(modulo_monoterm) {
        int newmult=1;
        for(unsigned int c=0; c<row_size(0); ++c) 
            newmult*=combin::factorial(column_size(c));
        for(unsigned int i=0; i<sym.size(); ++i) 
            sym.set_multiplicity(i, sym.signature(i)*newmult);
        }
    else {
        sym.permute_blocks.clear();
        for(unsigned int c=0; c<row_size(0); ++c) {
            unsigned int r=0;
            sym.value_permute.clear();
            sym.permutation_sign=-1;
            sym.input_asym.clear();
            while(r<number_of_rows() && c<row_size(r))
                sym.value_permute.push_back(rows[r++][c]);
            if(sym.value_permute.size()>1)
                sym.apply_symmetry();
            }
        }
//  sym.collect();
    }

template<class T>
void filled_tableau<T>::projector(combin::symmetriser<T>& sym, 
                                             combin::range_vector_t& sublengths_scattered) const
    {
    for(unsigned int r=0; r<number_of_rows(); ++r) 
        for(unsigned int c=0; c<row_size(r); ++c)
            sym.original.push_back(rows[r][c]);

    unsigned int offset=0;
    // symmetrise over boxes in rows
    for(unsigned int r=0; r<number_of_rows(); ++r) {
        sym.permutation_sign=1;
        sym.permute_blocks.clear();
        sym.block_length=1;
        sym.input_asym.clear();
        for(unsigned int c=0; c<row_size(r); ++c) 
            sym.permute_blocks.push_back(offset++);
        sym.apply_symmetry();
        }
    sym.permute_blocks.clear();
    for(unsigned int c=0; c<row_size(0); ++c) {
        unsigned int r=0;
        sym.value_permute.clear();
        sym.permutation_sign=-1;
        while(r<number_of_rows() && c<row_size(r)) 
            sym.value_permute.push_back(rows[r++][c]);

        sym.sublengths_scattered=sublengths_scattered;

//      // Now setup sublengths_scattered to take into account
//      // asym_ranges.  These asym_ranges refer to values stored in the
//      // boxes of the full tableau. We need to find the locations of
//      // these values inside the full original, as that is what goes
//      // into sublengths_scattered.
//
//      sym.input_asym.clear();
//      sym.sublengths.clear();
//      sym.sublengths_scattered.clear();
//      for(unsigned int m=0; m<sym.value_permute.size(); ++m) {
//          // Try to find this value in an asym range.
//          unsigned int overlap=0;
//          for(unsigned int n=0; n<asym_ranges.size(); ++n) {
//              for(unsigned int nn=0; nn<asym_ranges[n].size(); ++nn) {
//                  if(sym.value_permute[m]==asym_ranges[n][nn]) {
//                      std::cout << "found " << sym.value_permute[m] << " in range" << std::endl;
//                      // FIXME: this assumes that even though asym_ranges[n] is a superset
//                      // of the current part of value_permute, elements are in the same order.
//                      ++m; ++nn;
//                      while(nn<asym_ranges[n].size()) {
//                          if(sym.value_permute[m]==asym_ranges[n][nn]) {
//                              std::cout << "same range: " << sym.value_permute[m] << std::endl;
//                              ++m;
//                              ++overlap;
//                              }
//                          ++nn;
//                          }
//                      break;
//                      }
//                  }
//              }
//          if(overlap>0) --m;
//          sym.sublengths.push_back(overlap+1);
//          }
//      unsigned int sum=0;
//      for(unsigned int m=0; m<sym.sublengths.size(); ++m)
//          sum+=sym.sublengths[m];
//
//      std::cout << sum << " " << sym.value_permute.size() << std::endl;
//      assert(sum==sym.value_permute.size());
        
        // All set to run...
        if(sym.value_permute.size()>1)
            sym.apply_symmetry();
        }
    }

template<class T>
filled_tableau<T>& filled_tableau<T>::operator=(const filled_tableau<T>& other)
    {
    rows=other.rows;
    tableau::operator=(other);
    return (*this);
    }

template<class T>
yngint_t tableaux<T>::total_dimension(unsigned int dim)
    {
    yngint_t totdim=0;
    typename tableau_container_t::const_iterator it=storage.begin();
    while(it!=storage.end()) {
        totdim+=(*it).dimension(dim);
        ++it;
        }
    return totdim;
    }
    
template<class T, class OutputIterator>
void LR_tensor(const tableaux<T>& tabs1, const T& tab2, unsigned int maxrows,
                   OutputIterator out, bool alltabs=false)
    {
    typename tableaux<T>::tableau_container_t::const_iterator it=tabs1.storage.begin();
    while(it!=tabs1.storage.end()) {
        LR_tensor((*it), tab2, maxrows, out, alltabs);
        ++it;
        }
    }

template<class T1, class T2>
void add_box(T1& tab1, unsigned int row1,
                 const T2& tab2, unsigned int row2, unsigned int col2)
    {
    tab1.add_box(row1, tab2(row2,col2));
    }

template<class T1>
void add_box(T1& tab1, unsigned int row1, 
                 const tableau& tab2, unsigned int row2, unsigned int col2) 
    {
    tab1.add_box(row1);
    }

typedef filled_tableau<std::pair<int, int> > keeptrack_tab_t;

template<class Tab, class OutputIterator>
void LR_add_box(const Tab& tab2, Tab& newtab,
                     unsigned int currow2, unsigned int curcol2, unsigned int startrow,
                     unsigned int maxrows, 
                     OutputIterator outit, 
                     keeptrack_tab_t& Ycurrent,
                     bool alltabs)
    {
    // Are we at the end of the current row of boxes in tab2 ?
    if((++curcol2)==tab2.row_size(currow2)) {
        // Are we at the end of tab2 altogether?
        if((++currow2)==tab2.number_of_rows()) { 
            *outit=newtab;  // Store the product tableau just created.
            return;
            }
        curcol2=0;
        startrow=0;
        }

    // Rule "row_by_row".
    for(unsigned int rowpos=startrow; rowpos<std::min(newtab.number_of_rows()+1,maxrows); ++rowpos) {
        // Rule "always_young".
        if(rowpos>0 && rowpos<newtab.number_of_rows())
            if(newtab.row_size(rowpos-1)==newtab.row_size(rowpos))
                continue; // No, would lead to non-Young tableau shape.

        // The column where the box will be added.
        unsigned int colpos=(rowpos==newtab.number_of_rows())?0:newtab.row_size(rowpos);
        
        // Rule "avoid_sym2asym".
        for(unsigned int rr=0; rr<rowpos; ++rr)
            if(Ycurrent(rr,colpos).first==(int)(currow2))
                goto rule_violated;

        // Rule "avoid_asym2sym".
        if(alltabs) // if not generating all tabs, ordered will take care of this already.
            for(unsigned int cc=0; cc<colpos; ++cc)
                if(Ycurrent(rowpos,cc).second==(int)(curcol2))
                    goto rule_violated;
        
        // Rule "ordered".
        if(!alltabs && currow2>0) {
            int numi=0, numimin1=0;
            if(rowpos>0) {
                for(unsigned int sr=0; sr<rowpos; ++sr)  // top to bottom
                    for(unsigned int sc=0; sc<Ycurrent.row_size(sr); ++sc) { // right to left
                        // Count all boxes from currow2 and from currow2-1.
                        if(Ycurrent(sr,sc).first==(int)(currow2))    ++numi;
                        if(Ycurrent(sr,sc).first==(int)(currow2)-1)  ++numimin1;
                        }
                }
            ++numi; // the box to be added
            if(numi>numimin1)
                goto rule_violated;

            // continue counting to see whether a previously valid box is now invalid
            for(unsigned int sr=rowpos; sr<Ycurrent.number_of_rows(); ++sr)  // top to bottom
                for(int sc=Ycurrent.row_size(sr)-1; sc>=0; --sc) { // right to left
                    if(Ycurrent(sr,sc).first==(int)(currow2))    ++numi;
                    if(Ycurrent(sr,sc).first==(int)(currow2)-1)  ++numimin1;
                    if(numi>numimin1)
                        goto rule_violated;
                    }
            }

        // Put the box at row 'rowpos' and call LR_add_box recursively
        // to add the other boxes. 
        Ycurrent.add_box(rowpos, std::pair<int,int>(currow2, curcol2));
        add_box(newtab, rowpos, tab2, currow2, curcol2);
        LR_add_box(tab2, newtab, currow2, curcol2, alltabs?0:rowpos, maxrows,
                      outit, Ycurrent, alltabs);
        
        // Remove the box again in preparation for trying to add it to other rows.
        newtab.remove_box(rowpos);
        Ycurrent.remove_box(rowpos);

       rule_violated: ;
        }
    }
          
template<class Tab, class OutputIterator>
void LR_tensor(const Tab& tab1, const Tab& tab2, unsigned int maxrows,
                    OutputIterator outit, bool alltabs=false) 
    {
    // Make a copy of tab1 because LR_add_box has to change it and
    // tab1 is const here.
    Tab newtab(tab1); 

    // Tableau which keeps track of the LR rules. It contains the
    // current (incomplete) shape of the tensor product, and for all boxes 
    // which come from tab2 we store the row and column of tab2 
    // from which they originated. Tab1 boxes have (-2,-2) stored.
    keeptrack_tab_t Ycurrent;
    Ycurrent.copy_shape(tab1);
    keeptrack_tab_t::iterator yi=Ycurrent.begin();
    while(yi!=Ycurrent.end()) {
        (*yi)=std::pair<int,int>(-2,-2);
        ++yi;
        }

    LR_add_box(tab2, newtab, 0, -1, 0, maxrows, outit, Ycurrent, alltabs);
    }

template<class T, class OutputIterator>
void LR_tensor(const tableaux<T>&, bool symmetric, unsigned int maxrows, OutputIterator outit)
    {
    assert(1==0);
    }



std::ostream& operator<<(std::ostream&, const tableau& );
    
template<class T>
std::ostream& operator<<(std::ostream&, const tableaux<T>& );

template<class T>
std::ostream& operator<<(std::ostream&, const filled_tableau<T>& );

template<class T> 
unsigned int filled_tableau<T>::number_of_rows() const
    {
    return rows.size();
    }

template<class T> 
unsigned int filled_tableau<T>::row_size(unsigned int num) const
    {
    assert(num<rows.size());
    return rows[num].size();
    }

template<class T>
T& filled_tableau<T>::operator()(unsigned int row, unsigned int col) 
    {
    assert(row<rows.size());
    assert(col<rows[row].size());
    return rows[row][col];
    }

template<class T>
const T& filled_tableau<T>::operator()(unsigned int row, unsigned int col)  const
    {
    assert(row<rows.size());
    assert(col<rows[row].size());
    return rows[row][col];
    }

template<class T>
const T& filled_tableau<T>::operator[](unsigned int boxnum) const
    {
    assert(1==0);
    assert(this->row<rows.size());
    assert(this->col<rows[this->row].size());
    return rows[this->row][this->col];
    }

template<class T>
filled_tableau<T>::~filled_tableau()
    {
    }

template<class T>
void filled_tableau<T>::add_box(unsigned int rownum) 
    {
    if(rownum>=rows.size())
        rows.resize(rownum+1);
    assert(rownum<rows.size());
    rows[rownum].push_back(T());
    }

template<class T>
void filled_tableau<T>::add_box(unsigned int rownum, T val) 
    {
    if(rownum>=rows.size())
        rows.resize(rownum+1);
    assert(rownum<rows.size());
    rows[rownum].push_back(val);
    }

template<class T>
void filled_tableau<T>::swap_columns(unsigned int c1, unsigned int c2)
    {
    assert(c1<row_size(0) && c2<row_size(0));
    assert(column_size(c1)==column_size(c2));
    for(unsigned int r=0; r<column_size(c1); ++r) {
        T tmp=(*this)(r,c1);
        (*this)(r,c1)=(*this)(r,c2);
        (*this)(r,c2)=tmp;
        }
    }

template<class T>
void filled_tableau<T>::remove_box(unsigned int rownum)
    {
    assert(rownum<rows.size());
    assert(rows[rownum].size()>0);
    rows[rownum].pop_back();
    if(rows[rownum].size()==0)
        rows.pop_back();
    }

template<class T>
void filled_tableau<T>::clear()
    {
    rows.clear();
    tableau::clear();
    }

template<class T>
std::ostream& operator<<(std::ostream& str, const tableaux<T>& tabs)
    {
    typename tableaux<T>::tableau_container_t::const_iterator it=tabs.storage.begin();
    while(it!=tabs.storage.end()) {
        str << (*it) << std::endl << std::endl;
        ++it;
        }
    return str;
    }

template<class T>
std::ostream& operator<<(std::ostream& str, const filled_tableau<T>& tab)
    {
    for(unsigned int i=0; i<tab.number_of_rows(); ++i) {
        for(unsigned int j=0; j<tab.row_size(i); ++j) {
//          str << "|" << tab(i,j) << "|";
            str << tab(i,j);
            }
        if(i==0) {
            str << "  " << tab.dimension(10);
            if(tab.has_nullifying_trace()) str << " null";
            }
        if(i!=tab.number_of_rows()-1) 
            str << std::endl;
        else
            str << " (" << tab.multiplicity << ")" << std::endl;
        }
    return str;
    }

};