synaps/mpol/MPoly_2_2.h

// Copyright (c) 2003  INRIA Sophia-Antipolis (France) and
//                     Department of Informatics and Telecommunications
//                     University of Athens (Greece).
// All rights reserved.
//
// Authors : Elias P. TSIGARIDAS <et@di.uoa.gr>
//           Athanasios Kakargias <grad0460@di.uoa.gr>

// Partially supported by INRIA's project "CALAMATA", a
// bilateral collaboration with National Kapodistrian University of
// Athens.
// Partially supported by the IST Programme of the EU as a Shared-cost
// RTD (FET Open) Project under Contract No  IST-2000-26473
// (ECG - Effective Computational Geometry for Curves and Surfaces)
//

#ifndef SYNAPS_ARITHM_ALGEBRAIC_MPOLY_2_2_H
#define SYNAPS_ARITHM_ALGEBRAIC_MPOLY_2_2_H


#include <synaps/init.h>
#include <synaps/usolve/Algebraic_config.h>
#include <synaps/arithm/let.h>
#include <synaps/mpol.h>


__ALGEBRAIC_BEGIN_NAMESPACE

//F := c[5]*X^2 + c[4]*Y^2 + c[3]*X*Y + c[2]*X + c[1]*Y + c[0];
template< typename RT > struct MPoly_2_2;

template < typename POLY,
           typename RT >
POLY
specialize_y(const MPoly_2_2<RT>& f, const typename NumberTraits<RT>::FT& a);


template< typename RT_ >
struct MPoly_2_2 {

  static unsigned long counter;

  typedef NumberTraits<RT_>    NTR;
  typedef typename NTR::RT      RT;
  typedef typename NTR::FT      FT;

  typedef std::vector<RT>       Vector;
  
  typedef RT                    value_type;
  typedef unsigned              size_type;

  std::vector<RT>  C_;
  unsigned long    id_;
  FT xc_;
  FT yc_;
  

  MPoly_2_2() : C_(6)
  {
    id_ = counter++;
  }

  MPoly_2_2(const Vector& v)
    : C_(v)
  {
    id_ = counter++;
  }

  MPoly_2_2(RT* C)
    : C_(6)
  {
    id_ = counter++;
    for (int i = 0; i < 6; ++i)
      C_[i] = C[i];
  }

  MPoly_2_2(
            const RT& c5, const RT& c4, const RT& c3, const RT& c2, const RT& c1, const RT& c0
            ) : C_(6)
  {
    id_ = counter++;
    C_[5] = c5;
    C_[4] = c4;
    C_[3] = c3;
    C_[2] = c2;
    C_[1] = c1;
    C_[0] = c0;
  }



  MPoly_2_2( const FT& c5, const FT& c4, const FT& c3,
             const FT& c2, const FT& c1, const FT& c0)
    : C_(6)
  {
    id_ = counter++;
    RT d = RT(1);

    if (numerator(c5) != 0) d *= denominator(c5);
    if (numerator(c4) != 0) d *= denominator(c4);
    if (numerator(c3) != 0) d *= denominator(c3);
    if (numerator(c2) != 0) d *= denominator(c2);
    if (numerator(c1) != 0) d *= denominator(c1);
    if (numerator(c0) != 0) d *= denominator(c0);

    C_[5] = numerator(c5) * (d / denominator(c5));
    C_[4] = numerator(c4) * (d / denominator(c4));
    C_[3] = numerator(c3) * (d / denominator(c3));
    C_[2] = numerator(c2) * (d / denominator(c2));
    C_[1] = numerator(c1) * (d / denominator(c1));
    C_[0] = numerator(c0) * (d / denominator(c0));
  }

  MPoly_2_2& make_lcoeff_positive()
  {
    int i = 5;
    while (C_[i] == RT(0)) --i;

    assert(i >= 1);
    if (C_[i] < RT(0)) {
      for (int j = i; j >= 0; --j)
        C_[j] = -C_[j];
    }
    return *this;
  }

  MPoly_2_2& make_lcoeff_X_positive()
  {
    int pos;
    if (C_[5] != 0) pos = 5;
    else if (C_[2] != 0) pos = 2;
    else pos = 3;

    if (pos != 3) {
      if (C_[pos] < 0) {
        for (int i = 0; i < 5; ++i) {
          C_[i] = -C_[i];
        }
      }
    }
    return *this;
  }

