#pragma once

#include <complex>
#include <random>

template <class T = double> class complex_normal_distribution {
private:
  std::normal_distribution<T> d;
  T δx;
  T δy;
  std::complex<T> μ;
  std::complex<T> eθ;

public:
  complex_normal_distribution(std::complex<T> μ, T σ, std::complex<T> κ) : μ(μ), d(0, σ / sqrt(2)) {
    δx = sqrt(1 + std::abs(κ));
    δy = sqrt(1 - std::abs(κ));
    eθ = std::polar((T)1, std::arg(κ));
  }

  template <class Generator> std::complex<T> operator()(Generator& g) {
    // First generate an axis-aligned complex number centered at the origin.
    std::complex<T> z(d(g) * δx, d(g) * δy);
    // Shift and rotate the number.
    return μ + eθ * z;
  }
};