[Boost-commit] svn:boost r51064 - in sandbox/monomials_horner_rule: . boost boost/math libs libs/math libs/math/doc libs/math/src libs/math/src/example

Subject: [Boost-commit] svn:boost r51064 - in sandbox/monomials_horner_rule: . boost boost/math libs libs/math libs/math/doc libs/math/src libs/math/src/example
From: erwann.rogard_at_[hidden]
Date: 2009-02-06 15:26:28

Author: e_r
Date: 2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
New Revision: 51064
URL: http://svn.boost.org/trac/boost/changeset/51064

Log:
Adding library monomials_horner_rule
Added:
   sandbox/monomials_horner_rule/
   sandbox/monomials_horner_rule/boost/
   sandbox/monomials_horner_rule/boost/math/
   sandbox/monomials_horner_rule/boost/math/monomials.hpp   (contents, props changed)
   sandbox/monomials_horner_rule/boost/math/monomials_properties.hpp   (contents, props changed)
   sandbox/monomials_horner_rule/boost/math/multi_factorial.hpp   (contents, props changed)
   sandbox/monomials_horner_rule/boost/math/multi_indexes.hpp   (contents, props changed)
   sandbox/monomials_horner_rule/boost/math/multi_indexes_derived.hpp   (contents, props changed)
   sandbox/monomials_horner_rule/libs/
   sandbox/monomials_horner_rule/libs/math/
   sandbox/monomials_horner_rule/libs/math/doc/
   sandbox/monomials_horner_rule/libs/math/doc/readme.txt   (contents, props changed)
   sandbox/monomials_horner_rule/libs/math/src/
   sandbox/monomials_horner_rule/libs/math/src/example/
   sandbox/monomials_horner_rule/libs/math/src/example/main.cpp   (contents, props changed)
   sandbox/monomials_horner_rule/libs/math/src/example/monomials.cpp   (contents, props changed)
   sandbox/monomials_horner_rule/libs/math/src/example/monomials.h   (contents, props changed)
   sandbox/monomials_horner_rule/libs/math/src/example/multi_factorial.cpp   (contents, props changed)
   sandbox/monomials_horner_rule/libs/math/src/example/multi_factorial.h   (contents, props changed)
   sandbox/monomials_horner_rule/libs/math/src/example/multi_indexes.cpp   (contents, props changed)
   sandbox/monomials_horner_rule/libs/math/src/example/multi_indexes.h   (contents, props changed)

