i ich pochodnych

ROCZNIKI POLSKIEGO TOWARZYSTWA MATEMATYCZNEGO SERIA III: MATEMATYKA STOSOWANA III (1974)

H.

WOŹNIAKOWSKI (Warszawa)

Analiza numeryczna algorytmów obliczania

wartości wielomianó~

i ich pochodnych

1. W ^poniższej pracy przedstawiamy jednoparametrową klasę algorytmów obliczania

wartości \vielomianu i jego znormalizowanych pochodnych rozpatrywaną przez Mary Shaw i J .F. Trauba w pracy [9].

Następnie dowodzimy, że numerycżna realizacja tych algorytmów w zmiennopozycyj- nej t-cyfrowej arytmetyce (fl) jest stabilna, co oznacza, ^{że każda} obliczona ^wartość jest do-

kładną wartością dla wielomianu o nieco zaburzonych współczynnikach.

Rozwa żarny również realizację jednego z tych algorytmów w arytmetyce fl2 oraz

fl

(por. [ 1 O], [ 5] ). Dla arytmetyki fl2 udowadniamy, ^że wszystkie obliczone ^wartości pocho- dnych odpowiadają, w pewnym przybliżeniu, temu samemu wielomianowi o nieco zaburzo- nych współczynnikach. Część pracy poświęcona analizie numerycznej jest polską wersją

pracy [ 11].

2. Wstęp. W literaturze światowej mozna ostatnio zaobserwować wiele prac poświęco­

nych działowi metod numerycznych: „computational complexity theory" (teoria procesów obliczeniowych (? )). Dziahen zajmuje się między innymi poszukiwaniem metod numerycz- nych minimalizujących ilość działań arytmetycznych potrzebnych do rozwiązania danego zagadnienia.

Wiele prac poświęconych jest problemowi obliczania ^wartości wielomianu i jego pochod- nych w możliwie najtańszy sposób (por. [l ], [ 2], [ 8], [ 9]).

I tak np. wiadomo ([ 1], [2] ), ^że algorytm Hornera (por. [ 4), [ 10]) jest jedynym algo- rytmem, który minimalizuje ilość dodawań i ^mnożeń przy obliczaniu ^wartości wielomianu P:

n

(1) P(y)

=Li

k=O

Jeśli jednak zastosujemy algorytm Hornera do obliczania wszystkich „znormalizowanych"

pochodnych p .,(x) (j) wielomianu stopnia n to

łatwo sprawdzić, że wykonamy~

n (n+ I) do-

dawań i ^mnożeń. l· Czy i w tym przypadku algorytm Hornera jest najtańszy?

Odpowiedź negatywną na to pytanie zawiera praca [ 9] , w której autorzy definiują mię-

dzy innymi algorytmy obliczania wszystkich pochodnych wielomianu kosztem

~

n (n+ I) do-

dawań i tylko 3n - 2 mnożeń i dzieleń. Ponadto, jeśli założymy, że stopień wielomianu n jest nieparzysty. to w pracy [ 9] znajdujemy algorytm

wymagający i

^{n (n+ I)}

^dodawań

ⁱ

3n - 3 mnożeń i dzieleń.

(79]

80 H. W o ź n i a k o w s k i

Tak więc widzimy, że w przypadku równoczesnego obliczania wszystkich (lub wielu) pochodnych znamy algorytmy efektywniejsze niż algorytm Hornera, afe czy te nowe algo- rytmy ^sąjuż optymalne? (tzn. czy minimalizują ilość działań?).

W ogólności, problem ten jest otwar~y chociaż w pracy [ 9] autorzy dowodzą, że równo- czesne obliczanie wszystkich wielkości

xl

