pade_approximants.hpp

Go to the documentation of this file.

// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2021 Simons Foundation
// Copyright (c) 2017 Igor Krivenko
//
// This program 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.
//
// This program 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 may obtain a copy of the License at
//     https://www.gnu.org/licenses/gpl-3.0.txt
//
// Authors: Michel Ferrero, Igor Krivenko, Olivier Parcollet, Hiroshi Shinaoka, Nils Wentzell

 
#pragma once
 
#include "./exceptions.hpp"
 
#include <gmpxx.h>
#include <triqs/arrays.hpp>
 
#include <complex>
#include <ostream>
 
namespace triqs::utility {
 
  // Alias for `std::complex<double>`.
  using dcomplex = std::complex<double>;


  struct gmp_complex {
    mpf_class re;

    mpf_class im;

    gmp_complex operator*(const gmp_complex &z) const { return {.re = z.re * re - z.im * im, .im = z.re * im + z.im * re}; }

    friend gmp_complex inverse(const gmp_complex &z) {
      mpf_class d = z.re * z.re + z.im * z.im;
      if (d == 0) TRIQS_RUNTIME_ERROR << "pade_approximant: GMP division by zero";
      return {.re = z.re / d, .im = -z.im / d};
    }

    gmp_complex operator/(const gmp_complex &z) const { return (*this) * inverse(z); }

    gmp_complex operator+(const gmp_complex &z) const { return {.re = z.re + re, .im = z.im + im}; }

    gmp_complex operator-(const gmp_complex &z) const { return {.re = re - z.re, .im = im - z.im}; }

    friend mpf_class real(const gmp_complex &z) { return z.re; }

    friend mpf_class imag(const gmp_complex &z) { return z.im; }

    [[nodiscard]] mpf_class norm() const { return real(*this) * real(*this) + imag(*this) * imag(*this); }

    gmp_complex &operator=(const std::complex<double> &z) {
      re = real(z);
      im = imag(z);
      return *this;
    }

    friend std::ostream &operator<<(std::ostream &out, gmp_complex const &z) {
      return out << " gmp_complex(" << z.re << "," << z.im << ")" << std::endl;
    }
  };

  class pade_approximant {
 
    nda::vector<dcomplex> z_in;
    nda::vector<dcomplex> a;
 
    public:
    static const int GMP_default_prec = 256;

    pade_approximant(const nda::vector<dcomplex> &z_in, const nda::vector<dcomplex> &u_in) : z_in(z_in), a(z_in.size()) {
 
      long N = z_in.size();
 
      // temporarily switch GMP's default precision
      unsigned long old_prec = mpf_get_default_prec();
      mpf_set_default_prec(GMP_default_prec);
 
      nda::array<gmp_complex, 2> g(N, N);
      gmp_complex MP_0 = {.re = 0.0, .im = 0.0};
      g()              = MP_0;
      for (long f = 0; f < N; ++f) { g(0, f) = u_in(f); }
 
      gmp_complex MP_1 = {.re = 1.0, .im = 0.0};
 
      for (long p = 1; p < N; ++p) {
 
        // truncate the continued fraction when |g| becomes very small
        if (g(p - 1, p - 1).norm() < 1.0e-20) break;
 
        for (long 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 (long 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);
    }

    dcomplex operator()(dcomplex z) const {
 
      dcomplex A1(0);
      dcomplex A2 = a(0);
      dcomplex B1(1.0);
 
      long N = a.size();
      for (long i = 0; i <= N - 2; ++i) {
        dcomplex Anew = A2 + (z - z_in(i)) * a(i + 1) * A1;
        dcomplex Bnew = 1.0 + (z - z_in(i)) * a(i + 1) * B1;
        A1            = A2 / Bnew;
        A2            = Anew / Bnew;
        B1            = 1.0 / Bnew;
      }
 
      return A2;
    }
  };

 
} // namespace triqs::utility