#include <iostream>
#include <complex>
 
// Smith's algorithm for robust complex division
std::complex<double> 
ComplexDiv(std::complex<double> const a, std::complex<double> b)
{
  if(std::abs(std::imag(b)) <= std::abs(std::real(b))) {
    double r   = std::imag(b) / std::real(b);
    double den = std::real(b) + (std::imag(b) * r); 
    return std::complex<double>
             (
             (std::real(a) + r * std::imag(a)) / den, 
             (std::imag(a) - r * std::real(a)) / den
             );
  } else {
    double r   = std::real(b) / std::imag(b);
    double den = std::real(b) * r + std::imag(b); 
    return std::complex<double>
             (
             (std::real(a) * r + std::imag(a)) / den,
             (std::imag(a) * r - std::real(a)) / den
             );  
  }
}
 
int main()
{
  std::complex<double> a = {1.0, 0.0};
  std::complex<double> b = {1.e-162, 0.0};
  std::cout << "Smith's Algorithm Result = " << ComplexDiv(a, b) << std::endl;
  std::cout << "Impementation Result     = " << a / b            << std::endl;
  return 0;
}