#include <Core/Symbolics/PROGRAM_COMPILER.h>
#include <Core/Arrays/INDIRECT_ARRAY.h>
#include <Core/Data_Structures/HASHTABLE.h>
#include <Core/Data_Structures/UNION_FIND.h>
#include <Core/Log/LOG.h>
#include <Core/Math_Tools/exchange_sort.h>
#include <Core/Symbolics/CODE_BLOCK.h>
#include <Core/Symbolics/COMPILER_DATA_OPS.h>
#include <Core/Symbolics/DOMINATORS.h>
#include <Core/Symbolics/INSTRUCTION.h>
#include <Core/Symbolics/PROGRAM.h>
#include <Core/Symbolics/PROGRAM_PARSE_NODE.h>
#include <Core/Utilities/PROCESS_UTILITIES.h>
#include <algorithm>
#include <cstdlib>
#include <cstring>
#include <set>
#ifdef USE_LEX_YACC
#include <Core/Symbolics/PROGRAM_LEX.ll.hpp>
#include <Core/Symbolics/PROGRAM_YACC.yy.hpp>
#endif

extern int yyparse();

namespace PhysBAM{

template<class T>
struct DISTRIBUTE_HELPER
{
    T coeff;
    ARRAY<REG> product;
};

template<class K,class V>
std::ostream& operator<<(std::ostream& o,const std::map<K,V>& m)
{
    o<<'(';
    for(auto&p:m) o<<p;
    return o<<')';
}

ARRAY<std::string> parse_identifiers;
ARRAY<double> parse_constants;
PROGRAM_PARSE_NODE* parse_root;

#define DEBUG if(debug) Test_Invariants();
#define VERBOSE if(verbose) {LOG::printf("\n%s\n",__FUNCTION__);Print();}

template<class MP,class FA,class FB,class F>
void Zip(const MP& a,const MP& b,FA fa,FB fb,F f)
{
    auto i=a.begin();
    auto j=b.begin();
    while(i!=a.end() && j!=b.end())
    {
        if(i->first<j->first) fa(*i++);
        else if(j->first<i->first) fb(*j++);
        else f(*i++,*j++);
    }
    while(i!=a.end()) fa(*i++);
    while(j!=b.end()) fb(*j++);
}

//#####################################################################
// Function For_Each_Instruction
//#####################################################################
template<class T,class F> void
For_Each_Instruction(const ARRAY<T*,MEM_LABEL>& code_blocks,F f)
{
    for(MEM_LABEL i{};i<code_blocks.m;i++)
        if(auto B=code_blocks(i))
            for(auto& o:*B)
                if(!f(o))
                    return;
}
//#####################################################################
// Function For_Each_Iterator
//#####################################################################
template<class T,class F> void
For_Each_Iterator(const ARRAY<T*,MEM_LABEL>& code_blocks,F f)
{
    typedef CODE_BLOCK::ITERATOR INST_IT;
    for(MEM_LABEL i{};i<code_blocks.m;i++)
        if(auto B=code_blocks(i))
            for(auto it=B->insts.begin();it!=B->insts.end();it++)
                if(!f(INST_IT{B,it}))
                    return;
}
//#####################################################################
// Function For_Each_Instruction_Reverse
//#####################################################################
template<class T,class F> void
For_Each_Instruction_Reverse(const ARRAY<T*,MEM_LABEL>& code_blocks,F f)
{
    for(MEM_LABEL i=code_blocks.m-1;i>=MEM_LABEL{};i--)
        if(auto B=code_blocks(i))
            for(auto it=B->insts.rbegin();it!=B->insts.rend();it++)
                if(!f(*it))
                    return;
}
//#####################################################################
// Function For_Each_Iterator_Reverse
//#####################################################################
template<class T,class F> void
For_Each_Iterator_Reverse(const ARRAY<T*,MEM_LABEL>& code_blocks,F f)
{
    typedef CODE_BLOCK::ITERATOR INST_IT;
    for(MEM_LABEL i=code_blocks.m-1;i>=MEM_LABEL{};i--)
        if(auto B=code_blocks(i))
            for(auto it=B->insts.end();it!=B->insts.begin();)
                if(!f(INST_IT{B,--it}))
                    return;
}
//#####################################################################
// Constructor
//#####################################################################
template<class T> PROGRAM_COMPILER<T>::
PROGRAM_COMPILER()
{
    Add_Constant(0);
    Add_Constant(1);
}
//#####################################################################
// Constructor
//#####################################################################
template<class T> PROGRAM_COMPILER<T>::
PROGRAM_COMPILER(const PROGRAM_COMPILER& pc)
    :num_in(pc.num_in),num_out(pc.num_out),num_reg(pc.num_reg),constants(pc.constants),
    constant_lookup(pc.constant_lookup),
    var_in(pc.var_in),var_out(pc.var_out),debug(pc.debug)
{
    if(debug) pc.Test_Invariants();
    for(auto& i:pc.dict)
    {
        assert(!i.second->ref);
        dict.Set(i.first,new DATA(*i.second));
    }
    for(auto B:pc.code_blocks)
    {
        if(!B)
        {
            code_blocks.Append(0);
            continue;
        }
        auto C=new CODE_BLOCK(*B);
        C->id=code_blocks.Append(C);
    }
    for(MEM_LABEL i{0};i<pc.code_blocks.m;i++)
    {
        auto A=pc.code_blocks(i);
        if(!A) continue;
        auto B=code_blocks(i);
        for(int j=0;j<2;j++)
        {
            if(A->next[j]) B->next[j]=code_blocks(A->next[j]->id);
            if(A->prev[j]) B->prev[j]=code_blocks(A->prev[j]->id);
        }
    }
    Compute_Use_Def();
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Destructor
//#####################################################################
template<class T> PROGRAM_COMPILER<T>::
~PROGRAM_COMPILER()
{
    for(auto&i:dict) delete i.second;
}
//#####################################################################
// Function Deduce_Op
//#####################################################################
template<class T> bool PROGRAM_COMPILER<T>::
Deduce_Op(OPERATION op,T& out,T* in0,T* in1,T* in2) const
{
    switch(op.op)
    {
        case op_mul: if((in0 && !*in0) || (in1 && !*in1)){out=0;return true;} return false;
        case op_div: if((in0 && !*in0)){out=0;return true;} return false;
        case op_pow: if(in0 && (*in0==0 || *in0==1)){out=*in0;return true;}
            if(in1 && *in1==0){out=1;return true;} return false;
        case op_or: if((in0 && *in0) || (in1 && *in1)){out=1;return true;} return false;
        case op_and: if((in0 && !*in0) || (in1 && !*in1)){out=0;return true;} return false;
        default:break;
    }
    return false;
}
//#####################################################################
// Function Before_Instruction
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Before_Instruction(const INST_IT& it,const OPERATION& type,REG src0,REG src1,REG src2)
{
    assert(type.Has_Dest());
    REG dest=New_Reg();
    auto ir=it->first->insts.insert(it->second,{type,dest,src0,src1,src2});
    Add_Use_Def({it->first,ir});
    return dest;
}
//#####################################################################
// Function Append_Instruction_Dest
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Append_Instruction_Dest(CODE_BLOCK* B,const OPERATION& type,REG dest,REG src0,REG src1,REG src2)
{
    B->insts.push_back({type,dest,src0,src1,src2});
    Add_Use_Def({B,std::prev(B->insts.end())});
    return dest;
}
//#####################################################################
// Function Append_Instruction_Dest
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Append_Instruction_Dest(const OPERATION& type,REG dest,REG src0,REG src1,REG src2)
{
    return Append_Instruction_Dest(code_blocks.Last(),type,dest,src0,src1,src2);
}
//#####################################################################
// Function Make_Adjoint
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Make_Adjoint(REG func_var) -> ADJOINT
{
    PHYSBAM_ASSERT(!func_var.Is_Const() && !func_var.Is_In());
    auto it=Get_Def(func_var);
    const INSTRUCTION& o=*it->second;
    PHYSBAM_ASSERT(o.type.Has_Dest());

    if(auto ptr=adjoints.Get_Pointer(func_var))
        return *ptr;

    ADJOINT a;
    REG zero=MEM_CONST(0);
    REG one=MEM_CONST(1);
    REG d,e,f,s,t;

    auto parent=func_var;
    if(auto p=diff_parent.Get_Pointer(func_var))
        parent=*p;
    if(auto p=diff_block_end.Get_Pointer(parent))
        it=*p;
    auto diff_last_it=it++;

    switch(o.type.op){
        case op_copy:a.val[0]={one};break;
        case op_neg:a.val[0]={one,1};break;
        case op_add:a.val[0]={one};a.val[1]={one};break;
        case op_sub:a.val[0]={one};a.val[1]={one,1};break;
        case op_phi:PHYSBAM_FATAL_ERROR("Cannot call adjoint for phi.");break;
        case op_mul:a.val[0]={o.src1};a.val[1]={o.src0};break;
        case op_div:a.val[0]={o.src1,0,1};
            a.val[1]={Before_Instruction(it,op_div,o.dest,o.src1),1};
            break;
        case op_inv:a.val[0]={Before_Instruction(it,op_mul,o.dest,o.dest),1};break;
        case op_sqrt:
            a.val[0]={Before_Instruction(it,op_div,Add_Constant((T).5),o.dest)};
            break;
        case op_exp:a.val[0]={o.dest};break;
        case op_ln:a.val[0]={o.src0,0,1};break;
        case op_pow:
            t=Before_Instruction(it,op_ln,o.src0);
            t=Before_Instruction(it,op_mul,o.dest,t);
            a.val[1]={t};
            t=Before_Instruction(it,op_sub,o.src1,one);
            t=Before_Instruction(it,op_pow,o.src0,t);
            t=Before_Instruction(it,op_mul,t,o.src1);
            a.val[0]={t};
            break;
        case op_sin:a.val[0]={Before_Instruction(it,op_cos,o.src0)};break;
        case op_cos:a.val[0]={Before_Instruction(it,op_sin,o.src0),1};break;
        case op_asin:
            t=Before_Instruction(it,op_mul,o.src0,o.src0);
            t=Before_Instruction(it,op_sub,one,t);
            t=Before_Instruction(it,op_sqrt,t);
            a.val[0]={t,0,1};
            break;
        case op_acos:
            t=Before_Instruction(it,op_mul,o.src0,o.src0);
            t=Before_Instruction(it,op_sub,one,t);
            t=Before_Instruction(it,op_sqrt,t);
            a.val[0]={t,1,1};
            break;
        case op_atan:
            t=Before_Instruction(it,op_mul,o.src0,o.src0);
            t=Before_Instruction(it,op_add,one,t);
            a.val[0]={t,0,1};
            break;
        case op_atan2:
            s=Before_Instruction(it,op_mul,o.src0,o.src0);
            t=Before_Instruction(it,op_mul,o.src1,o.src1);
            t=Before_Instruction(it,op_add,s,t);
            a.val[0]={Before_Instruction(it,op_div,o.src1,t)};
            a.val[1]={Before_Instruction(it,op_div,o.src0,t),1};
            break;
        case op_lt:case op_le:case op_gt:case op_ge:case op_eq:case op_ne:
        case op_not:case op_or:case op_and:
        case op_floor:case op_ceil:case op_rint:
            a.val[0]={zero};a.val[1]={zero}; break;
        case op_if:PHYSBAM_FATAL_ERROR("Cannot call adjoint for If.");break;
        case op_min:PHYSBAM_FATAL_ERROR("Cannot call adjoint for Min.");break;
        case op_max:PHYSBAM_FATAL_ERROR("Cannot call adjoint for Max.");break;
        case op_abs:PHYSBAM_FATAL_ERROR("Cannot call adjoint for Abs.");break;
        default:PHYSBAM_FATAL_ERROR("Missing diff instruction");
    }
    while(++diff_last_it!=it)
        diff_parent.Set(diff_last_it->second->dest,parent);
    diff_block_end.Set(parent,--it);
    adjoints.Set(func_var,a);
    return a;
}
//#####################################################################
// Function Diff_Internal
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Diff_Internal(REG func_var,REG diff_var) -> REG_AUG
{
    REG zero=MEM_CONST(0);
    if(func_var.Is_Const()) return {zero};
    if(func_var.Is_Dont_Care()) return {func_var.Dont_Care()};

    if(auto ptr=diffs.Get_Pointer({func_var,diff_var}))
        return *ptr;

    if(func_var.Is_In()) return {zero};
    auto it=Get_Def(func_var);
    
    const INSTRUCTION& o=*it++->second;
    PHYSBAM_ASSERT(o.type.Has_Dest());

    REG_AUG s0,s1,s2;
    if(o.type.Has_Src0()) s0=Diff_Internal(o.src0,diff_var);
    if(o.type.Has_Src1()) s1=Diff_Internal(o.src1,diff_var);
    if(o.type.Has_Src2()) s2=Diff_Internal(o.src2,diff_var);

    if(o.type.Use_Adjoint())
    {
        PHYSBAM_ASSERT(o.type.Has_Src0());
        auto adj=Make_Adjoint(func_var);

        auto parent=func_var;
        if(auto p=diff_parent.Get_Pointer(func_var))
            parent=*p;
        if(auto p=diff_block_end.Get_Pointer(parent))
            it=*p;
        auto diff_last_it=it++;

        if(debug) Test_Invariants();
        REG_AUG d={zero};
        if(o.type.Has_Src0()) Add_Mul(it,d,s0,adj.val[0]);
        if(o.type.Has_Src1()) Add_Mul(it,d,s1,adj.val[1]);
        if(o.type.Has_Src2()) Add_Mul(it,d,s2,adj.val[2]);

        while(++diff_last_it!=it)
            diff_parent.Set(diff_last_it->second->dest,parent);
        diff_block_end.Set(parent,--it);

        diffs.Set({func_var,diff_var},d);
        return d;
    }

    REG d,e,t;
    // First do ops without good adjoints.
    switch(o.type.op){
        case op_phi:
            s0=Prepend_Phi(it->first,s0,s1);
            break;
        case op_if:
            s0=If_Aug(it,o.src0,s1,s2);
            break;
        case op_min:
            t=Before_Instruction(it,op_lt,o.src0,o.src1);
            s0=If_Aug(it,t,s0,s1);
            break;
        case op_max:
            t=Before_Instruction(it,op_gt,o.src0,o.src1);
            s0=If_Aug(it,t,s0,s1);
            break;
        case op_abs:
            t=Before_Instruction(it,op_lt,o.src0,zero);
            s0=If_Aug(it,t,Neg(s0),s0);
            break;
        default:PHYSBAM_FATAL_ERROR("failed to handle non-adjoint case");
    }
    diffs.Set({func_var,diff_var},s0);
    return s0;
}
//#####################################################################
// Function Diff
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Diff(REG func_var,REG diff_var)
{
    diffs.Set({diff_var,diff_var},{MEM_CONST(1)});
    return Clear_Aug(Diff_Internal(func_var,diff_var));
}
//#####################################################################
// Function Diff
//#####################################################################
template<class T> ARRAY<REG,MEM_OUT> PROGRAM_COMPILER<T>::
Diff(REG diff_var)
{
    diffs.Set({diff_var,diff_var},{MEM_CONST(1)});
    ARRAY<REG,MEM_OUT> r(num_out);
    for(MEM_OUT i{};i<r.m;i++)
        r(i)=Clear_Aug(Diff_Internal(i,diff_var));
    return r;
}
//#####################################################################
// Function Diff
//#####################################################################
template<class T> ARRAY<REG> PROGRAM_COMPILER<T>::
Diff(const ARRAY_VIEW<const REG> func_vars,REG diff_var)
{
    diffs.Set({diff_var,diff_var},{MEM_CONST(1)});
    ARRAY<REG> r(func_vars.m);
    for(int i=0;i<func_vars.m;i++)
        r(i)=Clear_Aug(Diff_Internal(func_vars(i),diff_var));
    return r;
}
//#####################################################################
// Function Merge_Conditions
//#####################################################################
static std::pair<REG,bool>
Merge_Conditions(std::map<REG,bool>& r,const std::map<REG,bool>& a,
    const std::map<REG,bool>& b)
{
    std::pair<REG,bool> p={{},false};
    Zip(a,b,[](auto&){},[](auto&){},[&r,&p](auto&ia,auto&ib)
    {
        if(ia.second==ib.second) r.insert(ia);
        else p=ib;
    });
    return p;
}
//#####################################################################
// Function Gradient
//#####################################################################
template<class T> ARRAY<REG,MEM_IN> PROGRAM_COMPILER<T>::
Gradient(REG func_var)
{
    DEBUG;

    MEM_CONST zero=MEM_CONST(0);
    MEM_CONST one=MEM_CONST(1);
    MEM_LABEL func_var_block=Get_Def(func_var)->first->id;
    MEM_LABEL reversal_point=Get_Reversal_Point(func_var_block);
    MEM_LABEL end_point(Value(reversal_point)*2-1);

    ARRAY<std::map<REG,REG_AUG>,MEM_LABEL> block_diff_lookup(reversal_point+1);
    block_diff_lookup(reversal_point-1)[func_var]={one,false,false};

    for(MEM_LABEL i=reversal_point-1,r=reversal_point;i>=MEM_LABEL{};i--,r++)
        if(auto B=code_blocks(i))
        {
            auto C=code_blocks(r);
            auto& m=block_diff_lookup(i);
            if(auto F=B->next[0])
            {
                PHYSBAM_ASSERT(F->id>i);
                auto& f=block_diff_lookup(F->id);
                if(auto S=B->next[1])
                {
                    PHYSBAM_ASSERT(S->id>i);
                    auto& s=block_diff_lookup(S->id);
                    Zip(f,s,
                        [this,&m,C,zero](auto& fi)
                        {m[fi.first]=Prepend_Phi(C,fi.second,{zero,false,false});},
                        [this,&m,C,zero](auto& si)
                        {m[si.first]=Prepend_Phi(C,{zero,false,false},si.second);},
                        [this,&m,C](auto& fi,auto& si)
                        {m[fi.first]=Prepend_Phi(C,fi.second,si.second);});
                }
                else if(i!=reversal_point-1) m=f;
            }

            auto C_it=C->End_It();
            if(C_it!=C->Begin_It())
            {
                auto P=C_it;
                if((--P)->second->type.Is_Branch())
                    C_it=P;
            }
            for(auto it=B->insts.rbegin();it!=B->insts.rend();it++)
            {
                const INSTRUCTION& o=*it;
                if(!o.type.Has_Dest()) continue;
                auto dit=m.find(o.dest);
                if(dit==m.end()) continue;
                const auto& jc=dit->second;
                if(jc.reg==zero) continue;

                auto add=[this,&jc,&m,one,zero,C_it](REG src,REG_AUG diff)
                {
                    if(src.Is_Const() || diff.reg==zero) return;
                    auto p=m.insert({src,{zero,false,false}});
                    Add_Mul(C_it,p.first->second,jc,diff);
                };

                if(o.type.Use_Adjoint())
                {
                    PHYSBAM_ASSERT(o.type.Has_Src0());
                    auto adj=Make_Adjoint(o.dest);
                    if(o.type.Has_Src0()) add(o.src0,Trace_Register(adj.val[0],r));
                    if(o.type.Has_Src1()) add(o.src1,Trace_Register(adj.val[1],r));
                    if(o.type.Has_Src2()) add(o.src2,Trace_Register(adj.val[2],r));
                    m.erase(dit);
                    continue;
                }

                auto cond=[this,&jc,&m,zero,C_it](REG src,REG cond,bool first,bool sub=false)
                {
                    if(src.Is_Const()) return;
                    auto p=m.insert({src,{zero,false,false}});
                    auto& a=p.first->second;
                    auto b=Add(C_it,a,jc);
                    Make_Same_Aug(a,b);
                    a.reg=Before_Instruction(C_it,op_if,cond,first?a.reg:b.reg,first?b.reg:a.reg);
                };

                REG p;
                switch(o.type.op){
                    case op_if:
                        cond(o.src1,o.src0,true);
                        cond(o.src2,o.src0,false);
                        break;
                    case op_min:
                        p=Before_Instruction(C_it,op_lt,o.src0,o.src1);
                        cond(o.src0,p,false);
                        cond(o.src1,p,true);
                        break;
                    case op_max:
                        p=Before_Instruction(C_it,op_gt,o.src0,o.src1);
                        cond(o.src0,p,false);
                        cond(o.src1,p,true);
                        break;
                    case op_abs:
                        {
                            p=Before_Instruction(C_it,op_lt,o.src0,zero);
                            auto i=m.insert({o.src0,{zero,false,false}});
                            auto& a=i.first->second;
                            auto neg=Add(C_it,a,Neg(jc));
                            a=Add(C_it,a,jc);
                            Make_Same_Neg(a,neg);
                            a.reg=Before_Instruction(C_it,op_if,p,neg.reg,a.reg);
                        }
                        break;
                    case op_phi:
                        p=Prepend_Phi(B,one,zero);
                        cond(o.src0,p,false);
                        cond(o.src1,p,true);
                        break;
                    default: PHYSBAM_FATAL_ERROR("Missing diff instruction");
                }
                m.erase(dit);
            }
        }

    ARRAY<REG,MEM_IN> ret(num_in,init_new,zero);
    for(auto p:block_diff_lookup({}))
        if(p.first.Is_In())
            ret(p.first.In())=Clear_Aug(p.second);
    VERBOSE;
    DEBUG;
    return ret;
}
//#####################################################################
// Function Gradient
//#####################################################################
template<class T> ARRAY<ARRAY<REG,MEM_IN>> PROGRAM_COMPILER<T>::
Gradient(ARRAY_VIEW<const REG> func_vars)
{
    ARRAY<ARRAY<REG,MEM_IN>> ret(func_vars.m);
    for(int i=0;i<func_vars.m;i++)
        ret(i)=Gradient(func_vars(i));
    return ret;
}
//#####################################################################
// Function Directional_Derivative
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Directional_Derivative(REG func_var,ARRAY_VIEW<const REG,MEM_IN> dir)
{
    HASHTABLE<std::pair<REG,REG>,REG_AUG> diffs_dir;
    std::swap(diffs,diffs_dir);
    for(MEM_IN i{};i<dir.m;i++) diffs.Set({i,REG{}},{dir(i)});
    auto ret=Diff_Internal(func_var,REG{});
    std::swap(diffs,diffs_dir);
    return Clear_Aug(ret);
}
//#####################################################################
// Function Directional_Derivative
//#####################################################################
template<class T> ARRAY<REG,MEM_OUT> PROGRAM_COMPILER<T>::
Directional_Derivative(ARRAY_VIEW<const REG,MEM_IN> dir)
{
    HASHTABLE<std::pair<REG,REG>,REG_AUG> diffs_dir;
    std::swap(diffs,diffs_dir);
    for(MEM_IN i{};i<dir.m;i++) diffs.Set({i,REG{}},{dir(i)});
    ARRAY<REG,MEM_OUT> r(num_out);
    for(MEM_OUT i{};i<r.m;i++)
        r(i)=Clear_Aug(Diff_Internal(i,REG{}));
    std::swap(diffs,diffs_dir);
    return r;
}
//#####################################################################
// Function Register_Default_Func
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Register_Default_Func()
{
    typedef COMPILER_DATA_OPS<T> cd;

    Register_Func("inv",1,&cd::inv);
    Register_Func("sqrt",1,&cd::sqrt);
    Register_Func("exp",1,&cd::exp);
    Register_Func("ln",1,&cd::ln);
    Register_Func("sin",1,&cd::sin);
    Register_Func("cos",1,&cd::cos);
    Register_Func("asin",1,&cd::asin);
    Register_Func("acos",1,&cd::acos);
    Register_Func("atan",1,&cd::atan);
    Register_Func("atan2",2,&cd::atan2);
    Register_Func("floor",1,&cd::floor);
    Register_Func("ceil",1,&cd::ceil);
    Register_Func("rint",1,&cd::rint);
    Register_Func("abs",1,&cd::abs);
    Register_Func("min",2,&cd::min);
    Register_Func("max",2,&cd::max);
    Register_Func("pow",2,&cd::pow);
    Register_Func("mod",2,&cd::mod);
    Register_Func("tan",1,&cd::tan);
    Register_Func("cot",1,&cd::cot);
    Register_Func("sec",1,&cd::sec);
    Register_Func("csc",1,&cd::csc);
    Register_Func("acot",1,&cd::acot);
    Register_Func("asec",1,&cd::asec);
    Register_Func("acsc",1,&cd::acsc);
    Register_Func("trans",1,&cd::trans);
    Register_Func("dot",2,&cd::dot);
    Register_Func("cross",2,&cd::cross);
    Register_Func("outer",2,&cd::outer);
    Register_Func("col",2,&cd::col);
    Register_Func("row",2,&cd::row);
    Register_Func("tr",1,&cd::tr);
    Register_Func("zeros",-1,&cd::zeros);
    Register_Func("ones",-1,&cd::ones);
    Register_Func("id",1,&cd::id);
    Register_Func("norm2",1,&cd::norm2);
    Register_Func("norm",1,&cd::norm);
    Register_Func("det",1,&cd::det);
    Register_Func("inverse",1,&cd::inverse);
    Register_Func("vec",-1,&cd::vec);
}
//#####################################################################
// Function Register_Func
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Register_Func(const std::string& name,int num_in,DATA(*func)(PROGRAM_COMPILER<T>*pc,const ARRAY<DATA>& in))
{
    Func_Lookup().Set(name,{func,num_in});
}
//#####################################################################
// Function Exec_Func
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Exec_Func(const std::string& name,const ARRAY<DATA>& in) -> DATA
{
    static bool registered=false;
    if(!registered)
    {
        registered=true;
        Register_Default_Func();
    }
    if(auto*f=Func_Lookup().Get_Pointer(name))
    {
        PHYSBAM_ASSERT(in.m==f->second || f->second<0);
        return f->first(this,in);
    }
    LOG::printf("Function not found: %s\n",name);
    PHYSBAM_FATAL_ERROR();
}
//#####################################################################
// Function Print_Node
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Print_Node(PROGRAM_PARSE_NODE* node)
{
    if(node->type<256) LOG::printf("(%c",(char) node->type);
    else LOG::printf("(%i",node->type);
    if(node->a){LOG::printf(" ");Print_Node(node->a);}
    if(node->b){LOG::printf(" ");Print_Node(node->b);}
    LOG::printf(")");
}
//#####################################################################
// Function Process_Assign
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Process_Assign(DATA& a,const DATA& b) -> DATA&
{
    if(a.type==DATA::TYPE::undef)
    {
        PHYSBAM_ASSERT(a.ref);
        a.ref->Set_Type(b);
        for(auto&i:a.ref->x) i=New_Reg();
        a=*a.ref;
    }
    if(a.Same_Type(b))
        for(int i=0;i<a.x.m;i++)
            Append_Instruction_Dest(op_copy,a.x(i),b.x(i));
    else
        PHYSBAM_FATAL_ERROR("No match for "+a.Type_Str()+" = "+b.Type_Str());
    return a;
}
//#####################################################################
// Function Process_Componentwise
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Process_Componentwise(const DATA& a,const DATA& b,const OPERATION& op,
    bool allow_s_mv,const char* s0,const char* s1,const char* s2) -> DATA
{
    if(a.Same_Type(b))
    {
        DATA r(a);
        for(int i=0;i<a.x.m;i++)
            r.x(i)=Append_Instruction(op,a.x(i),b.x(i));
        return r;
    }
    if(b.type==DATA::TYPE::scalar)
    {
        DATA r(a);
        for(int i=0;i<a.x.m;i++)
            r.x(i)=Append_Instruction(op,a.x(i),b.x(0));
        return r;
    }
    if(allow_s_mv && a.type==DATA::TYPE::scalar)
    {
        DATA r(b);
        for(int i=0;i<b.x.m;i++)
            r.x(i)=Append_Instruction(op,a.x(0),b.x(i));
        return r;
    }
    PHYSBAM_FATAL_ERROR((std::string)"No match for "+s0+a.Type_Str()+s1+b.Type_Str()+s2+".");
}
//#####################################################################
// Function Process_Componentwise
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Process_Componentwise(const DATA& a,const OPERATION& op) -> DATA
{
    DATA r(a);
    for(int i=0;i<a.x.m;i++)
        r.x(i)=Append_Instruction(op,a.x(i));
    return r;
}
//#####################################################################
// Function Process_Componentwise
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Process_Mod(const DATA& a,const DATA& b,bool allow_s_mv,const char* s0,const char* s1,const char* s2) -> DATA
{
    auto c=Process_Componentwise(a,b,op_div,allow_s_mv,s0,s1,s2);
    auto d=Process_Componentwise(c,op_floor);
    auto e=Process_Componentwise(d,b,op_mul,true,s0,s1,s2);
    return Process_Componentwise(a,e,op_sub,allow_s_mv,s0,s1,s2);
}
//#####################################################################
// Function Process_Componentwise
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Process_Mul(const DATA& a,const DATA& b) -> DATA
{
    if(a.type==DATA::TYPE::scalar || b.type==DATA::TYPE::scalar)
        return Process_Componentwise(a,b,op_mul,true,"","*","");
    if(a.type==DATA::TYPE::mat && b.type==DATA::TYPE::vec && a.n==b.m)
    {
        DATA r(a.m);
        for(int i=0;i<a.m;i++)
        {
            REG x=MEM_CONST(0);
            for(int j=0;j<a.n;j++)
                x=Append_Instruction(op_add,x,Append_Instruction(op_mul,a(i,j),b.x(j)));
            r.x(i)=x;
        }
        return r;
    }
    if(a.type==DATA::TYPE::mat && b.type==DATA::TYPE::mat && a.n==b.m)
    {
        DATA r(a.m,b.n);
        for(int i=0;i<a.m;i++)
            for(int j=0;j<b.n;j++)
            {
                REG x=MEM_CONST(0);
                for(int k=0;k<a.n;k++)
                    x=Append_Instruction(op_add,x,Append_Instruction(op_mul,a(i,k),b(k,j)));
                r(i,j)=x;
            }
        return r;
    }
    PHYSBAM_FATAL_ERROR((std::string)"No match for "+a.Type_Str()+"*"+b.Type_Str()+".");
}
//#####################################################################
// Function Process_Node
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Process_Node(PROGRAM_PARSE_NODE* node) -> DATA
{
#ifdef USE_LEX_YACC
    if(node->type==TOKEN_FUNC)
    {
        PHYSBAM_ASSERT(node->a && node->a->type==TOKEN_IDENT);
        ARRAY<DATA> args;
        for(PROGRAM_PARSE_NODE* n=node->b;n;n=n->b)
            args.Append(Process_Node(n->a));
        return Exec_Func(parse_identifiers(node->a->val),args);
    }
    if(node->type==TOKEN_MATRIX)
    {
        int n=0,m=0,k=0;
        ARRAY<DATA> args;
        for(PROGRAM_PARSE_NODE* N=node->a;N;N=N->b)
        {
            k++;
            args.Append(Process_Node(N->a));
            if(N->type==';')
            {
                if(m>0){PHYSBAM_ASSERT(n==k);}
                else n=k;
                m++;
                k=0;
            }
        }
        if(m>0){PHYSBAM_ASSERT(n==k);}
        else n=k;
        if(n>0) m++;
        for(auto&i:args){PHYSBAM_ASSERT(i.type==DATA::TYPE::scalar);}
        DATA r(m,n);
        for(int i=0;i<m;i++)
            for(int j=0;j<n;j++)
                r.x(i+m*j)=args(j+n*i).x(0);
        return r;
    }

    DATA aa,bb;
    if(node->a) aa=Process_Node(node->a);
    if(node->b && node->type!='?' && node->type!=TOKEN_INDEX) bb=Process_Node(node->b);
    switch(node->type)
    {
        case TOKEN_LIST:
        case ',': return bb;
        case '=': return Process_Assign(aa,bb);
        case TOKEN_COMPADDEQ:
            return Process_Assign(aa,Process_Componentwise(aa,bb,op_add,false,"",".+=",""));
        case TOKEN_ADDEQ:
            return Process_Assign(aa,Process_Componentwise(aa,bb,op_add,false,"","+=",""));
        case TOKEN_COMPSUBEQ:
            return Process_Assign(aa,Process_Componentwise(aa,bb,op_sub,false,"",".-=",""));
        case TOKEN_SUBEQ:
            return Process_Assign(aa,Process_Componentwise(aa,bb,op_sub,false,"","-=",""));
        case TOKEN_COMPMULEQ:
            return Process_Assign(aa,Process_Componentwise(aa,bb,op_mul,false,"",".*=",""));
        case TOKEN_MULEQ:
            if(bb.type==DATA::TYPE::scalar)
                return Process_Assign(aa,Process_Componentwise(aa,bb,op_mul,false,"","*=",""));
            PHYSBAM_FATAL_ERROR("No match for "+aa.Type_Str()+"*="+bb.Type_Str());
        case TOKEN_COMPDIVEQ:
            return Process_Assign(aa,Process_Componentwise(aa,bb,op_div,false,"","./=",""));
        case TOKEN_DIVEQ:
            if(bb.type==DATA::TYPE::scalar)
                return Process_Assign(aa,Process_Componentwise(aa,bb,op_div,false,"","/=",""));
            PHYSBAM_FATAL_ERROR("No match for "+aa.Type_Str()+"/="+bb.Type_Str());
        case TOKEN_COMPMODEQ:
            return Process_Assign(aa,Process_Mod(aa,bb,false,"",".%=",""));
        case TOKEN_MODEQ:
            return Process_Assign(aa,Process_Mod(aa,bb,false,"","%=",""));
        case TOKEN_COMPPOWEQ:
            return Process_Assign(aa,Process_Componentwise(aa,bb,op_pow,false,"",".^=",""));
        case TOKEN_POWEQ:
            return Process_Assign(aa,Process_Componentwise(aa,bb,op_pow,false,"","^=",""));
        case '<':
            return Process_Componentwise(aa,bb,op_lt,true,"","<","");
        case '>':
            return Process_Componentwise(aa,bb,op_gt,true,"",">","");
        case TOKEN_EQ:
            return Process_Componentwise(aa,bb,op_eq,true,"","==","");
        case TOKEN_NE:
            return Process_Componentwise(aa,bb,op_ne,true,"","!=","");
        case TOKEN_LE:
            return Process_Componentwise(aa,bb,op_le,true,"","<=","");
        case TOKEN_GE:
            return Process_Componentwise(aa,bb,op_ge,true,"",">=","");
        case TOKEN_AND:
            return Process_Componentwise(aa,bb,op_and,true,"","&&","");
        case TOKEN_OR:
            return Process_Componentwise(aa,bb,op_or,true,"","||","");
        case TOKEN_COMPADD:
            return Process_Componentwise(aa,bb,op_add,true,"",".+","");
        case '+':
            if(node->b) return Process_Componentwise(aa,bb,op_add,true,"","+","");
            return Process_Componentwise(aa,op_copy);
        case TOKEN_COMPSUB:
            return Process_Componentwise(aa,bb,op_sub,true,"",".-","");
        case '-':
            if(node->b) return Process_Componentwise(aa,bb,op_sub,true,"","-","");
            return Process_Componentwise(aa,op_neg);
        case TOKEN_COMPMUL:
            return Process_Componentwise(aa,bb,op_mul,true,"",".*","");
        case '*': return Process_Mul(aa,bb);
        case TOKEN_COMPDIV:
            return Process_Componentwise(aa,bb,op_div,true,"","./","");
        case '/':
            if(bb.type==DATA::TYPE::scalar)
                return Process_Componentwise(aa,bb,op_div,true,"","/","");
            PHYSBAM_FATAL_ERROR("No match for "+aa.Type_Str()+"/"+bb.Type_Str());
        case TOKEN_COMPMOD:
            return Process_Mod(aa,bb,true,"",".%","");
        case '%':
            return Process_Mod(aa,bb,true,"","%","");
        case TOKEN_COMPPOW:
            return Process_Componentwise(aa,bb,op_pow,true,"",".^","");
        case '^':
            return Process_Componentwise(aa,bb,op_pow,true,"","^","");
        case '!':
            return Process_Componentwise(aa,op_not);
        case TOKEN_NUMBER:
            return DATA(Add_Constant(parse_constants(node->val)));
        case TOKEN_IDENT:{
            auto p=dict.Insert(parse_identifiers(node->val),nullptr);
            if(p.second)
            {
                *p.first=new DATA;
                auto d=**p.first;
                d.ref=*p.first;
                return d;
            }
            return **p.first;
        }
        case TOKEN_INDEX:
            {
                PHYSBAM_ASSERT(node->b->type==TOKEN_LIST);
                PHYSBAM_ASSERT(node->b->a);
                PHYSBAM_ASSERT(node->b->a->type==TOKEN_NUMBER);
                int idx0=parse_constants(node->b->a->val);
                PHYSBAM_ASSERT(idx0>=0 && idx0<aa.m);
                if(aa.type==DATA::TYPE::vec)
                {
                    PHYSBAM_ASSERT(!node->b->b);
                    return DATA(aa.x(idx0));
                }
                PHYSBAM_ASSERT(aa.type==DATA::TYPE::mat);
                PHYSBAM_ASSERT(node->b->b->type==TOKEN_NUMBER);
                int idx1=parse_constants(node->b->b->val);
                PHYSBAM_ASSERT(idx1>=0 && idx1<aa.n);
                return DATA(aa(idx0,idx1));
            }

        case '?':{
            PHYSBAM_ASSERT(node->b && node->b->type==':');
            DATA result;
            result.type=DATA::TYPE::undef;
            result.ref=&result;
            CODE_BLOCK* A=code_blocks.Last();
            CODE_BLOCK* B=new CODE_BLOCK;
            B->id=code_blocks.Append(B);
            if(node->b->a) Process_Assign(result,Process_Node(node->b->a));
            CODE_BLOCK* BB=code_blocks.Last();
            CODE_BLOCK* C=new CODE_BLOCK;
            C->id=code_blocks.Append(C);
            if(node->b->b) Process_Assign(result,Process_Node(node->b->b));
            CODE_BLOCK* CC=code_blocks.Last();
            CODE_BLOCK* D=new CODE_BLOCK;
            D->id=code_blocks.Append(D);
            PHYSBAM_ASSERT(aa.type==DATA::TYPE::scalar);
            INSTRUCTION in_A={op_brz,C->id,aa.x(0),B->id};
            A->insts.push_back(in_A);

            A->next[0]=B;
            B->prev[0]=A;
            A->next[1]=C;
            C->prev[0]=A;

            D->prev[0]=BB;
            D->prev[1]=CC;
            BB->next[0]=D;
            CC->next[0]=D;
            return result;}

        default: PHYSBAM_FATAL_ERROR("Unexpected node");
    }
#else
    PHYSBAM_FATAL_ERROR("Parsing support disabled; compile with USE_LEX_YACC.");
#endif
}
//#####################################################################
// Function Parse
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Parse(const std::string& str)
{
    code_blocks.Delete_Pointers_And_Clean_Memory();

    CODE_BLOCK* A=new CODE_BLOCK;
    A->id=code_blocks.Append(A);

    auto make_reg=[this](REG r)
    {
        REG x=New_Reg();
        INSTRUCTION in={op_copy,x,r};
        code_blocks({})->insts.push_back(in);
        return x;
    };

    dict.Set("pi",new DATA(Add_Constant(pi)));
    for(auto&p:var_in)
    {
        auto d=new DATA(p.second);
        for(auto&r:d->x) r=make_reg(r);
        dict.Set(p.first,d);
    }
    for(auto&p:var_out)
    {
        if(auto ptr=dict.Get_Pointer(p.first))
        {
            PHYSBAM_ASSERT((*ptr)->Same_Type(p.second));
        }
        else
        {
            auto d=new DATA(p.second);
            for(auto&r:d->x) r=New_Reg();
            dict.Set(p.first,d);
        }
    }

    Parse_Command(str);

    ARRAY<INSTRUCTION> zero_inst;
    for(auto&p:var_out)
    {
        auto&a=p.second.x;
        auto&b=dict.Get(p.first)->x;
        for(int i=0;i<a.m;i++)
        {
            INSTRUCTION in0={op_copy,b(i),MEM_CONST(0)};
            code_blocks({})->insts.push_front(in0);
            INSTRUCTION in1={op_copy,a(i),b(i)};
            code_blocks.Last()->insts.push_back(in1);
        }
    }

    VERBOSE;
    DEBUG;
    Make_SSA();
    Copy_Propagation();
    Remove_Dead_Code();
    Compress_Registers();
}
//#####################################################################
// Function Parse_Command
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Parse_Command(const std::string& str)
{
#ifdef USE_LEX_YACC
    YY_BUFFER_STATE state=yy_scan_string(str.c_str());
    yyparse();
    Process_Node(parse_root);
    parse_identifiers.Remove_All();
    parse_constants.Remove_All();
    yy_delete_buffer(state);
    yylex_destroy();
#else
    PHYSBAM_FATAL_ERROR("Parsing support disabled; compile with USE_LEX_YACC.");
#endif
}
//#####################################################################
// Function Print
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Print(const CODE_BLOCK* B) const
{
    for(const auto& o:*B) LOG::cout<<o<<"\n";
}
//#####################################################################
// Function Print
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Print() const
{
    for(MEM_CONST i{};i<constants.m;i++)
        LOG::printf("const c%d = %g\n", i, constants(i));

    LOG::printf("\n==== BLOCK CODE ====\n");
    for(MEM_LABEL i{};i<code_blocks.m;i++)
    {
        auto B=code_blocks(i);
        if(!B)
        {
            LOG::printf("== [block L%i] ==\n",i);
            continue;
        }
        LOG::printf("== block L%i ==\n",i);
        if(B->prev[0] && B->prev[1]) LOG::printf("(prev: %i %i)\n",B->prev[0]->id,B->prev[1]->id);
        Print(B);
        if(B->next[0] && B->next[0]->id!=i+1)
        {
            INSTRUCTION jmp(op_jmp,B->next[0]->id);
            LOG::cout<<jmp<<std::endl;
        }
    }
}
//#####################################################################
// Function Finalize
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Finalize(PROGRAM<T>& prog)
{
    DEBUG;
    auto& mem_defs=defs.a(mem_out);
    for(int j=0;j<mem_defs.m;j++){
        auto it=mem_defs(j);
        if(!it->first) continue;
        INSTRUCTION o=*it->second;
        if(o.type==op_copy && o.src0.Is_Reg()){
            Change_Instruction(it,{});
            Rename_Variable(o.src0,o.dest);
            DEBUG;
        }}

    Leave_SSA();

    int off_in=Value(num_reg),off_out=off_in+Value(num_in),off_const=off_out+Value(num_out);
    int offs[5]={0,off_const,off_in,off_out,0};
    ARRAY<int> labels(Value(code_blocks.m+1));
    for(MEM_LABEL i{};i<code_blocks.m;i++)
    {
        if(prog.code.m && prog.code.Last().type==op_jmp && prog.code.Last().dest==Value(i))
            prog.code.Pop();
        CODE_BLOCK* B=code_blocks(i);
        if(!B) continue;
        labels(Value(B->id))=prog.code.m;
        for(auto o:*B)
        {
            if(o.type==op_nop) continue;
            if(o.type==op_copy && o.dest==o.src0) continue;
            RAW_INSTRUCTION r;
            r.type=o.type;
            if(o.type.Has_Dest() || o.type.Is_Jump()){auto [t,v]=o.dest.Parts();r.dest=offs[t]+v;}
            if(o.type.Has_Src0()){auto [t,v]=o.src0.Parts();r.src0=offs[t]+v;}
            if(o.type.Has_Src1() || o.type.Is_Branch()){auto [t,v]=o.src1.Parts();r.src1=offs[t]+v;}
            if(o.type.Has_Src2()){auto [t,v]=o.src2.Parts();r.src2=offs[t]+v;}
            prog.code.Append(r);
        }

        if(B->next[0]) prog.code.Append({op_jmp,Value(B->next[0]->id),0,0});
    }
    if(prog.code.m && prog.code.Last().type==op_jmp && prog.code.Last().dest==Value(code_blocks.m))
        prog.code.Pop();
    labels.Last()=prog.code.m;
    for(int ip=0;ip<prog.code.m;ip++)
    {
        RAW_INSTRUCTION& o=prog.code(ip);
        if(o.type.Is_Jump()) o.dest=labels(o.dest);
        if(o.type.Is_Branch()) o.src1=labels(o.src1);
    }

    prog.num_out=num_out;
    prog.num_reg=Value(num_reg);
    prog.num_in=num_in;
    prog.constants=constants;
    VERBOSE;
    DEBUG;
    if(verbose) prog.Print();
}
//#####################################################################
// Function Optimize
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Optimize()
{
    Clear_Diff_Structures();
    Copy_Propagation();
    Remove_Dead_Code();
    Compress_Registers();

    Sparse_Conditional_Constant_Propagation();
    Simplify_Phis();
    Copy_Propagation();
    Simplify_With_Distributive_Law();
    Local_Common_Expresssion_Elimination();
    Remove_Dead_Code();
    Compress_Registers();
}
//#####################################################################
// Function Make_SSA
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Make_SSA()
{
    DEBUG;
    ARRAY<ARRAY<MEM_LABEL>,MEM_LABEL> graph_in(code_blocks.m),graph_out(code_blocks.m);
    for(MEM_LABEL i{};i<code_blocks.m;i++){
        if(!code_blocks(i)) continue;
        for(int s=0;s<2;s++){
            if(code_blocks(i)->prev[s])
                graph_in(i).Append(code_blocks(i)->prev[s]->id);
            if(code_blocks(i)->next[s])
                graph_out(i).Append(code_blocks(i)->next[s]->id);}}

    DOMINATORS dom(graph_in,graph_out);
    dom.Compute_Frontier();

    ARRAY<REG,MEM_LABEL> has_already(code_blocks.m),work(code_blocks.m);
    ARRAY<MEM_REG,MEM_REG> V;
    HASHTABLE<REG,HASHTABLE<MEM_LABEL> > HA;
    Set_Counting(V,num_reg);
    HASHTABLE<REG> V_used;
    for(auto v:V) V_used.Set(v);
    for(MEM_LABEL bl{};bl<code_blocks.m;bl++)
        if(CODE_BLOCK* B=code_blocks(bl))
            for(const auto& o:*B){
                if(o.type.Has_Dest())
                    HA.Get_Or_Insert(o.dest).Set(bl);
                if(o.type.Has_Src0())
                    HA.Get_Or_Insert(o.src0).Set(bl);
                if(o.type.Has_Src1())
                    HA.Get_Or_Insert(o.src1).Set(bl);
                if(o.type.Has_Src2())
                    HA.Get_Or_Insert(o.src2).Set(bl);}

    for(MEM_REG vi{};vi<V.m;vi++){

        MEM_REG v=V(vi);
        ARRAY<MEM_LABEL> W;
        HA.Get_Or_Insert(v).Get_Keys(W);
        work.Subset(W).Fill(vi);
        while(W.m){
            MEM_LABEL x=W.Pop_Value();
            for(int i=0;i<dom.frontier(x).m;i++){
                MEM_LABEL y=dom.frontier(x)(i);
                if(!has_already(y).Is_Value() || has_already(y)<vi){
                    if(code_blocks(y)->prev[0] && code_blocks(y)->prev[1]){
                        INSTRUCTION in={op_phi,v,v,v};
                        code_blocks(y)->insts.push_front(in);}
                    has_already(y)=vi;
                    if(!work(y).Is_Value() || work(y)<vi){
                        work(y)=vi;
                        W.Append(y);}}}}}

    ARRAY<ARRAY<REG>,MEM_REG> S(V.m);
    REG undef=MEM_CONST(0);
    for(MEM_REG i{};i<S.m;i++)
        S(i).Append(undef);
    Make_SSA_Relabel(num_reg,S,V_used,dom.children,{});
    for(MEM_REG i{};i<S.m;i++){
        assert(S(i).m==1 && S(i)(0)==undef);}

    Compute_Use_Def();
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Make_SSA_Relabel
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Make_SSA_Relabel(MEM_REG& count,ARRAY<ARRAY<REG>,MEM_REG>& S,HASHTABLE<REG>& V_used,const ARRAY<ARRAY<MEM_LABEL>,MEM_LABEL>& children,MEM_LABEL bl)
{
    ARRAY<MEM_REG> S_pop;
    CODE_BLOCK* B=code_blocks(bl);
    for(auto& o:*B){
        if(o.type!=op_phi){
            if(o.type.Has_Src0())
                if(V_used.Contains(o.src0))
                    o.src0=S(o.src0.Reg()).Last();
            if(o.type.Has_Src1())
                if(V_used.Contains(o.src1))
                    o.src1=S(o.src1.Reg()).Last();
            if(o.type.Has_Src2())
                if(V_used.Contains(o.src2))
                    o.src2=S(o.src2.Reg()).Last();}
        if(o.type.Has_Dest()){
            if(V_used.Contains(o.dest)){
                MEM_REG n=count++;
                S_pop.Append(o.dest.Reg());
                S(o.dest.Reg()).Append(n);
                o.dest=n;}}}

    for(int s=0;s<2;s++)
        if(CODE_BLOCK* nb=code_blocks(bl)->next[s]){
            int j=nb->prev[1]==code_blocks(bl);
            for(auto& o:*nb){
                if(o.type!=op_phi) break;
                REG* reg=j?&o.src1:&o.src0;
                *reg=S(reg->Reg()).Last();}}

    for(int i=0;i<children(bl).m;i++)
        Make_SSA_Relabel(count,S,V_used,children,children(bl)(i));

    for(int i=0;i<S_pop.m;i++)
        S(S_pop(i)).Pop();
}

struct LIVE_INTERVAL
{
    MEM_LABEL block;
    int b,e;

    int cmp(const LIVE_INTERVAL& l) const
    {
        if(block<l.block) return -1;
        if(l.block<block) return 1;
        if(e<=l.b) return -1;
        if(l.e<=b) return 1;
        return 0;
    }
};

std::ostream& operator<<(std::ostream& o,const LIVE_INTERVAL& l)
{
    LOG::printf("(%P: %P-%P)",l.block,l.b,l.e);
    return o;
}

bool Can_Merge(const ARRAY<LIVE_INTERVAL>& a,const ARRAY<LIVE_INTERVAL>& b)
{
    int i=0,j=0;
    while(i<a.m && j<b.m)
    {
        int c=a(i).cmp(b(j));
        if(c<0) i++;
        else if(c>0) j++;
        else return false;
    }
    return true;
}

ARRAY<LIVE_INTERVAL> Merge(const ARRAY<LIVE_INTERVAL>& a,const ARRAY<LIVE_INTERVAL>& b)
{
    ARRAY<LIVE_INTERVAL> o;
    auto append=[&o](const LIVE_INTERVAL& l)
    {
        if(o.m && o.Last().block==l.block && o.Last().e==l.b)
            o.Last().e=l.e;
        else o.Append(l);
    };
    
    int i=0,j=0;
    while(i<a.m && j<b.m)
    {
        int c=a(i).cmp(b(j));
        if(c<0) append(a(i++));
        else append(b(j++));
    }
    while(i<a.m) append(a(i++));
    while(j<b.m) append(b(j++));
    return o;
}

template<class T> void PROGRAM_COMPILER<T>::
Compute_Live(HASHTABLE<REG,ARRAY<LIVE_INTERVAL>>& live_ranges) const
{
    ARRAY<HASHTABLE<REG>,MEM_LABEL> live_at_end(code_blocks.m);
    for(MEM_OUT i{};i<num_out;i++) live_at_end.Last().Set(i);

    for(MEM_LABEL b=code_blocks.m-1;b>=MEM_LABEL{};b--)
    {
        auto* B=code_blocks(b);
        if(!B) continue;
        HASHTABLE<REG,LIVE_INTERVAL> hash;
        int i=B->insts.size();
        for(auto reg:live_at_end(b))
            hash.Set(reg,LIVE_INTERVAL{b,0,i});

        auto& insts=B->insts;
        for(auto it=insts.rbegin();it!=insts.rend();it++)
        {
            i--;
            if(it->type.Has_Dest())
            {
                auto li=hash.Get(it->dest);
                li.b=i;
                live_ranges.Get_Or_Insert(it->dest).Append(li);
                hash.Delete(it->dest);
            }
            if(it->type==op_phi)
            {
                if(!it->src0.Is_Dont_Care() && !it->src0.Is_Const())
                    live_at_end(B->prev[0]->id).Set(it->src0);
                if(!it->src1.Is_Dont_Care() && !it->src1.Is_Const())
                    live_at_end(B->prev[1]->id).Set(it->src1);
                continue;
            }
            if(it->type.Has_Src0() && !it->src0.Is_Const())
                hash.Insert(it->src0,{b,0,i});
            if(it->type.Has_Src1() && !it->src1.Is_Const())
                hash.Insert(it->src1,{b,0,i});
            if(it->type.Has_Src2() && !it->src2.Is_Const())
                hash.Insert(it->src2,{b,0,i});
        }
        for(auto&p:hash)
        {
            live_ranges.Get_Or_Insert(p.first).Append(p.second);
            for(int j=0;j<2;j++)
                if(B->prev[j])
                    live_at_end(B->prev[j]->id).Set(p.first);
        }
    }
    for(auto&lr:live_ranges)
        lr.second.Reverse();
}

REG Get_Reg(ARRAY<REG,MEM_REG>& remap,REG r)
{
    if(!r.Is_Reg()) return r;
    auto s=r.Reg();
    auto& u=remap(s);
    if(u==r) return r;
    u=Get_Reg(remap,u);
    return u;
};

//#####################################################################
// Function Leave_SSA
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Leave_SSA()
{
    DEBUG;
    Invalidate_Use_Def();
    
    HASHTABLE<REG,ARRAY<LIVE_INTERVAL>> live_ranges;
    Compute_Live(live_ranges);
    ARRAY<REG,MEM_REG> remap(num_reg);
    for(MEM_REG i{};i<remap.m;i++) remap(i)=i;
    
    auto try_merge=[&live_ranges,&remap](REG d,REG r)
    {
        if(d==r) return false;
        if(!r.Is_Reg()) std::swap(r,d);
        if(!r.Is_Reg()) return false;
        if(d.Is_Const()) return false;
        auto& a=live_ranges.Get(d);
        auto& b=live_ranges.Get(r);
        if(!Can_Merge(a,b)) return false;
        live_ranges.Set(d,Merge(a,b));
        live_ranges.Delete(r);
        remap(r.Reg())=d;
        return true;
    };

    For_Each_Instruction_Reverse(code_blocks,
        [&live_ranges,try_merge,&remap](auto& o)
        {
            if(!o.type.Has_Dest()) return true;
            REG dest=Get_Reg(remap,o.dest);
            if(o.type==op_phi)
            {
                REG src0=Get_Reg(remap,o.src0);
                REG src1=Get_Reg(remap,o.src1);
                if(!src0.Is_Dont_Care()) try_merge(dest,src0);
                if(!src1.Is_Dont_Care()) try_merge(dest,src1);
                return true;
            }
            if(o.type.Has_Src0() && try_merge(dest,Get_Reg(remap,o.src0))) return true;
            if(o.type.Has_Src1() && try_merge(dest,Get_Reg(remap,o.src1))) return true;
            if(o.type.Has_Src2() && try_merge(dest,Get_Reg(remap,o.src2))) return true;
            return true;
        });

    {
        std::map<std::pair<MEM_LABEL,int>,MEM_REG> find_end;
        for(auto&p:live_ranges)
            if(p.first.Is_Reg())
            {
                auto& li=p.second.Last();
                find_end[{li.block,li.e}]=p.first.Reg();
            }

        for(auto&p:live_ranges)
            if(p.first.Is_Reg() || p.first.Is_Out())
            {
                auto& li=p.second(0);
                std::pair<MEM_LABEL,int> k{li.block,li.b};
                auto it=find_end.upper_bound(k);
                if(it==find_end.begin()) continue;
                it--;
                try_merge(p.first,it->second);
                find_end.erase(it);
            }
    }

    {
        std::map<std::pair<MEM_LABEL,int>,MEM_REG> find_beg;
        for(auto&p:live_ranges)
            if(p.first.Is_Reg())
            {
                auto& li=p.second(0);
                find_beg[{li.block,li.b}]=p.first.Reg();
            }

        for(auto&p:live_ranges)
            if(p.first.Is_Reg() || p.first.Is_In())
            {
                auto& li=p.second.Last();
                std::pair<MEM_LABEL,int> k{li.block,li.e};
                auto it=find_beg.lower_bound(k);
                if(it==find_beg.end()) continue;
                try_merge(p.first,it->second);
                find_beg.erase(it);
            }
    }
    
    Relabel_Registers(remap);

    ARRAY<REG,MEM_REG> var_map(num_reg);
    for(MEM_REG i{};i<num_reg;i++) var_map(i)=i;
    // Be naive for now.
    for(MEM_LABEL i{};i<code_blocks.m;i++)
        if(CODE_BLOCK* B=code_blocks(i))
            for(auto& o:*B){
                if(o.type!=op_phi) break;
                if(o.src0==o.src1){
                    if(o.dest.Is_Reg()) var_map(o.dest.Reg())=o.src0;
                    else{
                        o.type=op_copy;
                        continue;}}
                else{
                    INSTRUCTION in0={op_copy,o.dest,o.src0};
                    INSTRUCTION in1={op_copy,o.dest,o.src1};
                    if(!o.src0.Is_Dont_Care())
                        code_blocks(i)->prev[0]->insts.push_back(in0);
                    if(!o.src1.Is_Dont_Care())
                        code_blocks(i)->prev[1]->insts.push_back(in1);}
                o.type=op_nop;}
    Relabel_Registers(var_map);
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Invalidate_Use_Def
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Invalidate_Use_Def()
{
    use_def_valid=false;
    for(int i=0;i<4;i++) uses.a(i).Remove_All();
    for(int i=0;i<4;i++) defs.a(i).Remove_All();
}
//#####################################################################
// Function Compute_Use_Def
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Compute_Use_Def()
{
    use_def_valid=true;
    for(int i=0;i<4;i++) uses.a(i).Remove_All();
    for(int i=0;i<4;i++) defs.a(i).Remove_All();
    uses.Resize(num_reg);
    uses.Resize(constants.m);
    uses.Resize(num_in);
    uses.Resize(num_out);
    defs.Resize(num_reg,{});
    defs.Resize(num_out,{});

    for(MEM_LABEL bl{};bl<code_blocks.m;bl++)
        if(CODE_BLOCK* B=code_blocks(bl))
            for(auto it=B->insts.begin();it!=B->insts.end();it++)
                Add_Use_Def({B,it});
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Relabel_Registers
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Relabel_Registers(ARRAY<REG,MEM_REG>& var_map)
{
    UNION_FIND<MEM_REG> uf(var_map.m);
    ARRAY<REG,MEM_REG> new_map(var_map.m);
    for(MEM_REG i{};i<var_map.m;i++)
        if(var_map(i).Is_Reg())
            uf.Union(var_map(i).Reg(),i);
    for(MEM_REG i{};i<var_map.m;i++)
        if(!var_map(i).Is_Reg())
            new_map(uf.Find(i))=var_map(i);

    MEM_REG k{};
    for(MEM_REG i{};i<var_map.m;i++){
        MEM_REG p=uf.Find(i);
        if(new_map(p)==REG()) new_map(p)=k++;
        new_map(i)=new_map(p);}

    For_Each_Iterator(code_blocks,
        [this,&new_map](const auto& it)
        {
            auto& o=*it->second;
            Remove_Use_Def(it);
            if(o.type.Has_Dest() && o.dest.Is_Reg())
                o.dest=new_map(o.dest.Reg());
            if(o.type.Has_Src0() && o.src0.Is_Reg())
                o.src0=new_map(o.src0.Reg());
            if(o.type.Has_Src1() && o.src1.Is_Reg())
                o.src1=new_map(o.src1.Reg());
            if(o.type.Has_Src2() && o.src2.Is_Reg())
                o.src2=new_map(o.src2.Reg());
            Add_Use_Def(it);
            return true;
        });

    num_reg=k;
    var_map.Exchange(new_map);
    // Not appropriate to update use-def, since not SSA anymore.  Invalidate
    // instead.  This should not be used anymore.
    for(int i=0;i<4;i++) uses.a(i).Remove_All();
    for(int i=0;i<4;i++) defs.a(i).Remove_All();
//    Update_Use_Def();
}
//#####################################################################
// Function SCCP_Visit_Instruction
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
SCCP_Visit_Instruction(INST_IT it,REG_ARRAY<VARIABLE_STATE>& variable_state,
    ARRAY<CODE_BLOCK*>& block_worklist,ARRAY<INST_IT>& op_worklist,ARRAY<bool,MEM_LABEL>& block_exec)
{
    const INSTRUCTION& o=*it->second;
    if(o.type==op_phi)
    {
        VARIABLE_STATE& d=variable_state(o.dest);
        if(d.state==VARIABLE_STATE::overdetermined) return;
        if(o.src0.Is_Dont_Care() || o.src1.Is_Dont_Care())
        {
            VARIABLE_STATE s=variable_state(o.src0.Is_Dont_Care()?o.src1:o.src0);
            if(s.state==VARIABLE_STATE::overdetermined || s.state==VARIABLE_STATE::constant)
            {
                if(d.state==s.state) return;
                d.state=s.state;
                op_worklist.Append(it);
            }
            return;
        }

        VARIABLE_STATE s0=variable_state(o.src0);
        VARIABLE_STATE s1=variable_state(o.src1);

        if(s0.state==VARIABLE_STATE::overdetermined || s1.state==VARIABLE_STATE::overdetermined){
            d.state=VARIABLE_STATE::overdetermined;
            op_worklist.Append(it);
            return;}

        if(s0.state==VARIABLE_STATE::constant && s1.state==VARIABLE_STATE::constant){
            if(s0.value!=s1.value) d.state=VARIABLE_STATE::overdetermined;
            else if(d.state==VARIABLE_STATE::constant) return;
            else{d.state=VARIABLE_STATE::constant;d.value=s0.value;}
            op_worklist.Append(it);
            return;}

        if(s0.state==VARIABLE_STATE::constant || s1.state==VARIABLE_STATE::constant){
            if(d.state==VARIABLE_STATE::constant) return;
            d.state=VARIABLE_STATE::constant;
            d.value=s0.state==VARIABLE_STATE::constant?s0.value:s1.value;
            op_worklist.Append(it);
            return;}
        return;}

    if(o.type.Has_Dest()){
        VARIABLE_STATE& d=variable_state(o.dest);
        if(d.state==VARIABLE_STATE::overdetermined) return;
        VARIABLE_STATE* s0=o.type.Has_Src0()?&variable_state(o.src0):0;
        VARIABLE_STATE* s1=o.type.Has_Src1()?&variable_state(o.src1):0;
        VARIABLE_STATE* s2=o.type.Has_Src2()?&variable_state(o.src2):0;

        bool sc0=(s0 && s0->state==VARIABLE_STATE::constant);
        bool sc1=(s1 && s1->state==VARIABLE_STATE::constant);
        bool sc2=(s2 && s2->state==VARIABLE_STATE::constant);
        bool wc0=(!s0 || s0->state==VARIABLE_STATE::constant);
        bool wc1=(!s1 || s1->state==VARIABLE_STATE::constant);
        bool wc2=(!s2 || s2->state==VARIABLE_STATE::constant);
        bool so0=(s0 && s0->state==VARIABLE_STATE::overdetermined);
        bool so1=(s1 && s1->state==VARIABLE_STATE::overdetermined);
        bool so2=(s2 && s2->state==VARIABLE_STATE::overdetermined);
        bool wo0=(!s0 || s0->state==VARIABLE_STATE::overdetermined);
        bool wo1=(!s1 || s1->state==VARIABLE_STATE::overdetermined);
        bool wo2=(!s2 || s2->state==VARIABLE_STATE::overdetermined);
        bool all_c=wc0 && wc1 && wc2;
        bool all_o=wo0 && wo1 && wo2;
        bool any_c=sc0 || sc1 || sc2;
        bool any_o=so0 || so1 || so2;
        
        // rhs is known; evaluate instruction
        if(all_c){
            if(d.state==VARIABLE_STATE::constant) return;
            d.state=VARIABLE_STATE::constant;
            d.value=Evaluate_Op(o.type,s0?s0->value:0,s1?s1->value:0,s2?s2->value:0);
            op_worklist.Append(it);
            return;}

        // Might be able to deduce result with partial information
        if(o.type.Can_Deduce()){
            // rhs is partially known; see if result can be deduced anyway
            if(any_c){
                T* a=(s0 && s0->state==VARIABLE_STATE::constant)?&s0->value:0;
                T* b=(s1 && s1->state==VARIABLE_STATE::constant)?&s1->value:0;
                T* c=(s2 && s2->state==VARIABLE_STATE::constant)?&s2->value:0;
                if(Deduce_Op(o.type,d.value,a,b,c)){
                    if(d.state==VARIABLE_STATE::constant) return;
                    d.state=VARIABLE_STATE::constant;
                    op_worklist.Append(it);
                    return;}
                else if(any_o){
                    d.state=VARIABLE_STATE::overdetermined;
                    op_worklist.Append(it);
                    return;}}
            if(all_o){
                d.state=VARIABLE_STATE::overdetermined;
                op_worklist.Append(it);
                return;}}
        else if(any_o){
            d.state=VARIABLE_STATE::overdetermined;
            op_worklist.Append(it);
            return;}
        return;}

    if(o.type==op_jmp){
        auto l=o.dest.Label();
        if(!block_exec(l)){
            block_exec(l)=true;
            block_worklist.Append(code_blocks(l));}
        return;}

    if(o.type==op_nop) return;
    // This leaves only op_brz and op_brnz
    auto d=o.dest.Label();
    auto s=o.src1.Label();

    VARIABLE_STATE s0=variable_state(o.src0);
    if(s0.state==VARIABLE_STATE::overdetermined){
        if(!block_exec(d)){
            block_exec(d)=true;
            block_worklist.Append(code_blocks(d));}
        if(!block_exec(s)){
            block_exec(s)=true;
            block_worklist.Append(code_blocks(s));}
        return;}

    if(s0.state!=VARIABLE_STATE::constant) return;

    if((o.type==op_brnz)!=(s0.value!=0)) d=s;
    if(!block_exec(d)){
        block_exec(d)=true;
        block_worklist.Append(code_blocks(d));}
}
//#####################################################################
// Function Sparse_Conditional_Constant_Propagation
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Sparse_Conditional_Constant_Propagation()
{
    DEBUG;
    REG_ARRAY<VARIABLE_STATE> variable_state;

    VARIABLE_STATE unknown={0,VARIABLE_STATE::unknown};
    VARIABLE_STATE overdetermined={0,VARIABLE_STATE::overdetermined};
    VARIABLE_STATE constant={0,VARIABLE_STATE::constant};
    variable_state.Resize(num_reg,unknown);
    variable_state.Resize(num_out,unknown);
    variable_state.Resize(num_in,overdetermined);
    variable_state.Resize(constants.m,constant);

    for(MEM_CONST i{};i<constants.m;i++)
        variable_state(i).value=constants(i);

    ARRAY<CODE_BLOCK*> block_worklist;
    ARRAY<INST_IT> op_worklist;
    ARRAY<bool,MEM_LABEL> block_exec(code_blocks.m);
    block_exec({})=true;
    block_worklist.Append(code_blocks({}));

    while(block_worklist.m || op_worklist.m){
        while(op_worklist.m){
            INST_IT op=op_worklist.Pop_Value();
            for(const auto& p:Get_Uses(*op->second)){
                if(block_exec(p.second->first->id))
                    SCCP_Visit_Instruction(p.second,variable_state,block_worklist,op_worklist,block_exec);}}

        while(block_worklist.m){
            CODE_BLOCK* block=block_worklist.Pop_Value();
            for(auto it=block->insts.begin();it!=block->insts.end();it++)
                SCCP_Visit_Instruction({block,it},variable_state,
                    block_worklist,op_worklist,block_exec);
            if(block->next[0] && (block->insts.empty() || !block->insts.back().type.Is_Jump()))
                if(!block_exec(block->next[0]->id)){
                    block_exec(block->next[0]->id)=true;
                    block_worklist.Append(block->next[0]);}}}

    for(int i=0;i<4;i++)
        for(int j=0;j<defs.a(i).m;j++)
            if(variable_state.a(i)(j).state==VARIABLE_STATE::constant)
                if(INST_IT it=defs.a(i)(j);it->first){
                    const INSTRUCTION& o=*it->second;
                    Change_Instruction(it,{op_copy,o.dest,Add_Constant(variable_state.a(i)(j).value)});
                    Propagate_Copy_And_Remove_Instruction(it);}

    for(MEM_LABEL i{};i<code_blocks.m;i++)
        if(CODE_BLOCK* B=code_blocks(i))
            if(!B->insts.empty()){
                INSTRUCTION& o=B->insts.back();
                if((o.type==op_brz || o.type==op_brz) && o.src0.Is_Const()){
                    if((o.type==op_brnz)==(constants(o.src0.Const())!=0)){
                        o.dest=o.src1;
                        exchange(B->next[0],B->next[1]);}
                    Change_Instruction(B->Last_It(),{});
                    Detach_Next(B,1);}}
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Propagate_Copy
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Propagate_Copy(REG old_var,REG new_var)
{
    assert(old_var!=new_var);
    HASHTABLE<const INSTRUCTION*,INST_IT>& ls=Get_Uses(old_var);
    while(ls.Size()){
        INST_IT it=ls.begin()->second;
        bool d0=it->second->type.Has_Src0() && it->second->src0==old_var;
        bool d1=it->second->type.Has_Src1() && it->second->src1==old_var;
        bool d2=it->second->type.Has_Src2() && it->second->src2==old_var;
        assert(d0 || d1 || d2);
        Remove_Use_Def(it);
        if(d0) it->second->src0=new_var;
        if(d1) it->second->src1=new_var;
        if(d2) it->second->src2=new_var;
        Add_Use_Def(it);}
}
//#####################################################################
// Function Copy_Propagation
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Copy_Propagation()
{
    DEBUG;
    // Copy propagation
    for(MEM_LABEL bl{};bl<code_blocks.m;bl++)
        if(CODE_BLOCK* B=code_blocks(bl))
            for(auto it=B->insts.begin();it!=B->insts.end();){
                auto it2=it++;
                if(it2->type==op_copy || (it2->type==op_phi && it2->src0==it2->src1))
                    Propagate_Copy_And_Remove_Instruction({B,it2});}
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Copy_Propagation
//#####################################################################
// Don't assume o.type==op_copy, but do assume it behaves as dest=src0
template<class T> void PROGRAM_COMPILER<T>::
Propagate_Copy_And_Remove_Instruction(INST_IT it)
{
    const INSTRUCTION& o=*it->second;
    Propagate_Copy(o.dest,o.src0);
    if(o.dest.Is_Reg()) Delete_Instruction(it);
    else Change_Instruction(it,{op_copy,o.dest,o.src0});
}
//#####################################################################
// Function Simplify_Phis
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Simplify_Phis()
{
    DEBUG;
    // Replace phi's with unreachble input by a copy.
    for(MEM_LABEL bl{};bl<code_blocks.m;bl++)
        if(CODE_BLOCK* B=code_blocks(bl)){
            if(!B->prev[0] || !B->prev[1]){
                int j=!B->prev[0];
                for(auto it=B->insts.begin();it!=B->insts.end();){
                    const INSTRUCTION& o=*it;
                    if(o.type!=op_phi) break;
                    auto it2=it++;
                    Change_Instruction({B,it2},{op_copy,o.dest,j?o.src1:o.src0});
                    Propagate_Copy_And_Remove_Instruction({B,it2});}
                B->prev[0]=B->prev[j];
                B->prev[1]=0;}
            else{
                for(auto it=B->insts.begin();it!=B->insts.end();){
                    const INSTRUCTION& o=*it;
                    auto it2=it++;
                    if(o.type!=op_phi) break;
                    if(o.src0!=o.src1) continue;
                    Propagate_Copy_And_Remove_Instruction({B,it2});}}}
    
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Remove_Dead_Code
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Remove_Dead_Code()
{
    DEBUG;
    // Prune unreachable basic blocks (bl=0 has no prev)
    ARRAY<CODE_BLOCK*> worklist;
    for(MEM_LABEL bl{1};bl<code_blocks.m;bl++)
        if(CODE_BLOCK* B=code_blocks(bl))
            if(!B->prev[0] && !B->prev[1])
                worklist.Append(B);
    for(int i=0;i<worklist.m;i++){
        CODE_BLOCK* B=worklist(i);
        for(int s=0;s<2;s++)
            if(CODE_BLOCK* C=B->next[s]){
                Detach_Next(B,s);
                if(!C->prev[0] && !C->prev[1])
                    worklist.Append(C);}}
    Simplify_Phis();
    for(int i=0;i<worklist.m;i++){
        CODE_BLOCK* B=worklist(i);
        while(!B->insts.empty()) Delete_Instruction({B,B->insts.begin()});
        code_blocks(B->id)=0;
        delete B;}

    // Merge basic blocks
    for(MEM_LABEL bl{};bl<code_blocks.m;bl++)
        if(CODE_BLOCK* B=code_blocks(bl))
            while(B->next[0] && !B->next[1] && !B->next[0]->prev[1]){
                CODE_BLOCK* nb=B->next[0];
                MEM_LABEL id=nb->id;
                for(int s=0;s<2;s++){
                    if(nb->next[s]){
                        int j=nb->next[s]->prev[1]==nb;
                        assert(nb->next[s]->prev[j]==nb);
                        nb->next[s]->prev[j]=B;}
                    B->next[s]=nb->next[s];
                    nb->next[s]=0;}
                assert(nb->prev[0]==B);
                nb->prev[0]=0;
                B->Merge_Code(nb);
                delete code_blocks(id);
                code_blocks(id)=0;}
    while(code_blocks.m && !code_blocks.Last())
        code_blocks.Pop();
    
    Compute_Use_Def();
    ARRAY<INST_IT> worklist2;
    for(int j=0;j<defs.a(mem_reg).m;j++)
        if(INST_IT it=defs.a(mem_reg)(j);it->first)
            if(!uses.a(0)(j).Size())
                worklist2.Append(it);

    while(worklist2.m){
        INST_IT it=worklist2.Pop_Value();
        REG r0=it->second->src0,r1=it->second->src1,r2=it->second->src2;
        if(it->second->type.Has_Src0() && r0.Is_Reg() && Get_Uses(r0).Size()==1) worklist2.Append(Get_Def(r0));
        if(it->second->type.Has_Src1() && r1.Is_Reg() && r1!=r0 && Get_Uses(r1).Size()==1) worklist2.Append(Get_Def(r1));
        if(it->second->type.Has_Src2() && r2.Is_Reg() && r2!=r0 && r2!=r1 && Get_Uses(r2).Size()==1) worklist2.Append(Get_Def(r2));
        Delete_Instruction(it);}

    for(MEM_LABEL bl{};bl<code_blocks.m;bl++)
        if(CODE_BLOCK* B=code_blocks(bl))
            for(auto it=B->insts.begin();it!=B->insts.end();){
                auto it2=it++;
                if(it2->type==op_nop)
                    Delete_Instruction({B,it2});}
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Detach_Next
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Detach_Next(CODE_BLOCK* B,int j)
{
    int phi_j=B->next[j]->prev[1]==B;
    assert(B->next[j]->prev[phi_j]==B);
    B->next[j]->prev[phi_j]=0;
    B->next[j]=0;
}
//#####################################################################
// Function Add_Constant
//#####################################################################
template<class T> MEM_CONST PROGRAM_COMPILER<T>::
Add_Constant(T x)
{
    typename std::map<T,MEM_CONST>::iterator it=constant_lookup.find(x);
    if(it!=constant_lookup.end()) return it->second;
    MEM_CONST i=constants.Append(x);
    constant_lookup[x]=i;
    Ensure_Valid_Use_Def(i);
    return i;
}
//#####################################################################
// Function Local_Common_Expresssion_Elimination
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Local_Common_Expresssion_Elimination(CODE_BLOCK* B)
{
    DEBUG;
    HASHTABLE<std::tuple<op_type,REG,REG,REG>,REG> expr;
    for(auto it=B->insts.begin();it!=B->insts.end();it++){
        auto& o=*it;
        if(!(o.type.Has_Dest())) continue;
        Reduce_In_Place({B,it});
        if(o.type==op_copy){
            Propagate_Copy(o.dest,o.src0);
            continue;}
        std::tuple<op_type,REG,REG,REG> in={o.type.op,o.src0,o.src1,o.src2};
        if(!(o.type.Has_Src1())) std::get<2>(in)={};
        if(!(o.type.Has_Src2())) std::get<3>(in)={};
        if(o.type.Commutes()) exchange_sort(std::get<1>(in),std::get<2>(in));
        REG e{};
        if(expr.Get(in,e)) Propagate_Copy(o.dest,e);
        else expr.Set(in,o.dest);}
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Local_Common_Expresssion_Elimination
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Local_Common_Expresssion_Elimination()
{
    DEBUG;
    for(MEM_LABEL bl{};bl<code_blocks.m;bl++){
        if(!code_blocks(bl)) continue;
        Local_Common_Expresssion_Elimination(code_blocks(bl));}
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Reduce_In_Place
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Reduce_In_Place(INST_IT it)
{
    REG const_0=MEM_CONST(0),const_1=MEM_CONST(1);
    auto& o=*it->second;
    switch(o.type.op){
        case op_add:
            if(o.src1==const_0){Change_Instruction(it,{op_copy,o.dest,o.src0});break;}
            if(o.src0==const_0){Change_Instruction(it,{op_copy,o.dest,o.src1});break;}
            break;
        case op_sub:
            if(o.src0==o.src1){Change_Instruction(it,{op_copy,o.dest,const_0});break;}
            if(o.src1==const_0){Change_Instruction(it,{op_copy,o.dest,o.src0});break;}
            if(o.src0==const_0){Change_Instruction(it,{op_neg,o.dest,o.src1});break;}
            break;
        case op_mul:
            if(o.src1==const_0 || o.src0==const_0){Change_Instruction(it,{op_copy,o.dest,const_0});break;}
            if(o.src1==const_1){Change_Instruction(it,{op_copy,o.dest,o.src0});break;}
            if(o.src0==const_1){Change_Instruction(it,{op_copy,o.dest,o.src1});break;}
            if(o.src1==Add_Constant(2)){Change_Instruction(it,{op_add,o.dest,o.src0,o.src0});break;}
            if(o.src0==Add_Constant(2)){Change_Instruction(it,{op_add,o.dest,o.src1,o.src1});break;}
            break;
        case op_div:
            if(o.src0==o.src1){Change_Instruction(it,{op_copy,o.dest,const_1});break;}
            if(o.src0==const_0 || o.src1==const_1){Change_Instruction(it,{op_copy,o.dest,o.src0});break;}
            if(o.src0==const_1){Change_Instruction(it,{op_inv,o.dest,o.src1});break;}
            break;
        case op_neg: break;
        case op_inv: break;
        case op_sqrt: break;
        case op_exp: break;
        case op_ln: break;
        case op_pow:
            if(o.src1==const_0){Change_Instruction(it,{op_copy,o.dest,const_1});break;}
            if(o.src1==const_1){Change_Instruction(it,{op_copy,o.dest,o.src0});break;}
            if(o.src1==Add_Constant(2)){Change_Instruction(it,{op_mul,o.dest,o.src0,o.src0});break;}
            if(o.src1==Add_Constant(-1)){Change_Instruction(it,{op_inv,o.dest,o.src0});break;}
            break;
        case op_sin: break;
        case op_cos: break;
        case op_asin: break;
        case op_acos: break;
        case op_atan: break;
        case op_atan2: break;
        case op_floor: break;
        case op_ceil: break;
        case op_rint: break;
        case op_min:
        case op_max:
            if(o.src0==o.src1) Change_Instruction(it,{op_copy,o.dest,o.src0});
            break;
        case op_abs: break;
        case op_if:
            if(o.src1==o.src2){Change_Instruction(it,{op_copy,o.dest,o.src1});break;}
            if(o.src0.Is_Const())
            {
                Change_Instruction(it,{op_copy,o.dest,o.src0!=const_0?o.src1:o.src2});
                break;
            }
            break;
        case op_lt:
        case op_gt:
        case op_ne:
            if(o.src0==o.src1){Change_Instruction(it,{op_copy,o.dest,const_0});break;}
            break;
        case op_le:
        case op_ge:
        case op_eq:
            if(o.src0==o.src1){Change_Instruction(it,{op_copy,o.dest,const_1});break;}
            break;
        case op_or:
        case op_and:
            if(o.src0==o.src1){Change_Instruction(it,{op_copy,o.dest,o.src0});break;}
            break;
        default: break;}
}
//#####################################################################
// Function Simplify_With_Distributive_Law_Helper_Prod
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Simplify_With_Distributive_Law_Helper_Prod(REG var,ARRAY<INST_IT>& extra_nodes,DISTRIBUTE_HELPER<T>& prod)
{
    if(var.Is_Const()){
        prod.coeff*=constants(var.Const());
        return;}

    if(var.Is_In() || var.Is_Out() || Get_Uses(var).Size()>1){
        prod.product.Append(var);
        return;}

    INST_IT it=Get_Def(var);
    assert(it->first);

    if(it->second->type==op_neg){
        extra_nodes.Append(it);
        prod.coeff=-prod.coeff;
        return Simplify_With_Distributive_Law_Helper_Prod(it->second->src0,extra_nodes,prod);}

    if(it->second->type!=op_mul){
        prod.product.Append(var);
        return;}

    extra_nodes.Append(it);
    Simplify_With_Distributive_Law_Helper_Prod(it->second->src0,extra_nodes,prod);
    Simplify_With_Distributive_Law_Helper_Prod(it->second->src1,extra_nodes,prod);
}
//#####################################################################
// Function Simplify_With_Distributive_Law_Helper
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Simplify_With_Distributive_Law_Helper(REG var,T coeff,ARRAY<INST_IT>& extra_nodes,ARRAY<DISTRIBUTE_HELPER<T> >& summands)
{
    if(var.Is_Const()){
        DISTRIBUTE_HELPER<T> prod={coeff*constants(var.Const())};
        summands.Append(prod);
        return;}

    if(extra_nodes.m && (var.Is_In() || var.Is_Out() || Get_Uses(var).Size()>1)){
        DISTRIBUTE_HELPER<T> prod={coeff};
        summands(summands.Append(prod)).product.Append(var);
        return;}

    INST_IT it=Get_Def(var);
    assert(it->first);

    if(it->second->type==op_mul){
        if(it->second->src0.Is_Const()){
            extra_nodes.Append(it);
            return Simplify_With_Distributive_Law_Helper(it->second->src1,coeff*constants(it->second->src0.Const()),extra_nodes,summands);}

        if(it->second->src1.Is_Const()){
            extra_nodes.Append(it);
            return Simplify_With_Distributive_Law_Helper(it->second->src0,coeff*constants(it->second->src1.Const()),extra_nodes,summands);}

        DISTRIBUTE_HELPER<T> prod={coeff};
        return Simplify_With_Distributive_Law_Helper_Prod(var,extra_nodes,summands(summands.Append(prod)));}

    if(it->second->type==op_neg){
        extra_nodes.Append(it);
        return Simplify_With_Distributive_Law_Helper(it->second->src0,-coeff,extra_nodes,summands);}

    if(it->second->type==op_copy){
        extra_nodes.Append(it);
        return Simplify_With_Distributive_Law_Helper(it->second->src0,coeff,extra_nodes,summands);}

    if(it->second->type!=op_add && it->second->type!=op_sub){
        DISTRIBUTE_HELPER<T> prod={coeff};
        summands(summands.Append(prod)).product.Append(var);
        return;}

    extra_nodes.Append(it);
    Simplify_With_Distributive_Law_Helper(it->second->src0,coeff,extra_nodes,summands);
    if(it->second->type==op_sub) coeff=-coeff;
    Simplify_With_Distributive_Law_Helper(it->second->src1,coeff,extra_nodes,summands);
}
//#####################################################################
// Function Simplify_With_Distributive_Law_Emit_Prod
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Simplify_With_Distributive_Law_Emit_Prod(DISTRIBUTE_HELPER<T>& prod,INST_IT before,bool& need_neg)
{
    need_neg=false;
    if(prod.coeff==0) return MEM_CONST(0);
    if(prod.product.m==0) return Add_Constant(prod.coeff);

    while(prod.product.m>1){
        REG src1=prod.product.Pop_Value();
        REG new_reg=Before_Instruction(before,op_mul,prod.product.Last(),src1);
        prod.product.Last()=new_reg;}

    if(prod.coeff==1) return prod.product(0);
    if(prod.coeff==-1){
        need_neg=true;
        return prod.product(0);}

    return Before_Instruction(before,op_mul,Add_Constant(prod.coeff),prod.product.Last());
}
//#####################################################################
// Function Simplify_With_Distributive_Law_Emit_Simple
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Simplify_With_Distributive_Law_Emit_Simple(ARRAY<DISTRIBUTE_HELPER<T> >& summands,INST_IT before,bool& need_neg)
{
    REG accum_reg{};
    need_neg=false;
    while(summands.m){
        bool neg_a=false;
        REG out_a=Simplify_With_Distributive_Law_Emit_Prod(summands.Last(),before,neg_a);
        summands.Pop();
        if(!accum_reg.Is_Value()){
            accum_reg=out_a;
            need_neg=neg_a;}
        else{
            REG new_reg;
            if(need_neg==neg_a) new_reg=Before_Instruction(before,op_add,accum_reg,out_a);
            else{
                if(need_neg) exchange(accum_reg,out_a);
                new_reg=Before_Instruction(before,op_sub,accum_reg,out_a);
                need_neg=false;}
            accum_reg=new_reg;}}
    return accum_reg;
}
//#####################################################################
// Function Simplify_With_Distributive_Law_Emit_Consts
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Simplify_With_Distributive_Law_Emit_Consts(ARRAY<DISTRIBUTE_HELPER<T> >& summands,INST_IT before,bool& need_neg)
{
    need_neg=false;
    T const_term=0;

    std::map<T,ARRAY<DISTRIBUTE_HELPER<T> > > consts;
    for(int i=0;i<summands.m;i++){
        if(summands(i).product.m==0) const_term+=summands(i).coeff;
        else consts[abs(summands(i).coeff)].Append(summands(i));}

    REG accum_reg{};

    ARRAY<DISTRIBUTE_HELPER<T> > no_coeff;
    typename std::map<T,ARRAY<DISTRIBUTE_HELPER<T> > >::iterator it;
    if((it=consts.find(0))!=consts.end()) consts.erase(it);
    if((it=consts.find(1))!=consts.end()){
        no_coeff.Exchange(it->second);
        consts.erase(it);}

    for(auto& it:consts){
        bool neg_a=false;
        for(int i=0;i<it.second.m;i++)
            it.second(i).coeff=sign(it.second(i).coeff);
        REG out_a=Simplify_With_Distributive_Law_Emit_Simple(it.second,before,neg_a);
        T c=it.first;
        if(neg_a) c=-c;
        REG new_reg=Before_Instruction(before,op_mul,Add_Constant(c),out_a);
        if(accum_reg.Is_Value())
            accum_reg=Before_Instruction(before,op_add,accum_reg,new_reg);
        else accum_reg=new_reg;}

    if(no_coeff.m){
        REG out_a=Simplify_With_Distributive_Law_Emit_Simple(no_coeff,before,need_neg);
        if(accum_reg.Is_Value()){
            accum_reg=Before_Instruction(before,need_neg?op_sub:op_add,accum_reg,out_a);
            need_neg=false;}
        else accum_reg=out_a;}

    if(!accum_reg.Is_Value()) return Add_Constant(const_term);

    if(const_term){
        accum_reg=Before_Instruction(before,need_neg?op_sub:op_add,Add_Constant(const_term),accum_reg);
        need_neg=false;}
    return accum_reg;
}
//#####################################################################
// Function Simplify_With_Distributive_Law_Emit
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Simplify_With_Distributive_Law_Emit(ARRAY<DISTRIBUTE_HELPER<T> >& summands,INST_IT before,bool& need_neg)
{
    HASHTABLE<REG,HASHTABLE<int> > has_var;
    for(int i=0;i<summands.m;i++)
        for(int j=0;j<summands(i).product.m;j++)
            has_var.Get_Or_Insert(summands(i).product(j)).Set(i);

    REG best_var{};
    int max_count=-1;
    for(const auto&p:has_var)
        if(p.second.Size()>max_count){
            max_count=p.second.Size();
            best_var=p.first;}

    if(max_count<=1){
        has_var.Remove_All();
        return Simplify_With_Distributive_Law_Emit_Consts(summands,before,need_neg);}

    int k=0;
    ARRAY<DISTRIBUTE_HELPER<T> > new_summands(max_count);
    HASHTABLE<int>& h=has_var.Get(best_var);
    for(auto n:h){
        new_summands(k).coeff=summands(n).coeff;
        new_summands(k++).product.Exchange(summands(n).product);
        summands(n).coeff=0;}
    for(int i=summands.m-1;i>=0;i--)
        if(summands(i).coeff==0){
            summands(i).coeff=summands.Last().coeff;
            summands(i).product.Exchange(summands.Last().product);
            summands.Pop();}
    for(int i=0;i<new_summands.m;i++)
        new_summands(i).product.Remove_Index_Lazy(new_summands(i).product.Find(best_var));
    has_var.Remove_All();
    bool neg_a=false,neg_b=false;
    REG out_a=Simplify_With_Distributive_Law_Emit(new_summands,before,neg_a);
    REG out_reg=Before_Instruction(before,op_mul,best_var,out_a);
    if(!summands.m){
        need_neg=neg_a;
        return out_reg;}
    REG out_b=Simplify_With_Distributive_Law_Emit(summands,before,neg_b);
    if(neg_a==neg_b){
        need_neg=neg_a;
        return Before_Instruction(before,op_add,out_reg,out_b);}
    need_neg=false;
    if(neg_b) return Before_Instruction(before,op_sub,out_reg,out_b);
    return Before_Instruction(before,op_sub,out_b,out_reg);
}
//#####################################################################
// Function Simplify_With_Distributive_Law
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Simplify_With_Distributive_Law(INST_IT it) -> INST_IT
{
    ARRAY<INST_IT> extra_nodes;
    ARRAY<DISTRIBUTE_HELPER<T> > summands;
    Simplify_With_Distributive_Law_Helper(it->second->dest,1,extra_nodes,summands);
    if(!extra_nodes.m) return --it;
    extra_nodes.Sort([](auto a,auto b){return &*a->second<&*b->second;});
    extra_nodes.Resize(std::unique(extra_nodes.begin(),extra_nodes.end())-extra_nodes.begin());
    for(int i=0;i<extra_nodes.m;i++)
        if(extra_nodes(i)!=it)
            Delete_Instruction(extra_nodes(i));
    bool need_neg=false;
    auto pit=it;
    --pit;
    REG out_reg=Simplify_With_Distributive_Law_Emit(summands,it,need_neg);
    if(need_neg) Change_Instruction(it,{op_neg,it->second->dest,out_reg});
    else{
        Change_Instruction(it,{op_copy,it->second->dest,out_reg});
        Propagate_Copy_And_Remove_Instruction(it);}
    return pit;
}
//#####################################################################
// Function Simplify_With_Distributive_Law
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Simplify_With_Distributive_Law()
{
    DEBUG;
    for(MEM_LABEL bl{};bl<code_blocks.m;bl++)
        if(CODE_BLOCK* B=code_blocks(bl))
            for(auto it=B->insts.end();it!=B->insts.begin();){
                it--;
                if(it->type!=op_neg && it->type!=op_add && it->type!=op_sub && it->type!=op_mul && it->type!=op_copy) continue;
                if(!it->dest.Is_In() && !it->dest.Is_Out() && Get_Uses(it->dest).Size()<=1) continue;
                it=Simplify_With_Distributive_Law({B,it})->second;
                it++;
            }
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Eliminate_Unused_Register
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Eliminate_Unused_Register(MEM_REG k)
{
    MEM_REG last{uses.a(mem_reg).m-1};
    if(k<last){
        if(defs(last)->first) defs(last)->second->dest=k;
        defs(k)=defs(last);
        Propagate_Copy(last,k);}
    defs.a(mem_reg).Pop();
    uses.a(mem_reg).Pop();
}
//#####################################################################
// Function Eliminate_Unused_Constant
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Eliminate_Unused_Constant(MEM_CONST k)
{
    if(Value(k)<=1) return;
    MEM_CONST last{uses.a(mem_const).m-1};
    constant_lookup.erase(constant_lookup.find(constants(k)));
    if(k<last){
        constant_lookup[constants(last)]=k;
        constants(k)=constants(last);
        Propagate_Copy(last,k);}
    constants.Pop();
    uses.a(mem_const).Pop();
}
//#####################################################################
// Function Compress_Registers
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Compress_Registers()
{
    DEBUG;
    HASHTABLE<REG,REG> h;
    num_reg={};
    h.Set(REG::Dont_Care(),REG::Dont_Care());
    auto remap=[this,&h](REG& r,bool dest=false)
    {
        if(!r.Is_Reg()) return;
        if(dest) h.Set(r,num_reg++);
        r=h.Get(r);
    };
    For_Each_Instruction(code_blocks,[this,remap](INSTRUCTION& in)
    {
        if(in.type.Has_Dest()) remap(in.dest,true);
        if(in.type.Has_Src0()) remap(in.src0);
        if(in.type.Has_Src1()) remap(in.src1);
        if(in.type.Has_Src2()) remap(in.src2);
        return true;
    });
    Compute_Use_Def();
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Remove_Use_Def
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Remove_Use_Def(INST_IT it)
{
    if(!use_def_valid) return;
    if(it->second->type.Has_Dest()) Set_Def(it->second->dest,{});
    if(it->second->type.Has_Src0()) Get_Uses(it->second->src0).Delete_If_Present(&*it->second);
    if(it->second->type.Has_Src1()) Get_Uses(it->second->src1).Delete_If_Present(&*it->second);
    if(it->second->type.Has_Src2()) Get_Uses(it->second->src2).Delete_If_Present(&*it->second);
}
//#####################################################################
// Function Add_Use_Def
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Add_Use_Def(INST_IT it)
{
    if(!use_def_valid) return;
    auto& in=*it->second;
    if(in.type.Has_Dest())
    {
        defs.Get_Resize(in.dest)=it;
        uses.Get_Resize(in.dest);
    }
    if(in.type.Has_Src0()) uses.Get_Resize(in.src0).Set(&*it->second,it);
    if(in.type.Has_Src1()) uses.Get_Resize(in.src1).Set(&*it->second,it);
    if(in.type.Has_Src2()) uses.Get_Resize(in.src2).Set(&*it->second,it);
}
//#####################################################################
// Function Delete_Instruction
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Delete_Instruction(INST_IT it)
{
    Remove_Use_Def(it);
    it->first->insts.erase(it->second);
}
//#####################################################################
// Function Change_Instruction
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Change_Instruction(INST_IT it,const INSTRUCTION& in)
{
    Remove_Use_Def(it);
    *it->second=in;
    Add_Use_Def(it);
}
//#####################################################################
// Function Shift_Registers
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Shift_Registers()
{
    DEBUG;
    MEM_REG shift=num_reg;
    num_reg+=Value(shift);
    auto sh=[shift](REG& v){if(v.Is_Reg()) v=v.Reg()+Value(shift);};

    for(auto B:code_blocks)
    {
        if(!B) continue;
        for(auto o:*B)
        {
            if(o.type.Has_Dest()) sh(o.dest);
            if(o.type.Has_Src0()) sh(o.src0);
            if(o.type.Has_Src1()) sh(o.src1);
            if(o.type.Has_Src2()) sh(o.src2);
        }
    }

    for(auto&h:dict)
        for(auto&a:h.second->x)
            sh(a);
    Compute_Use_Def();
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Rename_Variable
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Rename_Variable(REG var_old,REG var_new)
{
    if(var_old.Is_Out() || var_old.Is_Reg())
    {
        auto it=Get_Def(var_old);
        Remove_Use_Def(it);
        assert(it->second->dest==var_old);
        it->second->dest=var_new;
        Add_Use_Def(it);
    }

    ARRAY<INST_IT> values;
    Get_Uses(var_old).Get_Data(values);
    for(auto it:values)
    {
        auto&o=*it->second;
        Remove_Use_Def(it);
        if((o.type.Has_Src0()) && o.src0==var_old)
            o.src0=var_new;
        if((o.type.Has_Src1()) && o.src1==var_old)
            o.src1=var_new;
        if((o.type.Has_Src2()) && o.src2==var_old)
            o.src2=var_new;
        Add_Use_Def(it);
    }
}
//#####################################################################
// Function Append
//#####################################################################
template<class T> MEM_OUT PROGRAM_COMPILER<T>::
Append(const PROGRAM_COMPILER<T>& pc)
{
    DEBUG;
    PHYSBAM_ASSERT(pc.num_in<=num_in);
    HASHTABLE<REG,REG> var_map;

    for(MEM_REG i{};i<pc.num_reg;i++)
        var_map.Set(i,i+Value(num_reg));
    num_reg+=Value(pc.num_reg);

    for(MEM_CONST i{};i<pc.constants.m;i++)
        var_map.Set(i,Add_Constant(pc.constants(i)));

    for(MEM_IN i{};i<pc.num_in;i++)
        var_map.Set(i,i);

    MEM_OUT pc_num_out=pc.num_out;
    MEM_OUT num_out2=num_out;
    for(MEM_OUT i{};i<pc_num_out;i++)
        var_map.Set(i,i+Value(num_out));
    num_out+=Value(pc_num_out);
    
    auto sh=[&var_map](REG& v) {v=var_map.Get(v);};

    MEM_LABEL num_blocks=code_blocks.m;
    for(auto B:pc.code_blocks)
    {
        if(!B)
        {
            code_blocks.Append(0);
            continue;
        }
        auto C=new CODE_BLOCK(*B);
        C->id=code_blocks.Append(C);
        for(auto it=C->insts.begin();it!=C->insts.end();it++)
        {
            if(it->type.Has_Dest()) sh(it->dest);
            if(it->type.Has_Src0()) sh(it->src0);
            if(it->type.Has_Src1()) sh(it->src1);
            if(it->type.Has_Src2()) sh(it->src2);
            Add_Use_Def({C,it});
        }
    }
    for(MEM_LABEL i{0};i<pc.code_blocks.m;i++)
    {
        auto A=pc.code_blocks(i);
        if(!A) continue;
        auto B=code_blocks(i+Value(num_blocks));
        for(int j=0;j<2;j++)
        {
            if(A->next[j]) B->next[j]=code_blocks(A->next[j]->id+Value(num_blocks));
            if(A->prev[j]) B->prev[j]=code_blocks(A->prev[j]->id+Value(num_blocks));
        }
    }
    if(num_blocks && pc.code_blocks.m)
    {
        code_blocks(num_blocks-1)->next[0]=code_blocks(num_blocks);
        code_blocks(num_blocks)->prev[0]=code_blocks(num_blocks-1);
    }

    VERBOSE;
    DEBUG;
    return num_out2;
}
//#####################################################################
// Function Make_Output
//#####################################################################
template<class T> MEM_OUT PROGRAM_COMPILER<T>::
Make_Output(REG var)
{
    if(var.Is_Out()) return var.Out();
    if(debug) Test_Invariants();
    auto dest=New_Out();
    Append_Instruction_Dest(op_copy,dest,var);
    if(debug) Test_Invariants();
    return dest;
}
//#####################################################################
// Function Ensure_Valid_Use_Def
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Ensure_Valid_Use_Def(REG var)
{
    uses.Get_Resize(var);
    if(var.Is_Reg() || var.Is_Out()) defs.Get_Resize(var);
}
//#####################################################################
// Function Move_Outputs_To_Reg
//#####################################################################
template<class T> std::pair<MEM_REG,int> PROGRAM_COMPILER<T>::
Move_Outputs_To_Reg()
{
    DEBUG;
    std::pair<MEM_REG,int> r(num_reg,Value(num_out));
    Ensure_Valid_Use_Def(r.first+r.second-1);
    for(MEM_OUT i{};i<num_out;i++)
        Rename_Variable(i,New_Reg());
    var_out.Remove_All();
    num_out={};
    uses.a(mem_out).Remove_All();
    defs.a(mem_out).Remove_All();
    VERBOSE;
    DEBUG;
    return r;
}
//#####################################################################
// Function Test_Invariants
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Test_Invariants() const
{
    for(MEM_LABEL i{};i<code_blocks.m;i++)
        if(auto B=code_blocks(i))
        {
            for(int s=0;s<2;s++)
            {
                if(B->next[s])
                {
                    PHYSBAM_ASSERT(B->next[s]->prev[0]==B||B->next[s]->prev[1]==B);
                }
                if(B->prev[s])
                {
                    PHYSBAM_ASSERT(B->prev[s]->next[0]==B||B->prev[s]->next[1]==B);
                }
            }
            if(B->next[1])
            {
                PHYSBAM_ASSERT(B->insts.size());
                auto& o=B->insts.back();
                PHYSBAM_ASSERT(o.type.Is_Branch());
                PHYSBAM_ASSERT(o.dest.Label()==B->next[1]->id);
                PHYSBAM_ASSERT(o.src1.Label()==B->next[0]->id);
            }
        }
    std::set<REG> already_computed;
    auto is_valid_reg=[this,&already_computed](REG r,bool dest=false)
    {
        if(dest) already_computed.insert(r);
        else if(r.Is_Reg() || r.Is_Out())
            PHYSBAM_ASSERT(already_computed.contains(r));
        if(r.Is_Reg()){PHYSBAM_ASSERT(r.Reg()>=MEM_REG{} && r.Reg()<num_reg);}
        else if(r.Is_Const()){PHYSBAM_ASSERT(r.Const()>=MEM_CONST{} && r.Const()<constants.m);}
        else if(r.Is_In()){PHYSBAM_ASSERT(r.In()>=MEM_IN{} && r.In()<num_in);}
        else if(r.Is_Out()){PHYSBAM_ASSERT(r.Out()>=MEM_OUT{} && r.Out()<num_out);}
        else
        {
            auto [t,v]=r.Parts();
            PHYSBAM_ASSERT(t==5);
            PHYSBAM_ASSERT(v>=0 && v<2);
        }
    };
    For_Each_Instruction(code_blocks,[this,is_valid_reg](const INSTRUCTION& in)
    {
        if(in.type.Has_Dest()) is_valid_reg(in.dest,true);
        if(in.type.Has_Src0()) is_valid_reg(in.src0);
        if(in.type.Has_Src1()) is_valid_reg(in.src1);
        if(in.type.Has_Src2()) is_valid_reg(in.src2);
        return true;
    });
    if(use_def_valid)
    {
        For_Each_Instruction(code_blocks,[this](const INSTRUCTION& in)
        {
            if(in.type.Has_Dest())
            {
                auto it=Get_Def(in.dest);
                PHYSBAM_ASSERT(it->first);
                PHYSBAM_ASSERT(&*it->second==&in);
            }
            if(in.type.Has_Src0())
            {
                auto& u=Get_Uses(in.src0);
                PHYSBAM_ASSERT(u.Contains(&in));
                if(in.src0.Is_Reg() || in.src0.Is_Out())
                    PHYSBAM_ASSERT(Get_Def(in.src0)->first);
            }
            if(in.type.Has_Src1())
            {
                auto& u=Get_Uses(in.src1);
                PHYSBAM_ASSERT(u.Contains(&in));
                if(in.src1.Is_Reg() || in.src1.Is_Out())
                    PHYSBAM_ASSERT(Get_Def(in.src1)->first);
            }
            if(in.type.Has_Src2())
            {
                auto& u=Get_Uses(in.src2);
                PHYSBAM_ASSERT(u.Contains(&in));
                if(in.src2.Is_Reg() || in.src2.Is_Out())
                    PHYSBAM_ASSERT(Get_Def(in.src2)->first);
            }
            return true;
        });
    }
    bool phi_section=true;
    CODE_BLOCK* prev_block=0;
    For_Each_Iterator(code_blocks,
        [this,&prev_block,&phi_section](const auto& it)
        {
            if(it->first!=prev_block)
            {
                prev_block=it->first;
                phi_section=true;
            }
            auto& o=*it->second;
            if(o.type==op_phi)
            {
                PHYSBAM_ASSERT(phi_section);
            }
            else phi_section=false;
            return true;
        });
}
//#####################################################################
// Function Func_Lookup
//#####################################################################
template<class T> HASHTABLE<std::string,std::pair<COMPILER_DATA<T>(*)(PROGRAM_COMPILER<T>*pc,const ARRAY<COMPILER_DATA<T>>& in),int>>& PROGRAM_COMPILER<T>::
Func_Lookup()
{
    static HASHTABLE<std::string,std::pair<COMPILER_DATA<T>(*)(PROGRAM_COMPILER<T>*pc,const ARRAY<COMPILER_DATA<T>>& in),int>> h;
    return h;
}
//#####################################################################
// Function Clear_Neg
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Clear_Neg(REG_AUG ra) -> REG_AUG
{
    if(!ra.neg) return ra;
    if(ra.reg.Is_Const()) return {Add_Constant(-constants(ra.reg.Const())),false,ra.inv};
    auto it=ra.reg.Is_In()?Begin_It():Get_Def(ra.reg);
    if(!ra.reg.Is_In()) ++it;
    Move_Past_Phis(it);
    ra.reg=Before_Instruction(it,op_neg,ra.reg);
    ra.neg=false;
    return ra;
}
//#####################################################################
// Function Clear_Inv
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Clear_Inv(REG_AUG ra) -> REG_AUG
{
    if(!ra.inv) return ra;
    if(ra.reg.Is_Const()) return {Add_Constant(1/constants(ra.reg.Const())),ra.neg};
    auto it=ra.reg.Is_In()?Begin_It():Get_Def(ra.reg);
    if(!ra.reg.Is_In()) ++it;
    Move_Past_Phis(it);
    ra.reg=Before_Instruction(it,op_inv,ra.reg);
    ra.inv=false;
    return ra;
}
//#####################################################################
// Function Mul
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Mul(INST_IT it,REG_AUG a,REG_AUG b) -> REG_AUG
{
    MEM_CONST zero=MEM_CONST(0);
    MEM_CONST one=MEM_CONST(1);
    if(a.reg==zero || b.reg==zero) return {zero};
    if(a.reg==one) return a.neg?Neg(b):b;
    if(b.reg==one) return b.neg?Neg(a):a;
    bool neg=(a.neg!=b.neg);
    if(a.inv==b.inv)
        return {Before_Instruction(it,op_mul,a.reg,b.reg),neg,a.inv};
    if(a.inv) std::swap(a,b);
    return {Before_Instruction(it,op_div,a.reg,b.reg),neg};
}
//#####################################################################
// Function Add
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Add_Mul(INST_IT it,REG_AUG& o,REG_AUG a,REG_AUG b)
{
    o=Add(it,o,Mul(it,a,b));
}
//#####################################################################
// Function Add
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Add(INST_IT it,REG_AUG a,REG_AUG b) -> REG_AUG
{
    MEM_CONST zero=MEM_CONST(0);
    if(a.reg==zero) return b;
    if(b.reg==zero) return a;
    a=Clear_Inv(a);
    b=Clear_Inv(b);
    if(a.neg==b.neg) return {Before_Instruction(it,op_add,a.reg,b.reg),a.neg};
    if(a.neg) std::swap(a,b);
    return {Before_Instruction(it,op_sub,a.reg,b.reg)};
}
//#####################################################################
// Function If_Aug
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
If_Aug(INST_IT it,REG a,REG_AUG b,REG_AUG c) -> REG_AUG
{
    Make_Same_Aug(b,c);
    b.reg=Before_Instruction(it,op_if,a,b.reg,c.reg);
    return b;
}
//#####################################################################
// Function Prepend_Instruction
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Prepend_Instruction(CODE_BLOCK* B,const OPERATION& type,REG src0,REG src1,REG src2)
{
    return Before_Instruction(B->Begin_It(),type,src0,src1,src2);
}
//#####################################################################
// Function Prepend_Phi
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Prepend_Phi(CODE_BLOCK* B,REG src0,REG src1)
{
    assert(B->prev[0] && B->prev[1]);
    return Before_Instruction(B->Begin_It(),op_phi,src0,src1);
}
//#####################################################################
// Function Prepend_Phi
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Prepend_Phi(CODE_BLOCK* B,REG_AUG src0,REG_AUG src1) -> REG_AUG
{
    Make_Same_Aug(src0,src1);
    src0.reg=Prepend_Phi(B,src0.reg,src1.reg);
    return src0;
}
//#####################################################################
// Function Append_Instruction
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Append_Instruction(CODE_BLOCK* B,const OPERATION& type,REG src0,REG src1,REG src2)
{
    return Before_Instruction(B->End_It(),type,src0,src1,src2);
}
//#####################################################################
// Function Clear_Aug
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Clear_Aug(REG_AUG ra)
{
    return Clear_Neg(Clear_Inv(ra)).reg;
}
//#####################################################################
// Function Begin_It
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Begin_It() const -> INST_IT
{
    return code_blocks({})->Begin_It();
}
//#####################################################################
// Function Last_It
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Last_It() const -> INST_IT
{
    return code_blocks.Last()->Last_It();
}
//#####################################################################
// Function End_It
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
End_It() const -> INST_IT
{
    return code_blocks.Last()->End_It();
}
//#####################################################################
// Function Get_Def
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Get_Def(REG var) const -> INST_IT
{
    return defs(var);
}
//#####################################################################
// Function Set_Def
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Set_Def(REG var,INST_IT def)
{
    defs(var)=def;
}
//#####################################################################
// Function Get_Uses
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Get_Uses(REG var) const
    -> const HASHTABLE<const INSTRUCTION*,INST_IT>&
{
    return uses(var);
}
//#####################################################################
// Function Get_Uses
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Get_Uses(REG var)
    -> HASHTABLE<const INSTRUCTION*,INST_IT>&
{
    return uses(var);
}
//#####################################################################
// Function Get_Uses
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Get_Uses(const INSTRUCTION& in) const
    -> const HASHTABLE<const INSTRUCTION*,INST_IT>&
{
    return Get_Uses(in.dest);
}
//#####################################################################
// Function Clear_Diff_Structures
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Clear_Diff_Structures()
{
    diffs.Remove_All();
    adjoints.Remove_All();
    diff_block_end.Remove_All();
    diff_parent.Remove_All();
    reversal_points.Remove_All();
    parent_reg.Remove_All();
    trace_reg.Remove_All();
}
//#####################################################################
// Function Add_Reversal_Point
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Add_Reversal_Point()
{
    DEBUG;

    MEM_CONST zero=MEM_CONST(0);
    MEM_CONST one=MEM_CONST(1);
    MEM_LABEL num_blocks=code_blocks.m;
    
    struct cond_info
    {
        std::map<REG,bool> path;
        REG cond{};
        bool cond_val=false;
    };
    ARRAY<cond_info,MEM_LABEL> conditions(num_blocks);
    ARRAY<CODE_BLOCK*,MEM_LABEL> new_blocks(num_blocks);
    for(MEM_LABEL i{};i<num_blocks;i++)
        if(auto B=code_blocks(i))
        {
            new_blocks(i)=new CODE_BLOCK;
            auto& ir=conditions(i);
            if(auto R=B->prev[0])
            {
                auto ia=conditions(R->id);
                if(R->next[1])
                {
                    auto& o=R->insts.back();
                    bool br=R->next[1]==B;
                    ia.path[o.src0]=(br==(o.type==op_brnz));
                }
                if(auto S=B->prev[1])
                {
                    auto ib=conditions(B->prev[1]->id);
                    if(S->next[1])
                    {
                        auto& o=S->insts.back();
                        bool br=R->next[1]==B;
                        ib.path[o.src0]=(br==(o.type==op_brnz));
                    }
                    auto r=Merge_Conditions(ir.path,ia.path,ib.path);
                    if(r.first.Is_Value())
                    {
                        ir.cond=r.first;
                        ir.cond_val=r.second;
                    }
                    else
                    {
                        ir.cond=Prepend_Phi(B,zero,one),LOG::printf("line %P\n",__LINE__);
                        ir.cond_val=false;
                    }
                }
                else ir=ia;
            }
        }

    for(MEM_LABEL i=num_blocks-1;i>=MEM_LABEL{};i--)
        if(auto B=code_blocks(i))
        {
            auto C=new_blocks(i);
            for(int j=0;j<2;j++)
            {
                if(auto N=B->next[j])
                    C->prev[j]=new_blocks(N->id);
                if(auto N=B->prev[j])
                    C->next[j]=new_blocks(N->id);
            }
            C->id=code_blocks.Append(C);
        }
    code_blocks(new_blocks.m)->prev[0]=code_blocks(new_blocks.m-1);
    code_blocks(new_blocks.m-1)->next[0]=code_blocks(new_blocks.m);

    for(MEM_LABEL i=num_blocks-1;i>=MEM_LABEL{};i--)
        if(auto B=code_blocks(i))
        {
            auto C=new_blocks(i);
            if(B->prev[0] && B->prev[1])
            {
                MEM_LABEL id_0=B->prev[0]->id,o0(2*Value(num_blocks)-Value(id_0)-1);
                MEM_LABEL id_1=B->prev[1]->id,o1(2*Value(num_blocks)-Value(id_1)-1);
                auto& c=conditions(B->id);
                auto cond_reg=Trace_Register(c.cond,num_blocks);
                C->insts.push_back({c.cond_val?op_brnz:op_brz,o1,cond_reg,o0});
                Add_Use_Def(C->Last_It());
            }
        }

    reversal_points.Append(num_blocks);
    VERBOSE;
    DEBUG;
}
//#####################################################################
// Function Get_Reversal_point
//#####################################################################
template<class T> MEM_LABEL PROGRAM_COMPILER<T>::
Get_Reversal_Point(MEM_LABEL id)
{
    int r=reversal_points.Binary_Search(id);
    if(r==reversal_points.m) Add_Reversal_Point();
    return reversal_points(r);
}
//#####################################################################
// Function Trace_Register
//#####################################################################
template<class T> REG PROGRAM_COMPILER<T>::
Trace_Register(REG reg,MEM_LABEL id)
{
    if(!reg.Is_Reg()) return reg;
    MEM_LABEL start=Get_Def(reg)->first->id;
    if(start==id) return reg;
    if(start>id) return REG{};
    if(auto* pa=parent_reg.Get_Pointer(reg))
        return Trace_Register(*pa,id);
    auto& tr=trace_reg.Get_Or_Insert(reg);

    auto get=[&tr](MEM_LABEL id)
    {
        auto i=std::partition_point(tr.mapping.begin(),tr.mapping.end(),
            [id](const auto& a){return a.first<=id;});
        if(i==tr.mapping.begin()) return REG{};
        return (--i)->second;
    };
    if(!tr.mapping.m)
    {
        tr.mapping.Append({start,reg});
        for(MEM_LABEL b=start+1;b<code_blocks.m;b++)
        {
            auto B=code_blocks(b);
            PHYSBAM_ASSERT(B->prev[0]);
            REG ra=get(B->prev[0]->id);
            if(B->prev[1])
            {
                REG rb=get(B->prev[1]->id);
                if(ra!=rb)
                {
                    ra=Prepend_Phi(B,ra==REG{}?REG::Dont_Care():ra,
                        rb==REG{}?REG::Dont_Care():rb);
                    parent_reg.Set(ra,reg);
                }
            }
            if(tr.mapping.Last().second!=ra)
                tr.mapping.Append({b,ra});
        }
    }
    return get(id);
}
//#####################################################################
// Function Trace_Register
//#####################################################################
template<class T> auto PROGRAM_COMPILER<T>::
Trace_Register(REG_AUG reg,MEM_LABEL id) -> REG_AUG
{
    reg.reg=Trace_Register(reg.reg,id);
    return reg;
}
//#####################################################################
// Function Make_Same_Neg
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Make_Same_Neg(REG_AUG& ra,REG_AUG& rb)
{
    if(ra.reg==MEM_CONST(0) || ra.reg.Is_Dont_Care()) ra.neg=rb.neg;
    else if(rb.reg==MEM_CONST(0) || rb.reg.Is_Dont_Care()) rb.neg=ra.neg;
    else if(ra.neg!=rb.neg)
    {
        ra=Clear_Neg(ra);
        rb=Clear_Neg(rb);
    }
}
//#####################################################################
// Function Make_Same_Inv
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Make_Same_Inv(REG_AUG& ra,REG_AUG& rb)
{
    if(ra.reg==MEM_CONST(0)) rb=Clear_Inv(rb);
    else if(rb.reg==MEM_CONST(0)) ra=Clear_Inv(ra);
    else if(ra.reg==MEM_CONST(1) || ra.reg.Is_Dont_Care()) ra.inv=rb.inv;
    else if(rb.reg==MEM_CONST(1) || rb.reg.Is_Dont_Care()) rb.inv=ra.inv;
    else if(ra.inv!=rb.inv)
    {
        ra=Clear_Inv(ra);
        rb=Clear_Inv(rb);
    }
}
//#####################################################################
// Function Make_Same_Aug
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Make_Same_Aug(REG_AUG& ra,REG_AUG& rb)
{
    Make_Same_Inv(ra,rb);
    Make_Same_Neg(ra,rb);
}
//#####################################################################
// Function Move_Past_Phis
//#####################################################################
template<class T> void PROGRAM_COMPILER<T>::
Move_Past_Phis(INST_IT& it) const
{
    auto end=it->first->End_It();
    while(it!=end && it->second->type==op_phi)
        it++;
}
//#####################################################################
// Function Append_If_Else
//#####################################################################
template<class T> ARRAY<REG> PROGRAM_COMPILER<T>::
Append_If_Else(REG cond,std::function<ARRAY<REG>()> make_true,
    std::function<ARRAY<REG>()> make_false)
{
    CODE_BLOCK* A=code_blocks.Last();

    CODE_BLOCK* B=new CODE_BLOCK;
    B->id=code_blocks.Append(B);
    A->next[0]=B;
    B->prev[0]=A;
    ARRAY<REG> b_regs=make_true();
    CODE_BLOCK* BB=code_blocks.Last();

    CODE_BLOCK* C=new CODE_BLOCK;
    C->id=code_blocks.Append(C);
    A->next[1]=C;
    C->prev[0]=A;
    ARRAY<REG> c_regs=make_false();
    CODE_BLOCK* CC=code_blocks.Last();

    CODE_BLOCK* D=new CODE_BLOCK;
    D->id=code_blocks.Append(D);
    D->prev[0]=BB;
    D->prev[1]=CC;
    BB->next[0]=D;
    CC->next[0]=D;

    Append_Instruction_Dest(A,op_brz,C->id,cond,B->id);

    PHYSBAM_ASSERT(b_regs.m==c_regs.m);
    for(int i=0;i<b_regs.m;i++)
        b_regs(i)=Prepend_Phi(D,b_regs(i),c_regs(i));
    return b_regs;
}
//#####################################################################
// Function Add_Input
//#####################################################################
template<class T> MEM_IN PROGRAM_COMPILER<T>::
Add_Input(const std::string& name,MEM_IN in)
{
    if(in==MEM_IN(-1)) in=num_in++;
    var_in.Append({name,REG(in)});
    return in;
}
//#####################################################################
// Function Add_Input
//#####################################################################
template<class T> MEM_IN PROGRAM_COMPILER<T>::
Add_Input(const std::string& name,int m,MEM_IN in)
{
    auto& p=var_in(var_in.Add_End());
    p.first=name;
    auto&d=p.second;
    d.type=DATA::TYPE::vec;
    d.m=m;
    if(in==MEM_IN(-1))
    {
        in=num_in;
        num_in+=m;
    }
    auto reg=in;
    if(!m) return in;
    d.x.Resize(m);
    for(int i=0;i<m;i++) d.x(i)=reg++;
    return in;
}
//#####################################################################
// Function Add_Input
//#####################################################################
template<class T> MEM_IN PROGRAM_COMPILER<T>::
Add_Input(const std::string& name,int m,int n,MEM_IN in)
{
    auto& p=var_in(var_in.Add_End());
    p.first=name;
    auto&d=p.second;
    d.type=DATA::TYPE::mat;
    d.m=m;
    d.n=n;
    if(in==MEM_IN(-1))
    {
        in=num_in;
        num_in+=m*n;
    }
    auto reg=in;
    if(!m || !n) return in;
    d.x.Resize(m*n);
    for(int i=0;i<m*n;i++) d.x(i)=reg++;
    return in;
}
//#####################################################################
// Function Add_Output
//#####################################################################
template<class T> MEM_OUT PROGRAM_COMPILER<T>::
Add_Output(const std::string& name,MEM_OUT out)
{
    if(out==MEM_OUT(-1)) out=num_out++;
    var_out.Append({name,REG(out)});
    return out;
}
//#####################################################################
// Function Add_Output
//#####################################################################
template<class T> MEM_OUT PROGRAM_COMPILER<T>::
Add_Output(const std::string& name,int m,MEM_OUT out)
{
    auto& p=var_out(var_out.Add_End());
    p.first=name;
    auto&d=p.second;
    d.type=DATA::TYPE::vec;
    d.m=m;
    if(out==MEM_OUT(-1))
    {
        out=num_out;
        num_out+=m;
    }
    auto reg=out;
    if(!m) return out;
    d.x.Resize(m);
    for(int i=0;i<m;i++) d.x(i)=reg++;
    return out;
}
//#####################################################################
// Function Add_Output
//#####################################################################
template<class T> MEM_OUT PROGRAM_COMPILER<T>::
Add_Output(const std::string& name,int m,int n,MEM_OUT out)
{
    auto& p=var_out(var_out.Add_End());
    p.first=name;
    auto&d=p.second;
    d.type=DATA::TYPE::mat;
    d.m=m;
    d.n=n;
    if(out==MEM_OUT(-1))
    {
        out=num_out;
        num_out+=m*n;
    }
    auto reg=out;
    if(!m || !n) return out;
    d.x.Resize(m*n);
    for(int i=0;i<m*n;i++) d.x(i)=reg++;
    return out;
}
//#####################################################################
// Function Get_Const
//#####################################################################
template<class T> T PROGRAM_COMPILER<T>::
Get_Const(const DATA& a) const
{
    PHYSBAM_ASSERT(a.type==COMPILER_DATA<T>::TYPE::scalar);
    PHYSBAM_ASSERT(a.x(0).Is_Const());
    return constants(a.x(0).Const());
}
template struct PROGRAM_COMPILER<float>;
template struct PROGRAM_COMPILER<double>;
}
