Свой натуральный логарифм - Помощь студентам

RIO · 10.02.2008, 22:53

Здравствуйте, уважаемые форумчане ! ! ! Кто подскажет как написать натуральный логарифм, для своего типа данных, разложение в ряд долго считается ! В нете много чего обыскал, на русских, английских сайтах так и ничего не нашел !

B_N · 10.02.2008, 23:08

А что есть "свой тип данных"? Да и рядом, вообще-то, мухой собирается.

RIO · 11.02.2008, 02:34

Свой модуль который реализирует арефметические операции и т.д.

Stilet · 11.02.2008, 08:34

Цитата:

разложение в ряд долго считается !

Как определил?

B_N · 11.02.2008, 16:59

Как-то не понял я, в чём проблема...

Применяем подстановку x = (1 + y)/(1 - y), решаем ее по y, получаем y = (x - 1) / (x + 1), раскладываем в ряд:

Код:

#include <stdio.h>
#include <math.h>

double ln(double x, double accuracy, long *piterations)
{
	double y;
	double ln, factor, dividend;
	double tmp;

	if(x <= 0 || accuracy <= 0) return -1;

	y = (x - 1) / (x + 1);
	ln = 0;
	factor = 1;
	dividend = 2 * y;
	long passes = 0;

	while(1)
	{
		tmp = (dividend / factor);
		if(abs(tmp) < abs(accuracy)) break;

		ln += tmp;
		factor += 2;
		dividend *= (y*y);
		passes ++;
	}
	if(piterations) *piterations = passes;
	return ln;
}

int main()
{
	double	arg = ******; // здесь указываем аргумент
	double	lnstd, lnX, accuracy = 0.1;
	long	i, iterations;

	lnstd = log(arg);

	for(i = 1; i <= 20; i++)
	{
		lnX = lnnew(arg, accuracy, &iterations);
		printf("Acc. = -%02d | Log(x) = %.16f | Log(x)Std = %.16f | Iterations : %2d | delta = %.4e\n", i, lnX, lnstd, iterations, abs(lnstd - lnX));
		accuracy /= 10;
	}
	getchar();
}

получаем

Код:

Acc. = -01 | Log(x) = 2.1481617624230838 | Log(x)Std = 2.3025850929940459 | Iterations :  3 | delta = 1.5442e-001
Acc. = -02 | Log(x) = 2.2861255067762709 | Log(x)Std = 2.3025850929940459 | Iterations :  7 | delta = 1.6460e-002
Acc. = -03 | Log(x) = 2.3003122649112631 | Log(x)Std = 2.3025850929940459 | Iterations : 11 | delta = 2.2728e-003
Acc. = -04 | Log(x) = 2.3023645999148750 | Log(x)Std = 2.3025850929940459 | Iterations : 16 | delta = 2.2049e-004
Acc. = -05 | Log(x) = 2.3025618750676657 | Log(x)Std = 2.3025850929940459 | Iterations : 21 | delta = 2.3218e-005
Acc. = -06 | Log(x) = 2.3025825262481305 | Log(x)Std = 2.3025850929940459 | Iterations : 26 | delta = 2.5667e-006
Acc. = -07 | Log(x) = 2.3025849025313314 | Log(x)Std = 2.3025850929940459 | Iterations : 32 | delta = 1.9046e-007
Acc. = -08 | Log(x) = 2.3025850706521371 | Log(x)Std = 2.3025850929940459 | Iterations : 37 | delta = 2.2342e-008
Acc. = -09 | Log(x) = 2.3025850903297980 | Log(x)Std = 2.3025850929940459 | Iterations : 42 | delta = 2.6642e-009
Acc. = -10 | Log(x) = 2.3025850927828877 | Log(x)Std = 2.3025850929940459 | Iterations : 48 | delta = 2.1116e-010
Acc. = -11 | Log(x) = 2.3025850929682252 | Log(x)Std = 2.3025850929940459 | Iterations : 53 | delta = 2.5821e-011
Acc. = -12 | Log(x) = 2.3025850929919494 | Log(x)Std = 2.3025850929940459 | Iterations : 59 | delta = 2.0965e-012
Acc. = -13 | Log(x) = 2.3025850929937852 | Log(x)Std = 2.3025850929940459 | Iterations : 64 | delta = 2.6068e-013
Acc. = -14 | Log(x) = 2.3025850929940246 | Log(x)Std = 2.3025850929940459 | Iterations : 70 | delta = 2.1316e-014
Acc. = -15 | Log(x) = 2.3025850929940432 | Log(x)Std = 2.3025850929940459 | Iterations : 75 | delta = 2.6645e-015
Acc. = -16 | Log(x) = 2.3025850929940455 | Log(x)Std = 2.3025850929940459 | Iterations : 81 | delta = 4.4409e-016
Acc. = -17 | Log(x) = 2.3025850929940455 | Log(x)Std = 2.3025850929940459 | Iterations : 86 | delta = 4.4409e-016
Acc. = -18 | Log(x) = 2.3025850929940455 | Log(x)Std = 2.3025850929940459 | Iterations : 92 | delta = 4.4409e-016
Acc. = -19 | Log(x) = 2.3025850929940455 | Log(x)Std = 2.3025850929940459 | Iterations : 98 | delta = 4.4409e-016
Acc. = -20 | Log(x) = 2.3025850929940455 | Log(x)Std = 2.3025850929940459 | Iterations : 103 | delta = 4.4409e-016

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

