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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_quotient.hpp
* DESCRIPTION: Multiply matrices by putting entries under a common denominator
* COPYRIGHT  : (C) 2008 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_MATRIX_QUOTIENT_HPP
#define __MMX_MATRIX_QUOTIENT_HPP
#include <numerix/rational.hpp>
#include <algebramix/matrix_integer.hpp>
namespace mmx {
#define TMPL template<typename C>
#define NC Numerator_type(C)
#define DC Denominator_type(C)

/******************************************************************************
* Quotient matrix variant
******************************************************************************/

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

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

DEFINE_VARIANT (matrix_rational,
                matrix_quotient<matrix_naive>)

STMPL
struct matrix_variant_helper<rational> {
  typedef matrix_rational MV;
};

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

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

TMPL void
denominator_mul (NC* dest, const C* src, nat str, nat n, const DC& d) {
  for (; n != 0; n--, dest += str, src += str)
    *dest= numerator (*src) * (d / denominator (*src));
}

TMPL void
denominator_div (C* dest, const NC* src,
                 nat row_str, nat rows,
                 nat col_str, nat cols,
                 const DC* row_d, const DC* col_d)
{ 
  for (nat i=0; i < rows; i++, dest += row_str, src += row_str, row_d++) {
    C*        Dest = dest;
    const NC* Src  = src;
    const DC* Col_d= col_d;
    for (nat j=0; j < cols; j++, Dest += col_str, Src += col_str, Col_d++)
      *Dest= C (*Src) / ((*row_d) * (*Col_d));
  }
}

/******************************************************************************
* Multiply matrices by putting the entries under a common denominator
******************************************************************************/

template<typename V, typename W>
struct implementation<matrix_multiply,V,matrix_quotient<W> >:
  public implementation<matrix_linear,V>
{
  static const nat thr= 7;

TMPL static void
mul (C* d, const C* s1, const C* s2,
     nat r, nat l, nat c) 
{
  typedef typename Matrix_variant(Numerator_type(C)) NV;
  typedef implementation<vector_linear,NV> Vec_N;
  typedef implementation<vector_linear,NV> Vec_D;
  typedef implementation<matrix_multiply,V> Mat;
  typedef implementation<matrix_multiply,NV> Mat_N;
  typedef implementation<matrix_multiply,NV> Mat_D;
  typedef implementation<matrix_multiply,typename V::Naive> Naive_Mat;

  if (r <= thr || l <= thr || c <= thr)
    Naive_Mat::template mul<mul_op> (d, s1, s2, r, r, l, l, c, c);
  else {
    nat sd1  = aligned_size<DC,NV> (r);
    nat sd2  = aligned_size<DC,NV> (c);
    DC* d1   = mmx_new<DC> (sd1 + sd2);
    DC* d2   = d1 + sd1;
    nat sSrc1= aligned_size<NC,NV> (r*l);
    nat sSrc2= aligned_size<NC,NV> (c*l);
    nat sDest= aligned_size<NC,NV> (r*c);
    NC* Src1 = mmx_new<NC> (sSrc1 + sSrc2 + sDest);
    NC* Src2 = Src1 + sSrc1;
    NC* Dest = Src2 + sSrc2;
    for (nat i=0; i<r; i++) 
      d1[i] = least_common_denominator (s1 + Mat::index (i, 0, r, l),
                                        Mat::index (0, 1, r, l), l);
    //mmout << "d1= "; Vec_D::print (mmout, d1, r); mmout << "\n";
    for (nat j=0; j<c; j++) 
      d2[j] = least_common_denominator (s2 + Mat::index (0, j, l, c),
                                        Mat::index (1, 0, l, c), l);
    //mmout << "d2= "; Vec_D::print (mmout, d2, c); mmout << "\n";

    for (nat i=0; i<r; i++) 
      denominator_mul 
        (Src1 + Mat::index (i, 0, r, l), s1 + Mat::index (i, 0, r, l),
         Mat::index (0, 1, r, l), l, d1[i]);
    //mmout << "Src1= "; Vec_N::print (mmout, Src1, r*l); mmout << "\n";
    for (nat j=0; j<c; j++) 
      denominator_mul 
        (Src2 + Mat::index (0, j, l, c), s2 + Mat::index (0, j, l, c),
         Mat::index (1, 0, l, c), l, d2[j]);
    //mmout << "Src2= "; Vec_N::print (mmout, Src1, l*c); mmout << "\n";

    Mat_N::mul (Dest, Src1, Src2, r, l, c);
    //mmout << "Dest= "; Vec_N::print (mmout, Dest, r*c); mmout << "\n";

    denominator_div (d, Dest,
                     Mat::index (1, 0, r, c), r,
                     Mat::index (0, 1, r, c), c,
                     d1, d2);

    mmx_delete<DC> (d1, sd1 + sd2);
    mmx_delete<NC> (Src1, sSrc1 + sSrc2 + sDest);
  }
}

}; // implementation<matrix_multiply,V,matrix_quotient<NV> >

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