pV)¹(xr /j! dlaj=O, l, ... ,n wymaga co najmniej 2n - 1 ^mnożeń. A ^więc nie można oczekiwać poważnego zmniejszenia obecnie ^osiąganej

ilości 3n - 2 mnożeń lub dzieleń.

Poszukiwanie metod numerycznych minimalizujących ilość działań potrzebnych do roz-

wiązania danego zagadnienia wydaje ^się naturalne. Równocześnie musimy pamiętać, że na- wet najlepszą metodę możemy w praktyce ^wykonywać tylko w sposób przybliżony, na sku- tek nieuchronnych na ogół błędów zaokrągleń.

Zatem każda metoda przeznaczona do realizacji numerycznej, musi mieć pewną „od-

porność" na błędy zaokrągleń (co nazywamy numeryczną poprawnością lub stabilnością).

W dalszej ^części pracy dowodzimy, że rozważana w pracy [ 9 J jednoparametrowa kla- sa algorytmów obliczania ^wartości wielomianów i ich pochodnych (zawierająca jako przy- padki szczególne algorytm Hornera oraz algorytmy o najmniejszej znanej ilości działań)

jest numerycznie stabilna.

3. Algorytm S-T (Shaw-Traub). Przedstawiamy ^klas~ algorytmów obliczenia wartoś­

ci wielomianu i jego m pierwszych znormalizowanych pochodnych w danym punkcie x, x-:/= O, m ^.r;;;;; n, zdefiniowaną przez M. Shaw i

J

.F. Trauba w pracy [ 9]. Niech (tak jak w (1) ):

n

P(y)

= L

^an-kl

k=O będzie wielomianem stopnia n.

Niech

n+l=p·q

dla pewnych liczb naturalnych p,q, a ^następnie zdefiniujmy dwie funkcje s(j) = (n - j) mod q,

dlaj mod q =f O, dlaj mod q =O, gdzie j mod q oznacza ^resztę z dzielenia j przez q, j=O, l, ... , n:

A 1 go rytm S-T.

(2) i=O,l, ... ,n-1,

(3) Tj=a xs(O)

j

o '

j= 0,1, ... ,m, (4) ^Tj

=

^{Tj- l} + Tj xr(i-j)

i i-1 i-1 , j=O,I, ... ,m, i=j+l,j+2, ... ,n.

Za autorami pracy [ 9 J udowodnimy, ^że i

(5) ^Tj= i x s(i-j) '\' C(k .) L_, ,J ai-kx k--j ' k=j

Algorytmy obliczania ^wartości wielomianów i i'ch pochodnych 81

(6)

Tj=~

_n ^(j) _.,

,jmodą,

1·

gdzie .p(k, j) oznacza współczyrtnik Newtona.

Rzeczywiście. zależność

(5) Jest prawdziwa dla j

=

^-1,

^gdyż

C(-1, -1)

=

1, C(k,-1)

=

O dla k

>

^{-1, a stąd}

i'

T.-1 _{i -} == s(i+l) x L ')' C(k

' -

l) ai-k xk+l

=

xs(i+l) ai+ 1

k=-1 co jest zgodne z (2).

Z kolei, przyjmując i= j, mamy j

~i_ s(O) ~ C(k ') k-j s(O) rj - x L._; ,1 ai'-kx

=

^{x a}₀

k=j co jest zgodne z (3 ).

Stosując indukcję do (4) otrzymujemy

~I ~1

T! = xs(i-j) ~ C(k 1'-l)a. xk+I-j + xr(i-j) + s(i-l-j) ~ C(k J⁰)a. xk-j.

i L..J

'

i-1-k

L...J '

i-1-k

k~-1 k~

Pamiętając, że

(7) r(i-j)

=

s(i-j) - s(i-j-1) + l,

mamy

i-1 i

j _ s(i-j)( ~ { ( · ) (k ')2 k+l-j) _ s(i-j) ~C(k ') k-j Ti - x ai-j + L._; C k,1-l + C ,1 f ai-l-kx - x L.J ,1 ai-kx ,

k~ k~

co kończy dowód zależności (5 ). . .

Przyjmując

i= n i

pamiętając, że

xs(n-1)

=

x1 mod q , ze

związku

(5) otrzymujemy (6 ).

Korzystając ze ^związków (2)-(4) łatwo sprawdzić, że ilość dodawań w algorytmie S-T nie zależy od parametru q i wynosi (m + 1) (n -ł m). Natomiast ilość mnożeń wynosi

{

(p + 1) (ą - 1)

(m + 1) (p - r - 1) +

~

qr (r

+

1)

do zrealizowania (2), (3 ), do zrealizowania ( 4),

gdzie, tak jak poprzednio, n+ 1

=

p · q, natomiast r = entier ^(m/ą ).

Ze

związku

( 6) otrzymujemy

wartości

znormalizowanych pochodnych pU) (x )/j!

przemnożone

przez

xi

mod q .

Chcąc obliczyć wartości

pU) (x )/ j! musimy

wykonać

dodat- kowo m - r dzieleń.

Tak ^więc ogólna ilość mnożeń i ^dzieleń wynosi

m(n+ 1) 1

f

_mn (q) =n - 1 + q + _q - (m + 2)r + -2 q (r + 1 ).

82 H. W o ^ź n i a k o w s k i

W dwóch przypadkach ilość działań może być nieznacznie zmniejszona. Jeśli m

=

n, to p(n) (x )In! = a₀ i

możemy

nie

obliczać

tej

wielkości

z xq-l

p(:~

(x),

wykonując dzięki

temu o jedno dzielenie mniej. Analogicznie, jeśli q =n+ 1, to nie musimy obliczać xn+ 1 i wyko- nujemy w rezultacie o jedno mnożenie mniej.

4. Szczególne przypadki algorytmu S-T. Podstawiając za parametr q różne wartości otrzymujemy warianty algorytmu S-T charakteryzujące się różną ilością mnożeń i ^dzieleń:

I.ą=l:

i=O,l, ... ,n-1,

j =O, l, ...

,n,

j=O,l, ... ,m, i=j+l, ...

,n,

j =O, 1, ... , m.

Otrzymaliśmy

zatem algorytm Hornera, który wymaga wykonania (m

+

1) (n -

~m)

do-

dawań i ^mnożeń. Dla m = O algorytm Hornera jest jednym algorytmem, który minimalizuje

ilość dodawań i mnożeń, jak udowodnił Borodin (por. [ 1], [2) ).

li. q =n+ 1:

-1 n-i-1

Ti

=

ai+ 1 x ' i= O, l, ... ,n-1, j =O, l,

„„

m,

j =O, 1, ... , m, i= j+ l, ... ,n,

j =O, 1, ... , m.

Ten wariant algorytmu S-T jest algorytmem o najmniejszej znanej ilości działań arytme- tycznych i wymaga (m + 1) (n -

~

m)

dodawań,

2n - 1

mnożeń

im

dzieleń.

Ponadto,

jeś­

lim =n, możemy zrealizować ten algorytm kosztem 2n - 1 mnożeń i n - 1 dzieleń.

Pamiętając o twierdzeniu udowodnionym w pracy [ 9 ], a mianowicie, że obliczanie

wi,elkości xi

pU)(x)/j!,j =O, 1, ..

„ n,

wymaga co najmniej 2n - 1

mnożeń,

widzimy,

że

al- gorytm S-T dla q = n + 1 minimalizuje ilość mnożeń potrzebnych do obliczania wielkości

x j pU) ( x) / j ! , j ⁼ O, 1 , ... , n.

III. Niech n

będzie

nieparzyste, m =n, q =

ł

(n+ 1)

Ten wariant algorytmu S-T

wymaga~

n (n + 1)

dodawań

i 3n - 3

mnożeń

i

dzieleń.

IV. m

=

^{1, q}

=..Jn°+I.

Algorytm oblicza wartość P(x) i P'(x) kosztem 2n - 1 dodawań oraz n - 1 + 2 ~

mnożeń i dzieleń. Jeśli n + 1 nie jest kwadratem liczby naturalnej to za q przyjmujemy

c11ticr (~ ). '

Algorytmy obliczania wartości wielomianów i ich pochodnych

Ten algorytm jak również algorytm Munro (por. [ 8 J) są algorytmami o najmniejszej znanej ilości pziałań potrzebnej do obliczania P(x) oraz P'(x).

83

5. Analiza numeryc:.ma algorytmu S-T, realizowanego w arytmetyce fl. Zakładamy , że

przedstawiony algorytm S-T ^będzie realizowany w t cyfrowej zmiennopozycyjnej arytmety·

ce fl z poprawnym zaokrągleniem wyników działań arytmetycznych ([ 10) ).

·Niech rd(x) oznacza zmiennopozycyjną reprezentację rzeczywistej liczby x, tak ^więc rd(x)=x(l-ć:), lt:l~2-t,

a ponadto, jeśli a= rd(a), b = rd(b ), to obliczony wynik (8) fl (a o b) = (a o b) (1 - c),

' gdzie o oznacza dowolny z operatorów arytmetycznych (+, -, *,/).

Zagadnienia dotyczące przekroczenia zakresu licib arytmetyki fl (t~w . nadmiar lub nic domiar) nie będą rozważane w tym artykule. Abv ^uprościć dalsze oszacowania wprowadza- my relację „~" (por. [ 6] ). Niech a(t), b (t) będą funkcjami zmiennej t, O~ to~ t

<

+ 00.

Piszemy a(t)

1

b(t) wtedy i tylko wtedy, gdy istnieje stała K niezależna od t taka, ^że dla każdego! ;;.i: t_{0 ,}

I

a(t) - b(t) I~ 2-t

I

b(t) IK. Ponadto przez a(t) ~ b(t) rozumiemy a(t) ~ b(t) lub a(t)

1

b (t).

Utywając tej relacji możemy w wygodny sposób wyrazić oszacowania błędów typowyr h

złożonych operacji arytmetycznych. Na przykład obliczona w fl ^wartość iloczynu a ₁

*

a₂

* ...

*

^an spełnia zależność (por. [ 10] , str. 28):

n n

(9) fl (

n

_{i= 1}

^ai) =n

_{i= 1} ^{ai (1 -}

^{c ),}

gdzie

(lOa) 1 - c

=n

_i=2 n ⁽¹

^~

^ci)'

Stąd

- (n - l)rt

f

1 - (1 + z-t)n- l ~

c

~ 1 - (1 - 2-t)n-l

1

(n - 1 )2-t, tak ^więc zapiszemy:

(lOb)

I c

I~ (n - 1 )2-t,

l

zamiast relacji

Ie

^I~ (n - l)Z-t¹ (t1 = t - log2 1.06), wprowadzonej przy założeniu (n - l)Z-t

<

0.1 w [ 10], str'. 31-32.

Przejdźmy do analizy numerycznej algorytmu S-T. Bez zmniejszania ogólności, możem ·;

założyć, że współczynniki wielomianu ai oraz argument x są dane w fl: • ai = rd(a;), x = rd(x ).

Udowodnimv

TWIERDZENIE 1 (fi).

jeśli r!są

_I obliczonymi

wielkościami r/

_I w arytmetyce fl, to (11)

r/

^{= xs(i-j)}

L,

i C(k,j)ai-k (1

-11~;

)xk-i, k""'j

84 H. W o ź n i a k o w s k i

przy czym

(12)

177-:-

_z,-¹₁

1

~s(i+I)2-t, ₁ i= 0,1, ... ,n-l,

(13) j=O,I, ... ,m,

(14)

I

¹¹

~k I ~ {

_{z 1} 2·k + 1 -

o

_t ^{.k - j + s (i - j) } 2-t,}

j

=

O, 1, ... , m, i·

=

^j + 1, j + 2, ... , n, k

=

^{j, j} + 1, ..

„

i, gdzie

o

^ik oznacza symbol Kroneckera.

Do wód. Zacznijmy od zależności (2), (3). Z zależności (9), (IO) wynika

"'r-1

_i ⁼ _ai+Ix ^s(i+l)

(1 -

_{11i,-I '} ^{-I )}

'T!

⁼ aoxs(O) ( 1 -

77~.),

J }}

co jest równoważne z \ (11) i (12), (13).

Co więcej, możemy przyjąć

(15) dla k

>-I.

Obecnie zanalizujemy numeryczną realizację zależności (4). Z (8), (9), {IO) otrzymujemy

(16) T"'j _ ("'Tj-I "'j r(i-j) _{i - i-I + Ti-I} (l _ K 1.·) ) (I _ ^c-.i),

X t V I

gdzie

I

^K _z

^/I~

₁ ^{r(i - j)Z-t,}

Dowód zależności (14) przeprcwadzimy indukcyjnie ze ^względu na z·= O, l, ... ,n. Przy- padek z'= O został już sprawdzony. ^Załóżmy zatem, że zależności (11)-(14) są prawdziwe

"'j-1 -i

dla i - 1, ^i~ 1. Podstawiając (11) w miejsce Ti-I-, Ti-I otrzymujemy

Tf

_z

^{= (}

^xs(i-j)

L..t

i-I

'°"'

C(k,j-I)a. ^1-I-k ( 1 - 17f-I ) xk-j+I ^{i-I ,k} + k=j-I

i-I

+ xs(i-l-j)

L

^C(k,j} ai-l--k ( 1 -

11{_

₁,k) xk-j+r(i-j) ( 1 -

K/))

^{(1 -}

e,h

k=j

Zdefiniujmy

wielkości 'P/~

,

iJ; 1 ~ następująco

(17)

(

1 - 'Pili -^{j - (} 1 -

77

^j-1 i-l,k-1) 1 -⁽

ei ·

^j)'

1 -

iJ; ~ = (

1 - T/ (_ l ,k -l ) ( 1 - K

! ) ^{(1 -}

^{t; /).}

Pamiętając, że

(18)

otrzymujemy

Algorytmy obliczania wartości wielomianów i ich pochodn)'ch

r(i - j) = s(i - j) - s(i - j - 1) + 1

"'. (. ·)

~

^{k .(} C(k-1,j-l)ipk

+C(k-1,j)iJ;~).

r/ ⁼

^{x5 i-1}

^L.J

^{C(k, j) ai-kx -1 1 -}

k=j ' C(k,j)

Porównując z (11) widzimy, ^że

. C(k-1,j-l)ipk + C(k-1,j)iJJ,{

TJ} _ik = - - - -C(k,j) (19)

85

Zauważmy, że Tik

^jest

średnią ważoną ip 1 ^~

^{oraz iJ;}

k .

Z (16) i (17) otrzymujemy zatem za·

leżność rekurencyjną

I

TJk

I~

max

(I ip 1 ^{~ I,}

^{(1 - ó kj)}

^{I I/I~ I)} 1f

f2-t

+max (1.,,t/,k-ll, (1-ók)

l(77/_

₁,k-ll+r(i-j)2-t)) dla j = O, 1,

.„,

min (i - l, m); k = j, j + l, ... ,i.

Z założenia indukcyjnego oraz z (12), (13), (15) i (18) wynika

. t (

117/k I f

2- +max (2k - óik - j + s(i - j))2-t,

(2k - I - Ó ik - j + s(i - I - j) + r(i - j) )2-t)

=

(2k + l -

o

^{ik - j} + s(i - j) )2-t

co ^kończv dowód. .

Tak

~ięc udowodniliśmy, że każda

.obliczona

wielkość 1'f

^jest

dokładną wielkością

dla nieznacznie zaburzonego wielomianu

'P!:

'

_n

(20) ""i( ) - ~ pi Y - ~ an-k ( 1 -1]ik Y · j) k k=O

Zaburzenia

względne 77 1 ^{~ są}

^na

ogół różne

dla

różnych

_{i,j, ale}

posiadają

wspólne oszacowa- nie zgodnie z twierdzeniem 1.

Zauważmy, że

I

17!k. _I

I~

₁ (2k + 1 - Ó _ł 'k - j + s(i - j))2-t

~

₁ (2k + 1 - f> 'k +n - i)2-t _I

~

₁ 2n 2-t, a ponadto najlepsze oszacowania dla

77~

otrzymujemy dla q

=

1, najgorsze dla q

=

n + 1.

Oszacowania

T/k

dane przez (12), (13), (14)

zależą

liniowo od n.

Pamiętając, że błędy

reprezentacji współczynników ai oraz argumentu x powodują zaburzenia wyjściowego \Vic-

lomianu P do postaci

~

n an-k (1 - Ek)xk k=O

86 H. W o ź n i a k o w s k i

z najlepszymi oszacowaniami a priori dla zaburzeń

I

Ek ^I~ (k + 1 )2-t ~ (n+ 1 )rt,

I 1

dochodzimy do wniosku, ^że algorytm S-T, dla dowolnej ^wartości q, jest algorytmem ^świet­

nie numerycznie stabilnym.

Dla i

=

n otrzymujemy WNIOSEK.

"'{j)

"'i P (x) j modą

Tn

=-.,- ..

_J· ^{x j}

⁼

0,1, ... ,m, gdzie

P~

jest dany

zależnością

(20) (por. (6)).

Przedyskutujemy twierdzenie 1 dla szczególnych ^wartości q (por. §4).

I. q = 1: Dla tej ^wartości q algorytm S-T sprowadza ^się do algorytmu Hornera. ^Ponieważ s(i) =O, więc

(21)

111ik1~(2k+l-j-8

_{n 1 n} k)rt.

Dlaj =O powyższe osz;icowanie możemy znaleźć w [4] i [10], str. 55.

Il. q =n + 1: Teraz s(i) =n· - i; oszacowania na

77~k są następujące

(22)

I

11 ik _n

I ~ (

₁ 2 k + 1 - 8 k ) 2-t. _n

Zauważmy, że dla j

=

O otrzymaliśmy dokładnie to samo oszacowanie co w algorytmie Hor- nera (por. (21)), natomiast dla j

>

O oszacowanie (22) jest nieznacznie gorsze od (21 ).

I .

III. Niech n będzie nieparzyste, q

= 2

(n + 1). Wielkości 71~k mają teraz oszacowanie jeślij

<

^q,

jeślij ~ q.

6. Zastosowa'łlie arytmetyki

fh

w algorytmie S-T dla q

=

ⁿ + 1. Algorytm S-T w przy- padku q =n + 1 ^{może być} realizowany z wykorzystaniem arytmetyki fl2 (por. [ 1 O], str.55 ):

(23) T_i ^{-1 . _} _{. - aj+ I X} ^n-i-I (fl)' i=O,l, ... ,n-1,

(24) (fl)' ^j

=

0,1, ... ,m,

(25) j . - j-I j

Ti . - Ti-I + Ti-I j

=

0,1, ... ,m, i= j + 1, ... ,n.

Z (25) i (5) wynika

(26)

Algorytmy obliczania wartości wielomianów i ich pochodnych 87

k

(27) Tkj-1 = L.J ~ C(p, 1-l ak-px . ) n-k-p p=j-1

dla wszystkich rozważanych i, j, k.

Jem obliczana suma

J:

n· ^{ai w fl}2 jest zapamiętana z mantysą w podwójnej precyzji, to

~l .

n n

(28) fl2 (

L,

^a

i) ⁼ L,

^{a i ( 1 - Di),}

i= 1 i= 1 gdzie

1.5.I ~~(n+

_{z 1} 1 - i)r2'

dla dowolnych liczb rzeczywistych ai = rd(ai) (por. [ 10], str. 37). Udowodnimy

TWIERDZENIE 2 (fl2 ).

jeśli 1/ ^są

obliczonymi

wielkościam.i

T/ przy

użyciu

(23 ), (24) w fl oraz (25) w f1_{2 ,} to

(29)

gdzie (30)

1'/

⁼

L,

i ^{C(k, j)ai-k (1 -} "1,·-k) (1 -

ri,{

)xn-i+k, k=j

lri.

_1-^kl~ ₁ (n-i+k)2-t,

I

_"lik ^{j '~}

l 2

³ (1 - uj-1 ^{t: )} k . 2 ^-2t dla wszystkich rozważanych i, j, k (por. twierdzenie 1 ).

Do wód. Obliczone z wzorów (23) i (24)

wartości T.-

^l,

f! spełpiają z~teżności:

I J ··'

"'-l_fl( n-i-1)- n-i-1( )

Ti - ai+lx -ai+lx 1- 77i+l, i=O,l, ... ,n-1,

j=O,l, ... ,m, gdzie

k = 0,1, ... ,n, na mocy (9), (10), co jest równoważne związkom (29), (30).

Wartości 'f/ ^są

obliczane i

zapamiętywane

w podwójnej precyzji_, tak

więc

z (25), (26) wynika

(31) i-1

"'i -

( )'

"'j-1) Ti - fl::t L; Tk

i-1

)' T"'kj-1 ( 1 - " j )

L "ik ¹

k=j-1 k=j-I

88 H. W o ^ź n i a k o w s k i

gdzie

na moc 'f (28 ).

Dowód twierdzenia 2 jest indukcyjny ze ^względu na j. Przypadek j = -1 jest ^już spraw- dzony. Z

założenia

indukcyjnego

możemy zastąpić Tj-l

^{w (31)}

wyrażeniem

(29) (por. (27)):

i-1 k

r( ⁼ L L

C(p, j-l)ak-p (1 - Tlk-p) (1

-111;

¹)xn·-k-P(1 -

e~).

k=j-1 p=j-1

Przyjmując r

=

i + p - k i zmieniając porządek sumowania mamy

i ~l

r!

i = "\'a. .L...J r=j i-r (1 - Tl. )xn-i+r i-r . k=i-

L

l +j-r ^{C(k+r-i, j-1) (1 -} Tlkj-kl+ - .) , r ⁱ (1 -

e

ⁱ !k. ) •

Dalej,

przyjmując

p = k + r - i oraz

wykorzystując postać r/ ^daną

przez (29), otrzymujemy

~

k-1 C(p, j-1 )(

rit~,~-k,p

+ s{p+i-k (1 -

rit~:=k,p))

j p=j-1

1 - Tlik = 1 - - - -_C(k, j) (32)

na mocy tożsamości

k-1

C(k, j) =

L

^{C(p, j-1 ).}

p=j-1,

Tak

więE:, Tl~

, podobnie jak w twierdzeniu l, jest

średnią ważoną

poprzednich

zaburzeń.

Stąd wynika oszacowanie (por. ( 31))

(33)

I

_{Tl·k ~. ł}

i I

~ _{1 ,_ 1 <;;} max _{p <: k-1}

(

1-r11p+'-k p + -

i-

_ł 1 _'

I

_' 3 2 ^{(k - p} ) - 2 2

t)

.

Ponieważ Tli~

1 =O,

więc łatwo sprawdzić, że rozwiązanie

(33) dlaj

~O

jest dane nierów-

nością

co ^kończy dowód twierdzenia 2.

Z twierdzenia 2 wynika, ^że dla

na mocy (29) mamy

]'i

_{i - ~} ^~

'\'

i C(k ·) "' '1 ai-kx n-i+k ' k=j

Algorytmy obliczania wartości wielomianów i ich pochodnych 89

co oznacza,

że

obliczone z wykorzystaniem fl2

wartości r/

^{{dane z}

^podwójną

^pre-

cyzją) są, praktycznie biorąc, dokładnymi wartościami zdefiniowanymi w algorytmie S-T

(ą

= n + 1) dla tego samego wielomianu

L

n ^{an-k tk' o}

współc~ynnikach

nieco zabu-

k =O rzonych w stosunku do wielomianu wyjściowego.

Ponieważ końcowe

wyniki

T~+

1 zazwyczaj

zaokrąglamy

do pojedynczej precyzji, za- tem, w pewnym przybliżeniu, możemy powiedzieć, że spełniona jest formuła Kahana ([ 3] ), tzn.: „otrzymane wyniki są nieco zaburzonym dokładnym rozwiązaniem zadania, o nieco zaburzonych danych". Ponadto widzimy,

że

oszacowania

zaburzeń względnych fi,~ są

mniejsze w arytmetyce fl2 . ^Mnożnik (n - 1 + 2k + 1 · - Ó ik) występujący w twierdzeniu l, jest teraz zastąpiony odpowiednio przez (n - i+ k ).

7. Zastosowanie arytmetyki

fi

^w algorytmie S-T dla q

=

n + 1. Omówimy teraz algorytm S-T dla q =n+ 1 w arytmetyce

fl.

Arytmetyka

fi

jest standardową arytmetyką fl z użyciem algorytmu Mr/>llera (por. [5], [7]) dodawania n liczb.

W przypadku q = n + 1 możemy wykorzystać M~llerowskie sumowanie w algorytmie S-T w sposób następujący:

ST 1 : for j : = O step 1 until m do begin

p:

=O;

s ·=Tj· . j'

(34) for i :

=

j + 1 step 1 until n do begin

a . = Ti-I.

. i-1 ' z : = s +a;

p : =

s - z +a+

p;

s : =z;

T! :

_J

=

s + p end end;

Udowodnimy

TWIERDZENIE 3 (fl).

Jeśli r! są, wielkościami

obliczonymi algorytmem ST 1 (34),to i

'

(35)

"'i_

Ti - L,_; ~ C k, J ai-k

( .) (

1 - f'l;1i i ) x n-le-i , k=j

gdzie

(36)

I

Tl !k.

I ~

(n - i + k + 2j + 2)

r

^t

' 1

dla wszystkich rozważanych i, j, k.

"' 1 ...

Do wód. Wartości T.- , Tl są obliczane w fl, tak więc

i J

T"'-1 _ fl ( n-i-1) _ n-i-I ( -1 )

i - ai+ 1 ^{x - a;+} 1 ^x 1 - fi;.-i , i=O,l„ .. ,n-1,

90

gdzie

H. W o ź n i a k o w s k i

j= 0,1,„.,m,

I

11.-¹₁

I ~

(n - i - l)r ^t,

'·=

¹

111!. _JJ

i ~n

₁ .rt,

na mocy (9), (10), co jest równoważne (35), (36). Z pracy [4] wynika

i-1 i-1

r/ ^=fl (z

_k=j-1

1t1) ⁼

_k=j-1

z ^{rtl{l - c/,),}

gdzie

I i I

c.k _i ^~2'2. ₁ ^-t

Przeprowadzając analogiczne rozumowanie jak w twierdzeniu 2, otrzymujemy ^wyrażenie na

T/~

dane przez (32), a

stąd

I

'f/.k ~ ^j

I

2 · 2- + t max

I

T/p+{-k,p j-

l ^I

•

I 1 j-I ~ p <; k-1 Zauważmy, że z (37) wynika

IT/.k

⁰

1~

(n-i+k+2)rt,

I 1

a stosując indukcję ze względu naj, łatwo otrzymujemy (36).

Twierdzenia 1,2 ,3 w przypadku i= n dają następujące oszacowania na ^względne zabu- rzenie współczynników ai-k:

I

T/!k.

I~

(2k + 1 - ó k)2-t, (fl)

I t n

I 1 - (1 -T/i-k) {l -

'f/~~)lf

k · rt, (fl_{2 )}

i

TJ _n ik

I~

₁ (k + 2j + 2)rt, (fl) k = j, j + 1 , . „, n.

Zauważmy, że

otrzymane oszacowania T/!k w

fi są

lepsze od

oszacowań

w fl jedynie dla małych wartościj, a oszacowania w fl₂ i

fi

^są praktycznie te same.

8. Zastosowania algorytmu S-T. Algorytm S-T obliczania znormalizowanych pochod- nych wielomianu w danym punkcie może być używany:

1° do rozwinięć potęgowych wielomianu. Mianowicie, jeśli n

P(y)

= L

^an-kl

k=O chcemy rozwinąć w punkcie x:

Algorytmy obliczania wartości wielomianów i"ch pochodnych

L

n ^{bn-k(y - x)k,} k=O

to szukane współczynniki bn-k są równe

b n-k

=

p(k)(x)/k!

=

Tk /xk n mod q dla k = O, l, ... ,n.

91

Współczynniki bn-k mogą być obliczone algorytmem S-T. Z przedstawionych twier-

dzeń wynika, że każdy obliczony współczynnik bn-k jest ^dokładny dla nieco zaburzone- go wielomianu.

Niech

P

^będzie wielomianem o obliczonych współczynnikach:

"' )'"' n k P(y)

=

L; bn-k (y - x) .

k=O

Porównajmy wartości wielomianów

P

^{oraz P.} Łatwo sprawdzić, że n

P(y) =

L

^{an-k(l - sk)l,}

k=O gdzie

(

k { 2k + 1 - ó

Is ^l~rt

lxl+ly-xl) X k nk

k i ly I

3k + 2

w fl, q dowolne, w fl_{2 ,} q = n + 1, w

ff,

q =n+ 1.

Tak więc, jeśli założymy, że

(I

x

I+ I

y - x

I)/ I

y

I

jest rzędu jedności, to wartości wielo- mianu

P

są dokładne dla wielomianu o nieco zaburzonych współczynnikach.

Zauważmy, że chociaż zaburzenie dla każdego

bn-k

^są na ogół różne, to dla ustalonego y potrafimy udowodnić, że wartości

P

^pochodzą od jednego, nieco zaburzonego wielomianu.

2° W wielu metodach iteracyjnych rozwiązywania równania P(y)

=

O, wykorzystujemy

wielkości

p(i) (x) /j!.

I tak np. w metodzie Newtona konstruujemy ciąg

P(xi)

X =x - - - i+l i p'(x.)'

i

gdzie P(xi) p' (x.)

'

Literatura

To n Tl

I

^xl mod q

n

[ l] A. B or od i n, Homer's Rule is Uniquely Optimal, w książce Theory of Machines and Computa- tions, edited by Kohavi and Paz, 1971, str. 45-58.

[2] - Computational Complexity - Theory and Practice, to apper in Currents in the Theory of Compu- ting, edited by A. Abo, Prentice - Hall, 1972. .

[3] W. K ah an, A survey oferroranalysis, IFIPCongress 71.

[ 4] A. K ie ł b as i ń s ki, Numeryczne obliczanie wartości wielomianu algorytmem Homera, Sprawoz- dania IMM i ZON UW, zeszyt 16, 1968.

[ 5] - Algorytm sumowania z poprawkami i niektóre jego zastosowania, Matematyka Stosowana 1 (1973), str. 23-41.

92 H. W o ź n i a k o w s k i

[6] A.Kiełbas iński, Analiza numeryczna algorytmu ortogonalizacji Grama-Schmz"dta, Matematyka Stosowana 2(1974), str. 15-35.

[7] O. M ~ 11 er, On a quasi double-precision in floating-poz"nt addition, BIT 5(1965), str. 37-50, 251-255.

[8) J.I. Mu nr o, Some results z"n the Study of Algorithms, Technical Report 32, Department of Com- puter Science, University of Toronto, 1971.

[ 9] Mary S h a w and J .F. T r a u b, On the number of multiplications for the evaluatz"on of a poly- nomial and some of its derivatives, Department of Computer Science, Carnegie-Mellon-University, 1972,J.S.C.M. Vol. 21, No. l,jaIIUary 1974, str. 161-167.

[ 1 O] J .H. W i 1 k i n s o n, Błędy zaokrągleń w procesach algebraicznych, Warszawa 196 7.

[ 11] H. W o ź n i a k o w s k i, Rounding error analysis for the evaluation of a polynomz"al and some of its derivatives. SIAMJ. Numer. Anal. Vol. 11, No. 4, September 1974.