Why these two snippets generate different value?

Question

When I run the following code:

#include <stdio.h>

int main()
{
    int i = 0;
    volatile long double sum = 0;
    for (i = 1; i < 50; ++i) /* first snippet */
    {
        sum += (long double)1 / i;
    }
    printf("%.20Lf\n", sum);
    sum = 0;
    for (i = 49; i > 0; --i) /* second snippet */
    {
        sum += (long double)1 / i;
    }
    printf("%.20Lf", sum);
    return 0;
}

The output is:

4.47920533832942346919
4.47920533832942524555

Shouldn't the two numbers be same? And more interestingly, the following code:

#include <stdio.h>

int main()
{
    int i = 0;
    volatile long double sum = 0;
    for (i = 1; i < 100; ++i) /* first snippet */
    {
        sum += (long double)1 / i;
    }
    printf("%.20Lf\n", sum);
    sum = 0;
    for (i = 99; i > 0; --i) /* second snippet */
    {
        sum += (long double)1 / i;
    }
    printf("%.20Lf", sum);
    return 0;
}

produces:

5.17737751763962084084
5.17737751763962084084

So why are they different then and same now?

Your code has undefined behaviour, since you promise `printf` a double but pass it a long double. — Kerrek SB, Nov 28 '15 at 17:24
The first loop in the second snippet runs twice as much, so I don't really understand the problem. — cadaniluk, Nov 28 '15 at 17:25
Please refer to http://stackoverflow.com/questions/588004/is-floating-point-math-broken — Weather Vane, Nov 28 '15 at 17:32
floating point addition is neither associative nor commutative (although it is *approximately*) — John Coleman, Nov 28 '15 at 17:37
@JohnColeman All floating-point operations that should be commutative are. — Pascal Cuoq, Nov 28 '15 at 18:17
@PascalCuoq a+b = b+a holds but a+ b + c need not be c + a + b, although perhaps that is a failure of associativity in disguise. — John Coleman, Nov 28 '15 at 19:10
@JohnColeman a + b = b + a **is** the definition of commutativity. https://en.wikipedia.org/wiki/Commutative_property#Mathematical_definitions . a+ b + c need not be c + a + b indeed for + to be commutative — Pascal Cuoq, Nov 28 '15 at 22:09
@PascalCuoq You are correct. In the absence of associativity certain identities which otherwise follow from the commutative law no longer hold. In particular -- ordinary addition obeys a generalized commutative law whereby any permutation of the terms in a multi-term sum yields the same result. This is what I was thinking of in my original comment. The relatively straightforward proof by induction that commutative implies generalized commutative is easily seen to depend on the associative law. — John Coleman, Nov 28 '15 at 22:35
@sunqingyao No -- that doesn't follow. In the presence of commutativity, a + b + c = c + b + a is equivalent to associativity. If a = -b and c is close to but not equal to a (hence close to the negative of b) then a + b + c = c but you could have c + b + a != c because of loss of precision in c + b — John Coleman, Nov 29 '15 at 19:06
@sun qingyao: No, it isn't in floating point arithmetic. For example, `100000000.0f + 4.0f + 4.0f` is still `100000000.0f`. However, `4.0f + 4.0f + 100000000.0f` is `100000008.0f` — AnT stands with Russia, Dec 03 '15 at 06:19
@PascalCuoq according to [this thread](http://stackoverflow.com/questions/24442725/is-floating-point-addition-commutative-in-c) `a+b == b+a` is not guaranteed to hold — M.M, Dec 03 '15 at 08:53
@M.M That question is C++. `a + b == b + a` has to hold (except for NaN results) in any C99 compiler that sets `FLT_EVAL_METHOD` ≥ 0. Although C++ incorporates the definition of `FLT_EVAL_METHOD` through cfloat since 2003, the implementors of C++ compilers do not care. Still, this question is tagged C. — Pascal Cuoq, Dec 03 '15 at 09:06

Netch · Accepted Answer · 2015-12-10T15:34:42.793

First, please correct your code. By C standard, %lf isn't principal for *printf ('l' is void, the data type remains double). To print long double, one should use %Lf. With your variant %lf, it's possible to get into a bug with improper format, cut-down value, etc. (You seem running 32-bit environment: in 64 bits, both Unix and Windows pass double in XMM registers, but long double otherwhere - stack for Unix, memory by pointer for Windows. On Windows/x86_64, you code will segfault because callee expects pointer. But, with Visual Studio, long double is AFAIK aliased to double, so you can remain ignorant of this change.)

Second, you can't be sure this code is not optimized by your C compiler to compile-time calculations (which can be done with more precision than default run-time one). To avoid such optimization, mark sum as volatile.

With these changes, your code shows:

At Linux/amd64, gcc4.8:

for 50:

4.47920533832942505776
4.47920533832942505820

for 100:

5.17737751763962026144
5.17737751763962025971

At FreeBSD/i386, gcc4.8, without precision setting or with explicit fpsetprec(FP_PD):

4.47920533832942346919
4.47920533832942524555

5.17737751763962084084
5.17737751763962084084

(the same as in your example);

but, the same test on FreeBSD with fpsetprec(FP_PE), which switches FPU to real long double operations:

4.47920533832942505776
4.47920533832942505820

5.17737751763962026144
5.17737751763962025971

identical to Linux case; so, in real long double, there is some real difference with 100 summands, and it is, in accordance with common sense, larger than for 50. But your platform defaults to rounding to double.

And, finally, in general, this is well-known effect of a finite precision and consequent rounding. For example, in this classical book, this misrounding of decreasing number series sum is explained in the very first chapters.

I am not really ready now to investigate source of results with 50 summands and rounding to double, why it shows such huge difference and why this difference is compensated with 100 summands. That needs much deeper investigation than I can afford now, but, I hope, this answer clearly shows you a next place to dig.

UPDATE: if it's Windows, you can manipulate FPU mode with _controlfp() and _controlfp_s(). In Linux, _FPU_SETCW does the same. This description elaborates some details and gives example code.

UPDATE2: using Kahan summation gives stable results in all cases. The following shows 4 values: ascending i, no KS; ascending i, KS; descending i, no KS; descending i, KS:

50 and FPU to double:

4.47920533832942346919 4.47920533832942524555
4.47920533832942524555 4.47920533832942524555

100 and FPU to double:

5.17737751763962084084 5.17737751763961995266
5.17737751763962084084 5.17737751763961995266

50 and FPU to long double:

4.47920533832942505776 4.47920533832942524555
4.47920533832942505820 4.47920533832942524555

100 and FPU to long double:

5.17737751763962026144 5.17737751763961995266
5.17737751763962025971 5.17737751763961995266

you can see difference disappeared, results are stable. I would assume this is nearly final point that can be added here :)

By the original C standard `%lf` in `printf` is undefined. By C99 and later `%lf` in `printf` is equivalent to `%f`. Yet it is a good idea to stick to `%lf` for `double` and `%f` for `float`, for consistency with `scanf`. So, `%lf` is useful. Just not in this case. — AnT stands with Russia, Dec 03 '15 at 06:17
I'm also talkinfg about `printf` group. What I mean is that `%lf` is useful in `printf` as part of a convention: "`%lf` for `double` and `%f` for `float`". — AnT stands with Russia, Dec 09 '15 at 08:20

Why these two snippets generate different value?

1 Answers1

Linked