0

I'm on Manjaro 64 bit, latest edition. HP pavilion g6, Codeblocks Release 13.12 rev 9501 (2013-12-25 18:25:45) gcc 5.2.0 Linux/unicode - 64 bit.

There was a discussion between students on why

  • sn = 1/n diverges
  • sn = 1/n^2 converges

So decided to write a program about it, just to show them what kind of output they can expect

#include <iostream>
#include <math.h>
#include <fstream>

using namespace std;

int main()
{
   long double sn =0, sn2=0; // sn2 is 1/n^2
   ofstream myfile;
   myfile.open("/home/Projects/c++/test/test.csv");

    for (double n =2; n<100000000;n++){
            sn += 1/n;
            sn2 += 1/pow(n,2);
            myfile << "For n = " << n << " Sn = " << sn << " and Sn2 = " << sn2 << endl;

    }

    myfile.close();
    return 0;
}

Starting from n=9944 I got sn2 = 0.644834, and kept getting it forever. I did expect that the compiler would round the number and ignore the 0s at some point, but this is just too early, no?

So at what theoretical point does 0s start to be ignored? And what to do if you care about all 0s in a number? If long double doesn't do it, then what does?

I know it seems like a silly question but I expected to see a longer number, since you can store big part of pi in long doubles. By the way same result for double too.

Lynob
  • 5,059
  • 15
  • 64
  • 114
  • It's all in how you are outputting the number. – David Schwartz Nov 24 '15 at 19:54
  • @DavidSchwartz so doesn't fit on line? doesnt make sense, Can't C++ write on next line if the first line doesnt fit? – Lynob Nov 24 '15 at 19:55
  • The default precision for output is 6 digits. You have to ask for more, if you want to see the rest, for example `cout << setprecision(15);` – Bo Persson Nov 24 '15 at 19:55
  • @BoPersson ahhhhh awesome and can I ask it to output till the end of the number, however long it might be? – Lynob Nov 24 '15 at 19:56
  • `n * n` can only be more accurate than `pow(n, 2)`. There is no reason to write the latter, although if you are “lucky” the compiler will transform it into the former (for some definition of “lucky”). – Pascal Cuoq Nov 24 '15 at 19:57
  • @PascalCuoq good to know!!! – Lynob Nov 24 '15 at 19:58
  • @Lynob - It doesn't get much "longer" than that. Perhaps 1 or 2 additional digits in best case. – Bo Persson Nov 24 '15 at 20:00
  • 1
    just a sidenote: If you really want to understand why one converges and the other not you should do some maths. Convergence/Divergence by definition cannot be proven by a program running for a finite time. – 463035818_is_not_an_ai Nov 24 '15 at 20:03
  • @tobi303 sure thing, I know, but just wanted to give an example – Lynob Nov 24 '15 at 20:11
  • @Lynob And exactly that is a dubious method here. The analytical proof is essentially a one-liner, and just requires knowledge of Riemann sums and two basic rules of integration. – decltype_auto Nov 24 '15 at 21:55
  • @decltype_auto yes yes, this is not a development project or anything, just trying to show students something so wrote it like that. not trying to do it the smart way or be a genius or anything – Lynob Nov 24 '15 at 23:31
  • @Lynob It is obvious that your "code" doesn't meet the standards of a serious project's code. But that's not the point. The point is, that numerics are not meant as a replacement for analytics, but as a supplement. – decltype_auto Nov 25 '15 at 00:51
  • @decltype_auto thats not the point of my program nor the point of my question. my program is just to show students what to expect, because seeing is believing, and my question is about why long double not outputting the correct value. As for SO users, they don't have to concern themselves with the purpose of my program, they should just answer the question, perhaps i'm writing a game that requires me to write that code, in order for the player to gain more coins or something. – Lynob Nov 25 '15 at 17:54
  • @decltype_auto if i were asking about my code, I would have asked it in chatroom, after all you are talking to a man who have 1k+ rep so clearly I've been here long enough to know what's a good question and what should be closed, and I know that my code does not meet the standards, but again, the question is not about my code. My only fault is trying to help some students by visualizing what they learned about show it here. Had I asked the same question but replaced my code with Nasa's Spike scheduling algorithm, then would have received more upvotes,next time ill lie and make my Q look sexy – Lynob Nov 25 '15 at 18:15
  • @Lynob Very well - sorry for having to put it more clearly then, but your "program" is just poor maths clad in poor code. Your "students" may benefit from it nevertheless if they understand it as an example for how **not** not to do maths and how **not** to code. – decltype_auto Nov 25 '15 at 22:30

2 Answers2

7

The code that you wrote suffers from a classic programming mistake: it sums a sequence of floating-point numbers by adding larger numbers to the sum first and smaller numbers later.

This will inevitably lead to precision loss during addition, since at some point in the sequence the sum will become relatively large, while the next member of the sequence will become relatively small. Adding a sufficiently small floating-point value to a sufficiently large floating-point sum does not affect the sum. Once you reach that point, it will look as if the addition operation is "ignored", even though the value you attempt to add is not zero.

You can observe the same effect if you try calculating 100000000.0f + 1 on a typical machine: it still evaluates to 100000000. This does not happen because 1 somehow gets rounded to zero. This happens because the mathematically-correct result 100000001 is rounded back to 100000000. In order to force 100000000.0f to change through addition, you need to add at least 5 (and the result will be "snapped" to 100000008).

So, the issue here is not that the compiler "rounds the number when it gets so small", as you seem to believe. Your 1/pow(n,2) number is probably fine and sufficiently precise (not rounded to 0). The issue here is that at some iteration of your cycle the small non-zero value of 1/pow(n,2) just cannot affect the sum anymore.

While it is true that adjusting output precision will help you to see better what is going on (as stated in the comments), the real issue is what is described above.

When calculating sums of floating-point sequences with large differences in member magnitudes, you should do it by adding smaller members of the sequence first. Using my 100000000.0f example again, you can easily see that 4.0f + 4.0f + 100000000.0f correctly produces 100000008, while 100000000.0f + 4.0f + 4.0f is still 100000000.

AnT stands with Russia
  • 312,472
  • 42
  • 525
  • 765
  • This isn't actually the problem. There's way more than enough precision available to compute this sum to 6 digits the naive way. The OP probably just didn't wait long enough. (Also, rather than starting with smaller numbers, I'd probably go for [compensated summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm).) – user2357112 Nov 24 '15 at 20:40
2

You're not running into precision issues here. The sum doesn't stop at 0.644834; it keeps going to roughly the correct value:

#include <iostream>
#include <math.h>
using namespace std;

int main() {
    long double d = 0;
    for (double n = 2; n < 100000000; n++) {
        d += 1/pow(n, 2);
    }
    std::cout << d << endl;
    return 0;
}

Result:

0.644934

Note the 9! That's not 0.644834 any more.

If you were expecting 1.644934, you should have started the sum at n=1. If you were expecting visible changes between successive partial sums, you didn't see those because C++ is truncating the representation of the sums to 6 significant digits. You can configure your output stream to display more digits with std::setprecision from the iomanip header:

myfile << std::setprecision(9);
user2357112
  • 260,549
  • 28
  • 431
  • 505