Added: sandbox/monomials_horner_rule/boost/math/monomials.hpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/boost/math/monomials.hpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,175 @@
+//////////////////////////////////////////////////////////////////////////////
+//  monomials.hpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MATH_IFGT_MONOMIALS_HPP_ER_2009
+#define BOOST_MATH_IFGT_MONOMIALS_HPP_ER_2009
+#include <vector>
+#include <algorithm>
+#include <iterator>
+#include <boost/assert.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/range.hpp>
+#include <boost/tuple/tuple.hpp>
+//#include <boost/math/special_functions/binomial.hpp>
+#include <boost/bind.hpp>
+#include <boost/math/monomials_properties.hpp>
+namespace boost{namespace math{namespace ifgt{
+// Terminology associated with monomials
+// \see http://mathworld.wolfram.com/Polynomial.html
+// A polynomial is a mathematical expression involving a sum of powers in one or
+// more variables  multiplied by coefficients. A polynomial in one variable
+// (i.e., a univariate polynomial) with constant coefficients  is given by
+// a_nx^n+...+a_2x^2+a_1x+a_0. (1)
+// The individual summands with the coefficients (usually) included are called
+// monomials (Becker and Weispfenning 1993, p. 191),
+// the products of the form x_1^(alpha_1)...x_n^(alpha_n)
+// in the multivariate case,
+// i.e., with the coefficients omitted, are called terms
+// (Becker and Weispfenning 1993, p. 188).
+// The highest power  in a univariate polynomial is called its order,
+// or sometimes its degree
+
+
+        // Example of monomials and their efficient computation
+        // variables = (a,b,c)
+        //
+        // 1
+        // >a >b >c
+        // >a^2 ab ac >b^2 bc >c^2
+        // >a^3 a^2b a^2 c ab^2 abc ac^2 >b^3 b^2c bc^2 > c^3
+        //
+        // each > above marks the begining of a "read range"
+        // the kth write range on the nth line is obtained by
+        // multiplying the kth read range on the (n-1)th line
+        // by the kth variable in (a,b,c)
+        // For example, (b^3, b^2c, bc^2) = b*(b^2, bc, c^2)
+
+    template<typename Cont = std::vector<double> >
+    class monomials{
+        typedef typename Cont::size_type                size_type;
+        typedef typename Cont::iterator                 Iter;
+        typedef std::vector<Iter>                       Iters;
+        typedef monomials                               self_type;
+        typedef monomials_properties<>                  properties;
+    public:
+        typedef typename Cont::value_type               value_type;
+        typedef const Cont&                             result_type;
+        monomials(){}
+        monomials(self_type& that)
+        :degree_(that.degree_),terms_(that.terms_){}
+        self_type& operator=(const self_type that){
+            if(&that!=this){
+                degree_ = that.degree_;
+                terms_ = that.terms_;
+            }
+            return *this;
+        }
+
+        //degree (also called total degree) is |alpha|=alpha_0+...+alpha_{n-1}
+        unsigned degree()const{return degree_;}
+        result_type operator()()const{return terms_;}
+
+        template<typename R>
+        void operator()(const R& variables,unsigned the_degree){
+            typedef std::vector<Iter> Iters;
+            static Iters read_begins;
+            {
+                size_type sz = properties::number_degree_less_than(
+                    the_degree,(unsigned)(size(variables)));
+                terms_.resize(sz);
+            }
+            Iter write_begin = begin(terms_);
+
+            *write_begin++ = 1;
+            read_begins.clear();
+            std::fill_n(
+                back_inserter(read_begins),
+                size(variables),
+                begin(terms_)
+            );
+
+            unsigned i = 1;
+            unsigned n = the_degree + 1;
+            while(i<n){
+                BOOST_ASSERT(size(read_begins)==size(variables));
+                write_begin = for_each(//e=
+                    make_zip_iterator(
+                        make_tuple(
+                            begin(read_begins),
+                            begin(variables)
+                        )
+                    ),
+                    make_zip_iterator(
+                        make_tuple(
+                            end(read_begins),
+                            end(variables)
+                        )
+                    ),
+                    zip_func<R>(write_begin)
+                ).write_begin;
+
+                ++i;
+            }
+            degree_ = the_degree;
+        }
+
+    private:
+        // keeping terms_ inside spares the client the need to clear it
+        // or keeping track of the number of terms.
+        // using the same Cont for improved memory management
+        unsigned degree_;
+        Cont terms_;
+
+        template<typename R>
+        struct zip_func
+        {
+            typedef typename range_iterator<R>::type              Iter1;
+            typedef typename Iter1::value_type                    value1_type;
+            typedef const tuple<Iter&,const value1_type&>& argument_type;
+            typedef void                                          result_type;
+
+            zip_func(Iter write_begin_)
+            :write_begin(write_begin_),read_end(write_begin){}
+
+            zip_func(const zip_func& that)
+            :write_begin(that.write_begin),read_end(that.read_end){}
+
+            zip_func& operator=(const zip_func& that){
+                if(&that!=this){
+                    write_begin = that.write_begin;
+                    read_end = that.read_end;
+                }
+                return *this;
+            }
+
+            result_type operator()(argument_type t){
+                typename Iter1::value_type x = (t.template get<1>());
+
+                Iter read_begin = t.get<0>();
+                t.template get<0>() = write_begin;
+
+                write_begin =
+                    transform(
+                        read_begin,
+                        read_end,
+                        write_begin,
+                        bind(
+                            std::multiplies<value_type>(),
+                            _1,
+                            x
+                        )
+                    );
+            }
+
+            Iter write_begin;
+            Iter read_end;
+
+        };
+    };
+}}}
+
+
+#endif

Added: sandbox/monomials_horner_rule/boost/math/monomials_properties.hpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/boost/math/monomials_properties.hpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,24 @@
+//////////////////////////////////////////////////////////////////////////////
+// monomials_properties.hpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MATH_MONOMIALS_PROPERTIES_HPP_ER_2009
+#define BOOST_MATH_MONOMIALS_PROPERTIES_HPP_ER_2009
+#include <boost/math/special_functions/binomial.hpp>
+namespace boost{namespace math{
+
+    template<typename T = double>
+    struct monomials_properties{
+        static unsigned long number_degree_less_than(
+            unsigned the_degree,unsigned dim){
+                return (unsigned long)(binomial_coefficient<T>(
+                    the_degree+dim,dim));
+        }
+    };
+
+}}
+
+
+#endif

