Рекуррентное соотношение Мюллера: проблемы с округлением чисел с плавающей точкой

Original author: Lincoln Atkinson
  • Translation
Некоторое время назад я натолкнулся на упражнение, которое выглядит не так уж и сложно:

Пусть последовательность xn определена так:

посчитайте x30.

Это не так уж и трудно закодировать, возможно реализовав xi как рекурсивную функцию. С обычными числами с плавающей запятой двойной точности, по мере увеличения i, результат красиво сходится к 100. Супер!

К сожалению, 100 даже близко не является правильным ответом. На самом деле последовательность сходится к 5.

Проблема


Эта последовательность известна под именем «Рекуррентное соотношение Мюллера» и синтезирована специально для демонстрации того, как быстро и драматически ошибки округления чисел с плавающей запятой приводят к катастрофическим результатам, при правильных (ну, неправильных) условиях. В этой работе (pdf) в деталях рассматриваются различные опасности округления, в частности и эта последовательность (стр. 14). Другое описание проблемы можно увидеть здесь (pdf).

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

Я больше ничего не могу добавить, но думаю, что этой проблемой стоит поделиться.

Вычисление корректного результата


Мы можем сравнить первые сто значений последовательности между обычными числами с плавающей запятой и арифметикой с произвольной точностью с помощью короткой программы на Mathematica. Мы используем удобный синтакс мемоизации Mathematica, так что вычисления происходят довольно быстро даже для 100 и более итераций.