  MPoly_2_2& make_lcoeff_Y_positive()
  {
    int pos;
    if (C_[4] != 0) pos = 4;
    else if (C_[1] != 0) pos = 1;
    else pos = 3;

    if (pos != 3) {
      if (C_[pos] < 0) {
        for (int i = 0; i < 5; ++i) {
          C_[i] = -C_[i];
        }
      }
    }
    return *this;
  }


  MPoly_2_2& revertXY()
  {
    std::swap(C_[5], C_[4]);
    std::swap(C_[2], C_[1]);
    return *this;
  }

  int degree(const int& idx) const
  {
    // degree wrt X if idx == 0
    // degree wrt Y otherwise
    if (idx == 0) {
      if (C_[5] != RT(0)) return 2;
      if ((C_[3] != RT(0))  || (C_[2] != RT(0))) return 1;
      return 0;
    }
    if (C_[4] != RT(0)) return 2;
    if ((C_[3] != RT(0))  || (C_[1] != RT(0))) return 1;
    return 0;
  }

  int degree() const
  {
    if ((C_[5] != RT(0)) || (C_[4] != RT(0)) || (C_[3] != RT(0))) return 2;
    if ((C_[2] != RT(0)) || (C_[1] != RT(0))) return 1;
    if (C_[0] != RT(0)) return 0;
    return -1;
  }

  void set_center( const FT& xc, const FT& yc)
  {
    xc_ = xc;
    yc_ = yc;
  }

  FT xc() const 
  {
    return xc_;
  }
  
  FT yc() const 
  {
    return yc_;
  }
  

  void set_coeff(const size_type& i, const RT& x)
  { C_[i] = x; }

  RT coeff(const size_type& i, const size_type& j) const
  {
    assert(i <= 2);
    assert(j <= 2);

    switch (i) {
    case 2:
      return C_[5];
      break;
    case 1:
      if (j == 1) return C_[3];
      if (j == 0) return C_[2];
      break;
    default:
      if (j == 2) return C_[4];
      if (j == 1) return C_[1];
      if (j == 0) return C_[0];
      break;
    }
  }

  RT operator[](const size_type& i) const
  {
    return C_[i];
  }

  unsigned long id() const { return id_; }

  template< typename POLY >
  int sign_at(const FT& x, const FT& y) const
  {
    return sign(eval(specialize_y<RT, POLY>(*this, y), x));
  }
};
template < typename RT >
unsigned long MPoly_2_2<RT>::counter = 0;


template < typename RT >
bool operator==( const MPoly_2_2<RT>& f1, const MPoly_2_2<RT>& f2)
{
  return (
          f1[0] == f2[0] &&
          f1[1] == f2[1] &&
          f1[2] == f2[2] &&
          f1[3] == f2[3] &&
          f1[4] == f2[4] &&
          f1[5] == f2[5] );
}

template < typename RT >
bool operator!=( const MPoly_2_2<RT>& f1, const MPoly_2_2<RT>& f2)
{
  return !(f1 == f2);
}



template < typename RT >
std::ostream&
print(std::ostream& os, const MPoly_2_2<RT>& f)
{
  if (f[5] != 0)
    os << "(" << f[5] << ")" << "*X^2";
  if (f[4] != 0)
    os << " + (" << f[4] << ")" << "*Y^2";
  if (f[3] != 0)
    os << "+ (" << f[3] << ")" << "*X*Y";
  if (f[2] != 0)
    os << " + (" << f[2] << ")" << "*X";
  if (f[1] != 0)
    os << " + (" << f[1] << ")" << "*Y";
  if (f[0] != 0)
    os << " + (" << f[0] << ")";
  return os;
}


template < typename RT >
std::ostream&
operator<<(std::ostream& output, const MPoly_2_2<RT>& f)
{
  return print(output, f);
}