Added: sandbox/monomials_horner_rule/boost/math/multi_factorial.hpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/boost/math/multi_factorial.hpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,91 @@
+//////////////////////////////////////////////////////////////////////////////
+// multi_factorial.hpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MATH_MULTI_FACTORIAL_HPP_ER_2009
+#define BOOST_MATH_MULTI_FACTORIAL_HPP_ER_2009
+#include <numeric>
+#include <boost/math/multi_indexes.hpp>
+#include <boost/math/special_functions/factorials.hpp>
+
+namespace boost{namespace math{
+
+
+    // This is deprecated
+    // Use multi_indexed_derived instead
+    template<unsigned int SourceSz>
+    class multi_factorial{
+      public:
+        typedef std::vector<unsigned>                   storage_type;
+        typedef typename storage_type::iterator         iter_type;
+        typedef iterator_range<iter_type>               iter_range_type;
+        typedef typename storage_type::value_type       value_type;
+        typedef typename storage_type::size_type        size_type;
+
+       static iter_range_type get(unsigned int degree){
+            static unsigned max_degree = 0;
+            static storage_type storage(1,(value_type)(0)); //initial value
+            size_type dist
+                = monomials_properties<>::number_degree_less_than(
+                    degree,SourceSz);
+            if(max_degree<degree){
+                typedef multi_indexes<SourceSz>                      mi_type;
+                typedef
+                    typename mi_type::iter_range_type      mi_iter_range_type;
+                mi_iter_range_type ir = mi_type::get(degree);
+                storage.clear();
+                transform(
+                    begin(ir),
+                    end(ir),
+                    back_inserter(storage),
+                    unary_op()
+                );
+                max_degree = degree;
+            }
+
+            iter_type b(begin(storage));
+            iter_type e(b); std::advance(e,dist);
+
+            return iter_range_type(b,e);
+        }
+      private:
+        multi_factorial(); //do not implement
+
+        typedef multi_indexes<SourceSz>  mi_type;
+        typedef typename mi_type::iter_range_type       mi_iter_range_type;
+        typedef typename mi_type::iter_type             mi_iter_type;
+        typedef typename mi_type::value_type            mi_value_type;
+        typedef typename mi_type::base_iter_type        mi_base_iter_type;
+        typedef typename mi_type::base_value_type       mi_base_value_type;
+
+        struct unary_op{
+            typedef mi_value_type                       argument_type;
+            typedef value_type                          result_type;
+            unary_op(){}
+            result_type operator()(argument_type ir){
+                return accumulate(
+                    begin(ir),
+                    end(ir),
+                    (value_type)(1), //initial value
+                    inner_op()
+                );
+            }
+            iter_type write_begin;
+            struct inner_op{
+                typedef value_type                      result_type;
+                typedef value_type                      first_argument_type;
+                typedef value_type                      second_argument_type;
+                inner_op(){}
+                result_type operator()(
+                    first_argument_type a1,
+                    second_argument_type a2)const{
+                        return a1 * (result_type)(factorial<double>(a2));
+                }
+            };
+        };
+    };
+}}
+
+#endif

