---

more/num/numeric.h

Go to the documentation of this file.

//  more/numeric.h -- Type independent numeric algorithms
//  Copyright (C) 1999--2001  Petter Urkedal (petter.urkedal@matfys.lth.se)

//  This file is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.

//  This file is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

//  As a special exception, you may use this file as part of a free
//  software library without restriction.  Specifically, if other files
//  instantiate templates or use macros or inline functions from this
//  file, or you compile this file and link it with other files to
//  produce an executable, this file does not by itself cause the
//  resulting executable to be covered by the GNU General Public
//  License.  This exception does not however invalidate any other
//  reasons why the executable file might be covered by the GNU General
//  Public License.

//  $Id: numeric.h,v 1.1 2002/05/30 18:01:40 petter_urkedal Exp $


#ifndef MORE_NUMERIC_H
#define MORE_NUMERIC_H

#include <numeric>
#include <list>
#include <vector>
#include <more/math/math.h>
#include <more/math/complex.h>

namespace more {
namespace num {

  //  Generic integral power for types with associative
  //  multiplication.
  template<typename T>
    inline T assoc_power(T x, unsigned int y) {
        T z = y&1? x : 1;
        while(y >>= 1) {
            x *= x;
            if(y & 1) z *= x;
        }
        return z;
    }
  template<typename T>
    inline T assoc_power(T x, int y) {
        if(y < 0) return T(1)/pow(x, (unsigned int)-y);
        else return pow(x, (unsigned int)y);
    }


  template<typename T> struct mutable_ { typedef T eval; };
  template<typename T> struct mutable_<const T> { typedef T eval; };


  template<typename Iterator, typename Iterator2>
    inline typename scalar_type_
           <typename std::iterator_traits<Iterator>::value_type>::eval
    dot(Iterator first, Iterator last, Iterator2 first2,
        typename scalar_type_
        <typename std::iterator_traits<Iterator>::value_type>::eval
        init = 0.0) {
        while(first != last) {
            init += adj(*first) * (*first2);
            ++first; ++first2;
        }
        return init;
    }

  template<typename Iterator, typename Iterator2, typename Rank2Iterator>
    inline typename scalar_type_
           <typename std::iterator_traits<Iterator>::value_type>::eval
    dot(Iterator first, Iterator last, Iterator2 first2,
        Rank2Iterator it_metric,
        typename scalar_type_
        <typename std::iterator_traits<Iterator>::value_type>::eval
        init = 0.0) {
        while(first != last) {
            init += adj(*first) *
                vector_vector_product(it_metric->begin(), it_metric->end(),
                                      first2);
            ++first;
            ++it_metric;
        }
        return init;
    }

  template<typename Iterator>
    inline typename norm_type_
           <typename std::iterator_traits<Iterator>::value_type>::eval
    norm(Iterator first, Iterator last,
         typename norm_type_
         <typename std::iterator_traits<Iterator>::value_type>::eval
         init = zero) {
        while(first != last) {
            init += norm(*first);
            ++first;
        }
        return init;
    }

  template<typename Iterator, typename Rank2Iterator>
    inline typename norm_type_
           <typename std::iterator_traits<Iterator>::value_type>::eval
    norm(Iterator first, Iterator last,
         Rank2Iterator it_metric,
         typename norm_type_
         <typename std::iterator_traits<Iterator>::value_type>::eval
         init = 0.0) {
        Iterator first2 = first;
        while(first != last) {
            init += real(adj(*first)*vector_vector_product
                         (it_metric->begin(), it_metric->end(), first2));
            ++first;
            ++it_metric;
        }
        return init;
    }

  template<typename Iterator>
    inline void normalize(Iterator first, Iterator last) {
        typename norm_type_
                 <std::iterator_traits<Iterator>::value_type>::eval
            x = 1.0/sqrt(norm(first, last));
        while(first != last) { *first *= x; ++first; }
    }

  template<typename Iterator, typename Rank2Iterator>
    inline void normalize(Iterator first, Iterator last,
                          Rank2Iterator it_metric) {
        typename norm_type_
                 <std::iterator_traits<Iterator>::value_type>::eval
            x = 1.0/sqrt(norm(first, last, it_metric));
        while(first != last) { *first *= x; ++first; }
    }

