#ifndef   MTT_FIELD_FM_TMPL_H
#define   MTT_FIELD_FM_TMPL_H

// ===========================================================================
//  field_fq_tmpl.h -- 素体F_m のテンプレートによる実装
// ---------------------------------------------------------------------------
// ◎数学的説明
//  mを素数とし、F_mを整数をmで割った余りの集合とする。このときF_mの中では
// 加減乗除の演算を自由にすることが出来る。このF_mを有限体(素体)と言う。
// 
// ◎使用上の説明
//  Tをmodulo mでの多倍長整数演算が可能なクラスであるとし、Tは以下のI/Fを持つ
// とする(各引数は重複指定しても良いものとする)：
//   1. 各種コンストラクタ
//       T(), T(int), T(T const&)
//       T(const char*)も追加(10進数文字列)
//   2. 各種アサイン関数
//       static void set(int, T&), static void set(T const&, T&)
//   3. 各種単項演算関数
//       static void {neg|inv}(T const&, T&)
//   4. 各種2項演算関数
//       static void {add|sub|mul}(int, T&), 
//       static void {add|sub|mul}(Tconst&, T&)
//   5. 各種3項演算関数
//       static void {add|sub|mul}(int, T const&, T&),
//       static void {add|sub|mul}(T const&, int, T&),
//       static void {add|sub|mul}(T const&, T const&, T&)
//   6. その他
//       void normalize(), bool isZero()
//       
// このとき、e_<T, T> が目的とする素体F_mを表すクラスとなる。ただし、このまま
// このクラスが用いられる訳ではない。
// 
// ◎実装上の説明
//  この実装ではtemplateを駆使し、「a = b + c」のような記述を add(b, c, a)の
// ような関数呼出に変更するテクニックを使用している。すべての式は 
//   e_<base_type, expression> 
// という型になり、代入もしくは型の変換時に評価され、関数呼出になる。
//   e = (a+b)*(c-d);
// という式を例にとって説明する。(a+b)*(c-d) という式だけでは計算は行なわれる
// ことはなく、
//   e_<b_<e_<T, T>, binary_plus_<T>, e_<T, T> >,
//      binary_star_<T>,
//      b_<e_<T, T>, binary_minus_<T>, e_<T, T> > >
// という型の変数に a, b, c, d の参照が保存される。実際の評価は変数eへの代入
// の際に行なわれ、まず関数e_<>.ev_()が呼ばれ、その中で binary_star_<>.ev_()
// が呼ばれ、その中で関数e_<>.ev_()、binary_{plus|minus}_<>.ev_()が呼ばれる
// という具合になるが、この呼出過程はコンパイラによりinline展開される。


namespace mttfm {
// ===========================================================================
//  "const&" comes after T iff T is e<> when one write "rr<T>::t_"
template<class T, class E> class e_;
template<class T> struct rr { typedef T t_; };
template<class T, class E> struct rr<e_<T, E> > { typedef e_<T, E>const& t_; };


// ===========================================================================
//  binary expressions
template<class LT, class Op, class RT>
struct b_ {
    typename rr<LT>::t_ vl_;
    typename rr<RT>::t_ vr_;
    b_(LT const& vl, RT const& vr) : vl_(vl), vr_(vr) {}
private:
    b_(); // to prevent a compiler calling the default constructor
};


// ===========================================================================
//  unary expressions
template<class Op, class T>
struct u_ {
    typename rr<T>::t_ v_;
    u_(T const& v) : v_(v) {}
private:
    u_(); // to prevent a compiler calling the default constructor
};


// ===========================================================================
//  type reduction : evaluates expression
template <class T, class E>
inline void se_(e_<T, E>const& ex, e_<T, T>& d) { ex.ev_(d); }


// ===========================================================================
//  the body of the class
template<class T>
class e_<T, T> {
public:
    // ----------------------------------------------------------------------
    //  type
    typedef T value_type;
    
    // ----------------------------------------------------------------------
    //  flags_
    enum {
        IS_ZERO = 1 << 0,
        ENUM_FORBIDS_EXTRA_COMMAS
    };
    