template < typename RT >
MPoly_2_2<RT>
diff(const MPoly_2_2<RT>& f, const int& idx)
{
  //F := c[5]*X^2 + c[4]*Y^2 + c[3]*X*Y + c[2]*X + c[1]*Y + c[0];
  // diff wrt X if idx == 0
  // diff wrt Y otherwise
  RT coeff[6] = {0, 0, 0, 0, 0, 0};
  if (idx == 0) {
    coeff[2] = RT(2) * f[5];
    coeff[1] = f[3];
    coeff[0] = f[2];
  } else {
    coeff[2] = f[3];
    coeff[1] = RT(2) * f[4];
    coeff[0] = f[1];
  }
  return MPoly_2_2<RT>(coeff);
}



template < typename RT >
MPoly_2_2<RT>
diff_x(const MPoly_2_2<RT>& f)
{
  //F := c[5]*X^2 + c[4]*Y^2 + c[3]*X*Y + c[2]*X + c[1]*Y + c[0];
  // diff wrt X
  RT coeff[6] = {0, 0, 0, 0, 0, 0};
  coeff[2] = RT(2) * f[5];
  coeff[1] = f[3];
  coeff[0] = f[2];
  return MPoly_2_2<RT>(coeff);
}




template < typename RT >
MPoly_2_2<RT>
diff_y(const MPoly_2_2<RT>& f)
{
  //F := c[5]*X^2 + c[4]*Y^2 + c[3]*X*Y + c[2]*X + c[1]*Y + c[0];
  // diff wrt Y
  RT coeff[6] = {0, 0, 0, 0, 0, 0};
  coeff[2] = f[3];
  coeff[1] = RT(2) * f[4];
  coeff[0] = f[1];
  return MPoly_2_2<RT>(coeff);
}


template < typename POLY,
           typename RT >
POLY
specialize_x(const MPoly_2_2<RT>& f, const typename NumberTraits<RT>::FT& a)
{
  POLY p(3, AsSize());

  RT n = numerator(a);
  RT d = denominator(a);
  RT d_2 = SQR(d);
  p[2] = f[4] * d_2;
  p[1] = d * ( f[1] * d + f[3] * n );
  p[0] = (f[5] * n + f[2] * d ) * n + f[0] * d_2;

  return p;
}



template < typename POLY,
           typename RT >
POLY
specialize_y(const MPoly_2_2<RT>& f, const typename NumberTraits<RT>::FT& a)
{

  POLY p(3, AsSize());

  RT n = numerator(a);
  RT d = denominator(a);
  RT d_2 = SQR(d);
  p[2] = f[5] * d_2;
  p[1] = d * ( f[2] * d + f[3] * n );
  p[0] = (f[4] * n + f[1] * d) * n + f[0] * d_2;

  return p;
}




template < typename POLY,
           typename RT >