Added: sandbox/monomials_horner_rule/boost/math/multi_indexes.hpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/boost/math/multi_indexes.hpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,276 @@
+//////////////////////////////////////////////////////////////////////////////
+// multi_indexes.hpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MULTI_INDEXES_HPP_ER_2009
+#define BOOST_MULTI_INDEXES_HPP_ER_2009
+#include <vector>
+#include <functional>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/range.hpp>
+#include <boost/tuple/tuple.hpp>
+#include <boost/bind.hpp>
+#include <boost/assert.hpp>
+#include <boost/iterator/vector2matrix_iterator.hpp>
+#include <boost/math/monomials_properties.hpp>
+
+
+namespace boost{namespace math{
+
+    // Example of multi indexes for SourceSz = 3, degree =2
+    // >(0,0,0)
+    // >(1,0,0) >(0,1,0) >(0,0,1)
+    // >(2,0,0) (1,1,0) (1,0,1) >(0,2,0) (0,1,1) >(0,0,2)
+    // Note that
+    // {(2,0,0) (1,1,0) (1,0,1)} = {(1,0,0) (0,1,0) (0,0,1)} + (1,0,0)
+    // {(0,2,0) (0,1,1)} = {(0,1,0) (0,0,1)} + (0,1,0) etc.
+
+    template<unsigned int SourceSz>
+    class multi_indexes{
+      public:
+        typedef std::vector<unsigned>                   storage_type;
+        typedef typename storage_type::iterator         base_iter_type;
+        typedef typename storage_type::size_type        size_type;
+        typedef storage_type::value_type                base_value_type;
+
+        typedef vector2matrix_iterator<base_iter_type>  iter_type;
+        typedef std::vector<iter_type>                  iters_type;
+        typedef typename iter_type::value_type          value_type;
+        typedef iterator_range<iter_type>               iter_range_type;
+        typedef typename iters_type::iterator           iters_iter_type;
+        typedef typename iter_type::difference_type     iter_diff_type;
+        typedef std::vector<iter_diff_type>             iters_diff_type;
+        typedef typename iters_diff_type::iterator      iters_diff_iter_type;
+
+        static iter_range_type get(unsigned int degree){
+            static unsigned max_degree = 0;
+            static storage_type storage(SourceSz,(base_value_type)(0));
+            static iters_diff_type read_begins_dists(SourceSz,0);
+
+// // keeping read_begins rather than read_begins_dists causes
+// // error 'gliblc detected' at runtime
+//            static iters_type   read_begins(
+//                SourceSz,
+//                make_vector2matrix_iterator(
+//                    begin(storage),
+//                    SourceSz
+//            ));
+
+            //e.g. SourceSz=3: (1,0,0, 0,1,0, 0,0,1)
+            static storage_type variables = variables_init();
+
+            if(degree>max_degree){
+                unsigned i = max_degree+1;
+                unsigned n = degree+1;
+                size_type sz_old = size(storage)/SourceSz;
+                size_type sz_new
+                    = properties::number_degree_less_than(degree,SourceSz);
+                storage.resize(sz_new*SourceSz);
+                iter_type write_begin = make_vector2matrix_iterator(
+                    begin(storage),
+                    SourceSz
+                );
+                std::advance(write_begin,sz_old);
+
+                iters_type read_begins;
+
+//old implem, do not uncomment
+//                iters_diff_iter_type b = begin(read_begins_dists);
+//                iters_diff_iter_type e = end(read_begins_dists);
+//                for(iters_diff_iter_type i = b; i<e; i++){
+//                    iter_type it(begin(storage),SourceSz);
+//                    std::advance(it,*i);
+//                    read_begins.push_back(it);
+//                }
+
+                transform(
+                    begin(read_begins_dists),
+                    end(read_begins_dists),
+                    back_inserter(read_begins),
+                    unary_op(
+                        make_vector2matrix_iterator(
+                            begin(storage),
+                            SourceSz
+                        )
+                    )
+                );
+
+                while(i<n){
+                    BOOST_ASSERT(size(read_begins_dists)==SourceSz);
+                    BOOST_ASSERT(std::distance(
+                        make_vector2matrix_iterator(
+                            begin(variables),
+                            SourceSz
+                        ),
+                        make_end_vector2matrix_iterator(
+                                begin(variables),
+                                end(variables),
+                                SourceSz
+                        )
+                    )==SourceSz);
+
+                    //(***) bug
+                    write_begin = for_each(
+                        make_zip_iterator(
+                            make_tuple(
+                                begin(read_begins),
+                                make_vector2matrix_iterator(
+                                    begin(variables),
+                                    SourceSz
+                                )
+                            )
+                        ),
+                        make_zip_iterator(
+                            make_tuple(
+                                end(read_begins),
+                                make_end_vector2matrix_iterator(
+                                    begin(variables),
+                                    end(variables),
+                                    SourceSz
+                                )
+                            )
+                        ),
+                        zip_func(write_begin)
+                    ).write_begin;
+                    ++i;
+                }
+                read_begins_dists.clear();
+                std::transform(
+                    begin(read_begins),
+                    end(read_begins),
+                    back_inserter(read_begins_dists),
+                    bind(
+                        std::distance<iter_type>,
+                        make_vector2matrix_iterator(
+                            begin(storage),
+                            SourceSz
+                        ),
+                        _1
+                    )
+                );
+                max_degree = degree;
+            }//if(degree>max_degree)
+
+            size_type dist
+                = properties::number_degree_less_than(degree,SourceSz);
+            iter_type b(begin(storage),SourceSz);
+            iter_type e = b;
+            std::advance(e,dist);
+            iter_range_type(b,e);
+            return iter_range_type(b,e);
+        }//get(unsigned degree)
+
+      private:
+        typedef monomials_properties<> properties;
+        multi_indexes(); //do not implement
+
+        static storage_type variables_init(){
+            storage_type variables(SourceSz*SourceSz,0);
+            for(unsigned i=0; i<SourceSz; i++){
+                variables[i*SourceSz+i]=1;
+            }
+            return variables;
+        }
+
+        struct unary_op{
+            iter_type b;
+            unary_op(iter_type b_):b(b_){}
+            unary_op(const unary_op& that):b(that.b){}
+            unary_op& operator=(const unary_op& that){
+                if(&that!=this){b=that.b;}
+                return *this;
+            }
+            typedef iter_type                   result_type;
+            typedef const iter_diff_type&       argument_type;
+            result_type operator()(argument_type diff)const{
+                result_type it = b;
+                std::advance(it,diff);
+                return it;
+            }
+        };
+
+        struct zip_func
+        {
+            typedef const tuple<
+                iter_type&,const value_type&>&                  argument_type;
+            typedef void                                        result_type;
+
+            iter_type write_begin;
+            iter_type read_end;
+            unsigned cnt;//TODO remove
+
+            zip_func(iter_type write_begin_)
+            :write_begin(write_begin_),read_end(write_begin),cnt(0){}
+
+            zip_func(const zip_func& that)
+            :write_begin(that.write_begin),
+            read_end(that.read_end),
+            cnt(that.cnt){}
+
+            zip_func& operator=(const zip_func& that){
+                if(&that!=this){
+                    write_begin = that.write_begin;
+                    read_end = that.read_end;
+                    cnt = that.cnt;
+                }
+                return *this;
+            }
+
+            result_type operator()(argument_type t){
+                value_type variable_ir = (t.template get<1>());
+                iter_type read_begin = (t.template get<0>());
+                t.template get<0>() = write_begin;
+                write_begin =
+                    for_each(
+                        read_begin,
+                        read_end,
+                        unary_op(write_begin,variable_ir)
+                    ).write_begin;
+            }
+
+            struct unary_op{
+                unary_op(iter_type write_begin_,value_type variable_ir)
+                :write_begin(write_begin_),variable(variable_ir){}
+
+                unary_op(const unary_op& that)
+                :write_begin(that.write_begin),variable(that.variable){}
+
+                unary_op& operator=(const unary_op& that)
+                {
+                    if(&that!=this){
+                        write_begin = that.write_begin;
+                        variable = that.variable;
+                    }
+                    return *this;
+                }
+
+                typedef typename iter_type::value_type     argument_type;
+                typedef void                               result_type;
+                result_type operator()(argument_type ir){
+
+                    BOOST_ASSERT(size(*write_begin)==size(ir));
+                    BOOST_ASSERT(size(ir)==size(variable));
+                    std::transform(
+                        begin(ir),
+                        end(ir),
+                        begin(variable),
+                        begin(*write_begin),
+                        bind(
+                            std::plus<base_value_type>(),
+                            _1,
+                            _2
+                        )
+                    );
+                    ++write_begin;
+                }
+                iter_type   write_begin;
+                value_type  variable;
+            };//unary_op
+        };//zip_func
+    };//multi_indexes
+}}
+
+
+#endif