    // ----------------------------------------------------------------------
    //  ctor, dtor
    e_() : flags_(IS_ZERO) {}
    explicit e_(int r) : flags_(0)
        { if(r) T::set(r, r_); else flags_ = IS_ZERO; }
    explicit e_(T const& r) : flags_(0), r_(r) {}
    e_(e_ const& fm) :flags_(fm.flags_)
        { if(!(flags_&IS_ZERO)) T::set(fm.r_, r_); }
    template<class E>
    e_(e_<T, E>const& fm) { se_(fm, *this); }
    ~e_() {}
    
    // ----------------------------------------------------------------------
    //  accessors
    T& getR() { return r_; }
    T const& getR() const { return r_; }
    int& getFlags() { return flags_; }
    int getFlags() const { return flags_; }
    bool isZero() const { return flags_ == IS_ZERO; }
    
    // ----------------------------------------------------------------------
    //  operations
    void inverse() { if(!(flags_&IS_ZERO)) T::inv(r_, r_); }
    void normalize() {
        if(!(flags_&IS_ZERO)) {
            r_.normalize();
            if(r_.isZero()) flags_ |= IS_ZERO;
        }
    }
    void negate() { if(!(flags_&REAL_IS_ZERO)) T::neg(r_, r_); }
    
    // ----------------------------------------------------------------------
    // operator=
    e_ const& operator=(int r) {
        if(r) {
            T::set(r, r_);
            flags_ = 0;
        } else flags_ = IS_ZERO;
        return *this;
    }
    e_ const& operator=(T const& r) {
        T::set(r, r_);
        flags_ = 0;
        return *this;
    }
    e_ const& operator=(e_ const& fm) {
        flags_ = fm.flags_;
        if(!(flags_&IS_ZERO)) T::set(fm.r_, r_);
        return *this;
    }
    template<class E>
    e_ const& operator=(e_<T, E>const& fm) { se_(fm, *this); return *this; }
    
    // ----------------------------------------------------------------------
    //  operator+=
    e_ const& operator+=(int r) {
        if(!r) return *this;
        if(flags_&IS_ZERO) {
            flags_ = 0;
            T::set(r, r_);
        } else T::add(r, r_);
        return *this;
    }
    e_ const& operator+=(T const& r) {
        if(flags_&IS_ZERO) {
            flags_ = 0;
            T::set(r, r_);
        } else T::add(r, r_);
        return *this;
    }
    e_ const& operator+=(e_ const& fm) {
        if(fm.flags_&IS_ZERO) return *this;
        if(flags_&IS_ZERO) {
            flags_ = 0;
            T::set(fm.r_, r_);
        } else T::add(fm.r_, r_);
        return *this;
    }
    template<class E> // with type reduction
    e_ const& operator+=(e_<T, E>const& fm) {
        e_<T, T> tmp(fm);
        return *this += tmp;
    }
    
    // ----------------------------------------------------------------------
    //  operator-=
    e_ const& operator-=(int r) {
        if(!r) return *this;
        if(flags_&IS_ZERO) {
            flags_ = 0;
            T::set(-r, r_);
        } else T::add(-r, r_);
        return *this;
    }
    e_ const& operator-=(T const& r) {
        if(flags_&IS_ZERO) {
            flags_ = 0;
            T::neg(r, r_);
        } else T::sub(r, r_);
        return *this;
    }
    e_ const& operator-=(e_ const& fm) {
        if(fm.flags_&IS_ZERO) return *this;
        if(flags_&IS_ZERO) {
            flags_ = 0;
            T::neg(fm.r_, r_);
        } else T::sub(fm.r_, r_);
        return *this;
    }
    template<class E> // with type reduction
    e_ const& operator-=(e_<T, E>const& fm) {
        e_<T, T> tmp(fm);
        return *this -= tmp;
    }
    
