Почему эти два фрагмента генерируют разное значение?

Question

Когда я запускаю следующий код:

#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; 
} 

Выход:

4.47920533832942346919 
4.47920533832942524555 

не должны два числа такой же? И еще интересно, следующий код:

#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; 
} 

производит:

5.17737751763962084084 
5.17737751763962084084 

Так почему же они разные, то и же сейчас?

Answer 1

Во-первых, исправьте свой код. По стандарту C %lf не является директором для * printf ('l' недействителен, тип данных остается двойным). Чтобы напечатать длинный двойной, следует использовать %Lf. С вашим вариантом %lf можно попасть в ошибку с неправильным форматом, сокращаемым значением и т. Д. (Кажется, вы работаете с 32-разрядной средой: в 64 бит, как Unix, так и Windows проходят дважды в регистры XMM, но длинный double otherwhere - стек для Unix, память по указателю для Windows. В Windows/x86_64 код будет segfault, потому что callee ожидает указатель. Но, с Visual Studio, длинный двойной AFAIK aliased удваивается, поэтому вы можете не знать об этом изменении.)

Во-вторых, вы не можете быть уверены, что этот код не оптимизирован вашим компилятором C для вычисления времени компиляции (что может быть сделано с большей точностью, чем по умолчанию). Чтобы избежать такой оптимизации, отметьте sum как изменчивый.

С учетом этих изменений, ваш код показывает:

В Linux/amd64, gcc4.8:

для 50:

4.47920533832942505776 
4.47920533832942505820 

для 100:

5.17737751763962026144 
5.17737751763962025971 

На FreeBSD/i386, gcc4.8, без настройки точности или с явным fpsetprec (FP_PD):

4.47920533832942346919 
4.47920533832942524555 

5.17737751763962084084 
5.17737751763962084084 

(то же, что и в вашем примере);

, но тот же самый тест на FreeBSD с fpsetprec (FP_PE), который переключает FPU в реальные длинные двойные операции:

4.47920533832942505776 
4.47920533832942505820 

5.17737751763962026144 
5.17737751763962025971 

идентичен случае Linux; так что в real длинный двойной, есть реальная разница с 100 слагаемыми, и это, по здравому смыслу, больше, чем на 50. Но ваша платформа по умолчанию имеет округление до двойного.

И, наконец, в общем, это хорошо известный эффект конечной точности и последующего округления. Например, в this classical book это обход уменьшения суммы номеров номеров объясняется в первых главах.

Я не готов сейчас исследовать источник результатов с 50 слагаемыми и округлять до двух, почему он показывает такую ​​огромную разницу и почему эта разница компенсируется 100 слагаемыми.Это требует гораздо более глубокого расследования, чем я могу себе позволить сейчас, но, надеюсь, этот ответ ясно показывает вам следующее место, чтобы копать.

UPDATE: если это Windows, вы можете управлять режимом FPU с помощью _controlfp() и _controlfp_s(). В Linux, _FPU_SETCW делает то же самое. This description детализирует некоторые детали и дает пример кода.

UPDATE2: использование Kahan summation дает стабильные результаты во всех случаях. Ниже показано 4 значения: по возрастанию i, без KS; возрастание i, KS; нисходящий i, no KS; по убыванию я, КС:

50 и ФПУ дважды:

4.47920533832942346919 4.47920533832942524555 
4.47920533832942524555 4.47920533832942524555 

100 и ФПУ дважды:

5.17737751763962084084 5.17737751763961995266 
5.17737751763962084084 5.17737751763961995266 

50 и FPU долго дважды:

4.47920533832942505776 4.47920533832942524555 
4.47920533832942505820 4.47920533832942524555 

100 и FPU для удвоения:

5.17737751763962026144 5.17737751763961995266 
5.17737751763962025971 5.17737751763961995266 

вы можете видеть, что разница исчезла, результаты стабильны. Я бы предположил, что это почти конечная точка, которую можно добавить здесь :)