void
resultant_y_2_2(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  POLY A(2, AsSize());
  POLY B(3, AsSize());
  POLY R(5, AsSize());

  RT t8 =  g[3];
  RT t5 =  g[1];
  RT t69 =  ( t8 *  t5);
  RT t28 = -2 * t69;
  RT t13 =  f[0];
  RT t27 = -2 * t13;
  RT t15 =  g[0];
  RT t22 = 2 * t15;
  RT t1 =  g[4];
  RT t2 = f[1];
  RT t3 =  ( t1 * t2);
  RT t10 =  f[3];
  RT t17 =  ( t3 *  t10);
  RT t18 =  f[2];
  RT t29 =  ( t1 *  t1);
  RT t12 =  ( t29 *  t18);
  RT t11 =  (t2 *  t8 +  t10 *  t5);
  RT t14 = t1 * t13;
  RT t32 =  (t2 *  t5);
  RT t9 = -2 * t14 - t32;
  RT t19 = t1 * t18;
  RT t4 =  f[4];
  RT t20 =  g[2];
  RT t21 =  ( t4 *  t20);
  RT t7 = t21 + t19;
  RT t6 = 2 * t21 + t11;
  A[0] = (-t3 + t4 * t5) * t1;
  A[1] = (t4 * t8 - t1 * t10) * t1;
  RT t16 = t4 * t15;
  B[0] = (-t14 + t16) * t1;
  B[1] = (-t19 + t21) * t1;
  RT t23 =  f[5];
  RT t24 = t1 * t23;
  RT t25 =  g[5];
  RT t26 = t4 * t25;
  B[2] = -(t24 - t26) * t1;
  RT t37 =  (t2 * t2);
  RT t38 = t1 * t37;
  RT t40 = t4 * t4;
  RT t43 = t13 * t4;
  RT t44 = t5 * t5;
  R[0] = -t1 * (t29 * t13 * t13 - t14 * t32 + t43 * t44 + t9 * t16 + (t38 + t15 * t40) * t15);
  RT t52 = t1 * t4;
  RT t72 = t18 * t4;
  R[1] = -t1 * (2 * t13 * t12 + t17 * t22 + 2 * t43 * t69 + t72 * t44 - t7 * t32 - t11 * t14 + (t38 + t52 * t27 + t40 * t22) * t20 + (-2 * t19 - t11) * t16);
  RT t94 = t8 * t10;
  RT t96 = t10 * t10;
  RT t110 = t8 * t8;
  R[2] = t1 * (t72 * t28 - t9 * t26 + (-t44 * t23 - t110 * t13) * t4 + t11 * t21 + (-t20 * t20 - 2 * t15 * t25) * t40 - 2 * t20 * t17 + (-t15 * t96 - t25 * t37) * t1 + t6 * t19 + (-t18 * t18 + t23 * t27) * t29 + (t16 + t14) * t94 + (t32 + 2 * t16) * t24);
  RT t133 = t1 * t96;
  RT t144 = t23 * t4;
  R[3] = t1 * (-2 * t23 * t12 + t28 * t144 - t133 * t20 - t72 * t110 + t11 * t26 + t7 * t94 + 2 * (-t17 + t52 * t18 - t40 * t20) * t25 + t6 * t24);
  R[4] = -t1 * (t29 * t23 * t23 + (-t26 - t24) * t94 + t144 * t110 - 2 * t24 * t26 + (t133 + t40 * t25) * t25);

  AA = A;
  BB = B;
  RR = R;
  return;
}


template < typename POLY,
           typename RT >
void
resultant_y_2_1(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  // Resultant wrt to Y
  // The result is poly in X
  //cout << __FUNCTION__ << endl;

  POLY A(2, AsSize());
  POLY B(3, AsSize());
  POLY R(5, AsSize());

  A[0] =  g[1];
  A[1] =  g[3];
  RT t11 = -2 * A[1] * A[0];
  RT t7 =  f[4];
  B[2] =  g[5];
  RT t9 = -2 * B[2] * t7;
  B[1] = g[2];
  RT t6 = B[2] * A[0] + B[1] * A[1];
  B[0] = g[0];
  RT t3 = B[1] * A[0] + B[0] * A[1];
  RT t1 = pow( A[0],  2);
  RT t2 =  f[0];
  RT t4 = B[0] * A[0];
  RT t5 =  f[1];
  R[0] = -t1 * t2 + t4 * t5 - t7 * pow( B[0],  2);
  RT t15 =  f[3];
  RT t20 =  f[2];
  R[1] = -t20 * t1 + t2 * t11 + t4 * t15 - 2 * B[1] * t7 * B[0] + t3 * t5;
  RT t36 =  f[5];
  RT t39 = pow(A[1],  2);
  R[2] = B[0] * t9 + t20 * t11 - t1 * t36 - pow( B[1],  2) * t7 - t39 * t2 + t3 * t15 + t6 * t5;
  RT t41 = B[2] * A[1];
  R[3] = t41 * t5 + B[1] * t9 - t39 * t20 + t36 * t11 + t6 * t15;
  R[4] = -t39 * t36 - pow( B[2],  2) * t7 + t41 * t15;

  AA = A;
  BB = B;
  RR = R;
  return;
}


template < typename POLY,
           typename RT >
void
resultant_y_1_2(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  //cout << __FUNCTION__ << endl;
  resultant_y_2_1<POLY>(g, f, AA, BB, RR);
  return;
}



template < typename POLY,
           typename RT >
