C++ Theta function implementation

Question

I am trying to implement this function:

but it's not working. A minimal, verifiable example looks like:

#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] {sin(PI / N * (i + 1 - micro))};
    double den[countN] {0.5};

    for (int m=1; m<4; ++m)
    {
        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));

    }
    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).

I won't claim that this is the "best" of questions, but if it's going to be down-voted, at least let me know why so I can correct it. What I want is not for someone else to do it for me, but to let me know my error in implementation. I should have added I am a beginner in C++. — a concerned citizen, Jun 22 '16 at 11:02
I didn't vote, but if you're going to show us code, then show only the code that you intend us to look at, not code that includes lines that should be *ignored*. — eerorika, Jun 22 '16 at 11:09
There are only two std::cout lines in there and they only print out the temporary steps, I didn't think they would get in the way. But I'll delete the lines, I'm sorry for the confusion. — a concerned citizen, Jun 22 '16 at 11:12
the problem statement should be in the question, not referenced in a slide of a pdf somewhere. — UmNyobe, Jun 22 '16 at 11:15
I recommend using temporary variables, and splitting that monster of a statement to multiple trivial statements (goes for both numerator and denumerator) to more easily verify correctness. — eerorika, Jun 22 '16 at 11:26
I split the inner for loop to have one `for(m=0..N)` for `num` and one `for(1..N)` for `den`, while leaving the initial variables uninitialized (0), and now it works! I don't understand, because (in the same pdf) there is also a `sinh/cosh` version and that one is implemented like in the OP, and it works! Only, that one is static, not a function of `[i]`. Well, the good news it works, now that I split it, thank you for the guides. Do you think `for(m) { a=...; b=...; }` is the same as `for(m) { a=...; }` and `for(n) { b=...; }` ? — a concerned citizen, Jun 22 '16 at 11:34
@aconcernedcitizen `double den[countN] {0.5};` means `double den[countN] {0.5,0.0,0.0,0.0};`, not that `num` or `den` need to be arrays if you're going to just use one value — PeterT, Jun 22 '16 at 11:37
Your code doesn't have a `main()`, so nobody can compile and compare results. You should [edit] your question to give it a [mcve]. — Toby Speight, Jun 22 '16 at 11:39
@PeterT Ah! I knew the initialization was like that, but I forgot to count it for greater `i`! which means the first term is right, but the rest have a `0.5` less, each! That's my mistake! Thank you! This will not be a lesson easily forgotten. — a concerned citizen, Jun 22 '16 at 11:40
@TobySpeight I'm sorry for the lack of a complete code. Now that I know my error, rather than edit a 4th time, could this be marked as solved? PeterT is the one that spotted the error. — a concerned citizen, Jun 22 '16 at 11:48
@TobySpeight I modified it, but it's as good as solved now. Thank you all that helped. — a concerned citizen, Jun 22 '16 at 12:33
I turned the comment by @Peter into an answer; accept it if it solves your problem. Note that you shouldn't change the erroneous code in your question, as that means the answer no longer directly addresses it - I'll change that back for you. — Toby Speight, Jun 22 '16 at 16:19

Toby Speight · Accepted Answer · 2016-06-23T10:28:16.367

1

You had

double den[countN] {0.5};

this initializes the first element of den to 0.5 and all the other elements to 0.0 (default initialization). In other words, the above is equivalent to

double den[countN] {0.5, 0.0, 0.0, 0.0};

with as many zeros as necessary to fill the array. You probably wanted to initialize all the elements to 0.5. In your case, the easiest way to do that is when you first use that element - or, since you access only the single element den[i] during the lifetime of den, make it a plain double rather than an array:

for (int i=0; i<countN; ++i) {
    double micro {N % 2 ? 0.0 : 0.5};  // 0 for odd N, 1/2 for even N
    double num{sin(PI / N * (i + 1 - micro))};
    double den{0.5};

    for (int m=1; m<4; ++m) {
        num += pow(-1, m) * pow(q, m*(m + 1)) * sin((2 * m + 1) * PI / N * (i + 1 - micro));
        den += pow(-1, m) * pow(q, m*m) * cos(2 * m * PI / N * (i + 1 - micro));
    }
    Omega[i] = fabs(pow(q, 0.25) * num / den);
}

edited Jun 23 '16 at 10:28

answered Jun 22 '16 at 16:18

Toby Speight

27,591
48
66
103

Thank you for the solution (and sorry for the mix up). It didn't occur to me that `num` and `den` are only local to `for(m)`, maybe because I had initialized them before. At any rate, there is also a `σ₀` which uses the same sums but with `sinh/cosh`, and with `Λ` instead of `PI/N*(i+1-µ)`, and there I had used fixed values (non-array), exactly the way you show up in there. But, my bad reasoning was that I should use an array since `Ω_i` is a function of `i`. Marked as solved. – a concerned citizen Jun 23 '16 at 07:40
Still, with this occasion, I also found out I can initialize a vector like this: `std::vector matrix(nrOfTimes, value)`. If I use `std::vector` as a fixed array, non-dynamically resized, is there a speed penalty? – a concerned citizen Jun 23 '16 at 07:41
That's a question of its own - see if it's already been answered here on StackOverflow, and if you can't find an answer, then please ask it as a new question. You're almost certainly not the only person interested to know! – Toby Speight Jun 23 '16 at 10:26
You're right, in fact, here's my answer: https://stackoverflow.com/questions/3664272/is-stdvector-so-much-slower-than-plain-arrays – a concerned citizen Jun 23 '16 at 18:37

C++ Theta function implementation

1 Answers1