Хабр Курсы для всех
РЕКЛАМА
Практикум, Хекслет, SkyPro, авторские курсы — собрали всех и попросили скидки. Осталось выбрать!
Хабр доступен 24/7 благодаря поддержке друзей
Комментарии 116
Как-то очень неожиданное поведение. особенно когда оно происходит не плавно, а вот так скачком.
Там нет скачка. У указанного рекуррентного соотношения есть три неподвижных точки: 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 уменьшится в полтора раза и это явление будет рекурсивно повторяться.
Из них первая и третья — устойчивы, а вторая — неустойчива. Это как в механике: вершина холма — точка неустойчивого равновесия, низ лунки — устойчивого.
Начальные значения подобраны специальным образом, чтобы при точных вычислениях сходиться к 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 уменьшится в полтора раза и это явление будет рекурсивно повторяться.
Как-то мне разные нехорошие мысли приходят в голову об оригинальном авторе.
Очевидно что у полинома три неподвижных точки: 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.
Ну и при чем здесь ошибки округления? А тем более разные языки программирования и матлабы?
Очевидно что у полинома три неподвижных точки: 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 она сходиться не может в принципе, это точка отталкивания а не притяжения.
При заданных начальных условиях xn = (3n+1+5n+1)/(3n+5n)
Очевидно, что предел равен 5, при вычислениях в рациональных числах для исходной рекуррентной формулы мы также получаем последовательные приближения к 5.
Т. е. ожидаемый и правильный ответ должен быть равен 5.
Другое дело, что действительно, точка 5 неустойчива и при вычислениях с привнесённой ошибкой мы сразу от неё уйдём.
Так в этом и есть смысл статьи: вычисления с плавающей точкой сложны и могут приводить к неожиданным результатам, даже когда числа не являются малыми/большими, надо внимательно за ними следить.
Очевидно, что предел равен 5, при вычислениях в рациональных числах для исходной рекуррентной формулы мы также получаем последовательные приближения к 5.
Т. е. ожидаемый и правильный ответ должен быть равен 5.
Другое дело, что действительно, точка 5 неустойчива и при вычислениях с привнесённой ошибкой мы сразу от неё уйдём.
Так в этом и есть смысл статьи: вычисления с плавающей точкой сложны и могут приводить к неожиданным результатам, даже когда числа не являются малыми/большими, надо внимательно за ними следить.
А как получилось, что точка 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, оба больше единицы.
Вы как всегда правы, я недоглядел.
В некоторых языках есть встроенная рациональная арифметика. Очень удобно, когда делить нужно.
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
При двух рекурсивных вызовах в середине функции? Не работает, конечно.
На самом деле, непонятно, почему конечно. Мне кажется, что любую рекурсивную функцию с несколькими вызовами рекурсии можно переписать в рекурсивную функцию с одним вызовом рекурсии, а затем развернуть ее в цикл (с аккумулятором). Чтобы это сделать, нужно преобразовать исходную функцию, чтобы она возвращала не только новое значение, а и 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
N-1 предыдущихГодится только для последовательностей вида Ak = f(Ak-1, ..., Ak-N), где f – какая-то «комбинирующая функция». Но последовательности могут быть и такие, что k-ый член определяется, например, с помощью k/2, k/3 и т.д.
Если взять совсем любую рекурсивную функцию вроде mergesort'а, то, кажется мне, никаким конечным N тут не обойтись.
Как раз mergesort без рекурсии пишется легко, потому что там сначала идут рекурсивные вызовы для половинок массива, а потом результаты сливаются. Проще всего — слить сначала пары соседних элементов, потом четвёрки, потом восьмёрки… С помощью небольшой магии можно добиться даже такого же порядка выполнения, как в рекурсивном mergesort.
А нет, причина там в том, что во-первых, дерево рекурсии известно сразу и полностью, а во-вторых, хранение результатов не требует памяти — они закодированы в сортируемом массиве. Поэтому можно выполнять вычисления прямо из глубины, с любого места и в любом порядке. Главное, знать, какие узлы дерева уже вычислены, а какие — нет.
Идея все равно остается той же — возвращаем из функции не скаляр очередного результата, а вектор, содержащий новый результат и все, от чего он зависит. На вход функции передается вектор, в котором все результаты рекурсивных вызовов, кроме одного, уже посчитаны, нужно только ими воспользоваться и передать дальше. Как именно функция зависит от предыдущих вызовов, неважно.
На самом деле, подход, конечно годится не для любых рекурсивных функций, а только для тех, которые представимы в виде функций такого вида:
Так что, например, функция, типа
именно данной процедурой не возьмется. Впрочем, мне кажется, что если еще немного обобщить, то формулу можно будет подогнать и под этот случай — каждый последующий pre может зависеть от любого предыдущего. Вроде не должно это особо усложнить алгоритм.
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 может зависеть от любого предыдущего. Вроде не должно это особо усложнить алгоритм.
На самом деле, существует общий алгоритм преобразования, работающий для любых функций. Сначала из функции делается конечный автомат — а затем вводится стек для хранения состояний. Все.
Проблема в том, что ни ваш алгоритм, ни мой не экономят память — а просто загоняют проблему в другое место. Это, конечно, полезно для безстековых платформ — но слабо применимо в обычных случаях.
Проблема в том, что ни ваш алгоритм, ни мой не экономят память — а просто загоняют проблему в другое место. Это, конечно, полезно для безстековых платформ — но слабо применимо в обычных случаях.
Как мне кажется, проблема как раз и состоит в том, чтобы строить этот конечный автомат — иначе бы поголовно все компиляторы гарантировали бы разворот рекурсии. Между тем очень обидно, что даже простейшая рекурсия, не требующая введения аккумуляторов, практически никаким компилятором не гарантируется, что будет оптимизирована. А ведь в свете смещения мейнстрима в сторону функциональщины, это должна быть приоритетная задача!
Нет никакой проблемы в построении конечного автомата на основе программы. Компилятор того же C# умеет делать это автоматически аж для двух классов методов (генераторы и асинхронные методы). А в новом C++/CLI обещают аж автоматические конечные автоматы для использования в любых целях.
Проблема в другом. При использование этого метода не появляется никаких преимуществ перед обычной рекурсией.
Проблема в другом. При использование этого метода не появляется никаких преимуществ перед обычной рекурсией.
Ещё есть вариант преобразования с помощью CPS:
вот обсуждение на StackOverflow.
Вкратце, идея в том, чтобы действия, которые надо выполнить с полученным результатом, передать в первый рекурсивный вызов. Тут тоже никаких преимуществ перед обычной рекурсией, разве что память, которая потребляется, обычно меняется со стека на кучу.
вот обсуждение на StackOverflow.
Вкратце, идея в том, чтобы действия, которые надо выполнить с полученным результатом, передать в первый рекурсивный вызов. Тут тоже никаких преимуществ перед обычной рекурсией, разве что память, которая потребляется, обычно меняется со стека на кучу.
При использование этого метода не появляется никаких преимуществ перед обычной рекурсией.
Трудно согласиться. Всё-таки, штатный стек программы ограничен. А преобразование рекурсивной функции в автомат со стековой памятью позволит выдержать эквивалент рекурсивного вызова глубиной в сотни миллионов, если такое вдруг понадобится.
Интересно, можно ли для организации такого автомата воспользоваться уже реализованной в компиляторе C# упаковкой контекста, которая используется в реализации лямбда-выражений и в конструкции yield return.
Да, можно — даже двумя способами (через IEnumerable и через Task). Но в каждом из них придется реализовывать рекурсивный вызов особым способом — и внимательно следить, чтобы он не свелся после компиляции к настоящему рекурсивному вызову.
Мда.
Было:
Стало:
Было:
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;
}
Но, в принципе, так можно реализовать перекрёстные вызовы любого числа функций…
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
}
}
}
}
И даже так:
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;
}
}
Вы слишком извратились с кодом. Ну, или недостаточно извратились…
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 };
}
А какие ещё есть интересные примеры нестандартного использования yield? Могу себе представить сопрограммы (две параллельно работающие функции, периодически передающие управление друг другу) и цепочку функций занимающихся обработкой большого потока данных, каждая из которых может в любой точке кода запросить промежуточный результат у предыдущей функции и отдать свой результат следующей. Хотя сейчас, конечно, это лучше делать на потоках.
Потоки занимают системные ресурсы, и переключаются в не самые подходящие моменты.
В Unity что-то yield-подобное используется для реализации то ли асинхронности, то ли сопрограмм (я видел в декомпиляторе много «нестандартных» конечный автоматов).
Но сейчас, конечно же, для тех же целей куда удобнее использовать задачи и async/await.
В Unity что-то yield-подобное используется для реализации то ли асинхронности, то ли сопрограмм (я видел в декомпиляторе много «нестандартных» конечный автоматов).
Но сейчас, конечно же, для тех же целей куда удобнее использовать задачи и async/await.
Так ведь задачи это те же потоки, только в другой обёртке? Или нет?
И реализовать сопрограммы на async/await… это как?
И реализовать сопрограммы на async/await… это как?
Нет. Задача — это отложенно вычисляемое значение — причем упор на самом значении, а не на вычислении.
Или еще можно назвать задачу однократным событием. Самый славный метод задачи — ContinueWith — это просто подписка на событие. «Вызови меня, когда значение станет определенным». Многопоточность задачами поддерживается (и эта поддержка и является причиной довольно высокой цены задач для простой экономии стека) — но не требуется.
Рекурсию на задачах можно переписать вот так:
Вся магия — в вызове await Task.Yield(). Этот вызов прерывает текущий процесс вычислений, ставя задачу в очередь в пул потоков (откуда она, естественно, сразу же забирается на исполнение обратно). В данном случае он используется для очистки стека вызовов.
Если хочется сделать вычисление совсем однопоточным — то можно сделать вот так:
Или еще можно назвать задачу однократным событием. Самый славный метод задачи — 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;
}
Ясно. Я посмотрел код, сгенерированный из этой 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();
}
но она ждать не хочет.
Разумеется, ждать не хочет. Вы же ее два раза в другом потоке запустили (сначала через async void, потом через Task.Run).
Ну или через ExecRunFunc.
Console.WriteLine(RecFunc(4).Result) — не пробовали?)Ну или через ExecRunFunc.
Ладно. А как её подождать, если она запущена в другом потоке? Только через Events, или есть более канонические механизмы?
И насчёт сопрограмм вопрос остаётся. Пока мне удалось свести его к ситуации «задача A не может выдать результата, пока не будет запущена задача B. Последняя сообщает первой, какой у той будет результат, после чего A заканчивается, а B остаётся ждать нового запуска задачи A, чтобы узнать от неё свой результат». Но возможно ли такое вообще? Или в задачах сразу есть какой-нибудь аналог yield return?
И насчёт сопрограмм вопрос остаётся. Пока мне удалось свести его к ситуации «задача A не может выдать результата, пока не будет запущена задача B. Последняя сообщает первой, какой у той будет результат, после чего A заканчивается, а B остаётся ждать нового запуска задачи A, чтобы узнать от неё свой результат». Но возможно ли такое вообще? Или в задачах сразу есть какой-нибудь аналог yield return?
Wait() и Return именно что ждут результатов задачи, которая, возможно, запущена в другом потоке (а возможно, ее еще никто даже не запускал).
По поводу же сопрограмм — посмотрите на «утиные интерфейсы» Awaitable и Awaiter — и как они реализованы в рамках Task.Yield().
По поводу же сопрограмм — посмотрите на «утиные интерфейсы» Awaitable и Awaiter — и как они реализованы в рамках Task.Yield().
Wait() и Return именно что ждут результатов задачи, которая, возможно, запущена в другом потоке
У меня не ждёт. t.Wait() заканчивается на 3 секунды раньше, чем await RecFunc(15). Или нужен другой Wait()?
Кстати, скорость работы RecFunc — 30000 вызовов в секунду :)
Вы замеряете время работы не RecFunc — а CallFunc. При этом время работы замеряется не целиком — а до первого оператора await. Не удивительно, что оно так быстро работает :)
1. Уберите Task.Run. У вас задача — вызвать RunFunc — или создать отдельный поток? Со второй задачей вы справились :)
2. Сделайте CallFunc возвращающим Task. Функция, объявленная как async void, переходит в фоновой режим работы без возможности дождаться завершения после первого же вызова await.
(Функция, объявленная как async Task, действует так же — но с Task можно дальше работать)
1. Уберите Task.Run. У вас задача — вызвать RunFunc — или создать отдельный поток? Со второй задачей вы справились :)
2. Сделайте CallFunc возвращающим Task. Функция, объявленная как async void, переходит в фоновой режим работы без возможности дождаться завершения после первого же вызова await.
(Функция, объявленная как async Task, действует так же — но с Task можно дальше работать)
Спасибо, второй способ помог.
Нет, это совсем не быстро :) После перехода на Release скорость возросла до 780000 вызовов в секунду (в 25 раз). Но простая рекурсивная реализация работает быстрее ещё в 400 раз… (300 миллионов вызовов).
Реализация через yield и стек — 4 миллиона вызовов в секунду.
Нет, это совсем не быстро :) После перехода на Release скорость возросла до 780000 вызовов в секунду (в 25 раз). Но простая рекурсивная реализация работает быстрее ещё в 400 раз… (300 миллионов вызовов).
Реализация через yield и стек — 4 миллиона вызовов в секунду.
Блин, это не два способа — это два исправления к одному и тому же коду!
Да, примерно такой оверхэд и должен был получиться. Еще бы: на каждый вызов создается TaskCompletionSource, который создает Task. Потом три вызова await. Если первый просто кладет продолжение в пул потоков — то последущие два создают новую задачу, привязанную к предыдущей — и это для того, что бы положить продолжение в пул потоков по совершению задачи. По завершению метода задача резолвится — с использованием атомарных операций, которые сами по себе не образец быстродействия (даже когда им никто не мешает).
В случае же с yield ограничивающим фактором наверняка стал GC. Возвращаемый IEnumerable, к сожалению, структурой не сделать ну никак.
Да, примерно такой оверхэд и должен был получиться. Еще бы: на каждый вызов создается TaskCompletionSource, который создает Task. Потом три вызова await. Если первый просто кладет продолжение в пул потоков — то последущие два создают новую задачу, привязанную к предыдущей — и это для того, что бы положить продолжение в пул потоков по совершению задачи. По завершению метода задача резолвится — с использованием атомарных операций, которые сами по себе не образец быстродействия (даже когда им никто не мешает).
В случае же с yield ограничивающим фактором наверняка стал GC. Возвращаемый IEnumerable, к сожалению, структурой не сделать ну никак.
Про сопрограммы нашёл непосредственную реализацию:
code.logos.com/blog/2010/10/coroutines_with_c_5s_await.html
Интересная штука. Правда, пока не понял, зачем там очереди — неужели только для того, чтобы чередовать два объекта?
code.logos.com/blog/2010/10/coroutines_with_c_5s_await.html
Интересная штука. Правда, пока не понял, зачем там очереди — неужели только для того, чтобы чередовать два объекта?
Проблема в другом. При использование этого метода не появляется никаких преимуществ перед обычной рекурсией.
Как тут уже заметили, разве одно то, что стек не растет бесполезным образом с каждым рекурсивным вызовом, само по себе не достаточно?
А в новом C++/CLI обещают аж автоматические конечные автоматы для использования в любых целях.
Что значит «автоматические конечные автоматы» и что это дает? Где почитать?
Слышу разные вещи, обещаний море — трогать самому не приходилось.
Вот нашел что-то похожее. habrahabr.ru/company/infopulse/blog/203724
Вот нашел что-то похожее. habrahabr.ru/company/infopulse/blog/203724
И как же это будет выглядеть на примере 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);
}
Если не ошибаюсь, то при представлении в виде обобщенной функции в 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});
Понятно, что её можно расписать как-то так:
Но честно ли это?
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;
}
}
}
Но честно ли это?
А почему не должно быть честно? Выведенная мной формуля описывает точное преобразование. Возможно, я слишком быстро перескочил от итерации 2 к итерации 3, что это оказалось непонятно, вот более подробные выкладки (здесь итерация 6 — это то, что раньше было на итерации 3):
На самом деле, не нужно было останавливаться на одном рекурсивном вызове, кажется, уже сейчас можно развернуть рекурсию в цикл. Нужно подумать. Кстати, рабочая функция genF (напомню, это функция, генерирующая исходную рекурсивную функцию с несколькими точками рекурсии по функциям условия выхода из рекурсии, значения выхода из рекурсии, преобразования аргументов для рекурсивного вызова и преобразования результатов после рекурсивного вызова) выглядит так:
-- Итерация 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
Деление близких (17485/17484) чисел даёт, наверное, даже большую ошибку. Так как ведущей оказывается единица (или девятка), а дальше куча нулей (девяток) перед первой осмысленной цифрой, на которые тратится вся точность типа. В принципе, как понимаю, оно тут и имеет место быть
Насколько я понял из первого слинкованного документа:
xn = (3n+1+5n+1)/(3n+5n) без ошибок округления.
При их наличии
xn ≈ (3n+1+5n+1+γn*100n+1)/(3n+5n+γn*100n) где n = 3,4,5,… а γn — малое число представляющее собой ошибку округления.
Далее получается xn ≈ 100 — (95+97*(3/5)n)/(20n*γn + 1 + (3/5)n)
Т. е. само наличие малой ошибки (как бы мала она не была) смещает предел с 5 к 100.
Тут лучше сами документы почитать, так как там есть ещё более дикие примеры типа гладкой функции G(x).
xn = (3n+1+5n+1)/(3n+5n) без ошибок округления.
При их наличии
xn ≈ (3n+1+5n+1+γn*100n+1)/(3n+5n+γn*100n) где n = 3,4,5,… а γn — малое число представляющее собой ошибку округления.
Далее получается xn ≈ 100 — (95+97*(3/5)n)/(20n*γn + 1 + (3/5)n)
Т. е. само наличие малой ошибки (как бы мала она не была) смещает предел с 5 к 100.
Тут лучше сами документы почитать, так как там есть ещё более дикие примеры типа гладкой функции G(x).
Самую большую ошибку дает вычитание близких значений. Деление близких чисел как раз работает довольно точно — до тех пор, пока из результата не вычитается единица.
С чего бы?
В случае вычитания близких значений разница в экспонентах минимальна или её нет вообще.
Соответственно точность вычисления мантиссы максимальна.
0.71
0/01111110/01101011100001010001111
минус
0.7
0/01111110/01100110011001100110011
=
0.00999999
0/01111000/01000111101011100000000
В результате можно обнаружить все биты разницы мантиссы 1.01000111101011100 (первая единица скрыта) и увеличившуюся на 6 разрядов точность представления числа
В случае вычитания близких значений разница в экспонентах минимальна или её нет вообще.
Соответственно точность вычисления мантиссы максимальна.
0.71
0/01111110/01101011100001010001111
минус
0.7
0/01111110/01100110011001100110011
=
0.00999999
0/01111000/01000111101011100000000
В результате можно обнаружить все биты разницы мантиссы 1.01000111101011100 (первая единица скрыта) и увеличившуюся на 6 разрядов точность представления числа
и увеличившуюся на 6 разрядов точность представления числаТочность не может вырасти магически. Если эти разряды появились из ниоткуда — значит, в них мусор!
Если учитывать погрешности, то вычисление получается вот таким:
0,71 (1 ± ε) — 0,7 (1 ± ε) = 0,01 ± 1,41 ε = 0,01 (1 ± 141 ε)
Относительная погрешность возросла в 141 раз!
Эти 6 разрядов стали доступны для последующих вычислений в районе 0.01 — в этом вся суть плавающей точки.
Ни 0.7, ни 0.71 невозможно точно представить на конечной разрядной сетке.
Сама по себе операция ieee754(0,71) — ieee754(0.7) не теряет ни бита информации.
Какая может быть погрешность у 0x35C28F — 0x333333?
Ни 0.7, ни 0.71 невозможно точно представить на конечной разрядной сетке.
Сама по себе операция ieee754(0,71) — ieee754(0.7) не теряет ни бита информации.
Какая может быть погрешность у 0x35C28F — 0x333333?
При чем тут потери информации? Операция вычитания теряет точность, а не информацию.
Потеря точности => потеря информации. С потерей точности теряются информация о нижних битах.
Точность теряется пропорционально разнице экспонент. Вследствии выравнивания разрядной сетки параметров.
В приведённом мной примере никакой потери точности нет.
В приведённом мной примере никакой потери точности нет.
Что вы понимаете под «точностью»? Лично я понимаю под этим словом величину, равную минус логарифму относительной погрешности.
На всякий случай поясню, что такое «относительная погрешность». Это — потенциальное относительное отклонение вычисленной величины от (неизвестного) истинного значения.
εx = (x — x`) / max(x, x`)
На всякий случай поясню, что такое «относительная погрешность». Это — потенциальное относительное отклонение вычисленной величины от (неизвестного) истинного значения.
εx = (x — x`) / max(x, x`)
Не согласен. При вычитании/сложении мы теряем часть информации (если не всю) из меньшего числа, но в памяти у нас попрежнему находятся значимые цифры. При делении же, всю точность типа может заспамить нулями и девятками. И если у вас получится не «1» — уже повезло.
Вообще, если при работе с числами с плавающей точкой есть возможность приблизить оперируемые числа к нулю — лучше это сделать.
Например, ((2+ 0.0001) + (2 + 0.00000001))/2 приведёт к потере данных с большей вероятностью, чем 2 + (0.0001 + 0.00000001)/2.
В любом случае, если вы ведёте какие-то серьёзные вычисления во float — надо понимать, зачем и что вы делаете. Яркий тому пример — статья про разработку SkyForge — там есть описание работы вывернутым буфером глубины.
Вообще, если при работе с числами с плавающей точкой есть возможность приблизить оперируемые числа к нулю — лучше это сделать.
Например, ((2+ 0.0001) + (2 + 0.00000001))/2 приведёт к потере данных с большей вероятностью, чем 2 + (0.0001 + 0.00000001)/2.
В любом случае, если вы ведёте какие-то серьёзные вычисления во float — надо понимать, зачем и что вы делаете. Яркий тому пример — статья про разработку SkyForge — там есть описание работы вывернутым буфером глубины.
Гм. Неудачно выразился:
*f(((2+ 0.0001) + (2 + 0.00000001))/2) приведёт к потере данных с большей вероятностью, чем 2 + f'((0.0001 + 0.00000001)/2)
*f(((2+ 0.0001) + (2 + 0.00000001))/2) приведёт к потере данных с большей вероятностью, чем 2 + f'((0.0001 + 0.00000001)/2)
Что значит «заспамить девятками»? Девятка — это значащая цифра. И даже если мы, округляя 0,999999999999999654, получаем 1,000000000000000 — это не проблема. Потому что число при таком округлении изменилось не сильно.
«Хвост» числа не может быть более значимым, чем старшая 1. За исключением одного случая. Когда мы эту единицу потом вычитаем.
В статье о разработке SkyForge говорилось про буфер глубины. Что такое буфер глубины? Это нечто, что сравнивает глубину точек. Что такое сравнение? Это вычитание и проверка знака. Вот именно наличие вычитания и является проблемой.
Если бы не требовалось сравнивать глубины точек — то и выворачивать ничего не понадобилось бы.
«Хвост» числа не может быть более значимым, чем старшая 1. За исключением одного случая. Когда мы эту единицу потом вычитаем.
В статье о разработке SkyForge говорилось про буфер глубины. Что такое буфер глубины? Это нечто, что сравнивает глубину точек. Что такое сравнение? Это вычитание и проверка знака. Вот именно наличие вычитания и является проблемой.
Если бы не требовалось сравнивать глубины точек — то и выворачивать ничего не понадобилось бы.
> Что значит «заспамить девятками»? Девятка — это значащая цифра.
Намного эффективнее будет представить значение в виде 1+ (a — b) / b. Так вы получите всю доступную мощность типа для хранения данных, так как (а-b)/b близок к нулю. Это я и называю хранением бессмысленных девяток.
Причиной ошибки сравнения в статье про SkyForge на больших дистанция являлась как раз потеря точности при хранении больших чисел. Остававшейся точности типа не хватало для корректного сравнения.
Намного эффективнее будет представить значение в виде 1+ (a — b) / b. Так вы получите всю доступную мощность типа для хранения данных, так как (а-b)/b близок к нулю. Это я и называю хранением бессмысленных девяток.
Причиной ошибки сравнения в статье про SkyForge на больших дистанция являлась как раз потеря точности при хранении больших чисел. Остававшейся точности типа не хватало для корректного сравнения.
Сама по себе величина числа никаких проблем не создает. Попробуйте умножить все числа на 1020: они от этого станут намного больше — но это не помешает их успешно сравнивать друг с другом. Или наоборот: поделите все числа на 1020 — они станут очень маленькими, но сравнению это ничуть не поможет.
Проблемы исключительно из-за того, что числа близки друг к другу: модуль их разности много меньше модулей самих чисел.
Проблемы исключительно из-за того, что числа близки друг к другу: модуль их разности много меньше модулей самих чисел.
Насколько я понял ситуацию — разработчики SkyForge столкнулись с ситуацией, когда надо было определять ближе или дальше 2 поверхности, находящиеся на расстоянии, допустим, 20 см друг от друга, но при этом — 20км от камеры:
20 км = 2*10^4 м
20 см = 2*10^-2 м.
Разница в 6 порядков, при доп. операциях — это наверняка оказывается за пределами точности типа.
То есть в больших числах важно некое довольно малое ɛ, которое тип хранить не способен.
В таких условиях произвести сравнение вообще не возможно. Мы будем получать какой-то шум вместо актуальной картинки.
20 км = 2*10^4 м
20 см = 2*10^-2 м.
Разница в 6 порядков, при доп. операциях — это наверняка оказывается за пределами точности типа.
То есть в больших числах важно некое довольно малое ɛ, которое тип хранить не способен.
В таких условиях произвести сравнение вообще не возможно. Мы будем получать какой-то шум вместо актуальной картинки.
Думается мне, что если покрутить параметры(которые сейчас константами заданы) у такого отображения, то можно получить и несходящиеся итерации, а то и вовсе хаотическую последовательность.
У вас тоже первая мысль возникла — а есть ли хаотический режим у этой системы? :)
Ну, я не знаю как проверить это аналитически, но численно это можно было бы посмотреть по бассейнам притяжения по начальным условиям — если при изменении параметров есть тенденция к появлению фрактальных структур, значит скорее всего детерминированный хаос там есть.
По-моему, там ничего нет. Вы получите кривую
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) обратимо.
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) обратимо.
Позвольте полюбопытствовать, а как вы получили эти кривые?
Честно говоря, я её не проверял. Просто поверил, что исходная последовательность имеет вид (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), но в окрестности этой кривой число шагов возрастает.
Кроме того, построил несколько картинок, показывающих, куда и как быстро сходятся точки (y,z) — получилось, что все сходятся к (100,100), но в окрестности этой кривой число шагов возрастает.
В 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. Ну а все одновременно нулю они быть равны не могут.
Основная идея — заменить 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. Ну а все одновременно нулю они быть равны не могут.
Уточнение. Кривая является гиперболой и имеет уравнение y=8-15/z. Все точки этой кривой, кроме (3,3), сходятся к (5,5).
в матлабе расчет сразу сваливается на -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
Любопытства ради попробовал сделать это на PHP с помощью BC Math Functions
Чтобы правильно досчитать до 30, нужно bcscale выставить в 42
Чтобы правильно досчитать до 100, нужно bcscale выставить в 149
Сделал bcscale ( 500 );
[330] => 4.999999999999999999999999999999999999999999999999999999999999999999999999885708
Далее до i=386 держится чуть более 5 и потом уходит к 100.
Чтобы правильно досчитать до 30, нужно bcscale выставить в 42
Чтобы правильно досчитать до 100, нужно bcscale выставить в 149
Сделал bcscale ( 500 );
[330] => 4.999999999999999999999999999999999999999999999999999999999999999999999999885708
Далее до i=386 держится чуть более 5 и потом уходит к 100.
Только недавно попадал в похожую ловушку. Обычная с виду функция поиска корней квадратного уравнения с проверкой дискриминанта на близость к нулю (D < EPS). После изменения точности EPS с 1.0E-10 до 1.0E-100 результаты вычислений местами стали, мягко говоря, отличаться от правильных — причём на довольно простых (в смысле коэффициентов) уравнениях. Проблема та же — деление близких друг к другу чисел (отличающихся только N-ым знаком после запятой).
Проблема описана хорошо. Но где предлагаемые методы решения? Не в wolfram'ах и в maple'ах, а на каком-нибудь промышленном ЯП?
Для питона использовать mpmath.org, для других языков тоже должны быть либы.
Ага, есть ещё Gnu MP, но взять подобную либу — это не интересно, ибо она скорее всего бьёт по производительности. Наверняка, есть какие-нибудь приемчики, которые позволяют организовывать порядок вычислений с плавающей точкой таким образом, чтобы отодвинуть эффекты от округления на более далёкий план.
На каждый хитрый приём найдётся контрпример, который поставит алгоритм на колени.
Утверждение более чем голословное.
Ну, посудите сами — мы хотим оперировать бесконечными множествами (легко показать, что «идеальное» рациональное число — счётное, хоть и бесконечное), но при этом хотим хранить их в нескольких десятках битов. Понятно, что какой бы мы трюк для уплотнения не использовали, всегда можно подать на вход столько информации, что потери будут значительны. Чтобы потери были значительны надо сделать так, чтобы значимыми были биты достаточно длительного участка, а это можно сделать выполнив набор операций, делающих значимыми всё более и более «разнесённые» друг от друга биты. В какой-то момент мы «вылетаем» за число доступных нам битов, и все последующие вычисления становятся всё более и более ошибочными, так как захватывают всё больше и больше «не существующих бит».
А не знаю кто такие «вы», кто хотите оперировать бесконечными множествами чисел, но есть вполне большой круг задач, (например, задачи, связанные с обработкой сигналов), где множества значений вполне себе ограничены (по абсолютной величине), но множественные операции с ними, если они организованы влоб, приведут к большому накоплению шума и потере информации. А главное очень критична производительность и использование multi precision библиотек не допустимо. И в таких задачах есть способы, как выполнить нужные вычисления, не внося много шума, а даже уменьшив его.
Таким образом, не стоит говорить за всех, если вы в ваших задачах это не применимо.
Таким образом, не стоит говорить за всех, если вы в ваших задачах это не применимо.
Число рациональных чисел, ограниченных по модулю, счётно (то есть, бесконечно), как и число рациональных чисел вообще. Ваш аргумент не принимается даже отчасти.
Например, если мы об обработке сигналов — можно построить фильтр БИХ, который будет неустойчив (т.е. полностью забывать исходную информацию, перестав быть фильтром с «бесконечной» характеристикой, и сваливаться в шум) с любой мыслимой арифметикой, какую бы вы ни придумали, кроме бесконечно точной.
Например, если мы об обработке сигналов — можно построить фильтр БИХ, который будет неустойчив (т.е. полностью забывать исходную информацию, перестав быть фильтром с «бесконечной» характеристикой, и сваливаться в шум) с любой мыслимой арифметикой, какую бы вы ни придумали, кроме бесконечно точной.
> Ваш аргумент не принимается даже отчасти.
Ну-ка поясните подробнее каким образом сказанное мною связано с мощностью описанного компакта?
> Например,
Например, можно построить фильтр БИХ, который и с бесконечной точностью будет входить в генерацию.
Извините, но я не рассматриваю теоретические придумки заранее неправильных решений.
Вдумайтесь в то, что вы написали. Вы своим каментом утверждаете буквально следующее: «Если с помощью какого-то инструмента можно сотворить какую-нибудь хрень (типа вошедшего в генерацию БИХ-фильтра), то он абсолютно непригоден для применения нигде».
Ну-ка поясните подробнее каким образом сказанное мною связано с мощностью описанного компакта?
> Например,
Например, можно построить фильтр БИХ, который и с бесконечной точностью будет входить в генерацию.
Извините, но я не рассматриваю теоретические придумки заранее неправильных решений.
Вдумайтесь в то, что вы написали. Вы своим каментом утверждаете буквально следующее: «Если с помощью какого-то инструмента можно сотворить какую-нибудь хрень (типа вошедшего в генерацию БИХ-фильтра), то он абсолютно непригоден для применения нигде».
Ещё раз: на каждый хитрый приём найдётся контрпример. Я не утверждаю, что не может быть локальных оптимизаций, я утверждаю, что под любую заранее заданную конечную оптимизацию (в т.ч. с «выбором» оптимизаций) можно придумать контр-пример, в котором эта оптимизация будет себя вести хуже, чем обычная математика на флоатах.
По аналогии с архиваторами. Архиваторы работают отлично, но под любой архиватор можно придумать такой алгоритм, который даст на выходе файл большего размера, чем на входе.
По аналогии с архиваторами. Архиваторы работают отлично, но под любой архиватор можно придумать такой алгоритм, который даст на выходе файл большего размера, чем на входе.
У питона есть встроенные дроби:
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])
]
Универсального способа, наверное, нет, но стоит посмотреть в сторону Алгоритма Кэхэна
Спасибо. Я не имел ввиду универсальный способ, а вот подборку таких вот трюков типа Алгоритма Кэхэна в догонку к данной статье было бы интересно изучить.
Те статьи, на которые автор сослался как раз за авторством Уильяма Кэхэна. Он занимается проблемами вычислений с плавающей точкой больше 30 лет и не видит универсального способа, только понимание, только тесты.
Можно ещё наверное посоветовать его довольно известную статью Mathematics Written in Sand (pdf).
Можно ещё наверное посоветовать его довольно известную статью Mathematics Written in Sand (pdf).
Спасибо. Я всего лишь хотел сказать, что автор выбрал интересную тему и совершенно её не развил. Статья в том виде, как она сейчас, должна была быть лишь затравкой к более большой статье (циклу статей), о том, что со всем этим делать. В частности автор мог бы привести Кэхэна.
То, что универсального способа нет, это понятно. Но различные примеры для базовых случаев существуют и являются хорошей отправной точкой, чтобы «почувствовать» направление решения более сложных задач. И с этого и стоит начинать, чем изобретать самому. Об этом и должна была быть статья, ИМХО.
То, что универсального способа нет, это понятно. Но различные примеры для базовых случаев существуют и являются хорошей отправной точкой, чтобы «почувствовать» направление решения более сложных задач. И с этого и стоит начинать, чем изобретать самому. Об этом и должна была быть статья, ИМХО.
А вообще я туплю ) Только сейчас понял, что автор поста — это Вы, а сам пост — это перевод. А раз так, то упрекать вас бессмыслено. Но, кажется, днём статья не была помечена как перевод, или мне мерещится?
А такое обычно решать не надо. Как наверху написали, решение неустойчиво и в реальных физических задачах никогда не реализуется. А контроль сходимости к нужной точке всегда надо делать в любом случае.
Не смешите насчет не реализуется. Сплошь и рядом надо гладить.
Речь не про конкретный пример и его применимость к физическим задачам. А про целый класс приемов, которые позволяют уменьшить накопляемость погрешности.
Речь про то, что статья получилась незаконченной. В стиле: «а смотрите как можно сопроцессор раком поставить» с совершенно закономерным возникающим вопросом: «И что дальше?»
Речь про то, что статья получилась незаконченной. В стиле: «а смотрите как можно сопроцессор раком поставить» с совершенно закономерным возникающим вопросом: «И что дальше?»
Проверять формулы на устойчивость.
В реальности не бывает чисел 4 и 4.25, если мы говорим о результатах измерений — то бывают лишь 4 и 4.25 плюс минус неизвестная погрешность. А если уже в исходных данных есть погрешность — то схождение к 100 становится, внезапно, правильным ответом.
В этом плане куда хуже «баг», что последовательность может неожиданно сойтись к 5…
В реальности не бывает чисел 4 и 4.25, если мы говорим о результатах измерений — то бывают лишь 4 и 4.25 плюс минус неизвестная погрешность. А если уже в исходных данных есть погрешность — то схождение к 100 становится, внезапно, правильным ответом.
В этом плане куда хуже «баг», что последовательность может неожиданно сойтись к 5…
Все просто, x30 = 2328306745375394431036 / 465661390253305305137 :)
Код на Haskell
Код на 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
PostgreSQL:
Для double precision:
> 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 строк)
-----+----------------------------------
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
-----+-------------------
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
в postgres все равно сваливается:
117 | 4.9926409617676932037645167439762957583942847894
118 | 4.8526022942835196983796692766367894215589402192
119 | 1.9625018334983213858164746977314448782473776537
120 | -149.7768320341942100936570873585969161142777208
121 | 108.33830001081775763607044289660311421120052592
122 | 100.38482697300886042382536684044495629547021586
117 | 4.9926409617676932037645167439762957583942847894
118 | 4.8526022942835196983796692766367894215589402192
119 | 1.9625018334983213858164746977314448782473776537
120 | -149.7768320341942100936570873585969161142777208
121 | 108.33830001081775763607044289660311421120052592
122 | 100.38482697300886042382536684044495629547021586
На питоне (кэш использовал, чтобы не ждать долго, пока просчитаются итерации для каждого 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)
Я проверил на Python 3 c float и c Decimal.
Float разваливается на 12
Decimal держиться до 21
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
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Рекуррентное соотношение Мюллера: проблемы с округлением чисел с плавающей точкой