Added: sandbox/monomials_horner_rule/boost/math/multi_indexes_derived.hpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/boost/math/multi_indexes_derived.hpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,125 @@
+//////////////////////////////////////////////////////////////////////////////
+// multi_indexes_derived.hpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MULTI_INDEXES_DERIVED_HPP_ER_2009
+#define BOOST_MULTI_INDEXES_DERIVED_HPP_ER_2009
+#include <numeric>
+#include <cmath>
+#include <boost/math/multi_indexes.hpp>
+#include <boost/math/special_functions/factorials.hpp>
+
+namespace boost{namespace math{
+
+// Multi indexes derived are functions of multi_indexes
+// The particular function is controled by a policy:
+// multi_factorial_op :    (a,b,c) --> a! b! c!
+// multi_degree_op:        (a,b,c) --> a + b + c (also called total degree)
+// multi_power2divfact_op  (a,b,c) --> 2^{a+b+c}/(a! b! c!)
+
+            struct multi_factorial_op{
+                typedef unsigned                      result_type;
+                typedef unsigned                      first_argument_type;
+                typedef unsigned                      second_argument_type;
+                static result_type initial_value(){
+                    static result_type iv = (result_type)(1);
+                    return iv;
+                }
+                multi_factorial_op(){}
+                result_type operator()(
+                    first_argument_type a1,
+                    second_argument_type a2)const{
+                        return a1 * (result_type)(factorial<double>(a2));
+                }
+            };
+            struct multi_degree_op:public std::plus<unsigned>{
+                static unsigned initial_value(){
+                    static unsigned iv = (unsigned)(0);
+                    return iv;
+                }
+            };
+            struct multi_power2divfact_op{
+                typedef double                      result_type;
+                typedef double                      first_argument_type;
+                typedef unsigned                    second_argument_type;
+                static result_type initial_value(){
+                    static result_type iv = (result_type)(1);
+                    return iv;
+                }
+                multi_power2divfact_op(){}
+                result_type operator()(
+                    first_argument_type a1,
+                    second_argument_type a2)const{
+                        result_type d = std::pow(2.0,(result_type)(a2));
+                        result_type f = factorial<result_type>(a2);
+                        result_type r = d/f;
+                        return a1 * ((result_type)(r));
+                }
+            };
+
+
+    template<unsigned int SourceSz,typename Policy>
+    class multi_indexes_derived{
+      public:
+        typedef typename Policy::result_type            value_type;
+        typedef std::vector<value_type>                 storage_type;
+        typedef typename storage_type::iterator         iter_type;
+        typedef iterator_range<iter_type>               iter_range_type;
+        typedef typename storage_type::size_type        size_type;
+
+       static iter_range_type get(unsigned int degree){
+            static unsigned max_degree = 0;
+            static storage_type storage(1,Policy::initial_value());
+            size_type dist
+                = monomials_properties<>::number_degree_less_than(
+                    degree,SourceSz);
+            if(max_degree<degree){
+                typedef multi_indexes<SourceSz>            mi_type;
+                typedef
+                    typename mi_type::iter_range_type      mi_iter_range_type;
+                mi_iter_range_type ir = mi_type::get(degree);
+                storage.clear();
+                transform(
+                    begin(ir),
+                    end(ir),
+                    back_inserter(storage),
+                    unary_op()
+                );
+                max_degree = degree;
+            }
+
+            iter_type b(begin(storage));
+            iter_type e(b); std::advance(e,dist);
+
+            return iter_range_type(b,e);
+        }
+      private:
+        multi_indexes_derived(); //do not implement
+
+        typedef multi_indexes<SourceSz>  mi_type;
+        typedef typename mi_type::iter_range_type       mi_iter_range_type;
+        typedef typename mi_type::iter_type             mi_iter_type;
+        typedef typename mi_type::value_type            mi_value_type;
+        typedef typename mi_type::base_iter_type        mi_base_iter_type;
+        typedef typename mi_type::base_value_type       mi_base_value_type;
+
+        struct unary_op{
+            typedef mi_value_type                       argument_type;
+            typedef value_type                          result_type;
+            unary_op(){}
+            result_type operator()(argument_type ir){
+                return accumulate(
+                    begin(ir),
+                    end(ir),
+                    Policy::initial_value(), //initial value
+                    Policy()
+                );
+            }
+            iter_type write_begin;
+        };
+    };
+}}
+
+#endif