    // -------------------------------------
    //  operator*=
    e_ const& operator*=(int r) {
        if(!r) flags_ = IS_ZERO;
        else if(!(flags_&IS_ZERO)) T::mul(r, r_);
        return *this;
    }
    e_ const& operator*=(T const& r) {
        if(!(flags_&IS_ZERO)) T::mul(r, r_);
        return *this;
    }
    e_ const& operator*=(e_ const& fm) {
        if(fm.flags_ == IS_ZERO) flags_ = IS_ZERO;
        if(flags_ == IS_ZERO) return *this;
	T::mul(fm.r_, r_);
        return *this;
    }
    template<class E> // with type reduction
    e_ const& operator*=(e_<T, E>const& fm) {
        e_<T, T> tmp(fm);
        return *this *= tmp;
    }
protected:
    int flags_;
    T r_;
};


// ===========================================================================
//  expression classes
// -----------------------------------
// expressions like "a+a"
template<class T, class Op>
class e_<T, b_<e_<T, T>, Op, e_<T, T> > > {
    b_<e_<T, T>, Op, e_<T,T> > ex_;
public:
    e_(e_<T, T>const& vl, e_<T, T>const& vr) : ex_(vl, vr) {}
    void ev_(e_<T, T>& d) const { Op::ev_(ex_.vl_, ex_.vr_, d); }
};
// -----------------------------------
// expressions like "a+1"
template<class T, class Op>
class e_<T, b_<e_<T, T>, Op, int> >{
    b_<e_<T, T>, Op, int> ex_;
public:
    e_(e_<T, T>const& vl, int vr) : ex_(vl, vr) {}
    void ev_(e_<T, T>& d) const { Op::ev_(ex_.vl_, ex_.vr_, d); }
};
// -----------------------------------
// expressions like "1+a"
template<class T, class Op>
class e_<T, b_<int, Op, e_<T, T> > >{
    b_<int, Op, e_<T, T> > ex_;
public:
    e_(int vl, e_<T, T>const& vr) : ex_(vl, vr) {}
    void ev_(e_<T, T>& d) const { Op::ev_(ex_.vl_, ex_.vr_, d); }
};
// -----------------------------------
// expressions like "a*(b+c)"
template<class T, class Op, class E>
class e_<T, b_<e_<T, T>, Op, e_<T, E> > >{
    b_<e_<T, T>, Op, e_<T, E> > ex_;
public:
    e_(e_<T, T>const& vl, e_<T, E>const& vr) : ex_(vl, vr) {}
    void ev_(e_<T, T>& d) const {
        e_<T, T> tmp(ex_.vr_);
        Op::ev_(ex_.vl_, tmp, d);
    }
};
// -----------------------------------
// expressions like "a*b+c"
template<class T, class E, class Op>
class e_<T, b_<e_<T, E>, Op, e_<T, T> > >{
    b_<e_<T, E>, Op, e_<T, T> > ex_;
public:
    e_(e_<T, E>const& vl, e_<T, T>const& vr) : ex_(vl, vr) {}
    void ev_(e_<T, T>& d) const {
        e_<T, T> tmp(ex_.vl_);
        Op::ev_(tmp, ex_.vr_, d);
    }
};
// -----------------------------------
// expressions like "3*(a+b)"
template<class T, class Op, class E>
class e_<T, b_<int, Op, e_<T, E> > >{
    b_<int, Op, e_<T, E> > ex_;
public:
    e_(int vl, e_<T, E>const& vr) : ex_(vl, vr) {}
    void ev_(e_<T, T>& d) const {
        e_<T, T> tmp(ex_.vr_);
        Op::ev_(ex_.vl_, tmp, d);
    }
};
// -----------------------------------
// expressions like "(a+3)*b"
template<class T, class E, class Op>
class e_<T, b_<e_<T, E>, Op, int> >{
    b_<e_<T, E>, Op, int> ex_;
public:
    e_(e_<T, E>const& vl, int vr) : ex_(vl, vr) {}
    void ev_(e_<T, T>& d) const {
        e_<T, T> tmp(ex_.vl_);
        Op::ev_(tmp, ex_.vr_, d);
    }
};
// -----------------------------------
// expressions like "(a+3)*(b+c)"
template<class T, class E1, class Op, class E2>
class e_<T, b_<e_<T, E1>, Op, e_<T, E2> > >{
    b_<e_<T, E1>, Op, e_<T, E2> > ex_;
public:
    e_(e_<T, E1>const& vl, e_<T, E2>const& vr) : ex_(vl, vr) {}
    void ev_(e_<T, T>& d) const {
        e_<T, T> tmpl(ex_.vl_), tmpr(ex_.vr_);
        Op::ev_(tmpl, tmpr, d);
    }
};
// -----------------------------------
// expressions like "-a"
template<class T, class Op>
class e_<T, u_<Op, e_<T, T> > > {
    u_<Op, e_<T, T> > ex_;
public:
    e_(e_<T, T>const& v) : ex_(v) {}
    void ev_(e_<T, T>& d) const { Op::ev_(ex_.v_, d); }
};
// -----------------------------------
// expressions like "-(a+a)"
template<class T, class Op, class E>
class e_<T, u_<Op, e_<T, E> > > {
    u_<Op, e_<T, E> > ex_;
public:
    e_(e_<T, E>const& v) : ex_(v) {}
    void ev_(e_<T, T>& d) const {
        e_<T, T> tmp(ex_.v_);
        Op::ev_(tmp, d);
    }
};


// ===========================================================================
//  macros for operators
#define IsZero(x) (x.getFlags()&e_<T,T>::IS_ZERO)


// ===========================================================================
//  operator+
template<class T>
struct binary_plus_ {
    static void ev_(e_<T, T>const& lt, e_<T, T>const& rt, e_<T, T>& d) {
        int f = 0;
        if(IsZero(lt)) {
            if(IsZero(rt)) f = e_<T, T>::IS_ZERO;
            else T::set(rt.getR(), d.getR());
        } else {
            if(IsZero(rt)) T::set(lt.getR(), d.getR());
            else T::add(lt.getR(), rt.getR(), d.getR());
        }
        d.getFlags() = f;
    }
    static void ev_(e_<T, T>const& lt, int rt, e_<T, T>& d) {
        if(!rt) {
            d.getFlags() = lt.getFlags();
            if(!IsZero(d)) T::set(lt.getR(), d.getR());
            return;
        }
        d.getFlags() = 0;
        if(IsZero(lt)) T::set(rt, d.getR());
        else T::add(lt.getR(), rt, d.getR());
    }
    static void ev_(int lt, e_<T, T>const& rt, e_<T, T>& d) {
        ev_(rt, lt, d);
    }
};
template<class T, class LT, class RT>
inline e_<T, b_<e_<T, LT>, binary_plus_<T>, e_<T, RT> > >
operator+(e_<T, LT>const& lt, e_<T, RT>const& rt) {
    return e_<T, b_<e_<T, LT>, binary_plus_<T>, e_<T, RT> > >(lt, rt);
}
template<class T, class RT>
inline e_<T, b_<int, binary_plus_<T>, e_<T, RT> > >
operator+(int lt, e_<T, RT>const& rt) {
    return e_<T, b_<int, binary_plus_<T>, e_<T, RT> > >(lt, rt);
}
template<class T, class LT>
inline e_<T, b_<e_<T, LT>, binary_plus_<T>, int> >
operator+(e_<T, LT>const& lt, int rt) {
    return e_<T, b_<e_<T, LT>, binary_plus_<T>, int> >(lt, rt);
}


// ===========================================================================
//  operator-
template<class T>
struct binary_minus_ {
    static void ev_(e_<T, T>const& lt, e_<T, T>const& rt, e_<T, T>& d) {
        int f = 0;
        if(IsZero(lt)) {
            if(IsZero(rt)) f = e_<T, T>::IS_ZERO;
            else T::neg(rt.getR(), d.getR());
        } else {
            if(IsZero(rt)) T::set(lt.getR(), d.getR());
            else T::sub(lt.getR(), rt.getR(), d.getR());
        }
        d.getFlags() = f;
    }
    static void ev_(e_<T, T>const& lt, int rt, e_<T, T>& d) {
        if(!rt) {
            d.getFlags() = lt.getFlags();
            if(!RealIsZero(d)) T::set(lt.getR(), d.getR());
            return;
        }
        d.getFlags() = 0;
        if(RealIsZero(lt)) T::set(-rt, d.getR());
        else T::sub(lt.getR(), rt, d.getR());
    }
    static void ev_(int lt, e_<T, T>const& rt, e_<T, T>& d) {
        if(!lt) {
            d.getFlags() = rt.getFlags();
            if(!RealIsZero(d)) T::neg(rt.getR(), d.getR());
            return;
        }
        d.getFlags() = 0;
        if(IsZero(rt)) T::set(lt, d.getR());
        else T::sub(lt, rt.getR(), d.getR());
    }
};
template<class T, class LT, class RT>
inline e_<T, b_<e_<T, LT>, binary_minus_<T>, e_<T, RT> > >
operator-(e_<T, LT>const& lt, e_<T, RT>const& rt) {
    return e_<T, b_<e_<T, LT>, binary_minus_<T>, e_<T, RT> > >(lt, rt);
}
template<class T, class RT>
inline e_<T, b_<int, binary_minus_<T>, e_<T, RT> > >
operator-(int lt, e_<T, RT>const& rt) {
    return e_<T, b_<int, binary_minus_<T>, e_<T, RT> > >(lt, rt);
}
template<class T, class LT>
inline e_<T, b_<e_<T, LT>, binary_minus_<T>, int> >
operator-(e_<T, LT>const& lt, int rt) {
    return e_<T, b_<e_<T, LT>, binary_minus_<T>, int> >(lt, rt);
}


// ===========================================================================
//  operator*
template<class T>
struct binary_star_ {
    static void ev_(e_<T, T>const& lt, e_<T, T>const& rt, e_<T, T>& d) {
        if(IsZero(lt)||IsZero(rt)) {
            d.getFlags() = e_<T, T>::IS_ZERO;
            return;
        }
        if(&lt == &rt) {
            if(&lt == &d) T::mul(d.getR(), d.getR());
            else {
                d.getFlags() = lt.getFlags();
                T::mul(lt.getR(), lt.getR(), d.getR());
            }
        } else {
            if(&lt == &d) T::mul(rt.getR(), d.getR());
            else if(&rt == &d) T::mul(lt.getR(), d.getR());
            else T::mul(lt.getR(), rt.getR(), d.getR());
        }
    }
    static void ev_(e_<T, T>const& lt, int rt, e_<T, T>& d) {
        if(!rt) {
            d.getFlags() = e_<T, T>::IS_ZERO;
            return;
        }
        d.getFlags() = lt.getFlags();
        if(&lt == &d) {
            if(!IsZero(d)) T::mul(rt, d.getR());
            return;
        }
        if(!IsZero(lt)) T::mul(lt.getR(), rt, d.getR());
    }
    static void ev_(int lt, e_<T, T>const& rt, e_<T, T>& d) {
        ev_(rt, lt, d);
    }
};
template<class T, class LT, class RT>
inline e_<T, b_<e_<T, LT>, binary_star_<T>, e_<T, RT> > >
operator*(e_<T, LT>const& lt, e_<T, RT>const& rt) {
    return e_<T, b_<e_<T, LT>, binary_star_<T>, e_<T, RT> > >(lt, rt);
}
template<class T, class RT>
inline e_<T, b_<int, binary_star_<T>, e_<T, RT> > >
operator*(int lt, e_<T, RT>const& rt) {
    return e_<T, b_<int, binary_star_<T>, e_<T, RT> > >(lt, rt);
}
template<class T, class LT>
inline e_<T, b_<e_<T, LT>, binary_star_<T>, int> >
operator*(e_<T, LT>const& lt, int rt) {
    return e_<T, b_<e_<T, LT>, binary_star_<T>, int> >(lt, rt);
}


// ===========================================================================
//  operator-
template<class T>
struct unary_minus_ {
    static void ev_(e_<T, T>const& t, e_<T, T>& d) {
        d.getFlags() = t.getFlags();
        if(!IsZero(d)) T::neg(t.getR(), d.getR());
    }
};
template<class T, class E>
inline e_<T, u_<unary_minus_<T>, e_<T, E> > >
operator-(e_<T, E>const& v) {
    return e_<T, u_<unary_minus_<T>, e_<T, E> > >(v);
}


// ===========================================================================
//  function inverse
template<class T>
struct unary_inverse_ {
    static void ev_(e_<T, T>const& t, e_<T, T>& d) {
        d.getFlags() = t.getFlags();
        if(!IsZero(d)) T::inv(t.getR(), d.getR());
    }
};
template<class T, class E>
inline e_<T, u_<unary_inverse_<T>, e_<T, E> > >
inverse(e_<T, E>const& v) {
    return e_<T, u_<unary_inverse_<T>, e_<T, E> > >(v);
}


// ===========================================================================
//  end of macros for operators
#undef RealIsZero
#undef ImagIsZero


// ===========================================================================
//  end of namespace
}


#endif /* MTT_FIELD_FM_TMPL_H */