include/lacunaryx/lpolynomial.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : lpolynomial.hpp
* DESCRIPTION: Univariate lacunary polynomials
* COPYRIGHT  : (C) 2014  Gregoire Lecerf
*******************************************************************************
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
******************************************************************************/

#ifndef __MMX_LPOLYNOMIAL_HPP
#define __MMX_LPOLYNOMIAL_HPP
#include <basix/vector_sort.hpp>
#include <lacunaryx/lpolynomial_naive.hpp>
namespace mmx {

#define LPolynomial_variant(C,E) lpolynomial_variant_helper<C,E>::PV
#define TMPL_DEF \
  template<typename C, typename E=integer,\
           typename V=typename LPolynomial_variant(C,E) >
#define TMPL template<typename C, typename E, typename V>
#define TMPLK template<typename C, typename E, typename V, typename K>
#define SE typename signed_of_helper<E>::type
#define Format format<C>
#define Table table<C,E>
#define LPolynomial lpolynomial<C,E,V>
#define Abs_polynomial lpolynomial<Abs_type(C),E>
#define Real_polynomial lpolynomial<Real_type(C),E>
#define Center_polynomial lpolynomial<Center_type(C),E>
#define Radius_polynomial lpolynomial<Radius_type(C),E>
#define Lifted_polynomial lpolynomial<Lift_type(C),E>
#define Projected_polynomial lpolynomial<Project_type(C),E>
#define Reconstructed_polynomial lpolynomial<Reconstruct_type(C),E>
#define Bnd Bound_type(C)
TMPL class lpolynomial;
TMPL inline SE deg (const LPolynomial& P);
TMPL inline int degree (const LPolynomial& P);
TMPL inline nat N (const LPolynomial& P);

/******************************************************************************
* Lacunary univariate polynomials
******************************************************************************/

TMPL_DEF
class lpolynomial {
  MMX_ALLOCATORS;
public:
  Table rep;
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<lpolynomial_defaults,V> Pol;
  typedef typename Pol::template global_variables<LPolynomial> S;
public:
  inline Table operator * () const { return rep; }

  static inline generic get_variable_name () {
    return S::get_variable_name (); }
  static inline void set_variable_name (const generic& x) {
    S::set_variable_name (x); }

  inline lpolynomial () : rep (as<C> (0)) {}
  inline lpolynomial (const Table& t) : rep (t) {}
  inline lpolynomial (const Format& fm) : rep (promote (0, fm)) {}

  lpolynomial (const C& x, const E n= 0) {
    rep= Table (promote (0, get_format (x)));
    rep[n]= x; }

