Я делаю * не * хочу правильное округление для функции exp

Ответ 1

Чтобы ответить на общий вопрос о том, почему библиотечные функции необходимы для корректного округления результатов:

Плавная точка тяжелая, и часто время не поддается интуиции. Не каждый программист прочитал то, что у них должно было быть. Когда библиотеки использовали несколько слегка неточное округление, люди жаловались на точность библиотечной функции, когда их неточные вычисления неизбежно ошибались и порождали бессмысленность. В ответ авторы библиотек сделали свои библиотеки точно округленными, поэтому теперь люди не могут переложить вину на них.

Во многих случаях конкретные знания о алгоритмах с плавающей запятой могут значительно улучшить точность и/или производительность, например, в тестовом сценарии:

Взятие exp() чисел, очень близких к 0 в числах с плавающей запятой, проблематично, так как результат - это число, близкое к 1, а вся точность находится в разнице с одним, поэтому наиболее значимые цифры потеряны. Точнее (и значительно быстрее в этом тестовом сценарии) вычислять exp(x) - 1 через библиотечную функцию C math 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 binary32 (float) или их комбинация.

Обратите внимание, что лучшее описание проблемы обеспечит лучший ответ.

Ответ 2

Относительно вашего комментария к ответе @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 +..., которое уже дает вам все двойные точность для таких небольших аргументов. Теперь вам придется отбросить часть 1, как описано выше, а затем просто сводится к expm1 (a) = a * (1. + a * (1. + a/3.)/2.) И это должно пройти довольно быстро! Просто убедитесь, что a остается маленьким. Если он становится немного больше, просто добавьте следующий термин a ^ 4/24 (вы видите, как это сделать?).

<сильное → > РЕДАКТИРОВАТЬ < <

Я изменил программу тестирования OP следующим образом, чтобы проверить немного больше материала (обсуждение следует за кодом)

/* https://stackoverflow.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() --- */

Скомпилируйте и запустите его как test для воспроизведения сценария OP 0x2000... или запустите его с (до трех) необязательных аргументов test   # tests     timeexp, где #trials по умолчанию используется OP 1000000 и multipler по умолчанию 2 для OP 2x16 ^ 13 (измените его на 4 и т.д. Для других ее тестов). Для последнего аргумента timeexp введите 0, чтобы выполнить только expm1() (и мой ненужный taylor-like) расчет. Дело в том, чтобы показать, что неудачные временные случаи, отображаемые OP, исчезают с expm1(), который принимает "нет времени вообще" независимо от множителя.

Итак, по умолчанию выполняется test и test 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

Итак, первое, что вы заметите, это то, что OP cn с использованием exp() существенно отличается от cm1 == tm1, используя expm1(), а мой taylor ок. Если вы уменьшите 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. Не знаете, почему это произойдет.

→ __ SECOND EDIT __ <

Хорошо, @EOF, "вы заставили меня посмотреть" с вашим комментарием "heirarchical accumulation". И это действительно помогает приблизить сумму exp() ближе (гораздо ближе) к (предположительно правильной) сумме expm1(). Модифицированный код сразу же после обсуждения. Но одно обсуждение здесь: вспомните множитель сверху. Это прошло, и там же expon, так что знаменатель теперь 2 ^ expon, где значение по умолчанию 53, сопоставление OP по умолчанию (и Я считаю, что лучше совместить то, как она думала об этом). Хорошо, и вот код...

