/Users/mourrain/Devel/mmx/algebramix/include/algebramix/kronecker_polynomial.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : kronecker_polynomial.hpp
* DESCRIPTION: Multiplication of integer polynomials
* COPYRIGHT  : (C) 2009  Joris van der Hoeven and Gregoire Lecerf
*******************************************************************************
* 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_KRONECKER_POLYNOMIAL_HPP
#define __MMX_KRONECKER_POLYNOMIAL_HPP
#include <algebramix/polynomial.hpp>

namespace mmx {
#define TMPL template<typename C, typename V>
#define Polynomial polynomial<C,V>

/******************************************************************************
* Low level Kronecker multiplication for bivariate polynomials
******************************************************************************/

TMPL nat
max_polynomial_size (const Polynomial* src, nat n) {
  nat m= 0;
  for (nat i= 0; i < n; i++) m= max (N(src[i]), m);
  return m;
}

TMPL void
encode_kronecker (Polynomial& dest, const Polynomial* src, nat n, nat m) {
  if (n == 0) return;
  nat p= n * m;
  nat l= aligned_size<C,V> (p);
  C* x= mmx_new<C> (l); C* y= x;
  for (nat i= 0; i < n; i++, x += m, src++)
    for (nat j= 0; j < m; j++)
      x[j]= (*src)[j];
  dest= Polynomial (y, p, l, CF(src[0]));
}

TMPL void
decode_kronecker (Polynomial* dest, const Polynomial& src, nat n, nat m) {
  typedef implementation<polynomial_linear,V> Pol;
  if (n == 0) return;
  nat l= aligned_size<C,V> (m), i, j; C* x;
  const C* y= seg (src);
  for (i= 0, j= 0; i + 1 < n && j + m < N(src); i++, j += m, y += m) {
    x= mmx_new<C> (l);
    Pol::copy (x, y, m);
    dest[i]= Polynomial (x, m, l, CF(src));
  }
  if (i < n && N(src) > i * m) {
    m= N(src) - i * m;
    l= aligned_size<C,V> (m);
    x= mmx_new<C> (l);
    Pol::copy (x, y, m);
    dest[i]= Polynomial (x, m, l, CF(src));
    i++;
  }
  for (; i < n; i++) dest[i]= Polynomial (C(0));
}

TMPL inline void
mul_kronecker (Polynomial* dest,
               const Polynomial* s1, nat n1,
               const Polynomial* s2, nat n2) {
  typedef implementation<polynomial_linear,V> Pol;
  if (n1 == 0 || n2 == 0) return;
  if (n1 == 1) { Pol::mul_sc (dest, s2, s1[0], n2); return; }
  if (n2 == 1) { Pol::mul_sc (dest, s1, s2[0], n1); return; }
  nat m1= max_polynomial_size (s1, n1);
  nat m2= max_polynomial_size (s2, n2);
  nat m = m1 + m2 - 1;
  Polynomial x1, x2, y;
  encode_kronecker (x1, s1, n1, m);
  encode_kronecker (x2, s2, n2, m);
  y= x1 * x2;
  decode_kronecker (dest, y, n1 + n2 - 1, m); }

TMPL inline void
square_kronecker (Polynomial* dest, const Polynomial* s, nat n) {
  typedef implementation<polynomial_linear,V> Pol;
  if (n == 0) return;
  if (n == 1) { dest[0]= square_op::op (s[0]); return; }
  nat m = (max_polynomial_size (s, n) << 1) - 1;
  Polynomial x;
  encode_kronecker (x, s, n, m);
  x= square_op::op (x);
  decode_kronecker (dest, x, (n << 1) - 1, m); }

TMPL inline void
div_kronecker (Polynomial* dest,
               const Polynomial* s1, nat n1,
               const Polynomial* s2, nat n2) {
  typedef implementation<polynomial_linear,V> Pol;
  ASSERT (n2 != 0, "division by zero");
  if (n1 == 0) return;
  nat m1= max_polynomial_size (s1, n1);
  nat m2= max_polynomial_size (s2, n2);
  nat n= n1 - n2 + 1;
  Polynomial x1, x2, y;
  encode_kronecker (x1, s1, n1, m1);
  encode_kronecker (x2, s2, n2, m1);
  y= x1 / x2;
  nat l= default_aligned_size<C> (n1);
  Polynomial* tmp= mmx_new<Polynomial> (l);
  decode_kronecker (tmp, y, n1, m1);
  while (n1 > 0 && is_exact_zero (tmp[n1-1])) n1--;
  nat m= max_polynomial_size (tmp, n1);
  if (n1 <= n && m1 != m + m2 - 1) {
    mmx_delete<Polynomial> (tmp, l);
    ERROR ("unexact division");
  }
  Pol::copy (dest, tmp, n);
  mmx_delete<Polynomial> (tmp, l); }

#undef Polynomial
#undef TMPL
} // namespace mmx
#endif // __MMX_KRONECKER_POLYNOMIAL_HPP