  inline friend LPolynomial copy (const LPolynomial& P) {
    return LPolynomial (copy (*P)); }
  inline C operator[] (const E& i) const {
    return rep[i]; }
  LPolynomial& operator += (const LPolynomial& P);
  LPolynomial& operator -= (const LPolynomial& P);
  LPolynomial& operator *= (const C& c);
  LPolynomial& operator /= (const C& c);
};

TMPL
struct binary_return_type_helper<access_op,LPolynomial,E> {
  typedef C RET; };

TMPL inline nat N (const LPolynomial& P) {
  return N(*P); }

TMPL inline Format CF (const LPolynomial& P) {
  return CF1(*P); }

TMPL nat hash (const LPolynomial& c) { return hash (*c); }
TMPL nat exact_hash (const LPolynomial& c) { return exact_hash (*c); }
TMPL nat hard_hash (const LPolynomial& c) { return hard_hash (*c); }
TMPL bool operator == (const LPolynomial& c1, const LPolynomial& c2) {
  return (*c1) == (*c2); }  
TMPL bool operator != (const LPolynomial& c1, const LPolynomial& c2) {
  return (*c1) != (*c2); }
TMPL bool exact_eq (const LPolynomial& c1, const LPolynomial& c2) {
  return exact_eq (*c1, *c2); }
TMPL bool exact_neq (const LPolynomial& c1, const LPolynomial& c2) {
  return exact_neq(*c1,*c2); }
TMPL bool hard_eq (const LPolynomial& c1, const LPolynomial& c2) {
  return hard_eq (*c1, *c2); }  
TMPL bool hard_neq (const LPolynomial& c1, const LPolynomial& c2) {
 return hard_neq (*c1, *c2); }

DEFINE_UNARY_FORMAT_3 (lpolynomial);
STYPE_TO_TYPE(TMPL,scalar_type,LPolynomial,C);
STYPE_TO_TYPE(TMPL,monomial_type,LPolynomial,E);
UNARY_RETURN_TYPE(TMPL,abs_op,LPolynomial,Abs_polynomial);
UNARY_RETURN_TYPE(TMPL,ground_op,LPolynomial,Ground_type(C));
UNARY_RETURN_TYPE(TMPL,Re_op,LPolynomial,Real_polynomial);
UNARY_RETURN_TYPE(TMPL,center_op,LPolynomial,Center_polynomial);
UNARY_RETURN_TYPE(TMPL,radius_op,LPolynomial,Radius_polynomial);

/******************************************************************************
* Degree, valuation, and access
******************************************************************************/

TMPL inline SE
deg (const LPolynomial& P) {
  SE ret= -1;
  for (iterator<pair<E,C> > it= iterate (*P); busy (it); ++it) {
    if (cdr (*it) == 0) continue;
    ret= max ((SE) car (*it), ret);
  }
  return ret; }

TMPL inline SE
val (const LPolynomial& P) {
  if (N(P) == 0) return -1;
  iterator<pair<E,C> > it= iterate (*P);
  while (cdr (*it) == 0 && busy (it)) ++it; 
  if (!busy (it)) return -1;
  SE ret= car (*it); ++it;
  for (; busy (it); ++it) {
    if (cdr (*it) == 0) continue;
    ret= min ((SE) car (*it), ret);
  }
  return ret; }

TMPL inline SE
degree (const LPolynomial& P) {
  return deg (P); }

TMPL inline vector<C>
coefficients (const LPolynomial& P) {
  vector<C> ret (CF(P));
  for (iterator<pair<E,C> > it= iterate (*P); busy (it); ++it)
    if (cdr (*it) != 0) ret << cdr (*it);
  return ret; }

TMPL inline Ground_type(C)
ground (const LPolynomial& P) {
  return ground (P[0]); }

/******************************************************************************
* Casts
******************************************************************************/

template<typename T, typename F, typename TE, typename FE, 
         typename TV, typename FV> inline void
set_as (lpolynomial<T,TE,TV>& r, const lpolynomial<F,FE,FV>& p) {
  set_as (r.rep, p.rep); }

template<typename C, typename E, typename V, typename T> inline void
set_as (LPolynomial& r, const T& x) {
  r.rep= Table (promote (0, CF(r)));
  set_as (r.rep[0], x); }

template<typename T, typename F,
         typename TE, typename FE, 
         typename TV, typename FV>
struct as_helper<lpolynomial<T,TE,TV>,lpolynomial<F,FE,FV> > {
  static inline lpolynomial<T,TE,TV>
  cv (const lpolynomial<F,FE,FV>& p) {
    table<F,FE> f= *p;
    table<T,TE> r (as<T> (I (*p)), CF2(*p));
    for (iterator<FE> it= entries (f); busy (it); ++it)
      r[as<TE> (*it)]= as<T> (f[*it]);
    return r; }
};

/******************************************************************************
* Iterate
******************************************************************************/

TMPL iterator<pair<E,C> >
iterate (const LPolynomial& P) {
  return iterate (*P); }

TMPL vector<E>
supp (const LPolynomial& P) {
  vector<E> ret;
  for (iterator<pair<E,C> > it= iterate (*P); busy (it); ++it)
    if (cdr (*it) != 0) ret << car (*it);
  return ret; }

/******************************************************************************
* Output routines for dense polynomials
******************************************************************************/

TMPL inline generic var (const LPolynomial& P) {
  (void) P; return LPolynomial::get_variable_name (); }

TMPL inline void set_variable_name (const LPolynomial& P, const generic& x) {
  (void) P; return LPolynomial::set_variable_name (x); }

TMPL syntactic
flatten (const LPolynomial& P, const syntactic& v) {
  vector<E> w (entries (*P));
  sort_leq<gtr_op> (w);
  syntactic s= 0;
  for (nat i= 0; i < N(w); i++)
    if (P[w[i]] != 0)
      s= s + flatten (P[w[i]]) * pow (v, flatten (w[i]));
  return s; }

TMPL syntactic
flatten (const LPolynomial& P) {
  return flatten (P, as_syntactic (var (P))); }

/******************************************************************************
* Comparisons
******************************************************************************/

TMPL int
sign (const LPolynomial& P) {
  vector<E> v (entries (*P));
  sort (v);
  for (nat i=0; i<N(v); i++) {
    int s= sign (P[v[i]]);
    if (s != 0) return s;
  }
  return 0; }

TMPL int
compare (const LPolynomial& P1, const LPolynomial& P2) {
  vector<E> v (entries (*P1 + *P2));
  sort (v);
  for (nat i=0; i<N(v); i++) {
    int s= sign (P1[v[i]] - P2[v[i]]);
    if (s != 0) return s;
  }
  return 0; }

template<typename Op,typename C,typename E,typename V> inline bool
unary_hash (const LPolynomial& P) {
  return unary_hash<Op> (*P); }

template<typename Op,typename C,typename E,typename V> inline bool
binary_test (const LPolynomial& P1, const LPolynomial& P2) {
  vector<E> v (entries (*P1 + *P2));
  sort (v);
  for (nat i=0; i<N(v); i++)
    if (!Op::op (P1[v[i]], P2[v[i]])) return false;
  return true; }

PR_IDENTITY_OP_SUGAR(TMPL,LPolynomial)
COMPARE_SUGAR(TMPL,LPolynomial)
EQUAL_SCALAR_SUGAR(TMPL,LPolynomial,C)
COMPARE_SCALAR_SUGAR(TMPL,LPolynomial,C)
EQUAL_INT_SUGAR(TMPL,LPolynomial)
COMPARE_INT_SUGAR(TMPL,LPolynomial)
EQUAL_SCALAR_SUGAR_BIS(TMPL,LPolynomial,C)
COMPARE_SCALAR_SUGAR_BIS(TMPL,LPolynomial,C)

/******************************************************************************
* In place arithmetic operations
******************************************************************************/

TMPL LPolynomial&
LPolynomial::operator += (const LPolynomial& P) {
  rep += *P; return *this; }

TMPL LPolynomial&
LPolynomial::operator -= (const LPolynomial& P) {
  rep -= *P; return *this; }

TMPL LPolynomial&
LPolynomial::operator *= (const C& c) {
  rep *= c; return *this; }

TMPL LPolynomial&
LPolynomial::operator /= (const C& c) {
  rep /= c; return *this; }

/******************************************************************************
* Basic arithmetic operations
******************************************************************************/

TMPL LPolynomial
operator - (const LPolynomial& P) {
  return - *P; }

TMPL LPolynomial
operator + (const LPolynomial& P, const LPolynomial& Q) {
  return *P + *Q; }

TMPL LPolynomial
operator + (const LPolynomial& P, const C& c) {
  Table ret (*P);
  ret[0] = ret[0] + c;
  return ret; }

TMPL LPolynomial
operator + (const C& c, const LPolynomial& P) {
  Table ret (*P);
  ret[0] = c + ret[0];
  return ret; }

TMPL LPolynomial
operator - (const LPolynomial& P, const LPolynomial& Q) {
  return *P - *Q; }

TMPL LPolynomial
operator - (const LPolynomial& P, const C& c) {
  Table ret (*P);
  ret[0] = ret[0] - c;
  return ret; }

TMPL LPolynomial
operator - (const C& c, const LPolynomial& P) {
  Table ret (*P);
  ret[0] = c - ret[0];
  return ret; }

TMPL LPolynomial
operator * (const LPolynomial& P, const LPolynomial& Q) {
  Table ret (0);
  for (iterator<pair<E,C> > itP= iterate (*P); busy (itP); ++itP)
    for (iterator<pair<E,C> > itQ= iterate (*Q); busy (itQ); ++itQ) {
      E i= car(*itP), j= car(*itQ);
      ret[i+j] += cdr (*itP) * cdr (*itQ);
    }
  return ret; }

TMPL LPolynomial
operator * (const LPolynomial& P, const C& c) {
  return (*P) * c; }

TMPL LPolynomial
operator * (const C& c, const LPolynomial& P) {
  return c * (*P); }

TMPL LPolynomial
square (const LPolynomial& P) {
  Table ret (0);
  for (iterator<pair<E,C> > itP= iterate (*P); busy (itP); ++itP)
    ret[2*car(*itP)] += square (cdr (*itP));
  nat l= 0;
  for (iterator<pair<E,C> > itP= iterate (*P); busy (itP); ++itP, l++) {
    nat k= 0;
    for (iterator<pair<E,C> > itQ= iterate (*P); busy (itQ) && k < l; ++itQ, k++) {
      E i= car(*itP), j= car(*itQ);
      ret[i+j] += 2 * cdr (*itP) * cdr (*itQ);
    }
  }
  return ret; }

TMPL LPolynomial
operator / (const LPolynomial& P, const LPolynomial& Q) {
  ASSERT (deg (Q) != -1, "division by zero");
  LPolynomial R= P;
  Table ret (0);
  SE d, dQ= deg (Q);
  C lcQ= Q[dQ];
  while ((d= deg (R)) >= dQ) {
    C c= R[d] / lcQ;
    ret[d - dQ]= c;
    R = R - LPolynomial (c, d - dQ) * Q;
  }
  return ret; }

TMPL LPolynomial
operator / (const LPolynomial& P, const C& c) {
  return (*P) / c; }

ARITH_SCALAR_INT_SUGAR(TMPL,LPolynomial)

/******************************************************************************
* Powering
******************************************************************************/

TMPL LPolynomial
__binpow (const LPolynomial& P, const integer& n) {
  if (n == 1) return P;
  if (n == 2) return square (P);
  if (n <= 16) {
    LPolynomial Q= P;
    for (integer i= 1; i < n; i+= 1) Q = Q * P;
    return Q;
  }
  LPolynomial ret (square (_binpow (P, quo (n, 2))));
  if (rem (n, 2) != 0) ret = P * ret;
  return ret; }

TMPL LPolynomial
_binpow (const LPolynomial& P, const integer& n) {
  integer m= as<integer> (deg (P)) * n;
  if (bit_size (m) <= 8*sizeof(unsigned long long int)) {
    lpolynomial<C,unsigned long long int> Q=
      as<lpolynomial<C,unsigned long long int> > (P);
    return as<LPolynomial> (__binpow (Q, n));
  }
#if defined(BASIX_HAVE_INT128)
  if (bit_size (m) <= 8*sizeof(uint128)) {
    lpolynomial<C,uint128> Q= as<lpolynomial<C,uint128> > (P);
    return as<LPolynomial> (__binpow (Q, n));
  }
#endif
  return __binpow (P, n); }

TMPL LPolynomial
binpow (const LPolynomial& Q, const integer& n) {
  LPolynomial P= Q; simplify (P.rep);
  SE d= deg (P);
  if (n == 0) return LPolynomial (promote (1, CF(Q)), 0);
  if (d < 0 || n == 1) return P;
  if (N(P) == 1) return LPolynomial (P[d], as<E> (d * n));
  return _binpow (P, n); }

/******************************************************************************
* Usual operations
******************************************************************************/

TMPL LPolynomial
derive (const LPolynomial& P) {
  Table ret (0);
  for (iterator<pair<E,C> > it= iterate (*P); busy (it); ++it) {
    if (car (*it) == 0) continue;
    ret[car (*it) - 1] = promote (car (*it), cdr (*it)) * cdr (*it);
  }
  return ret; }

TMPL LPolynomial
xderive (const LPolynomial& P) {
  Table ret (0);
  for (iterator<pair<E,C> > it= iterate (*P); busy (it); ++it) {
    if (car (*it) == 0) continue;
    ret[car (*it)] = promote (car (*it), cdr (*it)) * cdr (*it);
  }
  return ret; }

TMPL C
eval (const LPolynomial& P, const C& c) {
  C ret= promote (0, CF (P));
  for (iterator<pair<E,C> > it= iterate (*P); busy (it); ++it) {
    ret += cdr (*it) * binpow (c, car (*it));
  }
  return ret; }

#undef TMPL_DEF
#undef TMPL
#undef TMPLK
#undef SE
#undef Format
#undef Table
#undef LPolynomial
#undef Abs_polynomial
#undef Real_polynomial
#undef Center_polynomial
#undef Radius_polynomial
#undef Lifted_polynomial
#undef Projected_polynomial
#undef Reconstructed_polynomial
#undef Bnd
} // mmx namespace
#endif // __MMX_LPOLYNOMIAL_HPP

---