ln(exp(x)+1)

есть функция ln(exp(x) + 1). нужно уметь вычислять её значения. причём, как для x ~0, так и для x >> 1. проблема в том, что при x~100 происходит переполнение при вычислении экспоненты. первое что приходит в голову:

Код (Text):

1. if (x > 3)
2. {
3.  return x;
4. }
5. else
6. {
7.  return ln(exp(x) + 1);
8. }

может ещё как-то можно?

 

хитрее придумал

y(x) = ln(exp(x) + 1);
y(x) = x + q(x);
q(x) = ln(exp(x) + 1) - x = ln(exp(x) + 1) - ln(exp(x)) = ln(exp(-x) + 1);

вобщем вот так:
y(x) = x + ln(exp(-x) + 1);

 

Ну ты герой, в Паскале по крайней мере exp(-x) также легко переполняется.

If x > ln(1e16) then return x else return ln(exp(x)+1)

 

с чего это, если аргумент только положителен?

 

Так она сначала щщитает exp(), а потом 1/x =)))

тождество у тебя интересное получилось, но я бы по своему считал, как привел выше

 

чего..? да ну нафиг. ни за что не поверю, что exp(-x) вычисляется как 1/exp(x). если только какими-нибудь извращенцами.

 

2010_07_14

загрузил новую версию, извините только, что без сорсов: их я выкосил.

Код (Text):

1. ; tbyte[logXY]=log(tbyte[x];tbyte[y]); logYX=log(y;x)=1/logXY
2. ; restore: fpu.
3. ; not err if(logXY=x)|(logXY=y)|(logYX=x)|(logYX=y)|(logXY=logYX)
4. proc log, x:dword, y:dword, logXY:dword, logYX:dword
5. ; log(x;y) = log(x;2)*log(2;y) = log(2;y)/log(2;x)
6.         push    eax
7.         lea     esp,[esp-108]
8.         fsave   [esp]
9.  
10.         fld1
11.         mov     eax,[x]
12.         fld     tbyte[eax]
13.         fabs
14.         fyl2x
15.         fld1
16.         mov     eax,[y]
17.         fld     tbyte[eax]
18.         fabs
19.         fyl2x
20.         fdivrp
21.         fld     st0
22.         mov     eax,[logXY]
23.         fstp    tbyte[eax]
24.         fld1
25.         fdivrp
26.         mov     eax,[logYX]
27.         fstp    tbyte[eax]
28.  
29.         frstor  [esp]
30.         lea     esp,[esp+108]
31.         pop     eax
32.         ret
33. endp

 

2010_07_14
извините за дурацкий комментарий

excuse me for a stupid comment :c

 

edemko это, конечно, хорошо, но вычислительные ресурсы у меня не x86.

 

эта функция забавно раскладывается в ряд Тейлора с a=0

 

ну теперь еще заметнее стало, что второй член просто поправка к тривиальному x=х.
Скока там твой арм дает значащих цифр? Боюсь, те 16 что я привел для даблов слишком жирно будет =))). Просто поставь порог.

 

ещё такой вопрос. может лучше степенной ряд посчитать? с точки зрения производительности.. один фиг арм эмулирует логарифм и экспоненту. только он их на удивление быстро эмулирует - вот эта формула занимает ~4000 тактов..

 

a ты смогешь прописать ряд вокруг ноля?
Пока пришло на ум вот что ln2 + 1/2x +...

 

еще один челен осилил +1/8x^2

 

Ряд Тейлора в нуле для функции ln(exp(x)+1) имеет вид

(B_n - числа Бернулли), только он тут мало чем поможет, ибо его радиус сходимости равен пи, так что ряд расходится при |x|>Pi.

 

diamond а как проц считает логарифм? он ведь вообще сходится в очень небольшом пределе.

 

Ну я через микроскоп процессоры не разглядывал, что на самом деле происходит, точно не знаю, но варианты могут быть следующие:
- свести вычисление логарифма к логарифму от числа из (0.5,1], представив аргумент в виде 2^k*y, k целое, y меньше единицы (собственно, именно так действительные числа обычно и хранятся) и применив формулу
ln(2^k*y) = k*ln(2)+ln(y), константу ln(2) можно предвычислить с любой нужной точностью;
- использовать не ряд Тейлора в единице, а разложение функции ln((1+x)/(1-x)), там тоже получается ряд с радиусом сходимости 1 (формула, кстати, есть на википедии в статье про логарифм), но когда x меняется от -1 до 1, аргумент логарифма пробегает все положительные значения;
- использовать комбинацию первых двух вариантов - ряд из второго сходится тем быстрее, чем ближе x к нулю, что соответствует аргументу логарифма, близкому к 1, как и в первом варианте.

 

cupuyc
Для больших значений аргумента твою функцию можно представить в виде
x + ln(1 + exp(-x)) = x + exp(-x) + exp(-2x)/2 + exp(-3x)/3 -...
Члены ряда быстро убывают, поэтому можно оставить первые несколько членов, а про вычисление экспоненты с большим отрицательным показателем написано в
http://www.cyberguru.ru/cpp-sources/algorithms/vytchislenie-eksponenty-sinusa-i-kosinusa-page2.html

 

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

И хотя ln(1+x) имеет радиус сходимости единица, а exp(x) имеет радиус сходимости бесконечность,
при эффективных вычислениях в обеих случаях аргумент сначала преобразуется к небольшой величине, например меньшей по абсолютному значению чем 0.5, а потом только подставляется в ряд. Конечная функция потом находится путем обратных преобразований. Для ln() это видимо явное прибавление порядка, для exp() это может быть последовательное возведение в квадрат.

Ничего похожего для ln(exp(x)+1) придумать нельзя, никаких преобразований большего входного аргумента к меньшему. Поэтому ряд для вычислений не пригоден.

Эта штука может тормозить если x примерно равен нулю

Так что остается то что дохтур сразу прописал на рецепте:
If x > ln(1e16) then return x else return ln(exp(x)+1)

 

А так тоже падает:

y(x) = x + ln((exp(-x/n))^n + 1); где n типа 100 (ну или можно оценку n сделать в зависимости от x).

?