//#####################################################################
// Copyright 2012.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYREGHT.txt.
//#####################################################################
#include <Core/Symbolics/PROGRAM_SCALAR.h>
#include <Core/Log/LOG.h>
#include <Core/Symbolics/CODE_BLOCK.h>
#include <Core/Symbolics/INSTRUCTION.h>
#include <Core/Symbolics/PROGRAM_COMPILER.h>
#include <Core/Vectors/VECTOR.h>

namespace PhysBAM{
namespace NS_PROGRAM_SCALAR
{
template<class T> PROGRAM_SCALAR<T> PROGRAM_SCALAR<T>::
Number(PROGRAM_COMPILER<T>& pc,T a)
{
    return {&pc,pc.Add_Constant(a)};
}

template<class T> PROGRAM_SCALAR<T>
operator+(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_copy,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator-(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_neg,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator+(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_add,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator-(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_sub,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator*(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_mul,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator/(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_div,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator<(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_lt,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator>(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_gt,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator<=(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_le,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator>=(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_ge,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator==(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_eq,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
operator!=(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_ne,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
cond(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b,const PROGRAM_SCALAR<T>& c)
{
    return {a.pc,a.pc->Append_Instruction(op_if,a.var,b.var,c.var)};
}

template<class T> ARRAY<PROGRAM_SCALAR<T>>
cond_lazy(const PROGRAM_SCALAR<T>& a,std::function<ARRAY<PROGRAM_SCALAR<ND<T>>>()> b,
    std::function<ARRAY<PROGRAM_SCALAR<ND<T>>>()> c)
{
    auto rb=[&b]()
    {
        auto s=b();
        ARRAY<REG> r(s.m);
        for(int i=0;i<s.m;i++) r(i)=s(i).var;
        return r;
    };
    auto rc=[&c]()
    {
        auto s=c();
        ARRAY<REG> r(s.m);
        for(int i=0;i<s.m;i++) r(i)=s(i).var;
        return r;
    };

    auto ar=a.pc->Append_If_Else(a.var,rb,rc);
    ARRAY<PROGRAM_SCALAR<T>> ap(ar.m);
    for(int i=0;i<ar.m;i++) ap(i)={a.pc,ar(i)};
    return ap;
}

template<class T> PROGRAM_SCALAR<T>
cond_lazy(const PROGRAM_SCALAR<T>& a,
    std::function<PROGRAM_SCALAR<ND<T>>()> b,
    std::function<PROGRAM_SCALAR<ND<T>>()> c)
{
    auto rb=[&b](){return ARRAY<REG>{b().var};};
    auto rc=[&c](){return ARRAY<REG>{c().var};};
    return {a.pc,a.pc->Append_If_Else(a.var,rb,rc)(0)};
}

template<class T> PROGRAM_SCALAR<T>
sqrt(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_sqrt,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
exp(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_exp,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
log(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_ln,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
sin(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_sin,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
cos(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_cos,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
tan(const PROGRAM_SCALAR<T>& a)
{
    return sin(a)/cos(a);
}

template<class T> PROGRAM_SCALAR<T>
asin(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_asin,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
acos(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_acos,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
atan(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_atan,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
floor(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_floor,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
ceil(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_ceil,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
rint(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_rint,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
min(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_min,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
max(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_max,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
abs(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_abs,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
Inverse(const PROGRAM_SCALAR<T>& a)
{
    return {a.pc,a.pc->Append_Instruction(op_inv,a.var)};
}

template<class T> PROGRAM_SCALAR<T>
atan2(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_atan2,a.var,b.var)};
}

template<class T> PROGRAM_SCALAR<T>
hypot(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return sqrt(a*a+b*b);
}

template<class T> PROGRAM_SCALAR<T>
pow(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return {a.pc,a.pc->Append_Instruction(op_pow,a.var,b.var)};
}

template<class T> std::ostream&
operator<<(std::ostream& o, const PROGRAM_SCALAR<T>& s)
{
    PHYSBAM_FATAL_ERROR();
}

template<class T> std::istream&
operator>>(std::istream& i, PROGRAM_SCALAR<T>& s)
{
    PHYSBAM_FATAL_ERROR();
}

template<class T> PROGRAM_SCALAR<T>
Magnitude_Squared(const PROGRAM_SCALAR<T>& s)
{
    return s*s;
}

template<class T> PROGRAM_SCALAR<T>
Dot_Product(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b)
{
    return a*b;
}

template<class T> PROGRAM_SCALAR<T>
wrap(const PROGRAM_SCALAR<T>& a,const PROGRAM_SCALAR<T>& b,const PROGRAM_SCALAR<T>& c)
{
    return max(b,min(a,c));
}

template<class T> PROGRAM_SCALAR<T>
Diff(const PROGRAM_SCALAR<T>& a,MEM_IN var)
{
    return {a.pc,a.pc->Diff(a.var,var)};
}

template<class T> void
Diff(PROGRAM_SCALAR<T>* r,const PROGRAM_SCALAR<T>* a,int m,MEM_IN var)
{
    if(!m) return;
    ARRAY<REG> o(m);
    for(int i=0;i<m;i++) o(i)=a[i].var;
    auto d=a->pc->Diff(o,var);
    for(int i=0;i<m;i++)
        r[i]={a->pc,d(i)};
}

template<class T> void
Diff(PROGRAM_SCALAR<T>* r,const PROGRAM_SCALAR<T>& a,
    const MEM_IN* var,int n)
{
    for(int j=0;j<n;j++)
        r[j]={a.pc,a.pc->Diff(a.var,var[j])};
}

template<class T> void
Diff(PROGRAM_SCALAR<T>* r,int si,int sj,
    const PROGRAM_SCALAR<T>* a,int m,const MEM_IN* var,int n)
{
    if(!m) return;
    ARRAY<REG> o(m);
    for(int i=0;i<m;i++) o(i)=a[i].var;
    for(int j=0;j<n;j++)
    {
        auto d=a->pc->Diff(o,var[j]);
        for(int i=0;i<m;i++)
            r[si*i+sj*j]={a->pc,d(i)};
    }
}

template<class T> void PROGRAM_SCALAR<T>::
Print() const
{
    LOG::printf("output reg: %P\n",var);
    pc->Print();
}

template PROGRAM_SCALAR<float> operator+(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> operator-(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> operator+(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator-(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator*(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator/(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator<(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator>(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator<=(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator>=(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator==(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> operator!=(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> cond(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b,const PROGRAM_SCALAR<float>& c);
template PROGRAM_SCALAR<float> abs(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> sqrt(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> exp(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> log(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> sin(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> cos(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> tan(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> asin(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> acos(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> atan(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> floor(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> ceil(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> rint(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> Inverse(const PROGRAM_SCALAR<float>& a);
template PROGRAM_SCALAR<float> atan2(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> hypot(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> pow(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template std::ostream& operator<<(std::ostream& o, const PROGRAM_SCALAR<float>& s);
template std::istream& operator>>(std::istream& i, PROGRAM_SCALAR<float>& s);
template PROGRAM_SCALAR<float> Magnitude_Squared(const PROGRAM_SCALAR<float>& s);
template PROGRAM_SCALAR<float> Dot_Product(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b);
template PROGRAM_SCALAR<float> wrap(const PROGRAM_SCALAR<float>& a,const PROGRAM_SCALAR<float>& b,const PROGRAM_SCALAR<float>& c);
template VECTOR<PROGRAM_SCALAR<float>,1> Gradient(const PROGRAM_SCALAR<float>& a,const VECTOR<MEM_IN,1>& vars);
template VECTOR<PROGRAM_SCALAR<float>,2> Gradient(const PROGRAM_SCALAR<float>& a,const VECTOR<MEM_IN,2>& vars);
template VECTOR<PROGRAM_SCALAR<float>,3> Gradient(const PROGRAM_SCALAR<float>& a,const VECTOR<MEM_IN,3>& vars);

template PROGRAM_SCALAR<double> operator+(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> operator-(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> operator+(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator-(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator*(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator/(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator<(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator>(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator<=(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator>=(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator==(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> operator!=(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> cond(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b,const PROGRAM_SCALAR<double>& c);
template PROGRAM_SCALAR<double> abs(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> sqrt(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> exp(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> log(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> sin(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> cos(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> tan(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> asin(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> acos(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> atan(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> floor(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> ceil(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> rint(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> Inverse(const PROGRAM_SCALAR<double>& a);
template PROGRAM_SCALAR<double> atan2(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> hypot(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> pow(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template std::ostream& operator<<(std::ostream& o, const PROGRAM_SCALAR<double>& s);
template std::istream& operator>>(std::istream& i, PROGRAM_SCALAR<double>& s);
template PROGRAM_SCALAR<double> Magnitude_Squared(const PROGRAM_SCALAR<double>& s);
template PROGRAM_SCALAR<double> Dot_Product(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b);
template PROGRAM_SCALAR<double> wrap(const PROGRAM_SCALAR<double>& a,const PROGRAM_SCALAR<double>& b,const PROGRAM_SCALAR<double>& c);
template PROGRAM_SCALAR<float> Diff(const PROGRAM_SCALAR<float>& a,MEM_IN var);
template PROGRAM_SCALAR<double> Diff(const PROGRAM_SCALAR<double>& a,MEM_IN var);
template PROGRAM_SCALAR<float> PROGRAM_SCALAR<float>::Number(PROGRAM_COMPILER<float>& pc,float a);
template PROGRAM_SCALAR<double> PROGRAM_SCALAR<double>::Number(PROGRAM_COMPILER<double>& pc,double a);
template void Diff(PROGRAM_SCALAR<double>*,PROGRAM_SCALAR<double> const&,MEM_IN const*,int);
template void Diff(PROGRAM_SCALAR<double>*,PROGRAM_SCALAR<double> const*,int,MEM_IN);
template void Diff(PROGRAM_SCALAR<double>*,int,int,PROGRAM_SCALAR<double> const*,int,MEM_IN const*,int);
template void Diff(PROGRAM_SCALAR<float>*,PROGRAM_SCALAR<float> const&,MEM_IN const*,int);
template void Diff(PROGRAM_SCALAR<float>*,PROGRAM_SCALAR<float> const*,int,MEM_IN);
template void Diff(PROGRAM_SCALAR<float>*,int,int,PROGRAM_SCALAR<float> const*,int,MEM_IN const*,int);
template void PROGRAM_SCALAR<double>::Print() const;
template void PROGRAM_SCALAR<float>::Print() const;
template PROGRAM_SCALAR<float> cond_lazy(PROGRAM_SCALAR<float> const&,std::function<PROGRAM_SCALAR<ND_HELPER<float>::TYPE> ()>,std::function<PROGRAM_SCALAR<ND_HELPER<float>::TYPE> ()>);
template PROGRAM_SCALAR<double> cond_lazy(PROGRAM_SCALAR<double> const&,std::function<PROGRAM_SCALAR<ND_HELPER<double>::TYPE> ()>,std::function<PROGRAM_SCALAR<ND_HELPER<double>::TYPE> ()>);
template ARRAY<PROGRAM_SCALAR<float>> cond_lazy(PROGRAM_SCALAR<float> const&,std::function<ARRAY<PROGRAM_SCALAR<ND_HELPER<float>::TYPE>> ()>,std::function<ARRAY<PROGRAM_SCALAR<ND_HELPER<float>::TYPE>> ()>);
template ARRAY<PROGRAM_SCALAR<double>> cond_lazy(PROGRAM_SCALAR<double> const&,std::function<ARRAY<PROGRAM_SCALAR<ND_HELPER<double>::TYPE>> ()>,std::function<ARRAY<PROGRAM_SCALAR<ND_HELPER<double>::TYPE>> ()>);
}
}
