/******************************************************************************* * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * * Copyright (C) 2012 by M. Ferrero, O. Parcollet, I. Krivenko * * TRIQS 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 3 of the License, or (at your option) any later * version. * * TRIQS 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 * TRIQS. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ #ifndef TRIQS_UTILITY_PADE_APPROXIMANTS_H #define TRIQS_UTILITY_PADE_APPROXIMANTS_H #include "pade_approximants.hpp" #include <triqs/arrays.hpp> #include <gmpxx.h> namespace triqs { namespace utility { typedef std::complex<double> dcomplex; // This implementation is based on a Fortran code written by // A. Poteryaev <Alexander.Poteryaev _at_ cpht.polytechnique.fr> // // The original algorithm is described in // H. J. Vidberg, J. W. Serene. J. Low Temp. Phys. 29, 3-4, 179 (1977) struct gmp_complex { mpf_class re, im; gmp_complex operator* (const gmp_complex &rhs){ return { rhs.re*re-rhs.im*im, rhs.re*im+rhs.im*re }; } friend gmp_complex inverse (const gmp_complex &rhs){ mpf_class d=rhs.re*rhs.re + rhs.im*rhs.im; return { rhs.re/d, -rhs.im/d}; } gmp_complex operator/ (const gmp_complex &rhs){ return (*this)*inverse(rhs); } gmp_complex operator+ (const gmp_complex &rhs){ return { rhs.re + re, rhs.im + im }; } gmp_complex operator- (const gmp_complex &rhs){ return { re - rhs.re, im - rhs.im }; } friend mpf_class real(const gmp_complex &rhs) { return rhs.re; } friend mpf_class imag(const gmp_complex &rhs) { return rhs.im; } gmp_complex& operator= (const std::complex<double> &rhs) { re = real(rhs); im = imag(rhs); return *this; } friend std::ostream & operator << (std::ostream & out,gmp_complex const & r) { return out << " gmp_complex("<<r.re<<","<<r.im<<")"<<std::endl ;} }; class pade_approximant { arrays::vector<dcomplex> z_in; // Input complex frequency points arrays::vector<dcomplex> a; // Pade coefficients public: static const int GMP_default_prec = 256; // Precision of GMP floats to use during a Pade coefficients calculation. pade_approximant(const arrays::vector<dcomplex> & z_in_, const arrays::vector<dcomplex> & u_in): z_in(z_in_), a(z_in.size()) { int N = z_in.size(); // Change the default precision of GMP floats. unsigned long old_prec = mpf_get_default_prec(); mpf_set_default_prec(GMP_default_prec); // How do we determine it? arrays::array<gmp_complex,2> g(N,N); gmp_complex MP_0 = {0.0, 0.0}; g() = MP_0; for (int f = 0; f<N; ++f) { g(0,f) = u_in(f); }; gmp_complex MP_1 = {1.0, 0.0}; for(int p=1; p<N; ++p) for(int j=p; j<N; ++j) { gmp_complex x = g(p-1,p-1)/g(p-1,j) - MP_1; gmp_complex y; y = z_in(j)-z_in(p-1); g(p,j) = x/y; } for(int j=0; j<N; ++j) { gmp_complex gj = g(j,j); a(j) = dcomplex(real(gj).get_d(), imag(gj).get_d()); } // Restore the precision. mpf_set_default_prec(old_prec); } // give the value of the pade continued fraction at complex number e dcomplex operator()(dcomplex e) const { dcomplex A1(0); dcomplex A2 = a(0); dcomplex B1(1.0), B2(1.0); int N = a.size(); for(int i=0; i<=N-2; ++i){ dcomplex Anew = A2 + (e - z_in(i))*a(i+1)*A1; dcomplex Bnew = B2 + (e - z_in(i))*a(i+1)*B1; A1 = A2; A2 = Anew; B1 = B2; B2 = Bnew; } return A2/B2; } }; }} #endif