Код:

double lnnew(double x, double accuracy, long *piterations)
{
	double	y;
	double	result, factor, dividend;
	double	tmp;
	double	newx;
	long	magn = 0;

	if(x <= 0 || accuracy <= 0) return -1;

	newx = x;
	while(newx > 4)
	{
		newx /= 2;
		magn ++;
	}

	y = (newx - 1) / (newx + 1);
	result = 0;
	factor = 1;
	dividend = 2 * y;
	long passes = 0;

	while(1)
	{
		tmp = (dividend / factor);
		if(abs(tmp) < abs(accuracy)) break;

		result += tmp;
		factor += 2;
		dividend *= (y*y);
		passes ++;
	}
	if(piterations) *piterations = passes;

	if(magn > 0) return result + magn * lnnew(2, accuracy, NULL); else return result;
}

B_N · 11.02.2008, 16:59

И теперь видим, что, например, нат. логарифм гугола вычисляется за..

Код:

Acc. = -01 | Log(x) = 221.4493718714588200 | Log(x)Std = 230.2585092994045800 | Iterations :  1 | delta = 8.8091e+000
Acc. = -02 | Log(x) = 229.6621702679662500 | Log(x)Std = 230.2585092994045800 | Iterations :  2 | delta = 5.9634e-001
Acc. = -03 | Log(x) = 230.2106982219337400 | Log(x)Std = 230.2585092994045800 | Iterations :  3 | delta = 4.7811e-002
Acc. = -04 | Log(x) = 230.2543424661317900 | Log(x)Std = 230.2585092994045800 | Iterations :  4 | delta = 4.1668e-003
Acc. = -05 | Log(x) = 230.2581273276654500 | Log(x)Std = 230.2585092994045800 | Iterations :  5 | delta = 3.8197e-004
Acc. = -06 | Log(x) = 230.2584730521519400 | Log(x)Std = 230.2585092994045800 | Iterations :  6 | delta = 3.6247e-005
Acc. = -07 | Log(x) = 230.2584739323942800 | Log(x)Std = 230.2585092994045800 | Iterations :  8 | delta = 3.5367e-005
Acc. = -08 | Log(x) = 230.2585058865846600 | Log(x)Std = 230.2585092994045800 | Iterations :  9 | delta = 3.4128e-006
Acc. = -09 | Log(x) = 230.2585089642227000 | Log(x)Std = 230.2585092994045800 | Iterations : 10 | delta = 3.3518e-007
Acc. = -10 | Log(x) = 230.2585092660293400 | Log(x)Std = 230.2585092994045800 | Iterations : 11 | delta = 3.3375e-008
Acc. = -11 | Log(x) = 230.2585092960442000 | Log(x)Std = 230.2585092994045800 | Iterations : 12 | delta = 3.3604e-009
Acc. = -12 | Log(x) = 230.2585092990630600 | Log(x)Std = 230.2585092994045800 | Iterations : 13 | delta = 3.4152e-010
Acc. = -13 | Log(x) = 230.2585092993696800 | Log(x)Std = 230.2585092994045800 | Iterations : 15 | delta = 3.4902e-011
Acc. = -14 | Log(x) = 230.2585092994009100 | Log(x)Std = 230.2585092994045800 | Iterations : 16 | delta = 3.6664e-012
Acc. = -15 | Log(x) = 230.2585092994044400 | Log(x)Std = 230.2585092994045800 | Iterations : 17 | delta = 1.4211e-013
Acc. = -16 | Log(x) = 230.2585092994044700 | Log(x)Std = 230.2585092994045800 | Iterations : 18 | delta = 1.1369e-013
Acc. = -17 | Log(x) = 230.2585092994044700 | Log(x)Std = 230.2585092994045800 | Iterations : 19 | delta = 1.1369e-013
Acc. = -18 | Log(x) = 230.2585092994044700 | Log(x)Std = 230.2585092994045800 | Iterations : 20 | delta = 1.1369e-013
Acc. = -19 | Log(x) = 230.2585092994044700 | Log(x)Std = 230.2585092994045800 | Iterations : 22 | delta = 1.1369e-013
Acc. = -20 | Log(x) = 230.2585092994044700 | Log(x)Std = 230.2585092994045800 | Iterations : 23 | delta = 1.1369e-013