/* https://stackoverflow.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. Таким образом, на самом низком уровне иерархии каждый бин должен (по модулю ошибок) иметь сумму # tests/ncparts. После этого аргумент binsize, который будет состоять из числа удвоений, суммированных в каждом бункере на каждом последующем уровне, до тех пор, пока на последнем уровне меньше (или равно) #bins binsize > . Итак, вот пример, делящий 1000000 проб в 50000 ящиков, что означает 20doubles/bin на самом низком уровне и 5doubles/bin после этого...

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

Обратите внимание, что c-n для 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...???

→ EDIT - пример, показывающий добавление пары "бинарное дерево", иерархия <

Пример ниже добавлен в соответствии с комментарием @EOF (Примечание: повторите копию предыдущего кода. Мне пришлось отредактировать вычисления nbins [ibin] для каждого следующего уровня до nbins [ibin] = (nbins [ibin-1] + (binsize-1))/binsize; от nbins [ibin] = (nbins [ibin-1] + 2 * binsize)/binsize;, который был "слишком консервативен" для создания ... 16,8,4, 2)

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

→ EDIT - показать элегантное решение @EOF в комментарии ниже <

"Добавление пара" может быть изящно выполнено рекурсивно, согласно приведенному ниже комментарию @EOF, который я воспроизвожу здесь. (Обратите внимание на случай 0/1 в конце рекурсии для обработки n четный/нечетный.)

  /* Quoting from EOF 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);
      } 

Ответ 3

Это "ответ" /продолжение в EOF, предшествующий комментариям его алгоритм trecu() и код для его предложения "суммирование двоичного дерева". "Предварительные требования" перед чтением этого читают эту дискуссию. Было бы неплохо собрать все это в одном организованном месте, но я еще этого не сделал...

... То, что я сделал, это собрать EOF trecu() в тестовую программу из предыдущего ответа, который я написал, изменив исходную тестовую программу OP. Но затем я обнаружил, что trecu() генерирует точно (и я имею в виду точно) тот же ответ, что и "простая сумма" c, используя exp(), а не сумма cm1, используя expm1(), которую мы ожидали от более точного суммирования двоичного дерева.

Но эта тестовая программа немного (может быть, две биты:) "свернуты" (или, как сказал EOF, "нечитабельно" ), поэтому я написал отдельную небольшую тестовую программу, приведенную ниже (с примерами прогонов и обсуждением ниже), отдельно тестировать/осуществлять trecu(). Более того, я также написал функцию bintreesum() в приведенный ниже код, который абстрагирует/инкапсулирует итеративный код для суммирования двоичного дерева, который я вложил в предыдущую тестовую программу. В этом предыдущем случае мой итеративный код действительно приблизился к ответу cm1, поэтому я ожидал, что EOF рекурсивный trecu() сделает то же самое. Долгий и короткий - это то, что ниже, то же самое происходит - bintreesum() остается близким к правильному ответу, в то время как trecu() уходит дальше, точно воспроизводя "равную сумму".

То, что мы суммируем ниже, это просто сумма (i), я = 1... n, что является просто известным n (n + 1)/2. Но это не совсем верно - чтобы воспроизвести проблему ОП, слагаемое не является суммой (i) в одиночку, а скорее суммой (1 + я * 10 ^ (- e)), где e можно указать в командной строке. Итак, для, скажем, n = 5, вы не получаете 15, а скорее 5.000... 00015, или для n = 6 вы получаете 6.000... 00021 и т.д. И чтобы избежать длинного длинного формата, я печатаю f ( ) sum-n, чтобы удалить эту целую часть. Хорошо??? Итак, вот код...

/* Quoting from EOF 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 И затем запустите с нулевым значением до трех необязательных команд- line args как trecu # tests e binsize По умолчанию: # tests = 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, все делают все правильно. То есть, верхняя строка, которая говорит "Sum", просто делает 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()!!! Когда я впервые посмотрел на него, я подумал, что это сработало. Но эта двойная рекурсия (для этого есть специальное имя?) В его случае по умолчанию:, по-видимому, более запутанным (по крайней мере для меня:), чем я думал. Как я уже сказал, он выполняет сумму, но не "двоичное дерево".

Хорошо, так кто хотел бы взять вызов и объяснить, что происходит в этой рекурсии trecu()? И, возможно, что более важно, исправить это, чтобы он сделал то, что предназначалось. Спасибо.