void
resultant_y_2_0(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  // Resultant wrt to Y
  // The result is poly in X
  //cout << __FUNCTION__ << endl;
  POLY A(2, AsSize());
  POLY B(3, AsSize());
  POLY R(5, AsSize());

  B[0] = g[0];
  B[1] = g[2];
  B[2] = g[5];
  RT t1 = f[4];
  R[0] = -f[4] * pow(B[0], 2);
  R[1] = -2 * f[4] * B[0] * B[1];
  R[2] = -t1 * (2 * B[0] * B[2] + pow(B[1], 2));
  R[3] = -2 * f[4] * B[1] * B[2];
  R[4] = -f[4] * pow(B[2], 2);

  AA = A;
  BB = B;
  RR = R;
  return;
}



template < typename POLY,
           typename RT >
void
resultant_y_0_2(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  return resultant_y_2_0<POLY>(g, f, AA, BB, RR);
}


template < typename POLY,
           typename RT >
void
resultant_y_1_1(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  // Resultant wrt to Y
  // The result is poly in X
  //cout << __FUNCTION__ << endl;
  POLY A(2, AsSize());
  POLY B(3, AsSize());
  POLY R(5, AsSize());

  A[1] =  g[3];
  A[0] =  g[1];
  RT t5 = -2 * A[1] * A[0];
  B[1] =  g[2];
  B[2] =  g[5];
  RT t4 = A[1] * B[1] + A[0] * B[2];
  B[0] = g[0];
  RT t2 = A[1] * B[0] + A[0] * B[1];
  RT t1 =  f[0];
  RT t3 =  f[1];
  R[0] = (-A[0] * t1 + B[0] * t3) * A[0];
  RT t11 = pow(A[0],  2);
  RT t12 =  f[2];
  RT t17 =  f[3];
  R[1] = t1 * t5 - t11 * t12 + A[0] * B[0] * t17 + t2 * t3;
  RT t21 = pow( A[1],  2);
  RT t27 =  f[5];
  R[2] = t12 * t5 - t21 * t1 - t11 * t27 + t4 * t3 + t2 * t17;
  R[3] = t27 * t5 - t21 * t12 + A[1] * B[2] * t3 + t4 * t17;
  R[4] = -(A[1] * t27 - B[2] * t17) * A[1];

  AA = A;
  BB = B;
  RR = R;
  return;
}



template < typename POLY,
           typename RT >
void
resultant_y_1_0(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  // Resultant wrt to Y
  // The result is poly in X
  //cout << __FUNCTION__ << endl;
  POLY A(2, AsSize());
  POLY B(3, AsSize());
  POLY R(5, AsSize());

  A[0] = f[1];
  A[1] = f[3];
  B[0] = f[0];
  B[1] = f[2];
  B[2] = f[5];
  R[0] = g[0];
  R[1] = g[2];
  R[2] = g[5];

  AA = A;
  BB = B;
  RR = R;
  return;
}



template < typename POLY,
           typename RT >
void
resultant_y_0_1(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  return resultant_y_1_0<POLY>(g, f, AA, BB, RR);
}


template < typename POLY,
           typename RT >
void
resultant_y_0_0(const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& AA, POLY& BB, POLY& RR)
{
  POLY A(2, AsSize());
  POLY B(3, AsSize());
  POLY R(5, AsSize());

  B[0] = f[0];
  B[1] = f[2];
  B[2] = f[5];

  R[0] = g[0];
  R[1] = g[2];
  R[2] = g[5];

  AA = A;
  BB = B;
  RR = R;
  return;
}


template < typename POLY,
           typename RT >
void
resultant_y(const MPoly_2_2<RT>& f1, const MPoly_2_2<RT>& f2,
            POLY& A, POLY& B, POLY& R)
{
  typedef void (* ResultantFunc)
    (const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& A, POLY& B, POLY& R);

  const ResultantFunc restab_[3][3] = {
    { resultant_y_0_0<POLY, RT>, resultant_y_0_1<POLY, RT>, resultant_y_0_2<POLY, RT> },
    { resultant_y_1_0<POLY, RT>, resultant_y_1_1<POLY, RT>, resultant_y_1_2<POLY, RT> },
    { resultant_y_2_0<POLY, RT>, resultant_y_2_1<POLY, RT>, resultant_y_2_2<POLY, RT> }
  };
  //cout << __FUNCTION__ << endl;

  restab_[f1.degree(1)][f2.degree(1)](f1, f2, A, B, R);
  return;
}




