Search code examples
c++complex-numbers

Complex numbers arithmetic goes wrong

Right now I am trying to implement the formula for the probabilities of a Poisson binomial distribution.

The formula is the last one in that section, with complex exponentials inside them: enter image description here

I have a very simple code that should do it, but it is outputting wrong probabilities.

#include <iostream>
#include <complex>
#include <cmath>
using namespace std;

int main(int argc, char const *argv[]) {
    int N=5; //number of coins
    double probabilities[N]={0.1,0.1,0.1,0.1,0.1}; //probabilities of coin landing head
    double distr[N+1];
    for (int j=0; j<N+1;j++){
        complex<double> temp1=0.0 + 0.0i;
        for (int k=0; k<N+1;k++){
           for(int l=0;l<N+1;l++){
           complex<double> temp2=1.0 + 0.0i;
               for (int m=0;m<N;m++){
                   temp2=temp2*(1.0+0.0i+(exp(l*2*M_PI*1.0i/(double(N)+1.0))-1.0+0.0i)*probabilities[m]);
               }
               temp2=temp2*exp(l*k*2*M_PI*1.0i/(double(N)+1.0));
               temp1=temp1+temp2/(double(N)+1.0);
           }
        }
        distr[j]=real(temp1);
    }
    for (int i=0;i<N+1;i++){
    cout<< distr[i] << ' ';
  }

The output of this code is [1,1,1,1,1,1], which is definitely not correct. I am thinking that maybe I have not been working correctly with complex numbers, but I do not see where I did something wrong. It is frustrating that such a simple program does not work :(.

Solution

  • From the code it is clear that temp1 does not depend on j, hence you get the same numbers, the sum over k. After you remove the outer loop over j and write distr[k] = real(temp1); and fix the sign in exp(l*k*...), you'll get the expected result:

    0.59049 0.32805 0.0729 0.0081 0.00045 1e-05

    Full code with some simplifications:

    int main() {
    using namespace std::literals::complex_literals;

    constexpr int N = 5;
    const double probabilities[N] = {.1, .1, .1, .1, .1};
    const auto c = std::exp(2i * M_PI / (N + 1.));

    double distr[N + 1];
    for (int k = 0; k <= N; ++k) {
        auto sum = std::complex<double>(0);
        for(int l = 0; l <= N; ++l) {
            auto prod = std::complex<double>(1);
            for (int m = 0; m < N; ++m)
                prod *= 1. + (std::pow(c, l) - 1.) * probabilities[m];
            sum += prod * std::pow(c, -l * k) / (N + 1.);
        }
        distr[k] = std::real(sum);
    }

    for (auto d : distr)
        std::cout << d << ' ';
}