f[y_z_:= 108 - (815 - 1500/z)/y;
 
xExact[0= 4;
xExact[1= 17/4;
xExact[n_:= xExact[n] = f[xExact[n-1], xExact[n-2]];
 
xFloat[0= 4;
xFloat[1= 4.25;
xFloat[n_:= xFloat[n] = f[xFloat[n-1], xFloat[n-2]];
 
TableForm[
 Table[{i, N[xExact[i], 20], N[xFloat[i], 20]}, {i, 0100}],
 TableHeadings ->
    {None, {"i""x[i] \"exact\"""x[i] floating point"}}
Результаты:
i     x[i] "exact"             x[i] floating point
----------------------------------------------------
0     4.0000000000000000000    4.0000000000000000000
1     4.2500000000000000000    4.25
2     4.4705882352941176471    4.47059
3     4.6447368421052631579    4.64474
4     4.7705382436260623229    4.77054
5     4.8557007125890736342    4.8557
6     4.9108474990827932004    4.91085
7     4.9455374041239167248    4.94554
8     4.9669625817627005987    4.96696
9     4.9800457013556311613    4.98004
10    4.9879794484783922601    4.98791
11    4.9927702880620680975    4.99136
12    4.9956558915066340266    4.96746
13    4.9973912683813441129    4.42971
14    4.9984339439448169190    -7.81676
15    4.9990600719708938678    168.943
16    4.9994359371468391480    102.04
17    4.9996615241037675378    100.1
18    4.9997969007134179127    100.005
19    4.9998781354779312492    100.
20    4.9999268795045999045    100.
21    4.9999561270611577381    100.
22    4.9999736760057124446    100.
23    4.9999842055202727079    100.
24    4.9999905232822276594    100.
25    4.9999943139585595936    100.
26    4.9999965883712560237    100.
27    4.9999979530213569080    100.
28    4.9999987718123113300    100.
29    4.9999992630872057846    100.
30    4.9999995578522583059    100.
31    4.9999997347113315242    100.
32    4.9999998408267904691    100.
33    4.9999999044960712411    100.
34    4.9999999426976416502    100.
35    4.9999999656185845961    100.
36    4.9999999793711506158    100.
37    4.9999999876226903184    100.
38    4.9999999925736141727    100.
39    4.9999999955441684970    100.
40    4.9999999973265010958    100.
41    4.9999999983959006566    100.
42    4.9999999990375403937    100.
43    4.9999999994225242361    100.
44    4.9999999996535145416    100.
45    4.9999999997921087250    100.
46    4.9999999998752652350    100.
47    4.9999999999251591410    100.
48    4.9999999999550954846    100.
49    4.9999999999730572908    100.
50    4.9999999999838343745    100.
51    4.9999999999903006247    100.
52    4.9999999999941803748    100.
53    4.9999999999965082249    100.
54    4.9999999999979049349    100.
55    4.9999999999987429610    100.
56    4.9999999999992457766    100.
57    4.9999999999995474659    100.
58    4.9999999999997284796    100.
59    4.9999999999998370877    100.
60    4.9999999999999022526    100.
61    4.9999999999999413516    100.
62    4.9999999999999648110    100.
63    4.9999999999999788866    100.
64    4.9999999999999873319    100.
65    4.9999999999999923992    100.
66    4.9999999999999954395    100.
67    4.9999999999999972637    100.
68    4.9999999999999983582    100.
69    4.9999999999999990149    100.
70    4.9999999999999994090    100.
71    4.9999999999999996454    100.
72    4.9999999999999997872    100.
73    4.9999999999999998723    100.
74    4.9999999999999999234    100.
75    4.9999999999999999540    100.
76    4.9999999999999999724    100.
77    4.9999999999999999835    100.
78    4.9999999999999999901    100.
79    4.9999999999999999940    100.
80    4.9999999999999999964    100.
81    4.9999999999999999979    100.
82    4.9999999999999999987    100.
83    4.9999999999999999992    100.
84    4.9999999999999999995    100.
85    4.9999999999999999997    100.
86    4.9999999999999999998    100.
87    4.9999999999999999999    100.
88    4.9999999999999999999    100.
89    5.0000000000000000000    100.
90    5.0000000000000000000    100.
91    5.0000000000000000000    100.
92    5.0000000000000000000    100.
93    5.0000000000000000000    100.
94    5.0000000000000000000    100.
95    5.0000000000000000000    100.
96    5.0000000000000000000    100.
97    5.0000000000000000000    100.
98    5.0000000000000000000    100.
99    5.0000000000000000000    100.
100   5.0000000000000000000    100.

Share post

Similar posts

AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 116

    +3
    Как-то очень неожиданное поведение. особенно когда оно происходит не плавно, а вот так скачком.
      +22
      Там нет скачка. У указанного рекуррентного соотношения есть три неподвижных точки: 3, 5 и 100.
      Из них первая и третья — устойчивы, а вторая — неустойчива. Это как в механике: вершина холма — точка неустойчивого равновесия, низ лунки — устойчивого.

      Начальные значения подобраны специальным образом, чтобы при точных вычислениях сходиться к 5. Однако малейшее отклонение (неизбежное при вычислениях с плавающей точкой) приведет к раскачке последовательности к 3 или 100.

      Например, как только результат вычислений выскочит за 5 (даже на одну гугольную), процесс пойдет в раскачку: если z и y больше 5, то 815-1500/z будет меньше 515 и после деления на y (который чуть-чуть больше 5) будет меньше 103. Эрго, 108-(815-1500/z)/y будет уже больше 5 уже почти на десять гугольных. На следующей итерации отклонение в положительную сторону от 5 вырастет еще в 10 раз и так далее, пока не достигнет следующей точки сгущения 100.

      Для любителей матана приведу грубое обоснование пассажа про десять гугольных. Пусть d — бесконечно малая величина и y = z = 5 + d + O(d^2). Тогда f(y, z) = 5 + 43/5 d + O(d^2). Вот этот множитель 43/5 (округленный до 10 в предыдущем абзаце) и раскачивает лодку. В то же время 3 и 100 — устойчивы, поскольку, например, при y = z = 100 + d + O(d^2) имеем f(y, z) = 100 + 157/2000 d + O(d^2), т. е. отклонение уменьшилось в 2000/157 раз.

      Существование специальных начальных значений, которые все же ведут к положению неустойчивого равновесия, объясняется тем фактом, что для k = sqrt(283/5)-6 = 1.52 и y = 5 + d + O(d^2), z = 5 + k*d + O(d^2) мы получим f(y, z) = 5 + d/k + O(d^2) — т. е. отклонение от 5 уменьшится в полтора раза и это явление будет рекурсивно повторяться.
        –2
        Как-то мне разные нехорошие мысли приходят в голову об оригинальном авторе.
        Очевидно что у полинома три неподвижных точки: 3, 5 и 100 — средняя школа, класс примерно восьмой (когда там учат полином на полином делить, не помю точно). При этом точка 5 неустойчива (опоздал с комментарием), а 3 и 100 вполне себе устойчивы. Если задать y=z=3 то итерации вполне себе устойчивы и никакие ошибки округления на результат не влияют.
        Если бы мы искали нули полинома последовательными приближениями то вполне бы получили и 3 и 100, однако здесь есть дополнительно условие на рекурсивность — на каждой следующей итерации x=f(y,z), z->y, y->x и надо выполнить условие чтобы не только начальные y и z попадали в область аттрактора 3, но из значения на каждой следующей итерации оставались в зоне аттрактора. Разлагая по y и z вблизи 3 (матан, первый курс, тут я уверен) получаем что х~=3+45*(y-3)-500/9*(y-3), то есть для схождения к 3 начальные значения должны быть примерно y=z/.63.
        В то же время, чтобы это соотношение выполнялось на следующем шаге (z->y, y->x), надо чтобы выполнялось примерно y~=.334*z, то есть условия схождения к 3 никогда не выполняются. Другими словами, в двумерном пространстве (y,z) мы так выбираем путь приближения к аттрактору 3, что он обязательно выйдет из зоны притяжения и уйдет к аттрактору 100.
        Ну и при чем здесь ошибки округления? А тем более разные языки программирования и матлабы?
        100 даже близко не является правильным ответом. На самом деле последовательность сходится к 5
        — что значит «на самом деле»? Последовательность сходится именно к 100, так вы ее задали, а к 5 она сходиться не может в принципе, это точка отталкивания а не притяжения.
          +6
          При заданных начальных условиях xn = (3n+1+5n+1)/(3n+5n)
          Очевидно, что предел равен 5, при вычислениях в рациональных числах для исходной рекуррентной формулы мы также получаем последовательные приближения к 5.
          Т. е. ожидаемый и правильный ответ должен быть равен 5.
          Другое дело, что действительно, точка 5 неустойчива и при вычислениях с привнесённой ошибкой мы сразу от неё уйдём.
          Так в этом и есть смысл статьи: вычисления с плавающей точкой сложны и могут приводить к неожиданным результатам, даже когда числа не являются малыми/большими, надо внимательно за ними следить.
          +1
          А как получилось, что точка 3 устойчива? Если взять z=3+a, y=3+b то получится f(y,z)=3+35*b-500/9*a+o(a)+o(b). Собственные значения у матрицы перехода от (a,b) к (b,f(y,z)-3) получаются 100/3 и 5/3, оба больше единицы.
            0
            Вы как всегда правы, я недоглядел.
              0
              А что посоветуете почитать, чтобы хотя бы понять, причём тут матрица перехода и её собственные значения?
                0
                Это на стыке матана и линейной алгебры. Надо смотреть такие вещи, как якобиан отображения, линейный оператор, норма матрицы… Конечно, есть и более специализированные курсы (какая-нибудь теория устойчивости), но для общего понимания ситуации можно обойтись без них.
                  0
                  Зря я линал недостаточно усердно ботал на первом курсе, эх.

                  Спасибо.
                    +1
                    Это не первый курс. Т.е. «матан» и «линал» действительно на первом, но вот элементы теории устойчивости у нас были в курсе дифференциальных уравнений, кажется, на втором. Начинать диффуры не освоив анализ — весьма сомнительная идея.
          0
          В некоторых языках есть встроенная рациональная арифметика. Очень удобно, когда делить нужно.

          def f(y,z)
            108 - (815 - 1500 / z) / y
          end
          
          @x = { 0 => Rational(4), 1 => Rational(17,4) }
          
          def x(n)
            @x[n] ||= f(x(n-1), x(n-2))
          end
          
          (2..30).each do |r|
            puts x(r).to_f
          end
          
          • UFO just landed and posted this here
              0
              При двух рекурсивных вызовах в середине функции? Не работает, конечно.
              • UFO just landed and posted this here
                  0
                  На самом деле, непонятно, почему конечно. Мне кажется, что любую рекурсивную функцию с несколькими вызовами рекурсии можно переписать в рекурсивную функцию с одним вызовом рекурсии, а затем развернуть ее в цикл (с аккумулятором). Чтобы это сделать, нужно преобразовать исходную функцию, чтобы она возвращала не только новое значение, а и N-1 предыдущих, где N — количество рекурсивных вызовов. Преобразование примерно такое (на примере функции вычисления последовательности Фибоначчи):
                  -- Пример на Lua
                  -- Исходная функция
                  function fib(x)
                    if x < 2 then return -x end;
                    return fib(x-1) + fib(x-2);
                  end
                  -- Итерация 1: явно выделяем аккумуляторную функцию (+(x,y))
                  function fib1(x)
                    if x < 2 then return -x end;
                  
                    local x1 = fib1(x-1);
                    local x2 = fib1(x-2);
                    return x1 + x2;
                  end
                  -- Итерация 2: возвращаем N-1 предыдущих результатов из функции
                  function fib2(x)
                    if x < 2 then return -x, nil end;
                  
                    local x1 = fib2(x-1);
                    local x2 = fib2(x-2);
                    return x1 + x2, x1;
                  end
                  -- Итерация 3: избавляемся от всех рекурсивных вызовов, кроме одного
                  function fib3(x)
                    local function f(x, stop)
                      if stop then return -x, nil end;
                      if x < 2 then return -x, f(x-1, true) end;
                  
                      local x1, x2 = f(x-1, stop);
                      return x1 + x2, x1;
                    end
                    return f(x, false);
                  end
                  

                  Если записывать в общем виде, то как-то так:
                  -- Общий вид рекурсивной функции с несколькими рекурсивными вызовами
                  -- @X Вектор аргументов функции
                  -- @check условие выхода из рекурсии, функция от вектора аргументов
                  -- @default действия, выполняемые при выходе из рекурсии, функция от вектора аргументов
                  -- @pre действия, выполняемые над аргументами при продолжении рекурсии, до вызова рекурсивной функции, функция от вектора аргументов
                  -- @post действия, выполняемые над результатом рекурсивного вызова при продолжении рекурсии, функция от вектора результатов
                  function genF(check, default, pre, post)
                    local function f(X)
                      if check(X) then return default(X) end;
                  
                      return post(
                        f(pre[1](X)),
                        f(pre[2](X)),
                        f(pre[3](X)),
                        ...
                      );
                    end
                    return f;
                  end
                  -- Преобразованная функция
                  function genG(check, default, pre, post)
                    local function f(X, stop)
                      if stop then return default(X), nil end;
                      if check(X) then return default(X), f(pre(X), true) end;
                  
                      return post(f(pre(X), stop)), x1;
                    end
                    return funtion(X) return f(X, false); end
                  end
                  
                    0
                    N-1 предыдущих
                    Годится только для последовательностей вида Ak = f(Ak-1, ..., Ak-N), где f – какая-то «комбинирующая функция». Но последовательности могут быть и такие, что k-ый член определяется, например, с помощью k/2, k/3 и т.д.

                    Если взять совсем любую рекурсивную функцию вроде mergesort'а, то, кажется мне, никаким конечным N тут не обойтись.
                      0
                      Как раз mergesort без рекурсии пишется легко, потому что там сначала идут рекурсивные вызовы для половинок массива, а потом результаты сливаются. Проще всего — слить сначала пары соседних элементов, потом четвёрки, потом восьмёрки… С помощью небольшой магии можно добиться даже такого же порядка выполнения, как в рекурсивном mergesort.
                        0
                        А нет, причина там в том, что во-первых, дерево рекурсии известно сразу и полностью, а во-вторых, хранение результатов не требует памяти — они закодированы в сортируемом массиве. Поэтому можно выполнять вычисления прямо из глубины, с любого места и в любом порядке. Главное, знать, какие узлы дерева уже вычислены, а какие — нет.
                        0
                        Идея все равно остается той же — возвращаем из функции не скаляр очередного результата, а вектор, содержащий новый результат и все, от чего он зависит. На вход функции передается вектор, в котором все результаты рекурсивных вызовов, кроме одного, уже посчитаны, нужно только ими воспользоваться и передать дальше. Как именно функция зависит от предыдущих вызовов, неважно.
                          0
                          На самом деле, подход, конечно годится не для любых рекурсивных функций, а только для тех, которые представимы в виде функций такого вида:
                          function f(X)
                            if check(X) then return default(X) end;
                          
                            return post(
                              X,
                              f(pre[1](X)),
                              f(pre[2](X)),
                              f(pre[3](X)),
                              ...
                            );
                          end
                          

                          Так что, например, функция, типа
                          function g(x)
                            if x % 2 == 0 return x / 2 end
                            return g(g(x-1))
                          end
                          

                          именно данной процедурой не возьмется. Впрочем, мне кажется, что если еще немного обобщить, то формулу можно будет подогнать и под этот случай — каждый последующий pre может зависеть от любого предыдущего. Вроде не должно это особо усложнить алгоритм.
                            +1
                            На самом деле, существует общий алгоритм преобразования, работающий для любых функций. Сначала из функции делается конечный автомат — а затем вводится стек для хранения состояний. Все.

                            Проблема в том, что ни ваш алгоритм, ни мой не экономят память — а просто загоняют проблему в другое место. Это, конечно, полезно для безстековых платформ — но слабо применимо в обычных случаях.
                              0
                              Как мне кажется, проблема как раз и состоит в том, чтобы строить этот конечный автомат — иначе бы поголовно все компиляторы гарантировали бы разворот рекурсии. Между тем очень обидно, что даже простейшая рекурсия, не требующая введения аккумуляторов, практически никаким компилятором не гарантируется, что будет оптимизирована. А ведь в свете смещения мейнстрима в сторону функциональщины, это должна быть приоритетная задача!
                                +1
                                Нет никакой проблемы в построении конечного автомата на основе программы. Компилятор того же C# умеет делать это автоматически аж для двух классов методов (генераторы и асинхронные методы). А в новом C++/CLI обещают аж автоматические конечные автоматы для использования в любых целях.

                                Проблема в другом. При использование этого метода не появляется никаких преимуществ перед обычной рекурсией.
                                  +1
                                  Ещё есть вариант преобразования с помощью CPS:
                                  вот обсуждение на StackOverflow.
                                  Вкратце, идея в том, чтобы действия, которые надо выполнить с полученным результатом, передать в первый рекурсивный вызов. Тут тоже никаких преимуществ перед обычной рекурсией, разве что память, которая потребляется, обычно меняется со стека на кучу.
                                    0
                                    Да, или так. Но это, на самом деле, тоже очень похоже на конечный автомат — меняется только способ перехода между состояниями, само же разделение исходного алгоритма на блоки остается в том же виде.
                                    +1
                                    При использование этого метода не появляется никаких преимуществ перед обычной рекурсией.

                                    Трудно согласиться. Всё-таки, штатный стек программы ограничен. А преобразование рекурсивной функции в автомат со стековой памятью позволит выдержать эквивалент рекурсивного вызова глубиной в сотни миллионов, если такое вдруг понадобится.
                                    Интересно, можно ли для организации такого автомата воспользоваться уже реализованной в компиляторе C# упаковкой контекста, которая используется в реализации лямбда-выражений и в конструкции yield return.
                                      0
                                      Да, можно — даже двумя способами (через IEnumerable и через Task). Но в каждом из них придется реализовывать рекурсивный вызов особым способом — и внимательно следить, чтобы он не свелся после компиляции к настоящему рекурсивному вызову.
                                        0
                                        Мда.
                                        Было:
                                        static int RecFunc(int k){
                                          int s=1;
                                          for(int i=0;i<k;i++){
                                            s+=RecFunc(i);
                                            s+=RecFunc(i/2);
                                          }
                                          return s;
                                        }
                                        

                                        Стало:
                                                static int ExecFunc(int k) {
                                                    Stack<IEnumerator<int[]>> stack=new Stack<IEnumerator<int[]>>();
                                                    stack.Push(RecFunc(k).GetEnumerator());
                                                    for(;;) {
                                                        stack.Peek().MoveNext();
                                                        int[] arr=stack.Peek().Current;
                                                        if(arr[0]==0) {  // it was return
                                                            int lastres=arr[2];
                                                            stack.Pop().MoveNext();  // end iterator
                                                            if(stack.Count==0) return lastres;
                                                            stack.Peek().Current[2]=lastres;
                                                        } else {  // call
                                                            stack.Push(RecFunc(arr[1]).GetEnumerator()); // recursive call
                                                        }
                                                    }
                                                }
                                        
                                                static IEnumerable<int[]> RecFunc(int k) {
                                                    int s=1;
                                                    int[] inf1=new int[3]{1,0,0}; // code, argument, placeholder for result
                                                    for(int i=0;i<k;i++) {
                                                        inf1[1]=i;   // RecFunc(i)
                                                        yield return inf1;
                                                        s+=inf1[2];
                                                        inf1[1]=i/2;  // RecFunc(i/2)
                                                        yield return inf1;
                                                        s+=inf1[2];
                                                    }
                                                    inf1[2]=s;
                                                    inf1[0]=0;
                                                    yield return inf1;
                                                }
                                        


                                          0
                                          Но, в принципе, так можно реализовать перекрёстные вызовы любого числа функций…

                                              class Program{
                                                  public static void Main(string[] args){
                                                      RRFunc.Execute(RecFunc,4);
                                                  }
                                          
                                                  public static IEnumerable<RFunc> RecFunc(object[] arg) {
                                                      int s=1;
                                                      int k=(int)arg[0];
                                                      for(int i=0;i<k;i++) {
                                                          arg[0]=i;   // RecFunc(i)
                                                          yield return RecFunc;
                                                          s+=(int)arg[0];
                                                          arg[0]=i/2;   // RecFunc(i/2)
                                                          yield return RecFunc;
                                                          s+=(int)arg[0];
                                                      }
                                                      Console.WriteLine("func({0})={1}",k,s);
                                                      arg[0]=s;
                                                      yield return null;
                                                  }
                                              }
                                          
                                          
                                              public delegate IEnumerable<RFunc> RFunc(object[] argres);
                                              class RRFunc {
                                                  public static object Execute(RFunc func,object obj) {
                                                      object[] arg=new object[1] { obj };
                                                      Stack<IEnumerator<RFunc>> stack=new Stack<IEnumerator<RFunc>>();
                                                      stack.Push(func(arg).GetEnumerator());
                                                      for(;;) {
                                                          stack.Peek().MoveNext();
                                                          RFunc next=stack.Peek().Current;
                                                          if(next==null) {  // it was return
                                                              stack.Pop().MoveNext();  // end iterator
                                                              if(stack.Count==0) return arg[0];
                                                          } else {  // call
                                                              stack.Push(next(arg).GetEnumerator()); // recursive call
                                                          }
                                                      }
                                                  }
                                              }
                                          
                                            0
                                            И даже так:
                                                public delegate IEnumerable<RFunc> RFunc(object[] argres);
                                                class RRFunc {
                                                    public static object Execute(RFunc func,object obj) {
                                                        object[] arg=new object[1] { obj };
                                                        Stack<IEnumerator<RFunc>> stack=new Stack<IEnumerator<RFunc>>();
                                                        for(;;) {
                                                            if(func!=null) stack.Push(func(arg).GetEnumerator()); // recursive call
                                                            else if(stack.Count==1) return arg[0];
                                                            else stack.Pop();
                                                            stack.Peek().MoveNext();
                                                            func=stack.Peek().Current;
                                                        }
                                                    }
                                                    public static IEnumerable<RFunc> Fact(object[] arg) {
                                                        int k=(int)arg[0];
                                                        if(k>2) {
                                                            arg[0]=k-1;
                                                            yield return Fact;
                                                            arg[0]=(int)arg[0]*k;
                                                        }
                                                        yield return null;
                                                    }
                                                }
                                            
                                            0
                                            Вы слишком извратились с кодом. Ну, или недостаточно извратились…
                                            abstract class Operation<T> {
                                              public abstract void Process(Stack<IEnumerator<Operation<T>>> stack);
                                              protected abstract void SetResult(T value);
                                            
                                              public static void DisposeFrame(IEnumerator<Operation<T>> en) {
                                                var d = en as IDisposable;
                                                if (d != null) d.Dispose();
                                              }
                                            }
                                            
                                            class Return<T> : Operation<T> { 
                                              public T Value { get; set; }
                                            
                                              public override void Process(Stack<IEnumerator<Operation<T>>> stack) {
                                                var top = stack.Pop();
                                                stack.Peek().Current.SetResult(Value);
                                                DisposeFrame(top);
                                              }
                                            
                                              protected override void SetResult(T value) {
                                                throw new NotSupportedException();
                                              }
                                            }
                                            
                                            class Call<T> : Operation<T> {
                                              public IEnumerable<Operation<T>> Sequence { get; private set; }
                                              public T Result { get; private set; }
                                            
                                              public void Init(IEnumerable<Operation<T>> seq) {
                                                Sequence  = seq;
                                                Result = default(T);
                                              }
                                            
                                              public override void Process(Stack<IEnumerator<Operation<T>>> stack) {
                                                stack.Push(Sequence.GetEnumerator());
                                              }
                                            
                                              protected override void SetResult(T value) {
                                                Result = value;
                                              }
                                            }
                                            
                                            static T Exec<T> (this IEnumerable<Operation<T>> sequence) {
                                              var stack = new Stack<IEnumerator<Operation<T>>>();
                                              var $call = new Call<T> { Sequence = sequence };
                                              stack.Push(new Operation<T>[]{ $call }).GetEnumerator());
                                            
                                              try {
                                                while (stack.Count > 0) {
                                                  if (stack.Peek().MoveNext()) 
                                                    stack.Peek().Current.Process(stack);
                                                  else if (stack.Count > 1)
                                                    throw new InvalidOperationException("Функция не вернула никакого значения");
                                                }
                                                return $call.Result;
                                              } catch (Exception ex) {
                                                Exception ex1 = null;
                                                while (stack.Count > 0) {
                                                  try { Operation<T>.DisposeFrame(stack.Pop()); }
                                                  catch (Exception ex2) { ex1 = ex2; }
                                                }
                                                /* Используется такой странный способ - чтобы по возможности сохранить место появления исходного исключения.
                                                    К сожалению, без рекурсии это невозможно в общем случае - но сохранить хотя бы верхний уровень надо попытаться. */
                                                if (ex1 != null)
                                                  throw ex1;
                                                else
                                                  throw;
                                              }
                                            }
                                            
                                            static IEnumerable<Operation<int>> RecFunc(int k) {
                                              var $call = new Call<int>();
                                            
                                              int s=1;
                                              for(int i=0;i<k;i++) {
                                                /* Здесь я нагло использую тот факт, что до первого вызова MoveNext() функция-генератор не получает управления */
                                                $call.Init(RecFunc(i)); 
                                                yield return $call;
                                                s+= $call.Result;
                                            
                                                $call.Init(RecFunc(i / 2));
                                                yield return $call;
                                                s+= $call.Result;
                                              }
                                              return new Result<T> { Value = s };
                                            }
                                            
                                              0
                                              А какие ещё есть интересные примеры нестандартного использования yield? Могу себе представить сопрограммы (две параллельно работающие функции, периодически передающие управление друг другу) и цепочку функций занимающихся обработкой большого потока данных, каждая из которых может в любой точке кода запросить промежуточный результат у предыдущей функции и отдать свой результат следующей. Хотя сейчас, конечно, это лучше делать на потоках.
                                                0
                                                Потоки занимают системные ресурсы, и переключаются в не самые подходящие моменты.

                                                В Unity что-то yield-подобное используется для реализации то ли асинхронности, то ли сопрограмм (я видел в декомпиляторе много «нестандартных» конечный автоматов).

                                                Но сейчас, конечно же, для тех же целей куда удобнее использовать задачи и async/await.
                                                  0
                                                  Так ведь задачи это те же потоки, только в другой обёртке? Или нет?
                                                  И реализовать сопрограммы на async/await… это как?
                                                    0
                                                    Нет. Задача — это отложенно вычисляемое значение — причем упор на самом значении, а не на вычислении.

                                                    Или еще можно назвать задачу однократным событием. Самый славный метод задачи — ContinueWith — это просто подписка на событие. «Вызови меня, когда значение станет определенным». Многопоточность задачами поддерживается (и эта поддержка и является причиной довольно высокой цены задач для простой экономии стека) — но не требуется.

                                                    Рекурсию на задачах можно переписать вот так:
                                                    static async Task<int> RecFunc(int k) {
                                                      int s=1;
                                                      for(int i=0;i<k;i++) {
                                                        await Task.Yield();
                                                    
                                                        s+= await RecFunc(i);
                                                        s+= await RecFunc(i / 2);
                                                      }
                                                      return s;
                                                    }
                                                    


                                                    Вся магия — в вызове await Task.Yield(). Этот вызов прерывает текущий процесс вычислений, ставя задачу в очередь в пул потоков (откуда она, естественно, сразу же забирается на исполнение обратно). В данном случае он используется для очистки стека вызовов.

                                                    Если хочется сделать вычисление совсем однопоточным — то можно сделать вот так:
                                                    class SequencedSynchronizationContext : SynchronizationContext {
                                                      private SendOrPostCallback callback;
                                                      private object state;
                                                      public override void Post(SendOrPostCallback d, Object state) {
                                                        // Сюда мы будем попадать при каждом await Task.Yield()
                                                        if (this.callback != null) throw new NotSupportedException();
                                                        this.callback = callback;
                                                        this.state = state;
                                                      }
                                                    
                                                      // Здесь надо бы перегрузить остальные методы, чтобы заставить их кидать NotSupportedException - но мне лень.
                                                    
                                                      public void Run() {
                                                        var old = SynchronizationContext.Current;
                                                        SetSynchronizationContext(this);
                                                        try {
                                                          while (callback != null) {
                                                            var cb = callback;
                                                            callback = null;
                                                            cb(state);
                                                          }
                                                        } finally {
                                                          SetSynchronizationContext(old);
                                                        }
                                                      }
                                                    }
                                                    
                                                    static int ExecRunFunc(int k) {
                                                      var ctx = new SequencedSynchronizationContext();
                                                      int result = 0;
                                                      ctx.Post(async () => result = await RecFunc(k));
                                                      ctx.Run();
                                                      return result;
                                                    }
                                                    
                                                      0
                                                      Ясно. Я посмотрел код, сгенерированный из этой RecFunc — страшная штука. Так же, как и в случае с yield, держит все переменные в отдельном классе и реализует state-машину с методом MoveNext, да ещё и выполняет кучу дополнительной работы, связанной с задачами.
                                                      Только один вопрос — как эту функцию вызывать, чтобы дождаться окончания работы? Я попробовал вот так:
                                                              static async void CallFunc() {
                                                                  int a=await RecFunc(4);
                                                                  Console.WriteLine(a);
                                                              }
                                                              static void Main(string[] args) {
                                                                  Task t=Task.Run(new Action(CallFunc));
                                                                  t.Wait();
                                                              }
                                                       

                                                      но она ждать не хочет.
                                                        0
                                                        Разумеется, ждать не хочет. Вы же ее два раза в другом потоке запустили (сначала через async void, потом через Task.Run).

                                                        Console.WriteLine(RecFunc(4).Result) — не пробовали?)

                                                        Ну или через ExecRunFunc.
                                                          0
                                                          Ладно. А как её подождать, если она запущена в другом потоке? Только через Events, или есть более канонические механизмы?

                                                          И насчёт сопрограмм вопрос остаётся. Пока мне удалось свести его к ситуации «задача A не может выдать результата, пока не будет запущена задача B. Последняя сообщает первой, какой у той будет результат, после чего A заканчивается, а B остаётся ждать нового запуска задачи A, чтобы узнать от неё свой результат». Но возможно ли такое вообще? Или в задачах сразу есть какой-нибудь аналог yield return?
                                                            0
                                                            Wait() и Return именно что ждут результатов задачи, которая, возможно, запущена в другом потоке (а возможно, ее еще никто даже не запускал).

                                                            По поводу же сопрограмм — посмотрите на «утиные интерфейсы» Awaitable и Awaiter — и как они реализованы в рамках Task.Yield().
                                                              0
                                                              Wait() и Return именно что ждут результатов задачи, которая, возможно, запущена в другом потоке

                                                              У меня не ждёт. t.Wait() заканчивается на 3 секунды раньше, чем await RecFunc(15). Или нужен другой Wait()?
                                                              Кстати, скорость работы RecFunc — 30000 вызовов в секунду :)
                                                                0
                                                                Вы замеряете время работы не RecFunc — а CallFunc. При этом время работы замеряется не целиком — а до первого оператора await. Не удивительно, что оно так быстро работает :)

                                                                1. Уберите Task.Run. У вас задача — вызвать RunFunc — или создать отдельный поток? Со второй задачей вы справились :)

                                                                2. Сделайте CallFunc возвращающим Task. Функция, объявленная как async void, переходит в фоновой режим работы без возможности дождаться завершения после первого же вызова await.
                                                                (Функция, объявленная как async Task, действует так же — но с Task можно дальше работать)
                                                                  0
                                                                  Спасибо, второй способ помог.
                                                                  Нет, это совсем не быстро :) После перехода на Release скорость возросла до 780000 вызовов в секунду (в 25 раз). Но простая рекурсивная реализация работает быстрее ещё в 400 раз… (300 миллионов вызовов).
                                                                  Реализация через yield и стек — 4 миллиона вызовов в секунду.
                                                                    0
                                                                    Блин, это не два способа — это два исправления к одному и тому же коду!

                                                                    Да, примерно такой оверхэд и должен был получиться. Еще бы: на каждый вызов создается TaskCompletionSource, который создает Task. Потом три вызова await. Если первый просто кладет продолжение в пул потоков — то последущие два создают новую задачу, привязанную к предыдущей — и это для того, что бы положить продолжение в пул потоков по совершению задачи. По завершению метода задача резолвится — с использованием атомарных операций, которые сами по себе не образец быстродействия (даже когда им никто не мешает).

                                                                    В случае же с yield ограничивающим фактором наверняка стал GC. Возвращаемый IEnumerable, к сожалению, структурой не сделать ну никак.
                                                                      0
                                                                      Там есть ещё один фактор, но не настолько драматичный. Даже если функции достаточно долгие, то код с yield выполняется примерно в 3-4 раза дольше, чем код без него. Просто потому, что локальные переменные приходится хранить в объекте, а не в регистрах.
                                                                0
                                                                Про сопрограммы нашёл непосредственную реализацию:
                                                                code.logos.com/blog/2010/10/coroutines_with_c_5s_await.html
                                                                Интересная штука. Правда, пока не понял, зачем там очереди — неужели только для того, чтобы чередовать два объекта?
                                                                  0
                                                                  Да, так и есть. Там используется очередь длины, не превышающей 2.
                                            0
                                            Проблема в другом. При использование этого метода не появляется никаких преимуществ перед обычной рекурсией.

                                            Как тут уже заметили, разве одно то, что стек не растет бесполезным образом с каждым рекурсивным вызовом, само по себе не достаточно?

                                            А в новом C++/CLI обещают аж автоматические конечные автоматы для использования в любых целях.

                                            Что значит «автоматические конечные автоматы» и что это дает? Где почитать?
                                  0
                                  И как же это будет выглядеть на примере quicksort? Сейчас он требует log(N) дополнительной памяти на вызовы. Удастся ли при разворачивании в цикл от неё избавиться?
                                  void QuickSort(T array[],int from,int to){
                                    if(to-from<=1) return;
                                    int mid=Split(array,from,to);
                                    QuickSort(array,from,mid);
                                    QuickSort(array,mid+1,to);
                                  } 
                                  
                                    0
                                    Если не ошибаюсь, то при представлении в виде обобщенной функции в gen<...> нужно подставить такие параметры:
                                    -- Генерируем функцию QuickSort. Принимает вектор аргументов. Все
                                    -- функции, из которых она состоит, также принимают вектор X = {array, from, to}.
                                    -- В реальном коде его можно было распаковать в аргументы функции.
                                    local QuickSort = getF({
                                      --Условие остановки рекурсии и значение при остановке
                                      check = function(X) return X[3] - X[2] <= 1 end,
                                      default = function() end,
                                      -- Преобразования аргументов для рекурсивных вызовов
                                      pre = {
                                        function(X) -- Для первого рекурсивного вызова
                                          local mid = Split(X[1], X[2], X[3])
                                          return {X[1], X[2], mid}
                                        end,
                                        function(X) -- Для второго рекурсивного вызова 
                                          local mid = Split(X[1], X[2], X[3])
                                          return {X[1], mid+1, X[3]}
                                        end,
                                      },
                                      -- Преобразование результата после рекурсивных вызовов
                                      post = function() end,
                                    });
                                    -- Вызов
                                    QuickSort({array, from, to});
                                    
                                      0
                                      Понятно, что её можно расписать как-то так:
                                      void QuickSort(T array[],int from,int to){
                                        var stack=null;
                                        while(to-from>1 || stack!=null){
                                          if(to-from>1){
                                            int mid=Split(array,from,to);
                                            stack=[mid+1,to,stack];
                                            to=mid;
                                          }else{
                                            (from,to,stack)=stack;
                                          }
                                        }
                                      }
                                      

                                      Но честно ли это?
                                        0
                                        А почему не должно быть честно? Выведенная мной формуля описывает точное преобразование. Возможно, я слишком быстро перескочил от итерации 2 к итерации 3, что это оказалось непонятно, вот более подробные выкладки (здесь итерация 6 — это то, что раньше было на итерации 3):
                                        -- Итерация 2: возвращаем N-1 предыдущих результатов из функции
                                        function fib2(x)
                                          if x < 2 then return -x, nil end;
                                        
                                          local x1 = fib2(x-1);
                                          local x2 = fib2(x-2);
                                          return x1 + x2, x1;
                                        end
                                        -- Итерация 3: используем все результаты, возвращаемые предыдущим вызовом
                                        function fib3(x)
                                          if x < 2 then return -x, nil end;
                                        
                                          local x1, p1 = fib3(x-1);-- p1 == x2 по определению
                                          local x2, p2 = fib3(x-2);-- p2 == x3 по определению
                                          return x1 + x2, x1;
                                        end
                                        -- Итерация 4: переименовываем результаты, чтобы увидеть истину
                                        function fib4(x)
                                          if x < 2 then return -x, nil end;
                                        
                                          local x1, x2 = fib4(x-1);
                                          local x2, x3 = fib4(x-2);
                                          return x1 + x2, x1;
                                        end
                                        -- Итерация 5: избавляемся от всех рекурсивных вызовов, кроме одного
                                        function fib5(x)
                                          if x < 2 then return -x, nil end;
                                        
                                          local x1, x2 = fib5(x-1);
                                          return x1 + x2, x1;
                                        end
                                        -- Итерация 6: вставляем защиту в начале итераций, пока не вычислим
                                        -- N (в данном случае, 2) значений, необходимых для рекурсии. Чтобы
                                        -- сохранить внешнюю сигнатуру, оборачиваем преобразованную функцию
                                        -- в функцию-хелпер
                                        function fib6(x)
                                          local function f(x, stop)
                                            if stop then return -x, nil end;
                                            if x < 2 then return -x, f(x-1, true) end;
                                        
                                            local x1, x2 = f(x-1, stop);
                                            return x1 + x2, x1;
                                          end
                                          return f(x, false);
                                        end
                                        

                                        На самом деле, не нужно было останавливаться на одном рекурсивном вызове, кажется, уже сейчас можно развернуть рекурсию в цикл. Нужно подумать. Кстати, рабочая функция genF (напомню, это функция, генерирующая исходную рекурсивную функцию с несколькими точками рекурсии по функциям условия выхода из рекурсии, значения выхода из рекурсии, преобразования аргументов для рекурсивного вызова и преобразования результатов после рекурсивного вызова) выглядит так:
                                        local function genF(t)
                                          local function none() end
                                        
                                          local check = t.check or none
                                          local default = t.default or none
                                          local pre = t.pre or {}
                                          local post = t.post or none
                                        
                                          local function f(...)
                                            -- Упаковываем все аргументы, принятые f в массив, иначе с
                                            -- ними работать нельзя
                                            local args = {...}
                                            -- unpack делает то же самое, что сделало бы
                                            -- return args[1], args[1], ..., args[#args]
                                            -- но только для произвольного количества элементов.
                                            -- Таким образом, во все функции передаются все значения,
                                            -- которые приняла f
                                            if check(unpack(args)) then return default(unpack(args)) end
                                        
                                            local results = {}
                                            for _, p in ipairs(pre) do
                                              results[#results + 1] = f(p(unpack(args)))
                                            end
                                        
                                            return post(unpack(results))
                                          end
                                          return f;
                                        end
                                        
                            +1
                            Деление близких (17485/17484) чисел даёт, наверное, даже большую ошибку. Так как ведущей оказывается единица (или девятка), а дальше куча нулей (девяток) перед первой осмысленной цифрой, на которые тратится вся точность типа. В принципе, как понимаю, оно тут и имеет место быть
                              +1
                              Насколько я понял из первого слинкованного документа:
                              xn = (3n+1+5n+1)/(3n+5n) без ошибок округления.
                              При их наличии
                              xn ≈ (3n+1+5n+1n*100n+1)/(3n+5nn*100n) где n = 3,4,5,… а γn — малое число представляющее собой ошибку округления.
                              Далее получается xn ≈ 100 — (95+97*(3/5)n)/(20nn + 1 + (3/5)n)
                              Т. е. само наличие малой ошибки (как бы мала она не была) смещает предел с 5 к 100.
                              Тут лучше сами документы почитать, так как там есть ещё более дикие примеры типа гладкой функции G(x).
                                0
                                Самую большую ошибку дает вычитание близких значений. Деление близких чисел как раз работает довольно точно — до тех пор, пока из результата не вычитается единица.
                                  0
                                  С чего бы?
                                  В случае вычитания близких значений разница в экспонентах минимальна или её нет вообще.
                                  Соответственно точность вычисления мантиссы максимальна.

                                  0.71
                                  0/01111110/01101011100001010001111

                                  минус

                                  0.7
                                  0/01111110/01100110011001100110011
                                  =

                                  0.00999999
                                  0/01111000/01000111101011100000000

                                  В результате можно обнаружить все биты разницы мантиссы 1.01000111101011100 (первая единица скрыта) и увеличившуюся на 6 разрядов точность представления числа

                                    0
                                    и увеличившуюся на 6 разрядов точность представления числа
                                    Точность не может вырасти магически. Если эти разряды появились из ниоткуда — значит, в них мусор!

                                    Если учитывать погрешности, то вычисление получается вот таким:

                                    0,71 (1 ± ε) — 0,7 (1 ± ε) = 0,01 ± 1,41 ε = 0,01 (1 ± 141 ε)

                                    Относительная погрешность возросла в 141 раз!
                                      0
                                      Эти 6 разрядов стали доступны для последующих вычислений в районе 0.01 — в этом вся суть плавающей точки.
                                      Ни 0.7, ни 0.71 невозможно точно представить на конечной разрядной сетке.
                                      Сама по себе операция ieee754(0,71) — ieee754(0.7) не теряет ни бита информации.
                                      Какая может быть погрешность у 0x35C28F — 0x333333?

                                        0
                                        При чем тут потери информации? Операция вычитания теряет точность, а не информацию.
                                          0
                                          Потеря точности => потеря информации. С потерей точности теряются информация о нижних битах.
                                            0
                                            Точность теряется пропорционально разнице экспонент. Вследствии выравнивания разрядной сетки параметров.
                                            В приведённом мной примере никакой потери точности нет.
                                              0
                                              Что вы понимаете под «точностью»? Лично я понимаю под этим словом величину, равную минус логарифму относительной погрешности.

                                              На всякий случай поясню, что такое «относительная погрешность». Это — потенциальное относительное отклонение вычисленной величины от (неизвестного) истинного значения.

                                              εx = (x — x`) / max(x, x`)
                                      0
                                      Не согласен. При вычитании/сложении мы теряем часть информации (если не всю) из меньшего числа, но в памяти у нас попрежнему находятся значимые цифры. При делении же, всю точность типа может заспамить нулями и девятками. И если у вас получится не «1» — уже повезло.

                                      Вообще, если при работе с числами с плавающей точкой есть возможность приблизить оперируемые числа к нулю — лучше это сделать.

                                      Например, ((2+ 0.0001) + (2 + 0.00000001))/2 приведёт к потере данных с большей вероятностью, чем 2 + (0.0001 + 0.00000001)/2.
                                      В любом случае, если вы ведёте какие-то серьёзные вычисления во float — надо понимать, зачем и что вы делаете. Яркий тому пример — статья про разработку SkyForge — там есть описание работы вывернутым буфером глубины.
                                        0
                                        Гм. Неудачно выразился:
                                        *f(((2+ 0.0001) + (2 + 0.00000001))/2) приведёт к потере данных с большей вероятностью, чем 2 + f'((0.0001 + 0.00000001)/2)
                                          +1
                                          Что значит «заспамить девятками»? Девятка — это значащая цифра. И даже если мы, округляя 0,999999999999999654, получаем 1,000000000000000 — это не проблема. Потому что число при таком округлении изменилось не сильно.

                                          «Хвост» числа не может быть более значимым, чем старшая 1. За исключением одного случая. Когда мы эту единицу потом вычитаем.

                                          В статье о разработке SkyForge говорилось про буфер глубины. Что такое буфер глубины? Это нечто, что сравнивает глубину точек. Что такое сравнение? Это вычитание и проверка знака. Вот именно наличие вычитания и является проблемой.

                                          Если бы не требовалось сравнивать глубины точек — то и выворачивать ничего не понадобилось бы.
                                            0
                                            > Что значит «заспамить девятками»? Девятка — это значащая цифра.
                                            Намного эффективнее будет представить значение в виде 1+ (a — b) / b. Так вы получите всю доступную мощность типа для хранения данных, так как (а-b)/b близок к нулю. Это я и называю хранением бессмысленных девяток.

                                            Причиной ошибки сравнения в статье про SkyForge на больших дистанция являлась как раз потеря точности при хранении больших чисел. Остававшейся точности типа не хватало для корректного сравнения.
                                              0
                                              Сама по себе величина числа никаких проблем не создает. Попробуйте умножить все числа на 1020: они от этого станут намного больше — но это не помешает их успешно сравнивать друг с другом. Или наоборот: поделите все числа на 1020 — они станут очень маленькими, но сравнению это ничуть не поможет.

                                              Проблемы исключительно из-за того, что числа близки друг к другу: модуль их разности много меньше модулей самих чисел.
                                                0
                                                Насколько я понял ситуацию — разработчики SkyForge столкнулись с ситуацией, когда надо было определять ближе или дальше 2 поверхности, находящиеся на расстоянии, допустим, 20 см друг от друга, но при этом — 20км от камеры:

                                                20 км = 2*10^4 м
                                                20 см = 2*10^-2 м.

                                                Разница в 6 порядков, при доп. операциях — это наверняка оказывается за пределами точности типа.
                                                То есть в больших числах важно некое довольно малое ɛ, которое тип хранить не способен.
                                                В таких условиях произвести сравнение вообще не возможно. Мы будем получать какой-то шум вместо актуальной картинки.
                                                  0
                                                  Совершенно верно, потому и надо менять формулу, избавляясь от таких сравнений.
                                      +1
                                      Думается мне, что если покрутить параметры(которые сейчас константами заданы) у такого отображения, то можно получить и несходящиеся итерации, а то и вовсе хаотическую последовательность.
                                        0
                                        У вас тоже первая мысль возникла — а есть ли хаотический режим у этой системы? :)
                                          0
                                          Ну, я не знаю как проверить это аналитически, но численно это можно было бы посмотреть по бассейнам притяжения по начальным условиям — если при изменении параметров есть тенденция к появлению фрактальных структур, значит скорее всего детерминированный хаос там есть.
                                            0
                                            По-моему, там ничего нет. Вы получите кривую
                                            z(p)= (3^p+5^p)/(3^(p-1)+5^(p-1))
                                            y(p)= (3^(p+1)+5^(p+1))/(3^p+5^p)
                                            точек, сходящихся к (5,5) — и всё. Остальные точки скатятся к (100,100). Фрактальности ожидать трудно, поскольку преобразование (y,z) -> (f(y,z),y) обратимо.
                                              0
                                              Позвольте полюбопытствовать, а как вы получили эти кривые?
                                                0
                                                Честно говоря, я её не проверял. Просто поверил, что исходная последовательность имеет вид (3^(n+1)+5^(n+1))/(3^n+5^n) — отсюда легко предположить, что все сходящиеся к 5 последовательности выглядят так же, только n смещены на дробную константу. И все точки (x_n,x_{n+1]) лягут на кривую (z(p),y(p)).
                                                Кроме того, построил несколько картинок, показывающих, куда и как быстро сходятся точки (y,z) — получилось, что все сходятся к (100,100), но в окрестности этой кривой число шагов возрастает.
                                                  +3
                                                  В PDF-ке все написано.

                                                  Основная идея — заменить xn на qn+1/qn, преобразовав «страшное» рекуррентное соотношение в обычное линейное. А все такие рекуррентные соотношения раскладываются в сумму геометрических прогрессий (ну, почти все) с основаниями, соответствующими корням характеристического многочлена.

                                                  Корни этого полинома — 3, 5 и 100 — потому qn = a 3n + b 5n + c 100n. Коэффициенты a, b и c зависят от начальных условий.

                                                  При этом, если c не равно 0 — то пределом исходной последовательности будет 100. Если c равно 0 — но b не равно 0 — пределом будет 5. Если b и c равны 0, а a не равно — пределом будет 3. Ну а все одновременно нулю они быть равны не могут.
                                                  0
                                                  Уточнение. Кривая является гиперболой и имеет уравнение y=8-15/z. Все точки этой кривой, кроме (3,3), сходятся к (5,5).
                                            0
                                            в матлабе расчет сразу сваливается на -7.
                                            4.0000
                                                4.2500
                                               -7.5147
                                             -130.7314
                                              217.9809
                                              114.1815
                                              104.3214
                                              100.9882
                                              100.3300
                                              100.0778
                                              100.0262
                                              100.0062
                                              100.0021
                                              100.0005
                                              100.0002
                                              100.0000
                                            
                                              0
                                              А вы точно exect'ом считали?
                                                0
                                                Это Вы просто попутали местами аргументы при вызове f.
                                                Вы считаете не
                                                x(i) = f (x(i-1), x(i-2)),
                                                a
                                                x(i) = f (x(i-2), x(i-1))

                                                Поменяйте местами и все у Вас сойдется. :)
                                                  0
                                                  И правда! А ведь я допустил ту же ошибку
                                                  0
                                                  На python и php получается подобный результат
                                                  +1
                                                  Любопытства ради попробовал сделать это на PHP с помощью BC Math Functions
                                                  Чтобы правильно досчитать до 30, нужно bcscale выставить в 42
                                                  Чтобы правильно досчитать до 100, нужно bcscale выставить в 149
                                                  Сделал bcscale ( 500 );
                                                  [330] => 4.999999999999999999999999999999999999999999999999999999999999999999999999885708
                                                  Далее до i=386 держится чуть более 5 и потом уходит к 100.
                                                    0
                                                    Только недавно попадал в похожую ловушку. Обычная с виду функция поиска корней квадратного уравнения с проверкой дискриминанта на близость к нулю (D < EPS). После изменения точности EPS с 1.0E-10 до 1.0E-100 результаты вычислений местами стали, мягко говоря, отличаться от правильных — причём на довольно простых (в смысле коэффициентов) уравнениях. Проблема та же — деление близких друг к другу чисел (отличающихся только N-ым знаком после запятой).
                                                      +4
                                                      Проблема описана хорошо. Но где предлагаемые методы решения? Не в wolfram'ах и в maple'ах, а на каком-нибудь промышленном ЯП?
                                                        0
                                                        Для питона использовать mpmath.org, для других языков тоже должны быть либы.
                                                          0
                                                          Ага, есть ещё Gnu MP, но взять подобную либу — это не интересно, ибо она скорее всего бьёт по производительности. Наверняка, есть какие-нибудь приемчики, которые позволяют организовывать порядок вычислений с плавающей точкой таким образом, чтобы отодвинуть эффекты от округления на более далёкий план.
                                                            0
                                                            На каждый хитрый приём найдётся контрпример, который поставит алгоритм на колени.
                                                              +2
                                                              Утверждение более чем голословное.
                                                                +3
                                                                Ну, посудите сами — мы хотим оперировать бесконечными множествами (легко показать, что «идеальное» рациональное число — счётное, хоть и бесконечное), но при этом хотим хранить их в нескольких десятках битов. Понятно, что какой бы мы трюк для уплотнения не использовали, всегда можно подать на вход столько информации, что потери будут значительны. Чтобы потери были значительны надо сделать так, чтобы значимыми были биты достаточно длительного участка, а это можно сделать выполнив набор операций, делающих значимыми всё более и более «разнесённые» друг от друга биты. В какой-то момент мы «вылетаем» за число доступных нам битов, и все последующие вычисления становятся всё более и более ошибочными, так как захватывают всё больше и больше «не существующих бит».
                                                                  0
                                                                  А не знаю кто такие «вы», кто хотите оперировать бесконечными множествами чисел, но есть вполне большой круг задач, (например, задачи, связанные с обработкой сигналов), где множества значений вполне себе ограничены (по абсолютной величине), но множественные операции с ними, если они организованы влоб, приведут к большому накоплению шума и потере информации. А главное очень критична производительность и использование multi precision библиотек не допустимо. И в таких задачах есть способы, как выполнить нужные вычисления, не внося много шума, а даже уменьшив его.

                                                                  Таким образом, не стоит говорить за всех, если вы в ваших задачах это не применимо.
                                                                    +1
                                                                    Число рациональных чисел, ограниченных по модулю, счётно (то есть, бесконечно), как и число рациональных чисел вообще. Ваш аргумент не принимается даже отчасти.

                                                                    Например, если мы об обработке сигналов — можно построить фильтр БИХ, который будет неустойчив (т.е. полностью забывать исходную информацию, перестав быть фильтром с «бесконечной» характеристикой, и сваливаться в шум) с любой мыслимой арифметикой, какую бы вы ни придумали, кроме бесконечно точной.
                                                                      0
                                                                      > Ваш аргумент не принимается даже отчасти.

                                                                      Ну-ка поясните подробнее каким образом сказанное мною связано с мощностью описанного компакта?

                                                                      > Например,

                                                                      Например, можно построить фильтр БИХ, который и с бесконечной точностью будет входить в генерацию.

                                                                      Извините, но я не рассматриваю теоретические придумки заранее неправильных решений.

                                                                      Вдумайтесь в то, что вы написали. Вы своим каментом утверждаете буквально следующее: «Если с помощью какого-то инструмента можно сотворить какую-нибудь хрень (типа вошедшего в генерацию БИХ-фильтра), то он абсолютно непригоден для применения нигде».
                                                                      +2
                                                                      Ещё раз: на каждый хитрый приём найдётся контрпример. Я не утверждаю, что не может быть локальных оптимизаций, я утверждаю, что под любую заранее заданную конечную оптимизацию (в т.ч. с «выбором» оптимизаций) можно придумать контр-пример, в котором эта оптимизация будет себя вести хуже, чем обычная математика на флоатах.

                                                                      По аналогии с архиваторами. Архиваторы работают отлично, но под любой архиватор можно придумать такой алгоритм, который даст на выходе файл большего размера, чем на входе.
                                                                        –2
                                                                        В реальных задачах входные данные имеют ограничения. И если, например, КИХ-фильтр начнет выдавать бред от того, что на вход ему подадут данные, которых в реальной жизни не будет никогда, то гори оно огнём.
                                                              0
                                                              У питона есть встроенные дроби:

                                                              from fractions import Fraction
                                                              
                                                              def f(y, z):
                                                                  return 108 - (815 - 1500/z)/y
                                                                  
                                                              X = [Fraction(4), Fraction(4.25)]
                                                              
                                                              
                                                              for x in xrange(2, 31):
                                                                  print float(X[0])
                                                                  X = [
                                                                      X[1],
                                                                      f(X[1], X[0])
                                                                  ]
                                                              

                                                              +1
                                                              Универсального способа, наверное, нет, но стоит посмотреть в сторону Алгоритма Кэхэна
                                                                0
                                                                Спасибо. Я не имел ввиду универсальный способ, а вот подборку таких вот трюков типа Алгоритма Кэхэна в догонку к данной статье было бы интересно изучить.
                                                                  +1
                                                                  Те статьи, на которые автор сослался как раз за авторством Уильяма Кэхэна. Он занимается проблемами вычислений с плавающей точкой больше 30 лет и не видит универсального способа, только понимание, только тесты.
                                                                  Можно ещё наверное посоветовать его довольно известную статью Mathematics Written in Sand (pdf).
                                                                    0
                                                                    Спасибо. Я всего лишь хотел сказать, что автор выбрал интересную тему и совершенно её не развил. Статья в том виде, как она сейчас, должна была быть лишь затравкой к более большой статье (циклу статей), о том, что со всем этим делать. В частности автор мог бы привести Кэхэна.

                                                                    То, что универсального способа нет, это понятно. Но различные примеры для базовых случаев существуют и являются хорошей отправной точкой, чтобы «почувствовать» направление решения более сложных задач. И с этого и стоит начинать, чем изобретать самому. Об этом и должна была быть статья, ИМХО.
                                                                      0
                                                                      А вообще я туплю ) Только сейчас понял, что автор поста — это Вы, а сам пост — это перевод. А раз так, то упрекать вас бессмыслено. Но, кажется, днём статья не была помечена как перевод, или мне мерещится?
                                                                    0
                                                                    А такое обычно решать не надо. Как наверху написали, решение неустойчиво и в реальных физических задачах никогда не реализуется. А контроль сходимости к нужной точке всегда надо делать в любом случае.
                                                                      0
                                                                      Не смешите насчет не реализуется. Сплошь и рядом надо гладить.
                                                                        +2
                                                                        Речь не про конкретный пример и его применимость к физическим задачам. А про целый класс приемов, которые позволяют уменьшить накопляемость погрешности.

                                                                        Речь про то, что статья получилась незаконченной. В стиле: «а смотрите как можно сопроцессор раком поставить» с совершенно закономерным возникающим вопросом: «И что дальше?»
                                                                        +3
                                                                        Проверять формулы на устойчивость.

                                                                        В реальности не бывает чисел 4 и 4.25, если мы говорим о результатах измерений — то бывают лишь 4 и 4.25 плюс минус неизвестная погрешность. А если уже в исходных данных есть погрешность — то схождение к 100 становится, внезапно, правильным ответом.

                                                                        В этом плане куда хуже «баг», что последовательность может неожиданно сойтись к 5…
                                                                          0
                                                                          Вы кэпствуете. Проверка формулы на устройчивость — это очевидно. Но вот вы увидили, что ваше решение неустойчиво. Что делать дальше? Выбросить его в помойку? Или сразу вместе с компьютером? )
                                                                            0
                                                                            Дальше надо придумывать другое решение. Также следует проверить на устойчивость саму задачу. Возможно, не стоит и браться за задачу, выбросив ее на помойку, рядом с задачей останова.
                                                                        +2
                                                                        Все просто, x30 = 2328306745375394431036 / 465661390253305305137 :)

                                                                        Код на Haskell
                                                                        import Data.Ratio
                                                                        
                                                                        main :: IO ()
                                                                        main = print $ xs !! 30
                                                                          where xs = x0 : x1 : zipWith f (tail xs) xs
                                                                                x0 = 4
                                                                                x1 = toRational 4.25
                                                                                f y z = 108 - (815 - 1500 / z) / y
                                                                          +2
                                                                          PostgreSQL:
                                                                          > WITH RECURSIVE t(n, y, z) AS (
                                                                              VALUES (1, 4.25::numeric(160,150), 4::numeric(160,150))
                                                                            UNION ALL
                                                                              SELECT n+1, (108 - (815.0 - 1500.0/z) / y)::numeric(160,150), y FROM t WHERE n < 100
                                                                          )
                                                                          SELECT n, trunc(z, 30) FROM t;
                                                                          


                                                                          4.999999999999999999999781665718
                                                                          n | trunc
                                                                          -----+----------------------------------
                                                                          1 | 4.000000000000000000000000000000
                                                                          2 | 4.250000000000000000000000000000
                                                                          3 | 4.470588235294117647058823529411
                                                                          4 | 4.644736842105263157894736842105
                                                                          5 | 4.770538243626062322946175637393
                                                                          6 | 4.855700712589073634204275534441
                                                                          7 | 4.910847499082793200440259263788
                                                                          8 | 4.945537404123916724773383803167
                                                                          9 | 4.966962581762700598711938487257
                                                                          10 | 4.980045701355631161268607994290
                                                                          11 | 4.987979448478392260140132893976
                                                                          12 | 4.992770288062068097489592548328
                                                                          13 | 4.995655891506634026624028261567
                                                                          14 | 4.997391268381344112893869021642
                                                                          15 | 4.998433943944816919013897178190
                                                                          16 | 4.999060071970893867816839606286
                                                                          17 | 4.999435937146839147997850886499
                                                                          18 | 4.999661524103767537786887985727
                                                                          19 | 4.999796900713417912662993074633
                                                                          20 | 4.999878135477931249231732030020
                                                                          21 | 4.999926879504599904466452981025
                                                                          22 | 4.999956127061157738119015282294
                                                                          23 | 4.999973676005712444579015149356
                                                                          24 | 4.999984205520272707924179788142
                                                                          25 | 4.999990523282227659407207470022
                                                                          26 | 4.999994313958559593649811689520
                                                                          27 | 4.999996588371256023706380879336
                                                                          28 | 4.999997953021356907988412868139
                                                                          29 | 4.999998771812311329999364514508
                                                                          30 | 4.999999263087205784555322944942
                                                                          31 | 4.999999557852258305867636192615
                                                                          32 | 4.999999734711331524163448988670
                                                                          33 | 4.999999840826790469128306699144
                                                                          34 | 4.999999904496071241143611348508
                                                                          35 | 4.999999942697641650166096897705
                                                                          36 | 4.999999965618584596072420928552
                                                                          37 | 4.999999979371150615793644560440
                                                                          38 | 4.999999987622690318410255295625
                                                                          39 | 4.999999992573614172662417737389
                                                                          40 | 4.999999995544168496979305857826
                                                                          41 | 4.999999997326501095805051386575
                                                                          42 | 4.999999998395900656625319264598
                                                                          43 | 4.999999999037540393666415394250
                                                                          44 | 4.999999999422524236088689817269
                                                                          45 | 4.999999999653514541613196499408
                                                                          46 | 4.999999999792108724953511638899
                                                                          47 | 4.999999999875265234966920729470
                                                                          48 | 4.999999999925159140978285386289
                                                                          49 | 4.999999999955095484586299093272
                                                                          50 | 4.999999999973057290751537486102
                                                                          51 | 4.999999999983834374450835382511
                                                                          52 | 4.999999999990300624670469870213
                                                                          53 | 4.999999999994180374802270632782
                                                                          54 | 4.999999999996508224881358315504
                                                                          55 | 4.999999999997904934928813526203
                                                                          56 | 4.999999999998742960957287589006
                                                                          57 | 4.999999999999245776574372363786
                                                                          58 | 4.999999999999547465944623350009
                                                                          59 | 4.999999999999728479566773985431
                                                                          60 | 4.999999999999837087740064382411
                                                                          61 | 4.999999999999902252644038626262
                                                                          62 | 4.999999999999941351586423174610
                                                                          63 | 4.999999999999964810951853904353
                                                                          64 | 4.999999999999978886571112342463
                                                                          65 | 4.999999999999987331942667405424
                                                                          66 | 4.999999999999992399165600443235
                                                                          67 | 4.999999999999995439499360265934
                                                                          68 | 4.999999999999997263699616159558
                                                                          69 | 4.999999999999998358219769695733
                                                                          70 | 4.999999999999999014931861817440
                                                                          71 | 4.999999999999999408959117090463
                                                                          72 | 4.999999999999999645375470254278
                                                                          73 | 4.999999999999999787225282152566
                                                                          74 | 4.999999999999999872335169291540
                                                                          75 | 4.999999999999999923401101574924
                                                                          76 | 4.999999999999999954040660944954
                                                                          77 | 4.999999999999999972424396566972
                                                                          78 | 4.999999999999999983454637940183
                                                                          79 | 4.999999999999999990072782764110
                                                                          80 | 4.999999999999999994043669658466
                                                                          81 | 4.999999999999999996426201795079
                                                                          82 | 4.999999999999999997855721077047
                                                                          83 | 4.999999999999999998713432646228
                                                                          84 | 4.999999999999999999228059587737
                                                                          85 | 4.999999999999999999536835752642
                                                                          86 | 4.999999999999999999722101451585
                                                                          87 | 4.999999999999999999833260870951
                                                                          88 | 4.999999999999999999899956522570
                                                                          89 | 4.999999999999999999939973913542
                                                                          90 | 4.999999999999999999963984348125
                                                                          91 | 4.999999999999999999978390608875
                                                                          92 | 4.999999999999999999987034365325
                                                                          93 | 4.999999999999999999992220619195
                                                                          94 | 4.999999999999999999995332371517
                                                                          95 | 4.999999999999999999997199422910
                                                                          96 | 4.999999999999999999998319653742
                                                                          97 | 4.999999999999999999998991792177
                                                                          98 | 4.999999999999999999999395073945
                                                                          99 | 4.999999999999999999999637017138
                                                                          100 | 4.999999999999999999999781665718
                                                                          (100 строк)


                                                                          Для double precision:
                                                                          100 на 27 итерации
                                                                          n | z
                                                                          -----+-------------------
                                                                          1 | 4
                                                                          2 | 4.25
                                                                          3 | 4.47058823529412
                                                                          4 | 4.64473684210522
                                                                          5 | 4.77053824362508
                                                                          6 | 4.85570071256856
                                                                          7 | 4.91084749866063
                                                                          8 | 4.94553739553051
                                                                          9 | 4.966962408041
                                                                          10 | 4.98004220429301
                                                                          11 | 4.98790923279579
                                                                          12 | 4.99136264131455
                                                                          13 | 4.96745509555227
                                                                          14 | 4.42969049830883
                                                                          15 | -7.81723657845932
                                                                          16 | 168.939167671065
                                                                          17 | 102.039963152059
                                                                          18 | 100.09994751625
                                                                          19 | 100.004992040972
                                                                          20 | 100.000249579237
                                                                          21 | 100.00001247862
                                                                          22 | 100.000000623922
                                                                          23 | 100.000000031196
                                                                          24 | 100.00000000156
                                                                          25 | 100.000000000078
                                                                          26 | 100.000000000004
                                                                          27 | 100

                                                                          100 | 100
                                                                            +2
                                                                            в postgres все равно сваливается:
                                                                            117 | 4.9926409617676932037645167439762957583942847894
                                                                            118 | 4.8526022942835196983796692766367894215589402192
                                                                            119 | 1.9625018334983213858164746977314448782473776537
                                                                            120 | -149.7768320341942100936570873585969161142777208
                                                                            121 | 108.33830001081775763607044289660311421120052592
                                                                            122 | 100.38482697300886042382536684044495629547021586
                                                                              0
                                                                              Просто добавьте точности сколько вам нужно ;)
                                                                                0
                                                                                так в том-то и дело, что не понятно сколько добавлять. numeric(160,150) и так не мало
                                                                                  +2
                                                                                  Не поможет. Любая погрешность, сколь бы малой не была, очень быстро растет.
                                                                                0
                                                                                На питоне (кэш использовал, чтобы не ждать долго, пока просчитаются итерации для каждого i заново):

                                                                                from fractions import Fraction
                                                                                
                                                                                CACHE = {}
                                                                                
                                                                                
                                                                                def func(y, z):
                                                                                    return 108 - Fraction(815 - Fraction(1500, z), y) 
                                                                                
                                                                                
                                                                                def x(i):
                                                                                    if not i in CACHE:
                                                                                        if i == 0:
                                                                                            result = 4
                                                                                        elif i == 1:
                                                                                            result = Fraction('4.25')
                                                                                        else:
                                                                                            result = func(x(i - 1), x(i - 2))
                                                                                        CACHE[i] = result
                                                                                    return CACHE[i]
                                                                                
                                                                                
                                                                                def main():
                                                                                    for i in range(100):
                                                                                        print(i, float(x(i)))
                                                                                
                                                                                
                                                                                if __name__ == '__main__':
                                                                                    main()
                                                                                


                                                                                Вывод:

                                                                                (0, 4.0)
                                                                                (1, 4.25)
                                                                                (2, 4.470588235294118)
                                                                                (3, 4.644736842105263)
                                                                                (4, 4.770538243626063)
                                                                                (5, 4.855700712589074)
                                                                                ...
                                                                                (94, 5.0)
                                                                                (95, 5.0)
                                                                                (96, 5.0)
                                                                                (97, 5.0)
                                                                                (98, 5.0)
                                                                                (99, 5.0)
                                                                                
                                                                                  –1
                                                                                  Я проверил на Python 3 c float и c Decimal.

                                                                                  Float разваливается на 12
                                                                                  Код и результат
                                                                                  f = lambda y, z: (108 - (815 - 1500 / z) / y)
                                                                                  @functools.lru_cache()
                                                                                  def xer(i):
                                                                                      if i == 0: return 4.
                                                                                      if i == 1: return 4.25
                                                                                      return f(xer(i-1), xer(i-2))
                                                                                  

                                                                                  0 4.0
                                                                                  1 4.25
                                                                                  2 4.470588235294116
                                                                                  3 4.6447368421052175
                                                                                  4 4.770538243625083
                                                                                  5 4.855700712568563
                                                                                  6 4.91084749866063
                                                                                  7 4.945537395530508
                                                                                  8 4.966962408040999
                                                                                  9 4.980042204293014
                                                                                  10 4.987909232795786
                                                                                  11 4.991362641314552
                                                                                  12 4.967455095552268
                                                                                  13 4.42969049830883
                                                                                  14 -7.817236578459315
                                                                                  


                                                                                  Decimal держиться до 21
                                                                                  Decimal
                                                                                  @functools.lru_cache()
                                                                                  def xer(i):
                                                                                      if i == 0: return Decimal(4)
                                                                                      if i == 1: return Decimal('4.25')
                                                                                      return f(xer(i-1), xer(i-2))
                                                                                  

                                                                                  0 4
                                                                                  1 4.25
                                                                                  2 4.4705882352941176470588235
                                                                                  3 4.6447368421052631578947362
                                                                                  4 4.7705382436260623229461618
                                                                                  5 4.8557007125890736342039857
                                                                                  6 4.9108474990827932004342938
                                                                                  7 4.9455374041239167246519529
                                                                                  8 4.9669625817627005962571288
                                                                                  9 4.9800457013556311118526582
                                                                                  10 4.9879794484783912679439415
                                                                                  11 4.9927702880620482067468253
                                                                                  12 4.9956558915062356478184985
                                                                                  13 4.9973912683733697540253088
                                                                                  14 4.9984339437852482376781601
                                                                                  15 4.9990600687785413938424188
                                                                                  16 4.9994358732880376990501184
                                                                                  17 4.9996602467866575821700634
                                                                                  18 4.9997713526716167817979714
                                                                                  19 4.9993671517118171375788238
                                                                                  20 4.9897059157620938291040004
                                                                                  21 4.7951151851630947311130380
                                                                                  22 0.7281074924258006736651754
                                                                                  23 -581.7081261405031229400219627
                                                                                  24 105.8595186892360167901632650
                                                                                  

                                                                                    0
                                                                                    На самом деле до 11 и 18 соответственно. Дальше решение начинает отдаляться от истинного, сначала по чуть-чуть, потом быстрее.

                                                                                Only users with full accounts can post comments. Log in, please.