Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

more/num/numeric.h

Go to the documentation of this file.
00001 //  more/numeric.h -- Type independent numeric algorithms
00002 //  Copyright (C) 1999--2001  Petter Urkedal (petter.urkedal@matfys.lth.se)
00003 
00004 //  This file is free software; you can redistribute it and/or modify
00005 //  it under the terms of the GNU General Public License as published by
00006 //  the Free Software Foundation; either version 2 of the License, or
00007 //  (at your option) any later version.
00008 
00009 //  This file is distributed in the hope that it will be useful,
00010 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 //  GNU General Public License for more details.
00013 
00014 //  You should have received a copy of the GNU General Public License
00015 //  along with this program; if not, write to the Free Software
00016 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00017 
00018 //  As a special exception, you may use this file as part of a free
00019 //  software library without restriction.  Specifically, if other files
00020 //  instantiate templates or use macros or inline functions from this
00021 //  file, or you compile this file and link it with other files to
00022 //  produce an executable, this file does not by itself cause the
00023 //  resulting executable to be covered by the GNU General Public
00024 //  License.  This exception does not however invalidate any other
00025 //  reasons why the executable file might be covered by the GNU General
00026 //  Public License.
00027 
00028 //  $Id: numeric.h,v 1.1 2002/05/30 18:01:40 petter_urkedal Exp $
00029 
00030 
00031 #ifndef MORE_NUMERIC_H
00032 #define MORE_NUMERIC_H
00033 
00034 #include <numeric>
00035 #include <list>
00036 #include <vector>
00037 #include <more/math/math.h>
00038 #include <more/math/complex.h>
00039 
00040 namespace more {
00041 namespace num {
00042 
00043   //  Generic integral power for types with associative
00044   //  multiplication.
00045   template<typename T>
00046     inline T assoc_power(T x, unsigned int y) {
00047         T z = y&1? x : 1;
00048         while(y >>= 1) {
00049             x *= x;
00050             if(y & 1) z *= x;
00051         }
00052         return z;
00053     }
00054   template<typename T>
00055     inline T assoc_power(T x, int y) {
00056         if(y < 0) return T(1)/pow(x, (unsigned int)-y);
00057         else return pow(x, (unsigned int)y);
00058     }
00059 
00060 
00061   template<typename T> struct mutable_ { typedef T eval; };
00062   template<typename T> struct mutable_<const T> { typedef T eval; };
00063 
00064 
00065   template<typename Iterator, typename Iterator2>
00066     inline typename scalar_type_
00067            <typename std::iterator_traits<Iterator>::value_type>::eval
00068     dot(Iterator first, Iterator last, Iterator2 first2,
00069         typename scalar_type_
00070         <typename std::iterator_traits<Iterator>::value_type>::eval
00071         init = 0.0) {
00072         while(first != last) {
00073             init += adj(*first) * (*first2);
00074             ++first; ++first2;
00075         }
00076         return init;
00077     }
00078 
00079   template<typename Iterator, typename Iterator2, typename Rank2Iterator>
00080     inline typename scalar_type_
00081            <typename std::iterator_traits<Iterator>::value_type>::eval
00082     dot(Iterator first, Iterator last, Iterator2 first2,
00083         Rank2Iterator it_metric,
00084         typename scalar_type_
00085         <typename std::iterator_traits<Iterator>::value_type>::eval
00086         init = 0.0) {
00087         while(first != last) {
00088             init += adj(*first) *
00089                 vector_vector_product(it_metric->begin(), it_metric->end(),
00090                                       first2);
00091             ++first;
00092             ++it_metric;
00093         }
00094         return init;
00095     }
00096 
00097   template<typename Iterator>
00098     inline typename norm_type_
00099            <typename std::iterator_traits<Iterator>::value_type>::eval
00100     norm(Iterator first, Iterator last,
00101          typename norm_type_
00102          <typename std::iterator_traits<Iterator>::value_type>::eval
00103          init = zero) {
00104         while(first != last) {
00105             init += norm(*first);
00106             ++first;
00107         }
00108         return init;
00109     }
00110 
00111   template<typename Iterator, typename Rank2Iterator>
00112     inline typename norm_type_
00113            <typename std::iterator_traits<Iterator>::value_type>::eval
00114     norm(Iterator first, Iterator last,
00115          Rank2Iterator it_metric,
00116          typename norm_type_
00117          <typename std::iterator_traits<Iterator>::value_type>::eval
00118          init = 0.0) {
00119         Iterator first2 = first;
00120         while(first != last) {
00121             init += real(adj(*first)*vector_vector_product
00122                          (it_metric->begin(), it_metric->end(), first2));
00123             ++first;
00124             ++it_metric;
00125         }
00126         return init;
00127     }
00128 
00129   template<typename Iterator>
00130     inline void normalize(Iterator first, Iterator last) {
00131         typename norm_type_
00132                  <std::iterator_traits<Iterator>::value_type>::eval
00133             x = 1.0/sqrt(norm(first, last));
00134         while(first != last) { *first *= x; ++first; }
00135     }
00136 
00137   template<typename Iterator, typename Rank2Iterator>
00138     inline void normalize(Iterator first, Iterator last,
00139                           Rank2Iterator it_metric) {
00140         typename norm_type_
00141                  <std::iterator_traits<Iterator>::value_type>::eval
00142             x = 1.0/sqrt(norm(first, last, it_metric));
00143         while(first != last) { *first *= x; ++first; }
00144     }
00145 
00146   template<typename Iterator1, typename Iterator2, typename Scale>
00147     inline void assign_add_scaled(Iterator1 it1, Iterator1 it1_end,
00148                                   Iterator2 it2, Scale scale) {
00149         while(it1 != it1_end) {
00150             *it1 += scale*(*it2);
00151             ++it1; ++it2;
00152         }
00153     }
00154 
00155   template<typename Rank2Iterator>
00156     inline void orthonormalize(Rank2Iterator first, Rank2Iterator last) {
00157         Rank2Iterator it = first;
00158         while(it != last) { normalize(it->begin(), it->end()); ++it; }
00159         it = first;
00160         while(it != last) {
00161             for(Rank2Iterator it2 = first; it2 != it; ++it2)
00162                 assign_add_scaled(
00163                     it->begin(), it->end(), it2->begin(),
00164                     -dot(it2->begin(), it2->end(), it->begin()));
00165             normalize(it->begin(), it->end());
00166             ++it;
00167         }
00168     }
00169 
00170   template<typename Rank2Iterator, typename Rank2Iterator2>
00171     inline void orthonormalize(Rank2Iterator first, Rank2Iterator last,
00172                                Rank2Iterator2 it_metric) {
00173         Rank2Iterator it = first;
00174         while(it != last) {
00175             normalize(it->begin(), it->end(), it_metric);
00176             ++it;
00177         }
00178         it = first;
00179         while(it != last) {
00180             for(Rank2Iterator it2 = first; it2 != it; ++it2)
00181                 assign_add_scaled(
00182                     it->begin(), it->end(), it2->begin(),
00183                     -dot(it2->begin(), it2->end(), it->begin(),
00184                          it_metric));
00185             normalize(it->begin(), it->end(), it_metric);
00186             ++it;
00187         }
00188     }
00189 
00190   template<typename Iterator, typename Rank2Iterator>
00191     inline void project(Iterator first, Iterator last,
00192                         Rank2Iterator it_space_begin,
00193                         Rank2Iterator it_space_end) {
00194         typedef scalar_type_
00195                 <std::iterator_traits<Iterator>::value_type>::eval
00196                 scalar_type;
00197         std::list<scalar_type> coeff;
00198         for(Rank2Iterator it_space = it_space_begin;
00199             it_space != it_space_end; ++it_space)
00200             coeff.push_back(dot(it_space->begin(), it_space->end(), first));
00201 
00202         normalize(coeff.begin(), coeff.end());
00203         std::list<scalar_type>::iterator it_coeff = coeff.begin();
00204         std::fill(first, last, 0.0);
00205         for(Rank2Iterator it_space = it_space_begin;
00206             it_space != it_space_end; ++it_space) {
00207             assign_add_scaled(first, last, it_space->begin(), *it_coeff);
00208             ++it_coeff;
00209         }
00210     }
00211 
00212   template<typename Rank2Iterator, typename Rank2Iterator2>
00213     inline void project_collection(Rank2Iterator first, Rank2Iterator last,
00214                                    Rank2Iterator2 it_space_begin,
00215                                    Rank2Iterator2 it_space_end) {
00216         while(first != last) {
00217             project(first->begin(), first->end(),
00218                     it_space_begin, it_space_end);
00219             ++first;
00220         }
00221     }
00222 
00223   template<typename Iterator, typename Iterator2>
00224     inline
00225     typename times_type_
00226         < typename std::iterator_traits<Iterator>::value_type,
00227           typename std::iterator_traits<Iterator2>::value_type >::eval
00228     vector_vector_product(Iterator it, Iterator it_end, Iterator2 it2) {
00229         typename mutable_
00230           < typename times_type_
00231             < typename std::iterator_traits<Iterator>::value_type,
00232               typename std::iterator_traits<Iterator2>::value_type
00233             >::eval >::eval
00234             res = zero;
00235         while(it != it_end) {
00236             res += (*it) * (*it2);
00237             ++it;
00238             ++it2;
00239         }
00240         return res;
00241     }
00242 
00243   template<typename Rank2Iterator, typename Iterator>
00244     inline void
00245     multiply_matrix_to_vector(Rank2Iterator it, Rank2Iterator it_end,
00246                               Iterator it_dest, Iterator it_dest_end) {
00247         typedef typename std::iterator_traits<Iterator>::value_type
00248                 value_type;
00249         std::vector<value_type> src(it_dest, it_dest_end);
00250         while(it != it_end) {
00251             *it_dest = vector_vector_product(it->begin(), it->end(),
00252                                              src.begin());
00253             ++it;
00254             ++it_dest;
00255         }
00256     }
00257 
00258   template<typename Rank2Iterator, typename Rank2Iterator2>
00259     inline void
00260     multiply_matrix_to_matrix(Rank2Iterator it, Rank2Iterator it_end,
00261                               Rank2Iterator2 it2, Rank2Iterator2 it2_end) {
00262         while(it2 != it2_end) {
00263             multiply_matrix_to_vector(it, it_end, it2->begin(), it2->end());
00264             ++it2;
00265         }
00266     }
00267 }}
00268 
00269 #endif

Generated on Sat Sep 7 19:11:20 2002 for more with Doxygen 1.2.13.1. Doxygen 1.2.13.1 is written and copyright 1997-2002 by Dimitri van Heesch.