23 суммирования с точностью до 13-го знака. Остается чуть подправить условие, по которому определяется точность, чтобы учесть последующее умножение.

RIO · 11.02.2008, 23:25

B_N, спасибо тебе, от души спасибо большое, просто респект тебе, но есть довольно таки большая проблема, я не знаю С, программирую на Делфи только, если тебе не составит проблемы перевести этот код на Делфи, я просто всю жизнь тебе буду благодарен

B_N · 12.02.2008, 00:53

Код:

uses math, sysutils;

function logn(x, accuracy : double; var iterations : longint) : double;
var
    y, factor, dividend, tmp : double;
begin
    if( (x > 0) and (accuracy > 0) ) then begin
        y := (x - 1) / (x + 1);
        result := 0;
        factor := 1;
        dividend := 2 * y;
        iterations := 0;

        while(true) do begin
            tmp := (dividend / factor);
            if(abs(tmp) < abs(accuracy)) then break;

            result := result + tmp;
            factor := factor + 2;
            dividend := dividend * y * y;
            iterations := iterations + 1;
        end;
     end
    else result := -1;
end;

function lnnew(x, accuracy : double; var iterations : longint) : double;
var
    y, factor, dividend, tmp, newx : double;
    magn, dummy : longint;
begin
    if( (x > 0) and (accuracy > 0) )then begin
        magn := 0;
        newx := x;
        while(newx > 4) do begin newx := newx / 2; magn := magn + 1; end;

        y := (newx - 1) / (newx + 1);
        result := 0;
        factor := 1;
        dividend := 2 * y;
        iterations := 0;

        while(true) do begin
            tmp := (dividend / factor);
            if(abs(tmp) < abs(accuracy)) then break;

            result := result + tmp;
            factor := factor + 2;
            dividend := dividend * y * y;
            iterations := iterations + 1;
        end;
        if(magn > 0) then result := result + magn * lnnew(2, accuracy, dummy);
    end
    else result := -1;
end;


var
    arg, lnstd, lnX, accuracy : double;
    i, iterations : longint;

begin
    arg := 1e100;
    accuracy := 0.1;

    lnstd := ln(arg);

    for i := 1 to 20 do begin
        lnX := lnnew(arg, accuracy, iterations);
        writeln(format('Acc. = -%02d | Log(x) = %.16f | Log(x)Std = %.16f | Iterations : %2d | delta = %.4e',
            [i, lnX, lnstd, iterations, abs(lnstd - lnX)]));
        accuracy := accuracy / 10;
    end;
end.

