Я * не * хочу правильное округление для функции exp
GCC- реализация математической библиотеки C в системах Debian, по-видимому, (IEEE 754-2008) совместима с реализацией функции exp, подразумевая, что округление всегда должно быть правильным:
( из Википедии) Стандарт IEEE с плавающей запятой гарантирует, что сложение, вычитание, умножение, деление, объединенное умножение-сложение, квадратный корень и остаток с плавающей запятой дадут правильно округленный результат операции с бесконечной точностью. В стандарте 1985 года такая гарантия не была предоставлена для более сложных функций, и они, как правило, точны только с точностью до последнего бита. Тем не менее, стандарт 2008 года гарантирует, что соответствующие реализации дадут правильно округленные результаты, которые относятся к активному режиму округления; реализация функций, однако, не является обязательной.
Оказывается, я сталкиваюсь со случаем, когда эта функция фактически мешает, потому что точный результат exp функция часто почти точно посередине между двумя последовательными double значения (1), а затем программа выполняет множество дополнительных вычислений, теряя до 400 (!) скорости: на самом деле это было объяснением моего (плохо заданного:-S) вопроса № 43530011.
(1) Точнее, это происходит, когда аргумент exp оказывается, имеет вид (2 k + 1) × 2 -53, где k - довольно маленькое целое число (например, 242). В частности, вычисления, используемые pow (1. + x, 0.5) иметь тенденцию звонить exp с таким аргументом, когда x имеет порядок 2 -44.
Поскольку реализация правильного округления может быть очень трудоемкой в определенных обстоятельствах, я предполагаю, что разработчики также разработают способ получить немного менее точный результат (скажем, только до 0,6 ULP или что-то подобное) за раз который (приблизительно) ограничен для каждого значения аргумента в данном диапазоне... (2)
... а как это сделать??
(2) Я имею в виду, что я просто не хочу, чтобы некоторые исключительные значения аргумента, такие как (2 k + 1) × 2 -53, занимали бы гораздо больше времени, чем большинство значений того же порядка; но, конечно, я не возражаю, если некоторые исключительные значения аргумента идут намного быстрее, или если большие аргументы (в абсолютном значении) требуют большего времени вычисления.
Вот минимальная программа, демонстрирующая это явление:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
int main (void)
{
int i;
double a, c;
c = 0;
clock_t start = clock ();
for (i = 0; i < 1e6; ++i) // Doing a large number of times the same type of computation with different values, to smoothen random fluctuations.
{
a = (double) (1 + 2 * (rand () % 0x400)) / 0x20000000000000; // "a" has only a few significant digits, and its last non-zero digit is at (fixed-point) position 53.
c += exp (a); // Just to be sure that the compiler will actually perform the computation of exp (a).
}
clock_t stop = clock ();
printf ("%e\n", c); // Just to be sure that the compiler will actually perform the computation.
printf ("Clock time spent: %d\n", stop - start);
return 0;
}
Теперь после gcc -std=c99 program53.c -lm -o program53:
$ ./program53
1.000000e+06
Clock time spent: 13470008
$ ./program53
1.000000e+06
Clock time spent: 13292721
$ ./program53
1.000000e+06
Clock time spent: 13201616
С другой стороны, с program52 а также program54 (получил заменой 0x20000000000000 соответственно 0x10000000000000 а также 0x40000000000000):
$ ./program52
1.000000e+06
Clock time spent: 83594
$ ./program52
1.000000e+06
Clock time spent: 69095
$ ./program52
1.000000e+06
Clock time spent: 54694
$ ./program54
1.000000e+06
Clock time spent: 86151
$ ./program54
1.000000e+06
Clock time spent: 74209
$ ./program54
1.000000e+06
Clock time spent: 78612
Осторожно, это явление зависит от реализации! По-видимому, среди распространенных реализаций только те из систем Debian (включая Ubuntu) демонстрируют это явление.
P.-S.: Я надеюсь, что мой вопрос не является дубликатом: я тщательно и безуспешно искал похожий вопрос, но, возможно, я заметил, что использовал соответствующие ключевые слова…:-/
3 ответа
Чтобы ответить на общий вопрос о том, почему функции библиотеки должны давать правильно округленные результаты:
С плавающей точкой трудно, и часто противоречит интуиции. Не каждый программист прочитал то, что они должны иметь. Когда библиотеки допускали немного неточное округление, люди жаловались на точность библиотечной функции, когда их неточные вычисления неизбежно давали сбой и приводили к бессмыслице. В ответ авторы библиотеки сделали свои библиотеки точно округленными, так что теперь люди не могут переложить вину на них.
Во многих случаях конкретные знания об алгоритмах с плавающей запятой могут привести к значительному улучшению точности и / или производительности, как в тестовом примере:
Принимая exp() чисел очень близко к 0 в числах с плавающей точкой проблематично, так как результатом является число, близкое к 1 в то время как вся точность находится в разнице с единицей, поэтому самые значимые цифры теряются. Более точное (и значительно более быстрое в этом тесте) вычисление exp(x) - 1 через математическую библиотеку C expm1(x), Если exp() само по себе действительно нужно, это еще намного быстрее сделать expm1(x) + 1,
Аналогичная проблема существует для вычислений log(1 + x), для которого есть функция log1p(x),
Быстрое исправление, ускоряющее предоставленный тестовый пример:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
int main (void)
{
int i;
double a, c;
c = 0;
clock_t start = clock ();
for (i = 0; i < 1e6; ++i) // Doing a large number of times the same type of computation with different values, to smoothen random fluctuations.
{
a = (double) (1 + 2 * (rand () % 0x400)) / 0x20000000000000; // "a" has only a few significant digits, and its last non-zero digit is at (fixed-point) position 53.
c += expm1 (a) + 1; // replace exp() with expm1() + 1
}
clock_t stop = clock ();
printf ("%e\n", c); // Just to be sure that the compiler will actually perform the computation.
printf ("Clock time spent: %d\n", stop - start);
return 0;
}
Для этого случая время на моей машине таково:
Оригинальный код
1.000000e+06
Время проведенное: 21543338
Модифицированный код
1.000000e+06
Время потрачено: 55076
Программисты с расширенными знаниями о сопутствующих компромиссах могут иногда рассмотреть возможность использования приблизительных результатов, когда точность не критична
Опытный программист может написать примерную реализацию медленной функции, используя такие методы, как полиномы Ньютона-Рафсона, Тейлора или Макларина, особенно неточно округленные специальные функции из библиотек, таких как Intel MKL, AMD AMCL, ослабляя соответствие стандарту с плавающей запятой компилятора, снижая точность до двоичного кода ieee754 (float) или их комбинация.
Обратите внимание, что лучшее описание проблемы позволит получить лучший ответ.
Что касается вашего комментария к ответу @EOF, то "напишите свое собственное" замечание от @NominalAnimal здесь кажется достаточно простым, даже тривиальным, следующим образом.
Ваш исходный код выше, кажется, имеет максимально возможный аргумент для exp () a = (1 + 2 * 0x400) / 0x2000... = 4.55e-13 (это действительно должно быть 2 * 0x3FF, и я считаю 13 нули после вашего 0x2000... что делает его 2x16 ^ 13). Так что максимальный аргумент 4.55e-13 очень и очень мал.
И тогда тривиальное разложение Тейлора является exp(a)=1+a+(a^2)/2+(a^3)/6+..., что уже дает вам всю точность double для таких маленьких аргументов. Теперь вам придется отказаться от 1-й части, как описано выше, а затем она просто сводится к expm1(a)=a*(1.+a*(1.+a/3.)/2.) И что должен идти чертовски быстро! Просто убедитесь, что он остается маленьким. Если он станет немного больше, просто добавьте следующий термин, ^ 4/24 (вы видите, как это сделать?).
>> EDIT <<
Я изменил тестовую программу OP следующим образом, чтобы протестировать немного больше (обсуждение следует за кодом)
/* https://stackru.com/questions/44346371/
i-do-not-want-correct-rounding-for-function-exp/44397261 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define BASE 16 /*denominator will be (multiplier)xBASE^EXPON*/
#define EXPON 13
#define taylorm1(a) (a*(1.+a*(1.+a/3.)/2.)) /*expm1() approx for small args*/
int main (int argc, char *argv[]) {
int N = (argc>1?atoi(argv[1]):1e6),
multiplier = (argc>2?atoi(argv[2]):2),
isexp = (argc>3?atoi(argv[3]):1); /* flags to turn on/off exp() */
int isexpm1 = 1; /* and expm1() for timing tests*/
int i, n=0;
double denom = ((double)multiplier)*pow((double)BASE,(double)EXPON);
double a, c=0.0, cm1=0.0, tm1=0.0;
clock_t start = clock();
n=0; c=cm1=tm1=0.0;
/* --- to smooth random fluctuations, do the same type of computation
a large number of (N) times with different values --- */
for (i=0; i<N; i++) {
n++;
a = (double)(1 + 2*(rand()%0x400)) / denom; /* "a" has only a few
significant digits, and its last non-zero
digit is at (fixed-point) position 53. */
if ( isexp ) c += exp(a); /* turn this off to time expm1() alone */
if ( isexpm1 ) { /* you can turn this off to time exp() alone, */
cm1 += expm1(a); /* but difference is negligible */
tm1 += taylorm1(a); }
} /* --- end-of-for(i) --- */
int nticks = (int)(clock()-start);
printf ("N=%d, denom=%dx%d^%d, Clock time: %d (%.2f secs)\n",
n, multiplier,BASE,EXPON,
nticks, ((double)nticks)/((double)CLOCKS_PER_SEC));
printf ("\t c=%.20e,\n\t c-n=%e, cm1=%e, tm1=%e\n",
c,c-(double)n,cm1,tm1);
return 0;
} /* --- end-of-function main() --- */
Скомпилируйте и запустите его в качестве теста для воспроизведения сценария OP 0x2000... или запустите его с (до трех) необязательным аргументом test #trials multiplier timeexp, где #trials по умолчанию принимает значение 1000000 OP, а множитель по умолчанию равен 2 для 2x16 оператора OP. 13 (замените его на 4 и т. Д. Для других ее тестов). Для последнего аргумента timeexp введите 0, чтобы выполнить только вычисления expm1() (и мой ненужный тейлор-подобный). Смысл этого в том, чтобы показать, что неудачные хронометражи, отображаемые OP, исчезают с помощью expm1(), что занимает "вообще никакого времени" независимо от множителя.
По умолчанию запускается, тестирует и тестирует 1000000 4, производит (ладно, я назвал программу округления)...
bash-4.3$ ./rounding
N=1000000, denom=2x16^13, Clock time: 11155070 (11.16 secs)
c=1.00000000000000023283e+06,
c-n=2.328306e-10, cm1=1.136017e-07, tm1=1.136017e-07
bash-4.3$ ./rounding 1000000 4
N=1000000, denom=4x16^13, Clock time: 200211 (0.20 secs)
c=1.00000000000000011642e+06,
c-n=1.164153e-10, cm1=5.680083e-08, tm1=5.680083e-08
Итак, первое, что вы заметите, это то, что ОП с использованием exp () существенно отличается от cm1==tm1 с использованием expm1() и моего тейлора ок. Если вы уменьшите N, они придут к согласию следующим образом...
N=10, denom=2x16^13, Clock time: 941 (0.00 secs)
c=1.00000000000007140954e+01,
c-n=7.140954e-13, cm1=7.127632e-13, tm1=7.127632e-13
bash-4.3$ ./rounding 100
N=100, denom=2x16^13, Clock time: 5506 (0.01 secs)
c=1.00000000000010103918e+02,
c-n=1.010392e-11, cm1=1.008393e-11, tm1=1.008393e-11
bash-4.3$ ./rounding 1000
N=1000, denom=2x16^13, Clock time: 44196 (0.04 secs)
c=1.00000000000011345946e+03,
c-n=1.134595e-10, cm1=1.140730e-10, tm1=1.140730e-10
bash-4.3$ ./rounding 10000
N=10000, denom=2x16^13, Clock time: 227215 (0.23 secs)
c=1.00000000000002328306e+04,
c-n=2.328306e-10, cm1=1.131288e-09, tm1=1.131288e-09
bash-4.3$ ./rounding 100000
N=100000, denom=2x16^13, Clock time: 1206348 (1.21 secs)
c=1.00000000000000232831e+05,
c-n=2.328306e-10, cm1=1.133611e-08, tm1=1.133611e-08
А что касается времени exp () и expm1(), то убедитесь сами...
bash-4.3$ ./rounding 1000000 2
N=1000000, denom=2x16^13, Clock time: 11168388 (11.17 secs)
c=1.00000000000000023283e+06,
c-n=2.328306e-10, cm1=1.136017e-07, tm1=1.136017e-07
bash-4.3$ ./rounding 1000000 2 0
N=1000000, denom=2x16^13, Clock time: 24064 (0.02 secs)
c=0.00000000000000000000e+00,
c-n=-1.000000e+06, cm1=1.136017e-07, tm1=1.136017e-07
Вопрос: вы заметите, что как только вычисление exp () достигает N=10000 испытаний, его сумма остается постоянной независимо от большего N. Не уверен, почему это произойдет.
>> __ ВТОРОЕ РЕДАКТИРОВАНИЕ __<<
Хорошо, @EOF, "ты заставил меня выглядеть" с твоим комментарием "иерархического накопления". И это действительно работает, чтобы приблизить сумму exp () (намного ближе) к (предположительно правильной) сумме expm1(). Модифицированный код сразу же после обсуждения. Но здесь есть одно замечание: вспомните множитель сверху. Это ушло, и на том же месте находится expon, так что знаменатель теперь равен 2^expon, где значение по умолчанию равно 53, что соответствует значению OP по умолчанию (и я считаю, что лучше соответствует тому, как она думала об этом). Хорошо, а вот код...
/* https://stackru.com/questions/44346371/
i-do-not-want-correct-rounding-for-function-exp/44397261 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define BASE 2 /*denominator=2^EXPON, 2^53=2x16^13 default */
#define EXPON 53
#define taylorm1(a) (a*(1.+a*(1.+a/3.)/2.)) /*expm1() approx for small args*/
int main (int argc, char *argv[]) {
int N = (argc>1?atoi(argv[1]):1e6),
expon = (argc>2?atoi(argv[2]):EXPON),
isexp = (argc>3?atoi(argv[3]):1), /* flags to turn on/off exp() */
ncparts = (argc>4?atoi(argv[4]):1), /* #partial sums for c */
binsize = (argc>5?atoi(argv[5]):10);/* #doubles to sum in each bin */
int isexpm1 = 1; /* and expm1() for timing tests*/
int i, n=0;
double denom = pow((double)BASE,(double)expon);
double a, c=0.0, cm1=0.0, tm1=0.0;
double csums[10], cbins[10][65537]; /* c partial sums and heirarchy */
int nbins[10], ibin=0; /* start at lowest level */
clock_t start = clock();
n=0; c=cm1=tm1=0.0;
if ( ncparts > 65536 ) ncparts=65536; /* array size check */
if ( ncparts > 1 ) for(i=0;i<ncparts;i++) cbins[0][i]=0.0; /*init bin#0*/
/* --- to smooth random fluctuations, do the same type of computation
a large number of (N) times with different values --- */
for (i=0; i<N; i++) {
n++;
a = (double)(1 + 2*(rand()%0x400)) / denom; /* "a" has only a few
significant digits, and its last non-zero
digit is at (fixed-point) position 53. */
if ( isexp ) { /* turn this off to time expm1() alone */
double expa = exp(a); /* exp(a) */
c += expa; /* just accumulate in a single "bin" */
if ( ncparts > 1 ) cbins[0][n%ncparts] += expa; } /* accum in ncparts */
if ( isexpm1 ) { /* you can turn this off to time exp() alone, */
cm1 += expm1(a); /* but difference is negligible */
tm1 += taylorm1(a); }
} /* --- end-of-for(i) --- */
int nticks = (int)(clock()-start);
if ( ncparts > 1 ) { /* need to sum the partial-sum bins */
nbins[ibin=0] = ncparts; /* lowest-level has everything */
while ( nbins[ibin] > binsize ) { /* need another heirarchy level */
if ( ibin >= 9 ) break; /* no more bins */
ibin++; /* next available heirarchy bin level */
nbins[ibin] = (nbins[ibin-1]+(binsize-1))/binsize; /*#bins this level*/
for(i=0;i<nbins[ibin];i++) cbins[ibin][i]=0.0; /* init bins */
for(i=0;i<nbins[ibin-1];i++) {
cbins[ibin][(i+1)%nbins[ibin]] += cbins[ibin-1][i]; /*accum in nbins*/
csums[ibin-1] += cbins[ibin-1][i]; } /* accumulate in "one bin" */
} /* --- end-of-while(nprevbins>binsize) --- */
for(i=0;i<nbins[ibin];i++) csums[ibin] += cbins[ibin][i]; /*highest level*/
} /* --- end-of-if(ncparts>1) --- */
printf ("N=%d, denom=%d^%d, Clock time: %d (%.2f secs)\n", n, BASE,expon,
nticks, ((double)nticks)/((double)CLOCKS_PER_SEC));
printf ("\t c=%.20e,\n\t c-n=%e, cm1=%e, tm1=%e\n",
c,c-(double)n,cm1,tm1);
if ( ncparts > 1 ) { printf("\t binsize=%d...\n",binsize);
for (i=0;i<=ibin;i++) /* display heirarchy */
printf("\t level#%d: #bins=%5d, c-n=%e\n",
i,nbins[i],csums[i]-(double)n); }
return 0;
} /* --- end-of-function main() --- */
Хорошо, и теперь вы можете заметить два дополнительных аргумента командной строки после старого timeexp. Они являются ncparts для начального количества бинов, в которые будут распределены все #trials. Таким образом, на самом низком уровне иерархии каждый бин должен (по модулю ошибок:) иметь сумму двойных значений #trials/ncparts. Аргументом после этого является binsize, который будет количеством двойников, суммируемых в каждом бине на каждом последующем уровне, пока на последнем уровне не будет меньше (или равно) # бинов в качестве размера бина. Итак, вот пример, разделяющий 1000000 попыток на 50000 корзин, что означает 20 удвоений / корзина на самом низком уровне и 5 удвоений / корзина после этого...
bash-4.3$ ./rounding 1000000 53 1 50000 5
N=1000000, denom=2^53, Clock time: 11129803 (11.13 secs)
c=1.00000000000000465661e+06,
c-n=4.656613e-09, cm1=1.136017e-07, tm1=1.136017e-07
binsize=5...
level#0: #bins=50000, c-n=4.656613e-09
level#1: #bins=10002, c-n=1.734588e-08
level#2: #bins= 2002, c-n=7.974450e-08
level#3: #bins= 402, c-n=1.059379e-07
level#4: #bins= 82, c-n=1.133885e-07
level#5: #bins= 18, c-n=1.136214e-07
level#6: #bins= 5, c-n=1.138542e-07
Обратите внимание, что cn для exp () довольно хорошо сходится к значению expm1(). Но обратите внимание, что это лучше всего на уровне #5, и не сходится равномерно на всех. И заметьте, если вы разбиваете #trials только на 5000 начальных корзин, вы получите такой же хороший результат,
bash-4.3$ ./rounding 1000000 53 1 5000 5
N=1000000, denom=2^53, Clock time: 11165924 (11.17 secs)
c=1.00000000000003527384e+06,
c-n=3.527384e-08, cm1=1.136017e-07, tm1=1.136017e-07
binsize=5...
level#0: #bins= 5000, c-n=3.527384e-08
level#1: #bins= 1002, c-n=1.164153e-07
level#2: #bins= 202, c-n=1.158332e-07
level#3: #bins= 42, c-n=1.136214e-07
level#4: #bins= 10, c-n=1.137378e-07
level#5: #bins= 4, c-n=1.136214e-07
На самом деле, игра с ncparts и binsize, похоже, не демонстрирует особой чувствительности, и не всегда "чем больше, тем лучше" (т.е. меньше для binsize). Так что я не уверен, что именно происходит. Может быть ошибка (или два), или может быть еще один вопрос для @EOF...???
>> РЕДАКТИРОВАТЬ - пример, показывающий сложение пары "двоичное дерево" heirarchy<<
Пример ниже добавлен в соответствии с комментарием @EOF (Примечание: повторно скопируйте предыдущий код. Мне пришлось редактировать вычисление nbins[ibin] для каждого следующего уровня в nbins[ibin]=(nbins[ibin-1]+(binsize-1))) / binsize; from nbins[ibin] = (nbins [ibin-1] + 2 * binsize) / binsize; который был "слишком консервативным" для создания ... последовательности 16,8, 16,8)
bash-4.3$ ./rounding 1024 53 1 512 2
N=1024, denom=2^53, Clock time: 36750 (0.04 secs)
c=1.02400000000011573320e+03,
c-n=1.157332e-10, cm1=1.164226e-10, tm1=1.164226e-10
binsize=2...
level#0: #bins= 512, c-n=1.159606e-10
level#1: #bins= 256, c-n=1.166427e-10
level#2: #bins= 128, c-n=1.166427e-10
level#3: #bins= 64, c-n=1.161879e-10
level#4: #bins= 32, c-n=1.166427e-10
level#5: #bins= 16, c-n=1.166427e-10
level#6: #bins= 8, c-n=1.166427e-10
level#7: #bins= 4, c-n=1.166427e-10
level#8: #bins= 2, c-n=1.164153e-10
>> РЕДАКТИРОВАТЬ - показать элегантное решение @ EOF в комментарии ниже <<
"Добавление пары" может быть элегантно выполнено рекурсивно, согласно приведенному ниже комментарию @ EOF, который я воспроизводлю здесь. (Обратите внимание на случай 0/1 в конце рекурсии, чтобы обработать n четное / нечетное.)
/* Quoting from EOF's comment...
What I (EOF) proposed is effectively a binary tree of additions:
a+b+c+d+e+f+g+h as ((a+b)+(c+d))+((e+f)+(g+h)).
Like this: Add adjacent pairs of elements, this produces
a new sequence of n/2 elements.
Recurse until only one element is left.
(Note that this will require n/2 elements of storage,
rather than a fixed number of bins like your implementation) */
double trecu(double *vals, double sum, int n) {
int midn = n/2;
switch (n) {
case 0: break;
case 1: sum += *vals; break;
default: sum = trecu(vals+midn, trecu(vals,sum,midn), n-midn); break; }
return(sum);
}
Это "ответ"/ продолжение предыдущих комментариев EOF относительно его алгоритма trecu () и кода для его предложения "суммирование двоичного дерева". "Предпосылки", прежде чем читать это, читают эту дискуссию. Было бы неплохо собрать все это в одном организованном месте, но я этого еще не сделал...
... Что я сделал, так это встроил функцию trecu () EOF в тестовую программу из предыдущего ответа, который я написал, изменив исходную тестовую программу OP. Но потом я обнаружил, что trecu () сгенерировал точно (и я имею в виду именно тот же) тот же ответ, что и "простая сумма" c с использованием exp (), а не сумма cm1 с использованием expm1(), которую мы ожидали от более точного двоичного дерева суммирование.
Но эта тестовая программа немного (может быть, два бита:) "запутанная" (или, как сказал EOF, "нечитаемый"), поэтому я написал отдельную меньшую тестовую программу, приведенную ниже (с примерами прогонов и обсуждением ниже), по отдельности тест / упражнение trecu (). Более того, я также написал функцию bintreesum () в приведенном ниже коде, который абстрагирует / инкапсулирует итерационный код для суммирования двоичного дерева, который я встроил в предыдущую тестовую программу. В этом предыдущем случае мой итеративный код действительно приблизился к ответу cm1, поэтому я ожидал, что рекурсивный trecu () EOF сделает то же самое. Коротко говоря, это то, что ниже происходит то же самое - bintreesum () остается близким к правильному ответу, а trecu () уходит дальше, точно воспроизводя "простую сумму".
То, что мы суммируем ниже, это просто сумма (i),i=1...n, которая является просто известным n(n+1)/2. Но это не совсем правильно - чтобы воспроизвести проблему OP, summand - это не сумма (i), а сумма (1+i*10^(-e)), где e можно указать в командной строке. Так, скажем, для n=5 вы получите не 15, а 5.000...00015, или для n=6 вы получите 6.000...00021 и т. Д. И чтобы избежать длинного, длинного формата, я печатаю f() sum-n, чтобы удалить эту целочисленную часть. Хорошо??? Так вот код...
/* Quoting from EOF's comment...
What I (EOF) proposed is effectively a binary tree of additions:
a+b+c+d+e+f+g+h as ((a+b)+(c+d))+((e+f)+(g+h)).
Like this: Add adjacent pairs of elements, this produces
a new sequence of n/2 elements.
Recurse until only one element is left. */
#include <stdio.h>
#include <stdlib.h>
double trecu(double *vals, double sum, int n) {
int midn = n/2;
switch (n) {
case 0: break;
case 1: sum += *vals; break;
default: sum = trecu(vals+midn, trecu(vals,sum,midn), n-midn); break; }
return(sum);
} /* --- end-of-function trecu() --- */
double bintreesum(double *vals, int n, int binsize) {
double binsum = 0.0;
int nbin0 = (n+(binsize-1))/binsize,
nbin1 = (nbin0+(binsize-1))/binsize,
nbins[2] = { nbin0, nbin1 };
double *vbins[2] = {
(double *)malloc(nbin0*sizeof(double)),
(double *)malloc(nbin1*sizeof(double)) },
*vbin0=vbins[0], *vbin1=vbins[1];
int ibin=0, i;
for ( i=0; i<nbin0; i++ ) vbin0[i] = 0.0;
for ( i=0; i<n; i++ ) vbin0[i%nbin0] += vals[i];
while ( nbins[ibin] > 1 ) {
int jbin = 1-ibin; /* other bin, 0<-->1 */
nbins[jbin] = (nbins[ibin]+(binsize-1))/binsize;
for ( i=0; i<nbins[jbin]; i++ ) vbins[jbin][i] = 0.0;
for ( i=0; i<nbins[ibin]; i++ )
vbins[jbin][i%nbins[jbin]] += vbins[ibin][i];
ibin = jbin; /* swap bins for next pass */
} /* --- end-of-while(nbins[ibin]>0) --- */
binsum = vbins[ibin][0];
free((void *)vbins[0]); free((void *)vbins[1]);
return ( binsum );
} /* --- end-of-function bintreesum() --- */
#if defined(TESTTRECU)
#include <math.h>
#define MAXN (2000000)
int main(int argc, char *argv[]) {
int N = (argc>1? atoi(argv[1]) : 1000000 ),
e = (argc>2? atoi(argv[2]) : -10 ),
binsize = (argc>3? atoi(argv[3]) : 2 );
double tens = pow(10.0,(double)e);
double *vals = (double *)malloc(sizeof(double)*MAXN),
sum = 0.0;
double trecu(), bintreesum();
int i;
if ( N > MAXN ) N=MAXN;
for ( i=0; i<N; i++ ) vals[i] = 1.0 + tens*(double)(i+1);
for ( i=0; i<N; i++ ) sum += vals[i];
printf(" N=%d, Sum_i=1^N {1.0 + i*%.1e} - N = %.8e,\n"
"\t plain_sum-N = %.8e,\n"
"\t trecu-N = %.8e,\n"
"\t bintreesum-N = %.8e \n",
N, tens, tens*((double)N)*((double)(N+1))/2.0,
sum-(double)N,
trecu(vals,0.0,N)-(double)N,
bintreesum(vals,N,binsize)-(double)N );
} /* --- end-of-function main() --- */
#endif
Поэтому, если вы сохраните это как trecu.c, то скомпилируете его как cc –DTESTTRECU trecu.c –lm –o trecu, а затем запустите от нуля до трех необязательных аргументов командной строки как trecu #trials e binsize По умолчанию #trials=1000000 (как в программе OP), e=–10 и binsize=2 (для моей функции bintreesum () нужно делать сумму двоичного дерева, а не бины большего размера).
И вот некоторые результаты испытаний, иллюстрирующие проблему, описанную выше,
bash-4.3$ ./trecu
N=1000000, Sum_i=1^N {1.0 + i*1.0e-10} - N = 5.00000500e+01,
plain_sum-N = 5.00000500e+01,
trecu-N = 5.00000500e+01,
bintreesum-N = 5.00000500e+01
bash-4.3$ ./trecu 1000000 -15
N=1000000, Sum_i=1^N {1.0 + i*1.0e-15} - N = 5.00000500e-04,
plain_sum-N = 5.01087168e-04,
trecu-N = 5.01087168e-04,
bintreesum-N = 5.00000548e-04
bash-4.3$
bash-4.3$ ./trecu 1000000 -16
N=1000000, Sum_i=1^N {1.0 + i*1.0e-16} - N = 5.00000500e-05,
plain_sum-N = 6.67552231e-05,
trecu-N = 6.67552231e-05,
bintreesum-N = 5.00001479e-05
bash-4.3$
bash-4.3$ ./trecu 1000000 -17
N=1000000, Sum_i=1^N {1.0 + i*1.0e-17} - N = 5.00000500e-06,
plain_sum-N = 0.00000000e+00,
trecu-N = 0.00000000e+00,
bintreesum-N = 4.99992166e-06
Итак, вы можете видеть, что для запуска по умолчанию, e=–10, все делают все правильно. То есть, верхняя строка с надписью "Сумма" просто выполняет n (n + 1) / 2, поэтому, по-видимому, отображает правильный ответ. И все нижеуказанные соглашаются на стандартный тест e = –10. Но для случаев e=–15 и e=–16 ниже, trecu () точно совпадает с plain_sum, в то время как bintreesum остается довольно близко к правильному ответу. И, наконец, для e=–17, plain_sum и trecu () "исчезли", в то время как bintreesum () все еще там висит.
Таким образом, trecu () правильно выполняет суммирование, но его рекурсия, по-видимому, не выполняет тот тип "двоичного дерева", который мой более простой итеративный bintreesum (), по-видимому, делает правильно. И это действительно демонстрирует, что предложение EOF о "суммировании двоичного дерева" реализует значительное улучшение по сравнению с plain_sum для этих случаев типа 1+epsilon. Так что нам бы очень хотелось увидеть его рекурсивную работу trecu ()!!! Когда я первоначально смотрел на это, я думал, что это работало. Но эта двойная рекурсия (есть ли специальное название для этого?) В его значении по умолчанию: case, по-видимому, более запутан (по крайней мере, для меня:), чем я думал. Как я уже сказал, он делает сумму, а не "двоичное дерево".
Итак, кто хотел бы принять вызов и объяснить, что происходит в этой рекурсии trecu ()? И, может быть, еще важнее, исправить это так, чтобы он делал то, что задумано. Благодарю.