//#####################################################################
// Copyright 2017, Ounan Ding, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Grids/FACE_RANGE_ITERATOR.h>
namespace PhysBAM{
// These tables are generated by this program:
// Projects/grid_iterators/gen_face_range_iterator_table.cpp
static const long default_table[6]=
{
    0766000676000667000L
};
static const long ghost_table[6]=
{
    0762000672000663000L,
    0766004676004667004L,
    0724002634002624003L,
    0764042674042664043L,
    0344022244032244023L,
    0744422644432644423L
};
static const long delay_corners_table[6]=
{
    0442320442230443220L,
    0446324446234447224L,
    0426300436200427200L,
    0466340476240467240L,
    0366000276000267000L,
    0766400676400667400L
};
static const long duplicate_corners_table[6]=
{
    0762000672000663000L,
    0766004676004667004L,
    0726000636000627000L,
    0766040676040667040L,
    0366000276000267000L,
    0766400676400667400L
};
static const long omit_corners_table[6]=
{
    0542220452220443220L,
    0546224456224447224L,
    0524202434202425202L,
    0564242474242465242L,
    0344022254022245022L,
    0744422654422645422L
};
static const long inner_outer_mask=0100100010010001001L;
static const long clamp_mask=0000001000001000001L;
//#####################################################################
// Constructor
//#####################################################################
template<int d> FACE_RANGE_ITERATOR<d>::
FACE_RANGE_ITERATOR(const RANGE<IV>& outer,const RANGE<IV>& inner,
    RF flags,int side_input,int axis)
    :vecs{outer.min_corner,inner.min_corner,inner.max_corner,outer.max_corner}
{
    Initialize(flags,side_input,axis);
}
//#####################################################################
// Constructor
//#####################################################################
template<int d> FACE_RANGE_ITERATOR<d>::
FACE_RANGE_ITERATOR(const RANGE<IV>& range,int outer_ghost,int inner_ghost,
    RF flags,int side_input,int axis)
    :vecs{range.min_corner-outer_ghost,range.min_corner-inner_ghost,
        range.max_corner+inner_ghost,range.max_corner+outer_ghost}
{
    Initialize(flags,side_input,axis);
}
//#####################################################################
// Constructor
//#####################################################################
template<int d> FACE_RANGE_ITERATOR<d>::
FACE_RANGE_ITERATOR(const IV& counts,int outer_ghost,int inner_ghost,
    RF flags,int side_input,int axis)
    :vecs{IV()-outer_ghost,IV()-inner_ghost,
        counts+inner_ghost,counts+outer_ghost}
{
    Initialize(flags,side_input,axis);
}
//#####################################################################
// Constructor
//#####################################################################
template<int d> FACE_RANGE_ITERATOR<d>::
FACE_RANGE_ITERATOR(const RANGE<IV>& range,RF flags,int axis)
    :vecs{range.min_corner,range.min_corner,range.max_corner,range.max_corner}
{
    Initialize(flags,-1,axis);
}
//#####################################################################
// Constructor
//#####################################################################
template<int d> FACE_RANGE_ITERATOR<d>::
FACE_RANGE_ITERATOR(const IV& counts,int outer_ghost,RF flags,int axis)
    :vecs{IV()-outer_ghost,IV()-outer_ghost,
        counts+outer_ghost,counts+outer_ghost}
{
    Initialize(flags,-1,axis);
}
//#####################################################################
// Function Set_Range
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Set_Range(const RANGE<IV>& outer,const RANGE<IV>& inner)
{
    vecs[0]=outer.min_corner;
    vecs[1]=inner.min_corner;
    vecs[2]=inner.max_corner;
    vecs[3]=outer.max_corner;
}
//#####################################################################
// Function Set_Range
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Set_Range(const RANGE<IV>& range,int outer_ghost,int inner_ghost)
{
    vecs[0]=range.min_corner-outer_ghost;
    vecs[1]=range.min_corner-inner_ghost;
    vecs[2]=range.max_corner+inner_ghost;
    vecs[3]=range.max_corner+outer_ghost;
}
//#####################################################################
// Function Set_Range
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Set_Range(const IV& counts,int outer_ghost,int inner_ghost)
{
    vecs[0]=IV()-outer_ghost;
    vecs[1]=IV()-inner_ghost;
    vecs[2]=counts+inner_ghost;
    vecs[3]=counts+outer_ghost;
}
//#####################################################################
// Function Initialize
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Initialize(RF flags,int side_input,int axis)
{
    Encode(flags,side_input,axis);
    Reset(flags);
}
//#####################################################################
// Function Reset
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Reset(RF flags)
{
    if(!side_mask || !axis_mask) // No faces to visit; set to not valid.
    {
        if(!(flags&RF::reverse))
        {
            current[1](d-1)=0;
            face.index(d-1)=0;
        }
        else
        {
            current[0](d-1)=0;
            face.index(d-1)=-1;
        }
        return;
    }
    if(!(flags&RF::reverse)){
        side=2*d;
        face.axis=-1;
        Next_Axis_Side();}
    else{
        side=-1;
        face.axis=d;
        Prev_Axis_Side();}
}

// Search for octal character x at offsets given by mask, returning a mask
// of those offsets.  Note that (mask&result)==result.
static long find(long in,long mask,int x)
{
    long y=~(in^(mask*x)); // if octal character matches, put 7 at that spot.
    return y&(y>>1)&(y>>2)&mask; // Test for octal 7.
}

// Search for x&~1 and x|1 together.
static long find_pair(long in,long mask,int x)
{
    long y=~(in^(mask*x)); // if octal character matches, put 7 at that spot.
    return (y>>1)&(y>>2)&mask; // Test for octal 7.
}

static void do_inner(long* m,int num,RF flags)
{
    // Change 3->2, 4->5 at inner_outer_mask locations.
    if(any(flags&RF::skip_inner))
        for(int i=0;i<num;i++)
            m[i]^=find(m[i],inner_outer_mask,3)|find(m[i],inner_outer_mask,4);
}

static void do_outer(long* m,int num,RF flags)
{
    // Change 0->1, 7->6 at inner_outer_mask locations.
    if(any(flags&RF::skip_outer))
        for(int i=0;i<num;i++)
            m[i]^=find(m[i],inner_outer_mask,0)|find(m[i],inner_outer_mask,7);
}

//#####################################################################
// Function Encode
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Encode(RF flags,int side_input,int axis)
{
    if(any(flags&RF::axis_mask)) axis_mask=axis&((1<<d)-1);
    else axis_mask=axis<0?(1<<d)-1:(1<<axis);
    if(any(flags&RF::side_mask)) side_mask=side_input&((1<<2*d)-1);
    else side_mask=side_input<0?(1<<2*d)-1:(1<<side_input);
    flags&=~(RF::reverse|RF::axis_mask|RF::side_mask);
    if(!(flags&RF::ghost)) side_mask=1;

    // Explicitly check for flag combinations that do not make sense.
    if(!(flags&RF::ghost)){
        assert(!(flags&~RF::skip_outer));
        assert(vecs[0]==vecs[1] && vecs[2]==vecs[3]);}
    if(any(flags&RF::duplicate_corners)){
        assert(!(flags&~(RF::duplicate_corners|RF::ghost|RF::skip_outer|RF::skip_inner)));}
    if(any(flags&RF::omit_corners)){
        assert(!(flags&~(RF::omit_corners|RF::ghost|RF::skip_outer|RF::skip_inner)));}
    if(any(flags&RF::delay_corners)){
        assert(!(flags&~(RF::delay_corners|RF::ghost|RF::skip_outer|RF::skip_inner)));}

    const long* table=0;
    if(!(flags&RF::ghost)) table=default_table;
    else if(any(flags&RF::delay_corners)) table=delay_corners_table;
    else if(any(flags&RF::duplicate_corners)) table=duplicate_corners_table;
    else if(any(flags&RF::omit_corners)) table=omit_corners_table;
    else table=ghost_table;
    for(int i=0;i<2*d;i++) inst[i]=table[i];

    do_outer(inst,2*d,flags);
    if(!(flags&RF::ghost)) return;
    do_inner(inst,2*d,flags);

    bool need_repeat=false;
    for(int a=0;a<d;a++)
    {
        if(vecs[0](a)>vecs[1](a))
        {
            need_repeat=true;
            for(auto&m:inst)
                m&=~(7*find_pair(m,clamp_mask<<3*a,2)); // 2,3 -> 0
        }
        else if(vecs[0](a)>=vecs[1](a))
        {
            need_repeat=true;
            for(auto&m:inst)
                m-=2*find_pair(m,clamp_mask<<3*a,2); // 2 -> 0, 3 -> 1
        }
        if(vecs[2](a)>vecs[3](a))
        {
            need_repeat=true;
            for(auto&m:inst)
            {
                long n=find_pair(m,clamp_mask<<3*(a+3),4);
                m|=7*n; // 4,5 -> 7
                m-=n&~inner_outer_mask;
            }
        }
        else if(vecs[2](a)>=vecs[3](a))
        {
            need_repeat=true;
            for(auto&m:inst)
                m+=2*find_pair(m,clamp_mask<<3*(a+3),4); // 4 -> 6, 5 -> 7
        }
    }

    if(need_repeat)
        do_outer(inst,2*d,flags);
}
//#####################################################################
// Function Next_Helper
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Next_Helper()
{
    // next row of range
    face.index(d-1)=current[0](d-1);
    for(int i=d-2;i>=0;i--){
        if(++face.index(i)<current[1](i)) return;
        face.index(i)=current[0](i);}
    Next_Axis_Side();
}
//#####################################################################
// Function Prev_Helper
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Prev_Helper()
{
    // next row of range
    face.index(d-1)=current[1](d-1)-1;
    for(int i=d-2;i>=0;i--){
        if(--face.index(i)>=current[0](i)) return;
        face.index(i)=current[1](i)-1;}
    Prev_Axis_Side();
}
//#####################################################################
// Function Next_Axis_Side
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Next_Axis_Side()
{
    // finished current block; next block
    side++;
    if(!(side_mask>>side)){
        side=0;
        face.axis++;
        if(!(axis_mask>>face.axis)){
            face.index(d-1)=current[1](d-1); // Make invalid
            return;} // Done
        while(!(axis_mask&(1<<face.axis))) face.axis++;}
    while(!(side_mask&(1<<side))) side++;
    Fill_Current();
    face.index=current[0];
    if(!current[0].All_Less(current[1])) Next_Axis_Side();
}
//#####################################################################
// Function Prev_Axis_Side
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Prev_Axis_Side()
{
    // finished current block; next block
    side--;
    if(side<0 || !(side_mask&((2<<side)-1))){
        side=2*d-1;
        face.axis--;
        if(face.axis<0 || !(axis_mask&((2<<face.axis)-1))){
            face.index(d-1)=current[0](d-1)-1; // Make invalid
            return;} // Done
        while(!(axis_mask&(1<<face.axis))) face.axis--;}
    while(!(side_mask&(1<<side))) side--;
    Fill_Current();
    face.index=current[1]-1;
    if(!current[0].All_Less(current[1])) Prev_Axis_Side();
}
//#####################################################################
// Function Fill_Current
//#####################################################################
template<int d> void FACE_RANGE_ITERATOR<d>::
Fill_Current()
{
    // Assemble next block from rules.
    long mask=inst[side]>>(18*face.axis);
    for(int k=0;k<2;k++)
        for(int a=0;a<d;a++)
        {
            int m=(mask>>(3*(3*k+a)))&7;
            current[k](a)=vecs[m>>1](a)+(m&1);
        }
}
//#####################################################################
template class FACE_RANGE_ITERATOR<1>;
template class FACE_RANGE_ITERATOR<2>;
template class FACE_RANGE_ITERATOR<3>;
}
