Re: do you support MTL?

Hi Peter,

Thank you for your email and resolving my problems.  I implemented the
operators for the polynomial type, and now I'm able to multiply
vectors and matrixes of polynomials.  I'm attaching a complete program
at the bottom of this email.

Because I have the basic functionality going (a product of polynomial
vectors and matrixes), I think I will plugh ahead harnessing MTL.

I'm looking forward to the release of the new MTL.


Best,
Irek

Peter Gottschling wrote:

 > Hi Irek,
 >
 > Thank you for your interest in MTL and sorry for the trouble you run
 > into.  The problem is that nobody developed MTL for some years.
 > When we restarted it (with different authors), we saw many things
 > principally different.  Therefore, we decided to restart from the
 > beginning (which unfortunately goes slower than expected).
 > Therefore I tend to keep the maintenance efforts on MTL 2.xx minimal
 > (to not delay the delivery of the new version even further).
 >
 > However, the idea of polynomials as scalar types is interesting and
 > I tried to figure out your problems.  There were severals issues and
 > not everything was MTL's fault ;-)
 > 1. print_all_matrix is implemented as std::cout << x ..., it better
 > should be 'using std::cout; cout << x; ...'
 > - there are several such functions and we will fix this one day
 > - one can work around by defining the << operator in std namespace
 > 2. print_vector used iterator -> I changed it to const_iterator in
 > the hope that all relevant vectors have this
 > 3. instead of defining * and + operator for map it works better to
 > define a new type that derives from map (which is needed by
 > 4. anyway)
 > 4. the utilization of scale(v, 0) requires that polynomial is
 > constructible from int or double
 >
 > I add some code that solved the syntactical problems.  The
 > implementation of the numerical part I leave to you although that is
 > the more interesting.  I would appreciate if you can keep me
 > up-to-date occasionally (better directly to me instead of the list).
 >
 > Best,
 > Peter

********************************************************************

#include <iostream>
#include <map>
#include <vector>
#include <mtl/matrix.h>
#include <mtl/mtl.h>
#include <mtl/utils.h>
#include <mtl/blais.h>
#include <mtl/sparse1D.h>

using namespace mtl;
using namespace std;

struct polynomial : map<int, double>
{
   typedef map<int, double> base;

   polynomial() : base() {}

   polynomial(double c0) : base()
   {
     (*this)[0]= c0;
   }
};


typedef matrix<polynomial,
                 rectangle<>,
                 array< dense<> >,
                 row_major>::type Mat;


polynomial
operator *(const polynomial &p1, const polynomial &p2)
{
   polynomial result;

   for(polynomial::const_iterator i = p1.begin(); i != p1.end(); ++i)
     for(polynomial::const_iterator j = p2.begin(); j != p2.end(); ++j)
       result[i->first + j->first] += i->second * j->second;

   return result;
}


polynomial&
operator +=(polynomial &p1, const polynomial &p2)
{
   for(polynomial::const_iterator i = p2.begin(); i != p2.end(); ++i)
     p1[i->first] += i->second;
   return p1;
}


polynomial
operator +(const polynomial &p1, const polynomial &p2)
{
   polynomial result = p1;
   result += p2;
   return result;
}


namespace std
{
   ostream &operator << (ostream &out, const polynomial &p)
   {
     for(polynomial::const_reverse_iterator i = p.rbegin(); i != 
p.rend();)
       {
	out << i->second;

	if (i->first)
	  out << "*x^" << i->first;;
	
	if (++i != p.rend())
	  cout << " + ";
       }

     return out;
   }
}

int
main()
{
   Mat A(3, 3);

   A(0, 0) = polynomial();
   A(1, 1) = polynomial();
   A(2, 2) = polynomial();

   A(0, 0)[2] = 1;
   A(1, 1)[1] = 1;
   A(2, 2)[0] = 1;

   dense1D<polynomial> c(3), result(3);

   c[0][0] = 1;
   c[1][1] = 3;
   c[2][2] = 2;

   mult(A, c, result);

   print_all_matrix(A);

   print_vector(result);

   return 0;
}
_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/