template < typename POLY,
           typename RT >
POLY
resultant_y(const MPoly_2_2<RT>& f1, const MPoly_2_2<RT>& f2)
{
  POLY A;
  POLY B;
  POLY R;
  resultant_y(f1, f2, A, B, R);
  return R;
}




template < typename POLY,
           typename RT >
void
resultant_x(const MPoly_2_2<RT>& f1, const MPoly_2_2<RT>& f2,
            POLY& A, POLY& B, POLY& R)
{
  typedef void (* ResultantFunc)
    (const MPoly_2_2<RT>& f, const MPoly_2_2<RT>& g, POLY& A, POLY& B, POLY& R);

  const ResultantFunc restab_[3][3] = {
    { resultant_y_0_0<POLY, RT>, resultant_y_0_1<POLY, RT>, resultant_y_0_2<POLY, RT> },
    { resultant_y_1_0<POLY, RT>, resultant_y_1_1<POLY, RT>, resultant_y_1_2<POLY, RT> },
    { resultant_y_2_0<POLY, RT>, resultant_y_2_1<POLY, RT>, resultant_y_2_2<POLY, RT> }
  };
  //cout << __FUNCTION__ << endl;

  MPoly_2_2<RT> f(f1[4], f1[5], f1[3], f1[1], f1[2], f1[0]);
  MPoly_2_2<RT> g(f2[4], f2[5], f2[3], f2[1], f2[2], f2[0]);

  restab_[f.degree(1)][g.degree(1)](f, g, A, B, R);
  return;
}



template < typename POLY,
           typename RT >
POLY
resultant_x(const MPoly_2_2<RT>& f1, const MPoly_2_2<RT>& f2)
{
  POLY A;
  POLY B;
  POLY R;
  resultant_x(f1, f2, A, B, R);
  return R;

}

__ALGEBRAIC_END_NAMESPACE




__BEGIN_NAMESPACE_SYNAPS

// Conversions
namespace let {

  // Naive implementation. We have to change it in the near future!
  template < typename RT, typename O, typename R >
  MPol<RT, O, R>
  convert(const ALGEBRAIC::MPoly_2_2<RT>& f, type::As<MPol<RT, O, R> >)
  {
    typedef typename MPol<RT, O, R>::monom_t      Monom;

    MPol<RT, O, R>  h;
    MPol<RT, O, R> h5 =  f[5] *  Monom("x", 2);
    MPol<RT, O, R> h4 = f[4] *  Monom("y", 2);
    MPol<RT, O, R> h3x = f[3] *  Monom("x", 1) * Monom("y", 1);
    MPol<RT, O, R> h3y = Monom("y", 1);
    MPol<RT, O, R> h2 = f[2] *  Monom("x", 1);
    MPol<RT, O, R> h1 = f[1] *  Monom("y", 1);
    MPol<RT, O, R> h0 = f[0] *  Monom("x", 0);

    h = h0 + h1 + h2 + (h3x) + h4 + h5;
    return h;
  }

  template < typename RT, typename O, typename R >
    ALGEBRAIC::MPoly_2_2<RT> convert(const MPol<RT, O, R>& f, type::As< ALGEBRAIC::MPoly_2_2<RT> > )
  {
    typedef typename MPol<RT, O, R>::monom_t      Monom;

    ALGEBRAIC::MPoly_2_2<RT> h;
    h.set_coeff(5, MPOLDST::coeffof(f, Monom("x", 2)));
    h.set_coeff(4, MPOLDST::coeffof(f, Monom("y", 2)));
    h.set_coeff(3, MPOLDST::coeffof(f, Monom("x", 1) * Monom("y", 1)));
    h.set_coeff(2, MPOLDST::coeffof(f, Monom("x", 1)));
    h.set_coeff(1, MPOLDST::coeffof(f, Monom("y", 1)));
    h.set_coeff(0, MPOLDST::coeffof(f, Monom("x", 0)));

    return h;
  }
}

__END_NAMESPACE_SYNAPS

#endif // SYNAPS_ARITHM_ALGEBRAIC_MPOLY_2_2_H