mercredi 22 juin 2016

C++ Theta function implementation


I am trying to implement this function: formula, but it's not working. [EDIT] After user2079303's suggestion, PeterT's pointing out my error, and Toby Speight's suggestion, a mini-code that can be compiled to verify the validity of the function looks like this:

#include <iostream>
#include <cmath>

int main()
{
  int N {8};        // change this for testing <1..inf>
  double q {0.1 / N};
  int countN {static_cast<int>(floor(N / 2))};
  static const double PI {3.1415926535897932384626433832795};

  // Omega[i] =  Theta1(u,m) / Theta4(u,m)
  double Omega[countN];
  for (int i=0; i<countN; ++i)
  {
    double micro {!(N % 2) * 0.5};  // 0 for odd N, 1/2 for even N
    double num[countN];
    double den[countN];

    for (int m=0; m<4; ++m)
    {
      double tmp {0.5 + 0.5 * (m > 0)};
      num[i] += pow(-1, m) * pow(q, m*(m + 1)) * sin((2 * m + 1) * PI / N * (i + 1 - micro));
      den[i] += pow(-1, m) * pow(q, m*m) * cos(2 * m * PI / N * (i + 1 - micro)) * tmp;
      std::cout << tmp << "t" << num[i] << "t" << den[i] << "n";
    }
    Omega[i] = fabs(pow(q, 0.25) * num[i] / den[i]);
    std::cout << "  " << Omega[i] << "n";
  }

  // testing the values, they should be increasing in value
  for (const auto &elem: Omega)
    std::cout << elem << " ";
  std::cout << "n";

  return 0;
}

There is a minor simplification compared to the original: I factored 2 in both numerator and denominator and I used only the q^0.25 outside of the fraction. Also, countN is the r from the original document, micro is only the 1/2 for even N or 0 for odd N, and i is 0 for array's index but i+1 for calculations, but these shouldn't matter overall.

I tried this with wxMaxima:

Theta[1](x,y):=2*y^0.25*sum( (-1)^k*y^(k*(k+1))*sin((2*k+1)*x),k,0,n );
Theta[4](x,y):=1+2*sum( (-1)^k*y^(k^2)*cos(2*k*x),k,1,n );
n:4$
N:8$
a:0.05$
b(i):=%pi/N*(i-(1-mod(N,2))/2)$
for N:8 thru 9 do for i:1 thru N/2 do print(["N=",N,"i=",i],Theta[1](b(i),a)/Theta[4](b(i),a)),numer;

And the results, in C++:

(q=0.05; N=8)
Omega[0]=0.2018370065366672
Omega[1]=0.06058232646142273
Omega[2]=0.01205653570636574
Omega[3]=0.02127667733703158
(q=0.05; N=9)
Omega[0]=0.348078726440638
Omega[1]=0.1178366281313341
Omega[2]=2.559808325080287e-07
Omega[3]=0.02178788541277828

and in wxMaxima:

["N=",8,"i=",1]" "0.2018370065366672" "
["N=",8,"i=",2]" "0.5439269564954693" "
["N=",8,"i=",3]" "0.7569342043740249" "
["N=",8,"i=",4]" "0.850913653939989" "
["N=",9,"i=",1]" "0.348078726440638" "
["N=",9,"i=",2]" "0.6165773889432575" "
["N=",9,"i=",3]" "0.7800391631077094" "
["N=",9,"i=",4]" "0.8532352152763631

To my surprise, the first term is good, for bith N, so I can't tell what in my code is not right. Could someone help me spot the error?

To be clear about it: I am a beginner in C++ and I am not looking for someone to do it for me, but to let me know of my erros in coding (translating math to C++ code).


Aucun commentaire:

Enregistrer un commentaire