Added: sandbox/monomials_horner_rule/libs/math/doc/readme.txt
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/libs/math/doc/readme.txt	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,119 @@
+//////////////////////////////////////////////////////////////////////////////
+// monomials_horner_rule
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+///////////
+/ Contact /
+///////////
+
+Please send questions or suggestions to erwann.rogard_at_[hidden]
+
+/////////////
+/ Overview  /
+/////////////
+
+This collection of C++ classes computes multivariate monomials and derived quantities using Horner's rule a.k.a. graded lexicographic ordering.
+Individual classes contain specific documentation.
+
+//////////////////
+/ Requirements   /
+//////////////////
+
+Compiles fine under
+
+gcc version i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1 (Apple Inc. build 5490)
+
+The compiler search path must include
+boost_1_37_0
+boost/sandbox/miscellanea_iterator_facade/
+
+///////////////////////
+/ Output from main.cpp /
+///////////////////////
+-> example_monomials
+
+ -> d=1
+  ->p=1
+1 0   <-
+  ->p=2
+1 0 0   <-
+  ->p=3
+1 0 0 0   <-
+  ->p=4
+1 0 0 0 0   <-
+ <-
+ -> d=2
+  ->p=1
+1 0 1   <-
+  ->p=2
+1 0 1 0 0 1   <-
+  ->p=3
[phone]   <-
+  ->p=4
[phone]   <-
+ <-
+ -> d=3
+  ->p=1
+1 0 1 2   <-
+  ->p=2
[phone]   <-
+  ->p=3
[phone] 0 1 2 4 8   <-
+  ->p=4
[phone] 0 1 2 4 8 0 0 0 0 0 0 0 0 0 0 1 2 4 8 16   <-
+ <-
+ -> d=4
+  ->p=1
+1 0 1 2 3   <-
+  ->p=2
[phone]   <-
+  ->p=3
[phone] 0 0 0 0 0 0 0 0 0 0 1 2 3 4 6 9 8 12 18 27   <-
+  ->p=4
[phone] 0 0 0 0 0 0 0 0 0 0 1 2 3 4 6 9 8 12 18 27 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 3 4 6 9 8 12 18 27 16 24 36 54 81   <-
+ <-
+<-
+-> example_multi_indexes
+0 0 0
+1 0 0
+0 1 0
+0 0 1
+2 0 0
+1 1 0
+1 0 1
+0 2 0
+0 1 1
+0 0 2
+
+0 0 0
+1 0 0
+0 1 0
+0 0 1
+2 0 0
+1 1 0
+1 0 1
+0 2 0
+0 1 1
+0 0 2
+3 0 0
+2 1 0
+2 0 1
+1 2 0
+1 1 1
+1 0 2
+0 3 0
+0 2 1
+0 1 2
+0 0 3
+<-
+-> example_multi_factorial
+ -> multi_factorial<3>::degree(5)
[phone] 2 6 2 2 6 24 6 6 4 2 4 6 2 2 6 24 6 4 6 24 120 24 24 12 6 12 12 4 4 12 24 6 4 6 24 120 24 12 12 24 120
+ -> multi_degree<3>::degree(5)
+0 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+ -> multi_p2df<3>::degree(5)
[phone] 1.33333 4 4 4 8 4 1.33333 4 4 1.33333 0.666667 2.66667 2.66667 4 8 4 2.66667 8 8 2.66667 0.666667 2.66667 4 2.66667 0.666667 0.266667 1.33333 1.33333 2.66667 5.33333 2.66667 2.66667 8 8 2.66667 1.33333 5.33333 8 5.33333 1.33333 0.266667 1.33333 2.66667 2.66667 1.33333 0.266667
+<-
+