  template<typename Iterator1, typename Iterator2, typename Scale>
    inline void assign_add_scaled(Iterator1 it1, Iterator1 it1_end,
                                  Iterator2 it2, Scale scale) {
        while(it1 != it1_end) {
            *it1 += scale*(*it2);
            ++it1; ++it2;
        }
    }

  template<typename Rank2Iterator>
    inline void orthonormalize(Rank2Iterator first, Rank2Iterator last) {
        Rank2Iterator it = first;
        while(it != last) { normalize(it->begin(), it->end()); ++it; }
        it = first;
        while(it != last) {
            for(Rank2Iterator it2 = first; it2 != it; ++it2)
                assign_add_scaled(
                    it->begin(), it->end(), it2->begin(),
                    -dot(it2->begin(), it2->end(), it->begin()));
            normalize(it->begin(), it->end());
            ++it;
        }
    }

  template<typename Rank2Iterator, typename Rank2Iterator2>
    inline void orthonormalize(Rank2Iterator first, Rank2Iterator last,
                               Rank2Iterator2 it_metric) {
        Rank2Iterator it = first;
        while(it != last) {
            normalize(it->begin(), it->end(), it_metric);
            ++it;
        }
        it = first;
        while(it != last) {
            for(Rank2Iterator it2 = first; it2 != it; ++it2)
                assign_add_scaled(
                    it->begin(), it->end(), it2->begin(),
                    -dot(it2->begin(), it2->end(), it->begin(),
                         it_metric));
            normalize(it->begin(), it->end(), it_metric);
            ++it;
        }
    }

  template<typename Iterator, typename Rank2Iterator>
    inline void project(Iterator first, Iterator last,
                        Rank2Iterator it_space_begin,
                        Rank2Iterator it_space_end) {
        typedef scalar_type_
                <std::iterator_traits<Iterator>::value_type>::eval
                scalar_type;
        std::list<scalar_type> coeff;
        for(Rank2Iterator it_space = it_space_begin;
            it_space != it_space_end; ++it_space)
            coeff.push_back(dot(it_space->begin(), it_space->end(), first));

        normalize(coeff.begin(), coeff.end());
        std::list<scalar_type>::iterator it_coeff = coeff.begin();
        std::fill(first, last, 0.0);
        for(Rank2Iterator it_space = it_space_begin;
            it_space != it_space_end; ++it_space) {
            assign_add_scaled(first, last, it_space->begin(), *it_coeff);
            ++it_coeff;
        }
    }

  template<typename Rank2Iterator, typename Rank2Iterator2>
    inline void project_collection(Rank2Iterator first, Rank2Iterator last,
                                   Rank2Iterator2 it_space_begin,
                                   Rank2Iterator2 it_space_end) {
        while(first != last) {
            project(first->begin(), first->end(),
                    it_space_begin, it_space_end);
            ++first;
        }
    }

  template<typename Iterator, typename Iterator2>
    inline
    typename times_type_
        < typename std::iterator_traits<Iterator>::value_type,
          typename std::iterator_traits<Iterator2>::value_type >::eval
    vector_vector_product(Iterator it, Iterator it_end, Iterator2 it2) {
        typename mutable_
          < typename times_type_
            < typename std::iterator_traits<Iterator>::value_type,
              typename std::iterator_traits<Iterator2>::value_type
            >::eval >::eval
            res = zero;
        while(it != it_end) {
            res += (*it) * (*it2);
            ++it;
            ++it2;
        }
        return res;
    }

  template<typename Rank2Iterator, typename Iterator>
    inline void
    multiply_matrix_to_vector(Rank2Iterator it, Rank2Iterator it_end,
                              Iterator it_dest, Iterator it_dest_end) {
        typedef typename std::iterator_traits<Iterator>::value_type
                value_type;
        std::vector<value_type> src(it_dest, it_dest_end);
        while(it != it_end) {
            *it_dest = vector_vector_product(it->begin(), it->end(),
                                             src.begin());
            ++it;
            ++it_dest;
        }
    }

  template<typename Rank2Iterator, typename Rank2Iterator2>
    inline void
    multiply_matrix_to_matrix(Rank2Iterator it, Rank2Iterator it_end,
                              Rank2Iterator2 it2, Rank2Iterator2 it2_end) {
        while(it2 != it2_end) {
            multiply_matrix_to_vector(it, it_end, it2->begin(), it2->end());
            ++it2;
        }
    }
}}

#endif

---