include/algebramix/polynomial_quotient.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_quotient.hpp
* DESCRIPTION: Multiply polynomials by reduction to a common denominator
* COPYRIGHT  : (C) 2004  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_POLYNOMIAL_QUOTIENT_HPP
#define __MMX_POLYNOMIAL_QUOTIENT_HPP
#include <algebramix/quotient.hpp>
#include <algebramix/polynomial_dicho.hpp>
namespace mmx {
#define TMPL template<typename C>
#define NC Numerator_type(C)
#define DC Denominator_type(C)

/******************************************************************************
* Variants for polynomials with fractional coefficients
******************************************************************************/

template<typename V>
struct polynomial_quotient: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef polynomial_quotient<typename V::Positive> Positive;
  typedef polynomial_quotient<typename V::No_simd> No_simd;
  typedef polynomial_quotient<typename V::No_thread> No_thread;
  typedef polynomial_quotient<typename V::No_scaled> No_scaled;
};

template<typename F, typename V, typename W>
struct implementation<F,V,polynomial_quotient<W> >:
    public implementation<F,V,W> {};

/******************************************************************************
* Default variant
******************************************************************************/

template<typename C, typename V>
struct polynomial_variant_helper<quotient<polynomial<C,V>,
                                          polynomial<C,V> > > {
  typedef polynomial_quotient<polynomial_naive> PV;
};

/******************************************************************************
* Variable names, default is C(x)[y]
******************************************************************************/

template<typename V,typename W>
struct implementation<polynomial_defaults,V,polynomial_quotient<W> > {

  template<typename P>
  class global_variables {
    static inline generic& dyn_name () {
      static generic name = "x";
      return name; }
  public:
    static inline void set_variable_name (const generic& x) { dyn_name () = x; }
    static inline generic get_variable_name () { return dyn_name (); }
  };

  template<typename C, typename V1, typename V2>
  class global_variables<polynomial<
                           quotient<polynomial<C,V1>,
                                    polynomial<C,V1> >, V2> > {
    static inline generic& dyn_name () {
      static generic name = "y";
      return name; }
  public:
    static inline void set_variable_name (const generic& x) {
      dyn_name () = x; }
    static inline generic get_variable_name () {
      return dyn_name (); }
  };
};

/******************************************************************************
* Low level subroutines
******************************************************************************/

TMPL DC
least_common_denominator (const C* src, nat n) {
  switch (n) {
  case 0: return 1;
  case 1: return denominator (src[0]);
  case 2: return lcm (denominator (src[0]), denominator (src[1]));
  case 3: return lcm (denominator (src[0]), lcm (denominator (src[1]),
                                                 denominator (src[2])));
  default: return lcm (least_common_denominator (src, n>>1),
                       least_common_denominator (src + (n>>1), n - (n>>1)));
  }
}

TMPL void
denominator_mul (NC* dest, const C* src, nat n, const DC& d) {
  for (nat i=0; i<n; i++)
    dest[i]= numerator (src[i]) * (d / denominator (src[i]));
}

TMPL void
denominator_div (C* dest, const NC* src, nat n, const DC& d) {
  for (nat i=0; i<n; i++)
    dest[i]= C (src[i]) / d;
}

/******************************************************************************
* Multiply polynomials by putting the coefficients under a common denominator
******************************************************************************/

template<typename V, typename W>
struct implementation<polynomial_multiply,V,polynomial_quotient<W> >:
  public implementation<polynomial_linear,V>
{
  typedef polynomial_multiply_threshold<polynomial_quotient<W> > Th;
  typedef implementation<polynomial_multiply,W> Fallback;

TMPL static void
mul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  typedef typename Polynomial_variant(NC) NV;
  typedef implementation<polynomial_multiply,NV> Pol_N;
  if (min (n1, n2) < Threshold(C,Th))
    Fallback::mul (dest, s1, s2, n1, n2);
  else {
    DC  d1  = least_common_denominator (s1, n1);
    DC  d2  = least_common_denominator (s2, n2);
    /*
      mmout << "d1= " << d1 << ", " << "d2= " << d2 << "\n";
      for (nat i=0; i<n1; i++)
      mmout << "  s" << i << "\t" << s1[i] << "\n";
      for (nat i=0; i<n2; i++)
      mmout << "  t" << i << "\t" << s2[i] << "\n";
    */
    NC* Src1= mmx_new<NC> (n1);
    NC* Src2= mmx_new<NC> (n2);
    NC* Dest= mmx_new<NC> (n1+n2-1);
    denominator_mul (Src1, s1, n1, d1);
    denominator_mul (Src2, s2, n2, d2);
    Pol_N::mul (Dest, Src1, Src2, n1, n2);
    denominator_div (dest, Dest, n1+n2-1, d1*d2);
    mmx_delete<NC> (Src1, n1);
    mmx_delete<NC> (Src2, n2);
    mmx_delete<NC> (Dest, n1+n2-1);
  }
}

TMPL static void
square (C* dest, const C* s1, nat n1) {
  typedef typename Polynomial_variant(NC) NV;
  typedef implementation<polynomial_multiply,NV> Pol_N;
  if (n1 < Threshold(C,Th))
    Fallback::square (dest, s1, n1);
  else {
    DC  d1  = least_common_denominator (s1, n1);
    NC* Src1= mmx_new<NC> (n1);
    NC* Dest= mmx_new<NC> (n1+n1-1);
    denominator_mul (Src1, s1, n1, d1);
    Pol_N::square (Dest, Src1, n1);
    denominator_div (dest, Dest, n1+n1-1, square_op::op (d1));
    mmx_delete<NC> (Src1, n1);
    mmx_delete<NC> (Dest, n1+n1-1);
  }
}

TMPL static inline void
tmul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  (void) dest; (void) s1; (void) s2; (void) n1; (void) n2;
  ERROR ("transposed quotient multiplication not implemented");
}

}; // implementation<polynomial_multiply,V,polynomial_quotient<NV> >

#undef TMPL
#undef NC
#undef DC
}
#endif // __MMX_POLYNOMIAL_QUOTIENT_HPP

---