-
Data: 2014-07-29 10:30:29
Temat: Re: szybki logarytm
Od: Borneq <b...@a...hidden.pl> szukaj wiadomości tego autora
[ pokaż wszystkie nagłówki ]W dniu 2014-07-29 09:52, Borneq pisze:
>> całość:
>> http://pastebin.com/WuYW6MTJ
> x*(38.7609531719041669376926170246289866685554061249
47+
> x*(-245.85943859848429206309024729452081797918448383
011+
> x*(1305.02133731409794430037909996049017213519358597
46+
> x*(-5483.3897186528672757505743075974516855722609750
103+
> x*(18373.3692859474882123764452828071424328947070318
93+
> x*(-49720.330484326774507017799362836307490489909887
034+
Jest metoda szybsza niż Hornera, zwłaszcza dla wielomianów wysokiego
stopnia, mnoży od razu przez współczynnik kwadratowy.
>>>>
wklejam z .doc, który 10 lat temu napisałem, trochę tu będzie
nieczytelne, nie wiem jaką ma nazwę ten algorytm
Najpierw wielomian doprowadzamy do postaci, dla której współczynniki
przy dwóch najwyższych potęgach są równe jedności. Następnie rozkładamy
wielomian czynnikami kwadratowymi postaci y^2-a_i
Załóżmy, że wielomian pn(x) ma przy najwyższej potędze współczynnik
równy jedności, czyli jest postaci:
pn(x) = xn+an-1xn-1+an-2xn-2+...a1x+a0
po podstawieniu y = x + (an-1-1)/n ===> x = y - (an-1-1)/n
otrzymamy qn(y)
qn(y) = yn+yn-1+bn-2yn-2+...b1y+b0
Jeżeli zadany wielomian ma przy najwyższej potędze współczynnik różny od
jedności, czyli jest to wielomian wn(z) postaci:
wn(z) = cnzn+cn-1zn-1+cn-2zn-2+...c1z+c0
wtedy, aby otrzymać pn(x) należy albo podzielić cały wielomian przez cn
, albo podstawić
x = Sqrtn(cn)*z ===> z = x/ Sqrtn(cn)
zaleznie od tego, gdy chcemy przeskalować rezultat czy argument funkcji.
qn(y) = yn+yn-1+bn-2yn-2+...b1y+b0
dzieląc qn(y) przez y2-a1 otrzymujemy dla pewnego a1
qn(y) = (y2-a1)*( cn-2*yn-2+ cn-3*yn-3+cn-4*yn-4+...+c1*y+c0)+?1y+?1
Chcemy dobrać takie a1, aby po podzieleniu wielomianu przez (y2-a1)
otrzymać wielomian, który również będzie miał przy najwyższych dwóch
potęgach współczynniki równe jedności oraz reszta z dzielenia wielomianu
przez (y2-a1) była pojedyńczą liczbą.
qn(y) = (y2-a1)*(yn-2+yn-3+cn-4*yn-4+...+c1*y+c0)+?1
gdzie
cn-2=cn-3=1 oraz ?1=0
wzór rekurencyjny na cj ma postać:
cj=bj+2+a1*cj+2 (j=n-4,n-5,....-2)
przy czym c-1 = ?1, c-2 = ?1
Uwzględniając, że ?1 = 0 przedstawiamy ?1 jako wyrażenie zależne od a1 i bi
?1 = c-1 = b+a1*c1=b1+a1*(b3+a1*c3) =
b1+a1*b3+a12(b3+a1*c5) = b-1+a1*b3+a12*b5+...+anr-1*b2r-1+a1r
(2r = n-2 dla n parzystego lub n-1 dla n nieparzystego)
podstawiając ?1=0, otrzymamy dla a1 równanie
b1+a1*b3+a12*b5+...+anr-1*b2r-1+a1r = 0
Załóżmy, że ma pierwiastek rzeczywisty, niech a1 będzie pierwiastkiem
wtedy otrzymamy rozkład
qn(y)=(yn-a1)*(yn-2+yn-3+cm-4*ym-4+...c1*y+c0)+?1
dzieląc wielomian w nawiasach, jeżeli znowu znajdziemy pierwiastek,
otrzymamy
qn(y)=(y2-a1)*((y2-a2)*(yn-4+yn-5+dn-6*yn-6+...+a0)+
?2)+?1
a kiedy nie ma pierwiastka
qn(y) = (y2-a1)*(y*(yn-3+yn-4+....c1)+c0)+?1
Gdy za każdym razem mamy znaleziony pierwiastek ai, wtedy qn(y) ma rozkład:
qn(y) = (...(((?ny2+y+as)*(y2-as-1)+?s-1)*(y2-as-2)+?s-2)...
)*(y2-a1)+?1
(*)
?n = 1 dla n parzystych; 0 dla n nieparzystych
s = 1/2 *n dla n parzystych; 1/2 *(n+1) dla n nieparzystych
Algorytm obliczania pn(x)
albo x = x*Sqrtn(an) gdy najpierw przeskalowujemy argument funkcji
y = x+ (an-1-1)/n
z = y2
ws = ?n*z+y+as
wj = wj+1*(z-aj)+?j (j=s-1,s-2,...1)
pn(x) = w1
albo pn(x) = an*w1 gdy na końcu przeskalowujemy rezultat funkcji
przykład
p4(x) = x4+5x3+3x2+2
a3=5 => y = x+1 => x := y-1
q4(y)= y4+y3-6y2+5y+1 podstawiając ?1=0
0 = 5+a1
stąd a1 = -5
i
q4(y) = (y2+5)*(y2+y-11)+56
algorytm obliczania:
y := x+1
z := y*y
w := z+y-11
w ;= w*(z+5)+56
1.>
0.000204700z^7+0.001439274z^6+0.008328596z^5+0.04163
5012z^4+0.166667986z^3+0.500006347z^2+0.999999901z+0
.999999801
w przedziale <-1;1> maksymalna różnica wynosi około 2e-7
2. z:=x/0.297178123
otrzymamy
3. x^7 + 2.0895004648ox^6 + 3.5932515591ox^5
+ 5.3381571582ox^4 + 6.3504088014ox^3
+ 5.6616347459ox^2 + 3.3649849203ox + 0.999999801
x = y-(2.08950044541-1)/7
x = y-0.155642921
4.
y^7 + y^6 + 2.1506749053oy^5 + 3.1691354977oy^4
+ 3.7604523590oy^3 + 3.3533330221oy^2 + 1.9930983727oy
+ 0.5923035272
wyliczone:
5.x^5+x4+1.4184012627x3+2.4368618551x2+2.7217944997x
+1.5688833150
a: -0.732273642600117 (pierwiastek
x^3+2.1506749053x^2+3.760452359x+1.9930983727)
reszta: -0.5565483727
ten nie ma pierwiastka, ponieważ x^2+1.4184012627x+2.7217944997 nie ma
pierwiastków
i rozkłada się na
x^4+x3+1.4184012627x2+2.4368618551x+2.7217944997
i resztę 1.5688833150
rozkłada się na y^2-a dla a=-2.4368618551 (pierwiastek równania
x+2.4368618551)
i wielomian stopnia drugiego
x^2+x-1.0184605924
i resztę 5.2036422682
Następne wpisy z tego wątku
- 29.07.14 12:51 Borneq
- 29.07.14 15:12 feldmarszałek tusk
- 29.07.14 16:10 Borneq
- 29.07.14 17:04 A.L.
- 29.07.14 18:56 bartekltg
- 29.07.14 19:01 bartekltg
- 29.07.14 19:09 bartekltg
- 29.07.14 19:58 A.L.
- 29.07.14 20:04 A.L.
- 29.07.14 20:10 feldmarszałek tusk
- 29.07.14 20:25 Borneq
- 30.07.14 16:41 bartekltg
- 01.08.14 21:28 Borneq
- 01.08.14 21:37 feldmarszałek tusk
- 01.08.14 23:10 R.e.m.e.K
Najnowsze wątki z tej grupy
- TCL - problem z escape ostatniego \ w nawiasach {}
- Nauka i Praca Programisty C++ w III Rzeczy (pospolitej)
- testy-wyd-sort - Podsumowanie
- Tworzenie Programów Nieuprzywilejowanych Opartych Na Wtyczkach
- Do czego nadaje się QDockWidget z bibl. Qt?
- Bibl. Qt jest sztucznie ograniczona - jest nieprzydatna do celów komercyjnych
- Co sciaga kretynow
- AEiC 2024 - Ada-Europe conference - Deadlines Approaching
- Jakie są dobre zasady programowania programów opartych na wtyczkach?
- sprawdzanie słów kluczowych dot. zła
- Re: W czym sie teraz pisze programy??
- Re: (PDF) Surgical Pathology of Non-neoplastic Gastrointestinal Diseases by Lizhi Zhang
- CfC 28th Ada-Europe Int. Conf. Reliable Software Technologies
- Młodzi programiści i tajna policja
- Ada 2022 Language Reference Manual to be Published by Springer
Najnowsze wątki
- 2024-11-14 Dobra zmiana
- 2024-11-14 Czy prezydent może ułaskawić od zadośćuczynienia? [A. Lepper odszkodowania]
- 2024-11-14 Gliwice => Network Systems Administrator (IT Expert) <=
- 2024-11-14 Gliwice => Administrator Systemów Sieciowych (Ekspert IT) <=
- 2024-11-13 Filtr do pompy ruskiej
- 2024-11-12 Gdzie kosz?
- 2024-11-13 elektrycznie
- 2024-11-12 Jebane kurwa, kurwy.
- 2024-11-13 karta parkingowa
- 2024-11-13 Wl/Wyl (On/Off) bialy/niebieski
- 2024-11-12 I3C
- 2024-11-13 Kraków => DevOps Engineer (Junior or Regular level) <=
- 2024-11-13 Łódź => Senior SAP HANA Developer <=
- 2024-11-13 Zabrze => Senior PHP Symfony Developer <=
- 2024-11-13 Karlino => Konsultant wewnętrzny SAP (FI/CO) <=