Added: sandbox/monomials_horner_rule/libs/math/src/example/main.cpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/libs/math/src/example/main.cpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,20 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/main.cpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#include <libs/math/src/example/monomials.h>
+#include <libs/math/src/example/multi_indexes.h>
+#include <libs/math/src/example/multi_factorial.h>
+
+int main()
+{
+
+
+    example_monomials();
+    example_multi_indexes();
+    example_multi_factorial();
+
+    return 0;
+}

Added: sandbox/monomials_horner_rule/libs/math/src/example/monomials.cpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/libs/math/src/example/monomials.cpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,55 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/monomials.cpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#include <algorithm>
+#include <boost/math/special_functions/binomial.hpp>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/math/monomials.hpp>
+#include <libs/math/src/example/monomials.h>
+
+void example_monomials(){
+    std::cout << "-> example_monomials" << std::endl;
+
+    typedef std::vector<double> vars_type;
+    namespace math = boost::math;
+    typedef math::ifgt::monomials<vars_type> monoms_type;
+
+    unsigned max_d = 5;
+    unsigned max_p = 5;
+
+    vars_type vars;
+    {using namespace boost::assign; vars+=0.0,1.0,2.0,3.0,4.0;}
+
+
+    for(unsigned d = 1; d<max_d; d++){
+        std::cout << " -> d=" << d << std::endl;
+        vars_type::iterator e = boost::begin(vars);
+        std::advance(e,d);
+        vars_type sub;
+        copy(boost::begin(vars),e,back_inserter(sub));
+        BOOST_ASSERT(boost::size(sub)-d==0);
+
+        for(unsigned int p = 1; p<max_p; p++){
+            std::cout << "  ->p=" << p << std::endl;
+            //std::cout << boost::math::binomial_coefficient<double>(p+d,d)
+            //<< std::endl;
+
+            monoms_type monoms;
+
+            monoms(sub,p);
+            copy(
+                boost::begin(monoms()),
+                boost::end(monoms()),
+                std::ostream_iterator<double>(std::cout," ")
+            );
+            std::cout << "  <-" << std::endl;
+        }
+            std::cout << " <-" << std::endl;
+    }
+    std::cout << "<-" << std::endl;
+}

Added: sandbox/monomials_horner_rule/libs/math/src/example/monomials.h
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/libs/math/src/example/monomials.h	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,11 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/monomial.h
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_MONOMIALS_H_ER_200901
+#define LIBS_MATH_IFGT_EXAMPLE_MONOMIALS_H_ER_200901
+void example_monomials();
+
+#endif

