/Users/mourrain/Devel/mmx/algebramix/src/matrix_integer.cpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_integer.cpp
* DESCRIPTION: Multi-modular multiplication of integer matrices
* 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.
******************************************************************************/

#include <numerix/integer.hpp>

namespace mmx {

xnat
matrix_product_bit_size (const integer* s1, nat s1_rs, nat s1_cs,
                         const integer* s2, nat s2_rs, nat s2_cs,
                         nat r, nat l, nat c)
{
  double count= 1.0;
  xnat   sz   = 0;
  // NOTE: the largest integer in the result will be bounded by
  // count * 2^sz at the end. The bound is computed carefully,
  // so as to be both fast and relatively sharp.
  for (nat k=0; k<l; k++) {
    xnat sz1= 0, sz2= 0;
    const integer* ss1= s1 + k * s1_cs;
    const integer* ss2= s2 + k * s2_rs;
    for (nat i=0; i<r; i++, ss1 += s1_rs)
      sz1= max (sz1, bit_size (*ss1));
    for (nat j=0; j<c; j++, ss2 += s2_cs)
      sz2= max (sz2, bit_size (*ss2));
    xnat szs= sz1 + sz2;
    if (szs > sz) {
      if (szs >= sz + 30) {
        count= (count / 1.0e9) + 1.000000001;
        sz   = szs;
      }
      else {
        count= (count / ((double) (1 << (szs - sz)))) + 1.000000001;
        sz   = szs;
      }
    }
    else {
      if (sz >= szs + 30) count += 1.0e-9;
      else count += 1.000000001 / ((double) (1 << (sz - szs)));
    }
  }
  return sz + ((xnat) logb (count)) + 1;
}

} // namespace mmx