//#####################################################################
// Copyright 2013.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class PRIMITIVE_TENSORS
//##################################################################### 
#ifndef __PRIMITIVE_TENSORS__
#define __PRIMITIVE_TENSORS__

#include <Core/Math_Tools/cube.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/PRIMITIVE_MATRICES.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Tensors/DIAGONAL_TENSOR.h>
#include <Core/Tensors/PERMUTATION_TENSOR.h>
#include <Core/Tensors/SYMMETRIC_TENSOR.h>
#include <Core/Tensors/TENSOR.h>
#include <Core/Tensors/TENSORS_FORWARD.h>
#include <Core/Tensors/VEC_ID_SYM_TENSOR.h>
#include <Core/Tensors/VEC_ID_TENSOR.h>
#include <Core/Tensors/ZERO_TENSOR.h>
#include <Core/Vectors/VECTOR.h>
#include <cmath>
namespace PhysBAM{

/*
ZERO_VECTOR
VECTOR

ZERO_MATRIX
IDENTITY_MATRIX
SCALE_MATRIX
DIAGONAL_MATRIX
SYMMETRIC_MATRIX
MATRIX

ZERO_TENSOR
TENSOR
SYMMETRIC_TENSOR
DIAGONAL_TENSOR
VEC_ID_TENSOR
VEC_ID_SYM_TENSOR
PERMUTATION_TENSOR
*/

constexpr int ch3(int r,int i,int j,int k)
{
    if(r==0) return i;
    if(r==1) return j;
    return k;
}

consteval int add_type_ten(int a,int b)
{
    if(a>b) return add_type_ten(b,a);
    if(a==b) return a;
    if(b==ten_zero) return a;
    if(a==ten_gen) return a;
    if(b==ten_perm) return ten_gen;
    if(b==ten_diag) return ten_sym+(a-1)%3;
    if((a-1)%3!=(b-1)%3) return ten_gen;
    return ten_sym+(a-1)%3;
}

// -tensor
template<class A,class E=std::enable_if_t<IS_TENSOR<A>::value>,class R=A>
R __attribute__((flatten))
operator-(const A& a)
{
    R r;
    int b=0;
#define FUNCTION(i,j,k) (-get<i,j,k>(a))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

// tensor * scalar
template<class A,class T=typename A::SCALAR,
    class E=std::enable_if_t<IS_TENSOR<A>::value>,class R=A>
R __attribute__((flatten))
operator*(const A& a,T b)
{
    R r;
#define FUNCTION(i,j,k) (get<i,j,k>(a)*b)
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

// tensor / scalar
template<class A,class T=typename A::SCALAR,
    class E=std::enable_if_t<IS_TENSOR<A>::value>,class R=A>
R __attribute__((flatten))
operator/(const A& a,T b)
{
    R r;
#define FUNCTION(i,j,k) (get<i,j,k>(a)/b)
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

// scalar * tensor
template<class A,class T=typename A::SCALAR,
    class E=std::enable_if_t<IS_TENSOR<A>::value>,class R=A>
A __attribute__((flatten))
operator*(T b,const A& a)
{
    return a*b;
}

// tensor + tensor
template<class A,class B,class T=typename A::SCALAR,
    class E=std::enable_if_t<IS_TENSOR<A>::value&&IS_TENSOR<B>::value>,
    class F=std::enable_if_t<A::m==B::m&&A::n==B::n&&A::p==B::p>,
    class R=TEN_TYPE<T,add_type_ten(A::type,B::type),A::m,A::n,A::p>>
R __attribute__((flatten))
operator+(const A& a,const B& b)
{
    R r;
#define FUNCTION(i,j,k) (get<i,j,k>(a)+get<i,j,k>(b))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

// tensor - tensor
template<class A,class B,class T=typename A::SCALAR,
    class E=std::enable_if_t<IS_TENSOR<A>::value&&IS_TENSOR<B>::value>,
    class F=std::enable_if_t<A::m==B::m&&A::n==B::n&&A::p==B::p>,
    class R=TEN_TYPE<T,add_type_ten(A::type,B::type),A::m,A::n,A::p>>
R __attribute__((flatten))
operator-(const A& a,const B& b)
{
    R r;
#define FUNCTION(i,j,k) (get<i,j,k>(a)-get<i,j,k>(b))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

// tensor += tensor
template<class A,class B,class R=A,
    class E=std::enable_if_t<IS_TENSOR<A>::value&&IS_TENSOR<B>::value&&is_same<A,decltype(A()+B())>::value>>
R& __attribute__((flatten))
operator+=(A& r,const B& b)
{
    int a=0;
#define FUNCTION(i,j,k) (get<i,j,k>(r)+get<i,j,k>(b))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

// tensor -= tensor
template<class A,class B,class R=A,
    class E=std::enable_if_t<IS_TENSOR<A>::value&&IS_TENSOR<B>::value&&is_same<A,decltype(A()-B())>::value>>
R& __attribute__((flatten))
operator-=(A& r,const B& b)
{
    int a=0;
#define FUNCTION(i,j,k) (get<i,j,k>(r)-get<i,j,k>(b))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

consteval int tensor_product_type_ten(int s,int a,int b)
{
    if(a==mat_zero || b==vec_zero) return ten_zero;
    if(a==mat_id || a==mat_scale) return ten_vec_id+s;
    if(a==mat_gen) return ten_gen;
    return ten_sym+s;
}

// Tensor_Product(matrix,scalar)
template<int s,class A,class B,class T=typename A::SCALAR,
    class R=TEN_TYPE<T,tensor_product_type_ten(s,A::type,B::type),
        (s==0?B::m:A::m),ch3(s,A::m,B::m,A::n),(s==2?B::m:A::n)>,
    class E=std::enable_if_t<IS_MATRIX<A>::value&&IS_VECTOR<B>::value>>
R __attribute__((flatten))
Tensor_Product(const A& a,const B& b)
{
    R r;
#define FUNCTION(i,j,k) (get<(s==0?j:i),(s==2?j:k)>(a)*b(ch3(s,i,j,k)))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

// Fill_From
template<class A,class B,class R=A,
    class E=std::enable_if_t<IS_TENSOR<A>::value&&IS_TENSOR<B>::value>>
R& __attribute__((flatten))
Fill_From(A& r,const B& b)
{
    typedef decltype(A()+B()) Z;
    static_assert(is_same<A,Z>::value);
    static_assert(is_same<A,decltype(A()+B())>::value);
    int a=0;
#define FUNCTION(i,j,k) (get<i,j,k>(b))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
    return r;
}

consteval int transpose_type_ten(int u,int a)
{
    if(a==ten_zero || a==ten_gen || a==ten_perm || a==ten_diag) return a;
    int t=(a-1)%3;
    if(t==u) return a;
    return ten_sym+(3-t-u)%3;
}

constexpr bool is_sym_index(int u,int a)
{
    if(a<ten_sym || a>ten_vec_id_sym+2) return false;
    return (a-1)%3==u;
}

// Transposed
template<int s,int t,class A,
    class E=std::enable_if_t<IS_TENSOR<A>::value>,
    class T=typename A::SCALAR,
    int u=3-s-t,
    class R0=TEN_TYPE<T,transpose_type_ten(u,A::type),ch3(u,A::m,A::p,A::n),ch3(u,A::p,A::n,A::m),ch3(u,A::n,A::m,A::p)>,
    class R=conditional_t<is_sym_index(u,A::type),const A&,R0>>
R __attribute__((flatten))
Transposed(const A& a)
{
    if constexpr(is_sym_index(u,A::type)) return a;
    else
    {
        R r;
        int b=0;
#define FUNCTION(i,j,k) (get<ch3(u,i,k,j),ch3(u,k,j,i),ch3(u,j,i,k)>(a))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
        return r;
    }
}

consteval int tsp_type_ten(int u,int a)
{
    if(a==ten_zero || a==ten_perm) return ten_zero;
    if(a==ten_diag) return a;
    if(a==ten_gen) return ten_sym+u;
    int t=(a-1)%3;
    if(t==u) return a;
    return ten_sym+u;
}

// Twice_Symmetric_Part
template<int s,int t,class A,
    class E=std::enable_if_t<IS_TENSOR<A>::value>,
    class T=typename A::SCALAR,
    int u=3-s-t,
    class F=std::enable_if_t<(u==0&&A::n==A::p)||(u==1&&A::m==A::p)||(u==2&&A::m==A::n)>,
    class R=TEN_TYPE<T,tsp_type_ten(u,A::type),A::m,A::n,A::p>>
R __attribute__((flatten))
Twice_Symmetric_Part(const A& a)
{
    if constexpr(is_sym_index(u,A::type)) return a*2;
    else
    {
        R r;
        int b=0;
#define FUNCTION(i,j,k) (get<i,j,k>(a)+get<ch3(u,i,k,j),ch3(u,k,j,i),ch3(u,j,i,k)>(a))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
        return r;
    }
}

template<int r,int s,class MAT,class VEC>
enable_if_t<(IS_MATRIX<MAT>::value && IS_VECTOR<VEC>::value),decltype(Twice_Symmetric_Part<r,s>(Tensor_Product<r>(MAT(),VEC())))> __attribute__((flatten))
Symmetric_Tensor_Product(const MAT& a,const VEC& b)
{
    return Twice_Symmetric_Part<r,s>(Tensor_Product<r>(a,b));
}

template<class TEN>
std::enable_if_t<IS_TENSOR<TEN>::value,typename TEN::SCALAR> __attribute__((flatten))
Sum_Of_Squares(const TEN& t)
{
    typename TEN::SCALAR r=0;
    if constexpr(TEN::type<=ten_sym+2)
        for(auto&a:t.x)
            r+=a.Frobenius_Norm_Squared();
    else if constexpr(TEN::type==ten_zero){}
    else if constexpr(TEN::type==ten_perm)
        r=6*t.x*t.x;
    else
    {
        r=t.v.Magnitude_Squared();
        if constexpr(TEN::type>=ten_vec_id && TEN::type<=ten_vec_id+2)
            r*=TEN::un;
        else if constexpr(TEN::type>=ten_vec_id_sym && TEN::type<=ten_vec_id_sym+2)
            r*=2*(TEN::un+1);
    }
    return r;
}

//#####################################################################
// Function Contract<r,s>(T)
//##################################################################### 

template<int r,int s,class A,
    class E=std::enable_if_t<IS_TENSOR<A>::value>,
    int u=3-r-s,
    class F=std::enable_if_t<(u==0 && A::n==A::p)||(u==1 && A::m==A::p)||(u==2 && A::m==A::n)>,
    class T=typename A::SCALAR,
    int m=ch3(u,A::m,A::n,A::p),int n=(u==0?A::n:A::m),
    class R=VEC_TYPE<T,(A::type==ten_zero||A::type==ten_perm?vec_zero:vec_gen),m>>
R __attribute__((flatten))
Contract(const A& a)
{
    if constexpr(A::type==ten_zero || A::type==ten_perm) return {};
    else if constexpr(A::type==ten_vec_id+u) return a.v*A::un;
    else if constexpr(A::type>=ten_vec_id && A::type<=ten_vec_id+2) return a.v;
    else if constexpr(A::type==ten_vec_id_sym+u) return a.v*2;
    else if constexpr(A::type>=ten_vec_id_sym && A::type<=ten_vec_id_sym+2) return a.v*(A::un+1);
    else if constexpr(A::type==ten_diag) return a.v;
    else
    {
        return []<int... i,int... k>(const auto& a,std::integer_sequence<int,i...>,std::integer_sequence<int,k...>)
            {
                int j{};(void)j;
                if constexpr(u==0) return R{(j=i,(...+(a(j,k,k))))...};
                if constexpr(u==1) return R{(j=i,(...+(a(k,j,k))))...};
                if constexpr(u==2) return R{(j=i,(...+(a(k,k,j))))...};
            }(a,std::make_integer_sequence<int,m>{},std::make_integer_sequence<int,n>{});
    }
}

//#####################################################################
// Function Contract<s>(T,V)
//##################################################################### 

consteval int contract_ten_vec(int s,int a,int b)
{
    if(a==ten_zero || b==vec_zero) return mat_zero;
    if(a==ten_diag) return mat_diag;
    if(a==ten_sym+s) return mat_sym;
    if(a==ten_vec_id_sym+s) return mat_sym;
    if(a==ten_vec_id+s) return mat_scale;
    return mat_gen;
}

template<int s,class A,class B,
    class E=std::enable_if_t<IS_TENSOR<A>::value&&IS_VECTOR<B>::value>,
    class F=std::enable_if_t<(s==0 && A::m==B::m)||(s==1 && A::n==B::m)||(s==2 && A::p==B::m)>,
    class T=typename A::SCALAR,
    int m=(s==0?A::n:A::m),int n=(s==2?A::n:A::p),
    class R=MAT_TYPE<T,contract_ten_vec(s,A::type,B::type),m,n>>
R __attribute__((flatten))
Contract(const A& a,const B& b)
{
    if constexpr(A::type==ten_zero || B::type==vec_zero) return {};
    else if constexpr(A::type==ten_perm) return Cross_Product_Matrix(s==1?a.x*b:-a.x*b);
    else if constexpr(A::type==ten_diag) return R(b*a.v);
    else if constexpr(A::type==ten_vec_id+s) return R(a.v.Dot(b));
    else if constexpr(A::type==ten_vec_id+(s==0?1:0)) return Outer_Product(a.v,b);
    else if constexpr(A::type==ten_vec_id+(s==2?1:2)) return Outer_Product(b,a.v);
    else if constexpr(A::type==ten_vec_id_sym+s) return Symmetric_Outer_Product(b,a.v);
    else if constexpr(A::type==ten_vec_id_sym+(s==0?1:0)) return Outer_Product(b,a.v)+b.Dot(a.v);
    else if constexpr(A::type==ten_vec_id_sym+(s==2?1:2)) return Outer_Product(a.v,b)+b.Dot(a.v);
    else if constexpr(A::type==ten_sym+s)
    {
        return []<int... i,int... k>(auto&& a,auto&& b,std::integer_sequence<int,i...>,std::integer_sequence<int,k...>)
            {
                int j{};(void)j;
                return R{(j=i,((a.x(k).x[j]*b(k))+...))...};
            }(a,b,std::make_integer_sequence<int,m*(m+1)/2>{},std::make_integer_sequence<int,B::m>{});
    }
    else
    {
        static_assert(A::type<=ten_sym+2);
        return []<int... i,int... k>(auto&& a,auto&& b,std::integer_sequence<int,i...>,std::integer_sequence<int,k...>)
            {
                int j{};
                if constexpr(s==0) return R{(j=i,((a(k,j%m,j/m)*b(k))+...))...};
                if constexpr(s==1) return R{(j=i,((a(j%m,k,j/m)*b(k))+...))...};
                if constexpr(s==2) return R{(j=i,((a(j%m,j/m,k)*b(k))+...))...};
            }(a,b,std::make_integer_sequence<int,m*n>{},std::make_integer_sequence<int,B::m>{});
    }
}

//#####################################################################
// Function Contract<r,s>(T,M)
//##################################################################### 

consteval int contract_ten_mat(int r,int s,int a,int b)
{
    if(a==ten_zero || b==mat_zero) return ten_zero;
    if(b==mat_id || b==mat_scale) return a;
    if(a==ten_sym+r) return a;
    if(a==ten_vec_id+r) return a;
    if(a==ten_vec_id_sym+r) return ten_sym+r;
    if(a==ten_diag) return b==mat_diag?ten_diag:ten_sym+r;
    if(a>=ten_vec_id && a<=ten_vec_id+2 && (b==mat_diag || b==mat_sym))
        return ten_sym+(a-ten_vec_id);
    return ten_gen;
}

template<int t,int s,class A,class B,
    class E=std::enable_if_t<IS_TENSOR<A>::value&&IS_MATRIX<B>::value>,
    int bm=s==0?B::m:B::n,int bn=s==0?B::n:B::m,
    int m=(t==0?bn:A::m),int n=(t==1?bn:A::n),int p=(t==2?bn:A::p),
    class F=std::enable_if_t<(t==0 && A::m==bm)||(t==1 && A::n==bm)||(t==2 && A::p==bm)>,
    class T=typename A::SCALAR,
    class R=TEN_TYPE<T,contract_ten_mat(t,s,A::type,B::type),m,n,p>>
R __attribute__((flatten))
Contract(const A& a,const B& b)
{
    if constexpr(A::type==ten_zero || B::type==mat_zero) return {};
    else if constexpr(B::type==mat_id) return a;
    else if constexpr(B::type==mat_scale) return a*b.x;
    else if constexpr(A::type==ten_vec_id+t && s==0) return R{b.Transpose_Times(a.v)};
    else if constexpr(A::type==ten_vec_id+t && s==1) return R{b*a.v};
    else if constexpr(A::type>=ten_vec_id && A::type<=ten_vec_id+2)
    {
        static constexpr int u=(A::type-1)%3;
        if constexpr((t<3-t-u)!=(s==0)) return Tensor_Product<u>(b,a.v);
        else return Tensor_Product<u>(b.Transposed(),a.v);
    }
    else if constexpr(A::type==ten_vec_id_sym+t)
    {
        static constexpr int u=(A::type-1)%3,i=(u==0),j=2-(u==2);
        if constexpr(s==j-1) return Twice_Symmetric_Part<i,j>(Tensor_Product<i>(b,a.v));
        else return Twice_Symmetric_Part<i,j>(Tensor_Product<i>(b.Transposed(),a.v));
    }
    else if constexpr(A::type>=ten_vec_id_sym && A::type<=ten_vec_id_sym+2)
    {
        static constexpr int u=(A::type-1)%3;
        return Contract<t,s>(VEC_ID_TENSOR<T,u==0?1:0,m,m>(a.v),b)+Contract<t,s>(VEC_ID_TENSOR<T,u==2?1:2,m,m>(a.v),b);
    }
    else if constexpr(B::type==mat_diag)
    {
        R r;
    #define FUNCTION(i,j,k) (get<i,j,k>(a)*b.x(ch3(t,i,j,k)))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
        return r;
    }
    else if constexpr(A::type==ten_diag)
    {
        R r;
    #define SWP(x,y) (s==0?x:y),(s==0?y:x)
#define FUNCTION(i,j,k) (get<(t==0?j:i),(t==1?i:j),(t==2?i:k)>(a)*get<SWP((t==0?j:i),ch3(t,i,j,k))>(b))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_VISIT.h>
#undef SWP
        return r;
    }
    else
    {
        static_assert(A::type<=ten_sym+2 || A::type==ten_perm);
        static_assert(B::type==mat_gen || B::type==mat_sym);
        constexpr int q=bm;
        R r;
#define FUNCTION(i,j,k,l) (a(t==0?l:i,t==1?l:j,t==2?l:k)*b(s==0?l:ch3(t,i,j,k),s==1?l:ch3(t,i,j,k)))
#include <Core/Tensors/TENSOR_OPERATIONS_DOF_CONTRACT.h>
        return r;
    }
}

template<class A,class B,class C,class E=std::enable_if_t<IS_TENSOR<C>::value && IS_MATRIX<A>::value && IS_MATRIX<B>::value>>
auto Symmetric_Double_Contract_12(const C& t,const A& a,const B& b)
{return Twice_Symmetric_Part<1,2>(Contract<2,0>(Contract<1,0>(t,a),b));}

// tensor == tensor
template<class A,class B,class T=typename A::SCALAR,
    class E=std::enable_if_t<IS_TENSOR<A>::value&&IS_TENSOR<B>::value>,
    class F=std::enable_if_t<A::m==B::m&&A::n==B::n&&A::p==B::p>,
    class R=TEN_TYPE<T,add_type_ten(A::type,B::type),A::m,A::n,A::p>>
bool __attribute__((flatten))
operator==(const A& a,const B& b)
{
    if constexpr(A::m==0 || A::n==0 || A::p==0) return true;
    else if constexpr(R::type==ten_zero) return true;
    else if constexpr(R::type==ten_perm) return get<0,1,2>(a)==get<0,1,2>(b);
    else if constexpr(R::type==ten_diag)
    {
        return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
            {
                return (...&&(get<i,i,i>(a)==get<i,i,i>(b)));
            }(a,b,std::make_integer_sequence<int,R::m>{});
    }
    else if constexpr(R::type>=ten_vec_id && R::type<=ten_vec_id+2)
    {
        static constexpr int di=(R::type==ten_vec_id+0);
        static constexpr int dj=(R::type==ten_vec_id+1);
        static constexpr int dk=(R::type==ten_vec_id+2);
        return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
            {
                return (...&&(get<i*di,i*dj,i*dk>(a)==get<i*di,i*dj,i*dk>(b)));
            }(a,b,std::make_integer_sequence<int,R::um>{});
    }
    else if constexpr(R::type>=ten_vec_id_sym && R::type<=ten_vec_id_sym+2)
    {
        static constexpr int di=(R::type!=ten_vec_id_sym+2);
        static constexpr int dj=(R::type!=ten_vec_id_sym+0);
        static constexpr int dk=(R::type!=ten_vec_id_sym+1);
        if constexpr(R::m==1) return get<0,0,0>(a)==get<0,0,0>(b);
        else return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
                 {
                     return (...&&(get<(i+di)%R::m,(i+dj)%R::m,(i+dk)%R::m>(a)==get<(i+di)%R::m,(i+dj)%R::m,(i+dk)%R::m>(b)));
                 }(a,b,std::make_integer_sequence<int,R::m>{});
    }
    else if constexpr(R::type>=ten_sym && R::type<=ten_sym+2)
    {
        static constexpr int ss=R::un*(R::un+1)/2;
        return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
            {
                if constexpr(R::u==0)
                    return (...&&(get<i/ss,sym_i(R::un,i%ss),sym_j(R::un,i%ss)>(a)==get<i/ss,sym_i(R::un,i%ss),sym_j(R::un,i%ss)>(b)));
                else if constexpr(R::u==1)
                    return (...&&(get<sym_i(R::un,i%ss),i/ss,sym_j(R::un,i%ss)>(a)==get<sym_i(R::un,i%ss),i/ss,sym_j(R::un,i%ss)>(b)));
                else if constexpr(R::u==2)
                    return (...&&(get<sym_i(R::un,i%ss),sym_j(R::un,i%ss),i/ss>(a)==get<sym_i(R::un,i%ss),sym_j(R::un,i%ss),i/ss>(b)));
            }(a,b,std::make_integer_sequence<int,R::um*ss>{});
    }
    else
    {
        static_assert(R::type==ten_gen);
        return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
            {
                return (...&&(get<i/R::n/R::p,i%R::n,i/R::n%R::p>(a)==get<i/R::n/R::p,i%R::n,i/R::n%R::p>(b)));
            }(a,b,std::make_integer_sequence<int,R::m*R::n*R::p>{});
    }
}



template<class A,class B> std::enable_if_t<IS_TENSOR<A>::value&&IS_TENSOR<B>::value,decltype(A()+B())> Choose(const A&,const B&);
template<class T> T Choose_Zero(const T& a);
template<class T,int d> SCALE_MATRIX<T,d> Choose_Zero(const IDENTITY_MATRIX<T,d>& a);
template<class T,int m> typename conditional<m==0,FIXED_NUMBER<T,m>,T>::type Choose_Zero(const FIXED_NUMBER<T,m>& a);
}
#endif