Added: sandbox/monomials_horner_rule/libs/math/src/example/multi_factorial.cpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/libs/math/src/example/multi_factorial.cpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,61 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/multi_factorial.cpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <boost/range.hpp>
+#include <boost/math/multi_factorial.hpp>
+#include <boost/math/multi_indexes_derived.hpp>
+#include <libs/math/src/example/multi_factorial.h>
+void example_multi_factorial(){
+    std::cout << "-> example_multi_factorial" << std::endl;
+
+    namespace math = boost::math;
+
+    const unsigned d = 3;
+    //typedef math::multi_factorial<d>          mf_type;
+    typedef math::multi_indexes_derived<d,math::multi_factorial_op>   mf_type;
+    typedef math::multi_indexes_derived<d,math::multi_degree_op>      md_type;
+    typedef math::multi_indexes_derived<d,math::multi_power2divfact_op> mp_type;
+    typedef mf_type::iter_range_type            mf_iter_range_type;
+    typedef md_type::iter_range_type            md_iter_range_type;
+    typedef mp_type::iter_range_type            mp_iter_range_type;
+
+    unsigned int degree1 = 5;
+
+    mf_iter_range_type mf_ir = mf_type::get(degree1);
+
+    std::cout << " -> multi_factorial<" << d << ">";
+    std::cout << "::degree(" << degree1 << ")" << std::endl;
+    copy(
+        boost::begin(mf_ir),
+        boost::end(mf_ir),
+        std::ostream_iterator<double>(std::cout," ")
+    ); std::cout << std::endl;
+
+    md_iter_range_type md_ir = md_type::get(degree1);
+
+    std::cout << " -> multi_degree<" << d << ">";
+    std::cout << "::degree(" << degree1 << ")" << std::endl;
+    copy(
+        boost::begin(md_ir),
+        boost::end(md_ir),
+        std::ostream_iterator<double>(std::cout," ")
+    ); std::cout << std::endl;
+    mp_iter_range_type mp_ir = mp_type::get(degree1);
+
+    std::cout << " -> multi_p2df<" << d << ">";
+    std::cout << "::degree(" << degree1 << ")" << std::endl;
+    copy(
+        boost::begin(mp_ir),
+        boost::end(mp_ir),
+        std::ostream_iterator<double>(std::cout," ")
+    ); std::cout << std::endl;
+
+    std::cout << "<-" << std::endl;
+
+};

Added: sandbox/monomials_horner_rule/libs/math/src/example/multi_factorial.h
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/libs/math/src/example/multi_factorial.h	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,13 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/multi_factorial.h
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_MULTI_FACTORIAL_H_ER_200901
+#define LIBS_MATH_IFGT_EXAMPLE_MULTI_FACTORIAL_H_ER_200901
+
+void example_multi_factorial();
+
+
+#endif

Added: sandbox/monomials_horner_rule/libs/math/src/example/multi_indexes.cpp
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/libs/math/src/example/multi_indexes.cpp	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,44 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/multi_indexes.cpp
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <boost/range.hpp>
+#include <boost/math/multi_indexes.hpp>
+void example_multi_indexes(){
+    std::cout << "-> example_multi_indexes" << std::endl;
+
+    namespace math = boost::math;
+
+    typedef math::multi_indexes<3>      mi_type;
+    typedef mi_type::iter_type          mi_iter_type;
+    typedef mi_type::iter_range_type    mi_iter_range_type;
+
+    unsigned int degree1 = 2;
+
+    mi_iter_range_type ir = mi_type::get(degree1);
+
+    for(mi_iter_type i = boost::begin(ir); i<boost::end(ir); i++){
+        copy(
+            boost::begin(*i),
+            boost::end(*i),
+            std::ostream_iterator<double>(std::cout," ")
+        ); std::cout << std::endl;
+    }std::cout << std::endl;
+
+    ir = mi_type::get(degree1+1);//TODO bug if say get(degree+1)
+
+    for(mi_iter_type i = boost::begin(ir); i<boost::end(ir); i++){
+        copy(
+            boost::begin(*i),
+            boost::end(*i),
+            std::ostream_iterator<double>(std::cout," ")
+        ); std::cout << std::endl;
+    }
+
+    std::cout << "<-" << std::endl;
+};

Added: sandbox/monomials_horner_rule/libs/math/src/example/multi_indexes.h
==============================================================================
--- (empty file)
+++ sandbox/monomials_horner_rule/libs/math/src/example/multi_indexes.h	2009-02-06 15:26:27 EST (Fri, 06 Feb 2009)
@@ -0,0 +1,11 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/multi_indexes.h
+//  (C) Copyright 2009 Erwann Rogard
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_MULTI_INDEXES_H_ER_200901
+#define LIBS_MATH_IFGT_EXAMPLE_MULTI_INDEXES_H_ER_200901
+
+void example_multi_indexes();
+#endif