//##################################################################### 
// Copyright 2021, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <array>

namespace PhysBAM {

template<size_t m,size_t n,size_t s> consteval auto
up(std::array<std::array<double,m>,n> E,std::array<double,s> S)
{
    std::array<std::array<double,m>,n+s*2-2> A={};
    for(int i=0;i<n;i++) A[i+s-1]=E[i];

    std::array<std::array<double,m>,n+s-1> B = {};
    for(int i=0;i<n+s-1;i++)
        for(int j=0;j<s;j++)
            for(int k=0;k<m;k++)
                B[i][k]+=A[i+j][k]*S[j];
    
    std::array<std::array<double,m+1>,n+s-1> R={};
    double acc=0;
    for(int j=0;j<n+s-1;j++)
    {
        for(int i=0;i<m;i++)
        {
            R[j][i+1]=-B[j][i]/(i+1);
            acc+=B[j][i]/(i+1);
        }
        R[j][0]=acc;
    }

    return R;
}

template<size_t m,size_t n> consteval auto
pad(std::array<std::array<double,m>,n> A)
{
    std::array<std::array<double,m>,n+2> R={};
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++)
            R[i+1][j]=A[i][j];
    return R;
}

constexpr std::array<double,2> ST2={-1,1};
constexpr std::array<double,4> ST4={1/24.,-27/24.,27/24.,-1/24.};
constexpr std::array<double,6> ST6={-3/640.,25/384.,-75/64.,75/64.,-25/384.,3/640.};
constexpr std::array<double,8> ST8={5/7168.,-49/5120.,245/3072.,-1225/1024.,1225/1024.,-245/3072.,49/5120.,-5/7168.};

// N:=6;f:=`+`(seq(a||i*x^i,i=0..N-1));EV:=`+`(seq(f*c||x,x=0..N-1));DF:=subs(x=(N-1)/2,diff(f,x));S:=solve([coeffs(expand(DF-EV),indets(f))]);lprint(S);

typedef std::array<double,1> A0;
typedef std::array<double,3> A2;
typedef std::array<double,4> A3;
typedef std::array<double,5> A4;

#define NEW_CHAIN(name,deg,width,init)                                  \
    constexpr std::array<std::array<double,deg+1>,width> name##_ init;  \
    extern const auto name=name##_
#define INC_CHAIN(old_name,new_name,stencil)            \
    constexpr auto new_name##_=up(old_name##_,stencil); \
    extern const auto new_name=new_name##_;
#define PAD_CHAIN(old_name,new_name)                    \
    constexpr auto new_name##_=pad(old_name##_);        \
    extern const auto new_name=new_name##_;

namespace IS {
NEW_CHAIN(B0,0,1,{A0{1}});
INC_CHAIN(B0,B1,ST2);
INC_CHAIN(B1,B2,ST2);
INC_CHAIN(B2,B3,ST2);

NEW_CHAIN(C2,2,2,({A2{1,-4,3},A2{0,-2,3}}));
INC_CHAIN(C2,C3,ST2);
INC_CHAIN(C3,C4,ST2);

NEW_CHAIN(D2,2,1,({A2{1,-6,6}}));
INC_CHAIN(D2,D3,ST2);
INC_CHAIN(D3,D4,ST2);
INC_CHAIN(D4,D5,ST2);

NEW_CHAIN(F4,4,3,({A4{-3,2,30,-54,25},A4{-3,-2,132,-260,130},A4{0,0,18,-46,25}}));
INC_CHAIN(F4,F5,ST2);

NEW_CHAIN(G4,4,3,({A4{27.5,-55,-172.5,400,-200},A4{27.5,55,-1095,2080,-1040},A4{0,0,-172.5,400,-200}}));
INC_CHAIN(G4,G5,ST2);

NEW_CHAIN(H3,3,3,({A3{18,-50,46,-14},A3{18,50,-50},A3{0,0,4,14}}));
INC_CHAIN(H3,H4,ST2);

NEW_CHAIN(I3,3,4,({A3{0,-2,3,-1},A3{3,0,-6,3},A3{0,3,3,-3},A3{0,-1,0,1}}));
NEW_CHAIN(J3,3,4,({A3{0,-2,3,-1},A3{2,1,-6,3},A3{0,2,3,-3},A3{0,-1,0,1}}));

PAD_CHAIN(B0,B0p);
PAD_CHAIN(B1,B1p);
PAD_CHAIN(B2,B2p);
PAD_CHAIN(B3,B3p);
PAD_CHAIN(B0p,B0pp);
PAD_CHAIN(B1p,B1pp);
PAD_CHAIN(B2p,B2pp);
PAD_CHAIN(B3p,B3pp);
PAD_CHAIN(B1pp,B1ppp);

PAD_CHAIN(C2,C2p);
PAD_CHAIN(C3,C3p);
PAD_CHAIN(C4,C4p);
PAD_CHAIN(C2p,C2pp);
PAD_CHAIN(C3p,C3pp);
PAD_CHAIN(C4p,C4pp);

PAD_CHAIN(I3,I3p);
PAD_CHAIN(I3p,I3pp);
PAD_CHAIN(J3,J3p);
PAD_CHAIN(J3p,J3pp);
}

namespace IS4 {
INC_CHAIN(::PhysBAM::IS::B0,B1,ST4);
INC_CHAIN(::PhysBAM::IS::B1,B2,ST4);
INC_CHAIN(::PhysBAM::IS::B2,B3,ST4);
PAD_CHAIN(B1,B1p);
PAD_CHAIN(B2,B2p);
PAD_CHAIN(B3,B3p);

INC_CHAIN(::PhysBAM::IS::C2,C3,ST4);
INC_CHAIN(::PhysBAM::IS::C3,C4,ST4);
PAD_CHAIN(C3,C3p);
PAD_CHAIN(C4,C4p);

INC_CHAIN(::PhysBAM::IS::I3,I4,ST4);
INC_CHAIN(::PhysBAM::IS::J3,J4,ST4);
}

namespace IS6 {
INC_CHAIN(::PhysBAM::IS::B0,B1,ST6);
INC_CHAIN(::PhysBAM::IS::B1,B2,ST6);
INC_CHAIN(::PhysBAM::IS::B2,B3,ST6);
}

namespace IS8 {
INC_CHAIN(::PhysBAM::IS::B0,B1,ST8);
INC_CHAIN(::PhysBAM::IS::B1,B2,ST8);
INC_CHAIN(::PhysBAM::IS::B2,B3,ST8);
}

}