И несколько комментариев.
1. Выше я написал, что ln(1e100) для такой-то точности получается за 23 прохода - это не совсем так, на самом деле раньше, а дальше идет уже бесполезная трата времени, связанная с представлением типа double (см. http://en.wikipedia.org/wiki/Double_precision ). У него в мантиссе 52 бита, что дает максимум 15-16 значащих десятичных знаков.
2. Тот логарифм, который получается у меня, самую малость (в последнем знаке) отличается от "системного". Видимо это связано с тем, что библиотечный использует типы данных сопроцессора (в паскале можно использовать extended вместо double) или с использованием другого разложения.
3. Логарифм двойки выбран не случайно - представление вещественных чисел дает возможность напрямую получать параметр, который я получал с помощью деления, из экспоненты числа (см. ссылку).

Оставляю Вам поэкспериментировать со всем этим.

RIO · 12.02.2008, 01:59

Огромное спасибо ! ! ! Только ещё одна проблема : Как тебе повысить репутацию ? Кнопка : Добавить отзыв вроде не работает

10.02.2008, 22:53	#1
RIO Форумчанин Регистрация: 15.12.2007 Сообщений: 422	Свой натуральный логарифм Здравствуйте, уважаемые форумчане ! ! ! Кто подскажет как написать натуральный логарифм, для своего типа данных, разложение в ряд долго считается ! В нете много чего обыскал, на русских, английских сайтах так и ничего не нашел !

12.02.2008, 00:53	#8
B_N Новичок Джуниор Регистрация: 18.01.2008 Сообщений: 1,720	Код: uses math, sysutils; function logn(x, accuracy : double; var iterations : longint) : double; var y, factor, dividend, tmp : double; begin if( (x > 0) and (accuracy > 0) ) then begin y := (x - 1) / (x + 1); result := 0; factor := 1; dividend := 2 * y; iterations := 0; while(true) do begin tmp := (dividend / factor); if(abs(tmp) < abs(accuracy)) then break; result := result + tmp; factor := factor + 2; dividend := dividend * y * y; iterations := iterations + 1; end; end else result := -1; end; function lnnew(x, accuracy : double; var iterations : longint) : double; var y, factor, dividend, tmp, newx : double; magn, dummy : longint; begin if( (x > 0) and (accuracy > 0) )then begin magn := 0; newx := x; while(newx > 4) do begin newx := newx / 2; magn := magn + 1; end; y := (newx - 1) / (newx + 1); result := 0; factor := 1; dividend := 2 * y; iterations := 0; while(true) do begin tmp := (dividend / factor); if(abs(tmp) < abs(accuracy)) then break; result := result + tmp; factor := factor + 2; dividend := dividend * y * y; iterations := iterations + 1; end; if(magn > 0) then result := result + magn * lnnew(2, accuracy, dummy); end else result := -1; end; var arg, lnstd, lnX, accuracy : double; i, iterations : longint; begin arg := 1e100; accuracy := 0.1; lnstd := ln(arg); for i := 1 to 20 do begin lnX := lnnew(arg, accuracy, iterations); writeln(format('Acc. = -%02d \| Log(x) = %.16f \| Log(x)Std = %.16f \| Iterations : %2d \| delta = %.4e', [i, lnX, lnstd, iterations, abs(lnstd - lnX)])); accuracy := accuracy / 10; end; end. И несколько комментариев. 1. Выше я написал, что ln(1e100) для такой-то точности получается за 23 прохода - это не совсем так, на самом деле раньше, а дальше идет уже бесполезная трата времени, связанная с представлением типа double (см. http://en.wikipedia.org/wiki/Double_precision ). У него в мантиссе 52 бита, что дает максимум 15-16 значащих десятичных знаков. 2. Тот логарифм, который получается у меня, самую малость (в последнем знаке) отличается от "системного". Видимо это связано с тем, что библиотечный использует типы данных сопроцессора (в паскале можно использовать extended вместо double) или с использованием другого разложения. 3. Логарифм двойки выбран не случайно - представление вещественных чисел дает возможность напрямую получать параметр, который я получал с помощью деления, из экспоненты числа (см. ссылку). Оставляю Вам поэкспериментировать со всем этим. Последний раз редактировалось B_N; 12.02.2008 в 01:33.

Похожие темы
Тема	Автор	Раздел	Ответов	Последнее сообщение
Есть ли у вас свой сайт?	Delpher	Свободное общение	43	23.05.2009 07:55
In как этот логарифм в Delphi записать	delphin100	Общие вопросы Delphi	1	02.06.2008 19:20
Свой hint	Altera	Общие вопросы Delphi	3	14.05.2008 15:01
Как создать свой формат?	Rahim1993	Общие вопросы Delphi	2	04.03.2008 17:04
Логарифм в паскале	gamer123	Помощь студентам	1	20.01.2008 15:15

10.02.2008, 23:08	#2
B_N Новичок Джуниор Регистрация: 18.01.2008 Сообщений: 1,720	А что есть "свой тип данных"? Да и рядом, вообще-то, мухой собирается.

11.02.2008, 02:34	#3
RIO Форумчанин Регистрация: 15.12.2007 Сообщений: 422	Свой модуль который реализирует арефметические операции и т.д.

11.02.2008, 23:25	#7
RIO Форумчанин Регистрация: 15.12.2007 Сообщений: 422	B_N, спасибо тебе, от души спасибо большое, просто респект тебе, но есть довольно таки большая проблема, я не знаю С, программирую на Делфи только, если тебе не составит проблемы перевести этот код на Делфи, я просто всю жизнь тебе буду благодарен

12.02.2008, 01:59	#9
RIO Форумчанин Регистрация: 15.12.2007 Сообщений: 422	Огромное спасибо ! ! ! Только ещё одна проблема : Как тебе повысить репутацию ? Кнопка : Добавить отзыв вроде не работает