بسط تیلور توابع Ln و exp :
function Ln_Taylor(X : Extended; Count : Integer = 5):Extended;
const
n_0_1 = 0.01;
ln_1_01 = 0.00995033085316808284821535754426; // Ln(1.01)
n_1_0_High = (1.0 + n_0_1);
n_1_0_Low = (1.0 - n_0_1);
arr : array[0..21] of Extended = (
1.01, // 1.01 ^ (2 ^ 0 = 1)
1.0201, // 1.01 ^ (2 ^ 1 = 2)
1.04060401, // 1.01 ^ (2 ^ 2 = 4)
1.0828567056280801, // 1.01 ^ (2 ^ 3 = 8)
1.172578644923698520518625612016, // 1.01 ^ (2 ^ 4 = 16)
1.3749406785310970541622913505711, // 1.01 ^ (2 ^ 5 = 32)
1.8904618694795535717494712641154, // 1.01 ^ (2 ^ 6 = 64)
3.5738460799561286443002337196205, // 1.01 ^ (2 ^ 7 = 128)
12.77237580321778745481813393975, // 1.01 ^ (2 ^ 8 = 512)
163.13358365862322124527961077082, // 1.01 ^ (2 ^ 9 = 1024)
26612.566117305021291272917047288, // 1.01 ^ (2 ^ 10 = 2048)
708228675.34793125625127947801295, // 1.01 ^ (2 ^ 11 = 4096)
501587856585085410.3329144226122, // 1.01 ^ (2 ^ 12 = 8192)
2.515903778736202094338585203235e+35, // 1.01 ^ (2 ^ 13 = 16384)
6.3297718238591005455779113716959e+70, // 1.01 ^ (2 ^ 14 = 32768)
4.0066011342120564182784035050917e+141, // 1.01 ^ (2 ^ 15 = 65536)
1.6052852648669336927937427546953e+283, // 1.01 ^ (2 ^ 16 = 131072)
2.5769407815989014605653704987392e+566, // 1.01 ^ (2 ^ 17 = 262144)
6.6406237918675571564214160686182e+1132, // 1.01 ^ (2 ^ 18 = 524288)
4.4097884345117453067721383354009e+2265, // 1.01 ^ (2 ^ 19 = 1048576)
1.9446234037153549426176998268801e+4531, // 1.01 ^ (2 ^ 20 = 2097152)
3.7815601822774923352439350246879e+9062 // 1.01 ^ (2 ^ 21 = 4194304)
);
var
x_local, b, bb, b_n : Extended;
scale_count : Int32;
i, j : Integer;
begin
x_local := X;
scale_count := 0;
//------------------------------------------
// Ln(A/B) = Ln(A) - Ln(B)
// Ln(A) = Ln(A/B) + Ln(B)
// .......... (if B == 1.01 then)
//
// Ln(A) = Ln(A / 1.01) + Ln(1.01)
//
// ..........
// for i = 1 to count do
// {
// x /= 1.01
// }
// return Ln(x) + count * Ln(1.01)
//
//------------------------------------------
j := 22;
while ((x_local > n_1_0_High) and (j > 0)) do // x > 1.0
begin
Dec(j);
while ((x_local > n_1_0_High) and (arr[j] < x_local)) do
begin
x_local := x_local / arr[j];
scale_count := scale_count + Integer(1 shl j);
end;
end;
j := 22;
while ((x_local < n_1_0_Low) and (j > 0)) do // x < 1.0
begin
Dec(j);
while ((x_local < n_1_0_Low) and ((1.0 / arr[j]) > x_local)) do
begin
x_local := x_local * arr[j];
scale_count := scale_count - Integer(1 shl j);
end;
end;
// ---------- Start Taylor block -----------------------------------------------
//
// b^3 b^5 b^6 b^7
// Ln(a) = 2 ( b + ---- + ---- + ---- + ---- + ... )
// 3 5 7 9
// where :
// a - 1
// b = --------
// a + 1
//------------------------------------------------------------------------------
b := (x_local - 1.0) / (x_local + 1.0);
Result := b;
b_n := b;
bb := b * b;
i := 3;
while (i < (Count * 2)) do
begin
i := i + 2;
b_n := b_n * bb;
Result := Result + b_n / i;
end;
Result := Result * 2.0;
// ---------- End Taylor block -------------------------------------------------
Result := result + (ln_1_01 * scale_count);
end;
این بسط حول نقطه ی 1 سریع به جواب میرسه بخاطر همین تا جایی که میشه اعدادمون رو با روابط ریاضی میبریم نزدیک این مقدار بعد که Ln رو محاسبه کردیم دوباره با استفاده از روابط مذکور برمیگردونیم به مقدار مورد نظرمون.
به همین خاطر اول یک عدد نسبتا نزدیک به یک مثلا 1.01 رو انتخاب کردم بعد یک آرایه از توان های اون ساختم به این شکل :
اولین عضو آرایه = 1.01
دومین عضو آرایه = عضو قبلی ضرب در خودش
سومین عضو آرایه = عضو قبلی ضرب در خودش
...
بعد عدد X رو تا جایی که میشه به این اعداد تقسیم کردم تا کلا به ضریبی از 1.01 تقسیم کرده باشیم (البته به مقدار کافی، نه زیاد نه کم). بعد Ln رو حساب میکنیم بعد از محاسبه به تعدادی که تقسیم بر 1.01 کردیم به همون تعداد با Ln(1.01) جمع میکنیم. جواب بدست میاد.
برای اعداد کوچک تر از 1 هم تقسیم به ضرب تبدیل میشه و جمع به منها.
---------------------------------------------------
exp :
function Exp_Taylor(X:Extended; Count : Integer = 10):Extended;
const
arr_pow2 : array[0..12] of Extended = (
2.00000000000000000000,
4.00000000000000000000,
16.0000000000000000000,
256.000000000000000000,
65536.0000000000000000,
4294967296.00000000000,
18446744073709551616.0,
3.4028236692093846346337460743177e+38,
1.1579208923731619542357098500869e+77,
1.3407807929942597099574024998206e+154,
1.797693134862315907729305190789e+308,
3.231700607131100730071487668867e+616,
1.0443888814131525066917527107166e+1233
);
e = 2.7182818284590452353602874713527;
var
i, j: Integer;
x_n, x_local, e_pow : Extended;
pow2_count : Int32;
begin
x_local := Abs(X);
pow2_count := 0;
j := 13;
while ((x_local > 0.01) and (j > 0)) do
begin
Dec(j);
while ((x_local > 0.01) and (0.001 < (x_local / arr_pow2[j]))) do
begin
x_local := x_local / arr_pow2[j];
pow2_count := pow2_count + Integer(1 shl (j));
end;
end;
(*
x_local := Abs(X);
j := Trunc(x_local);
x_local := x_local - j;
*)
//---------------------------
Result := 1.0;
x_n := 1.0;
for i := 1 to Count do
begin
x_n := x_n * (x_local / i);
Result := Result + x_n;
end;
//---------------------------
(*
e_pow := 1.0;
for i := 1 to j do
e_pow := e_pow * e;
Result := Result * e_pow;
if X < 0 then
Result := 1.0 / Result;
*)
for i := 1 to pow2_count do
Result := Result * Result;
if X < 0 then
Result := 1.0 / Result;
end;
برای exp هم روشی مشابه داریم. به تعداد کافی عدد رو بر 2 تقسیم میکنیم تا کوچک بشه برسه نزدیک صفر. بعد exp رو حساب میکنیم بعد از محاسبه عدمون رو به همون تعداد که X رو بر 2 تقسیم کردیم، به توان 2 میرسونیم. یعنی اگر n بار به دو تقسیم کردیم، عدد به دست اومده از exp رو به توان 2 به توان n میرسونیم. یعنی یک حلقه ی n مرحله ای ه هر بار عدد رو در خودش ضرب میکنیم.
روش دیگه ای که اون رو تو متن کد غیر فعال کردم هم روشی مشابه هست که میایم بخش اعشاری عدد رو نگه میداریم و بخش صحیح رو ازش یگیریم و exp رو برای بخش اعشاری بین 0 تا 1 حساب میکنیم بعد به تعدادی که در بخش صحیح ، 1 داریم (خود بخش صحیح) عددمون رو در e که exp عدد 1 هست ضرب میکنیم . این روش هم خوبه ولی روش تقسیمات متوالی اجازه میده که عددمون بیشتر به صفر نزدیک بشه. مثلا عدد 0.9 رو که نزدیک صفر هست میتونیم با تقسیمات متوالی بر 2 باز هم ک.چیکتر کنیم. پس من روش تقسیم به 2 رو ترجیح میدم.
