ERKALE
ERKALE - DFT from Hel
 All Classes Functions Variables Friends Pages
mathf.h
1 /*
2  * This source code is part of
3  *
4  * E R K A L E
5  * -
6  * DFT from Hel
7  *
8  * Written by Susi Lehtola, 2010-2011
9  * Copyright (c) 2010-2011, Susi Lehtola
10  *
11  * This program is free software; you can redistribute it and/or
12  * modify it under the terms of the GNU General Public License
13  * as published by the Free Software Foundation; either version 2
14  * of the License, or (at your option) any later version.
15  */
16 
17 
18 
19 #include "global.h"
20 
21 #ifndef ERKALE_MATH
22 #define ERKALE_MATH
23 
24 #include <armadillo>
25 #include <vector>
26 #include <cfloat>
27 #include <string>
28 
30 double doublefact(int n);
32 double fact(int n);
34 double fact_ratio(int i, int r);
35 
37 double fgamma(double x);
38 
40 double sinc(double x);
41 
43 double bessel_jl(int l, double x);
44 
46 double boysF(int m, double x);
48 void boysF_arr(int mmax, double x, arma::vec & bf);
49 
51 double hyperg_1F1(double a, double b, double x);
52 
54 double choose(int m, int n);
55 
57 int getind(int l, int m, int n);
58 
60 template <class T> T max(const std::vector<T> & x) {
61  if(x.size()==0) {
62  ERROR_INFO();
63  throw std::runtime_error("Trying to get maximum value of empty array!\n");
64  }
65 
66  T m=x[0];
67  for(size_t i=1;i<x.size();i++)
68  if(m<x[i])
69  m=x[i];
70  return m;
71 }
72 
74 template <class T> T max4(const T & a, const T & b, const T & c, const T & d) {
75  return std::max(std::max(a,b),std::max(c,d));
76 }
77 
79 template <typename T> size_t find(const std::vector<T> & list, const T & val) {
80  for(size_t i=0;i<list.size();i++)
81  if(list[i]==val)
82  return i;
83 
84  return std::string::npos;
85 }
86 
88 double max_abs(const arma::mat & R);
90 double max_cabs(const arma::cx_mat & R);
92 double rms_norm(const arma::mat & R);
94 double rms_cnorm(const arma::cx_mat & R);
95 
97 template <class T> void reverse(std::vector<T> & a) {
98  size_t i;
99  size_t n=a.size();
100  T temp;
101  for(i=0;i<n/2;i++) {
102  temp=a[i];
103  a[i]=a[n-1-i];
104  a[n-1-i]=temp;
105  }
106 }
107 
109 template <class T> T sum(const std::vector<T> & a) {
110  T sum=0;
111  for(size_t i=0;i<a.size();i++) {
112  sum+=a[i];
113  }
114  return sum;
115 }
116 
118 template <class T> bool abscomp(const std::complex<T> & a, const std::complex<T> & b) {
119  return std::abs(a) < std::abs(b);
120 }
121 
123 std::vector<double> spline_interpolation(const std::vector<double> & xt, const std::vector<double> & yt, const std::vector<double> & x);
125 arma::vec spline_interpolation(const arma::vec & xt, const arma::vec & yt, const arma::vec & x);
127 double spline_interpolation(const std::vector<double> & xt, const std::vector<double> & yt, double x);
128 
130 arma::mat randu_mat(size_t M, size_t N, unsigned long int seed=0);
132 arma::mat randn_mat(size_t M, size_t N, unsigned long int seed=0);
133 
135 arma::mat real_orthogonal(size_t N, unsigned long int seed=0);
137 arma::cx_mat complex_unitary(size_t N, unsigned long int seed=0);
138 
140 double round(double x, unsigned n);
141 
149 arma::vec find_minima(const arma::vec & x, const arma::vec & y, size_t runave=0, double thr=DBL_MAX);
150 
152 template<typename T> size_t sorted_insertion(std::vector<T> & v, T t) {
153  // Get upper bound
154  typename std::vector<T>::iterator high;
155  high=std::upper_bound(v.begin(),v.end(),t);
156 
157  // Corresponding index is
158  size_t ind=high-v.begin();
159 
160  if(ind>0 && v[ind-1]==t) {
161  // Value already exists in vector - return position
162  return ind-1;
163  } else {
164  // Value doesn't exist in vector - add it
165  typename std::vector<T>::iterator pos=v.insert(high,t);
166  // and return the position
167  return pos-v.begin();
168  }
169 }
170 
179 void stencil(double z, const arma::vec & x, arma::mat & w);
180 
181 #endif