Sekos riba

3. Monotoninės sekos

1. Monotoninių sekų apibrėžimas.

Apibrėžimas. Seka $\{x_{n}\}$ vadinama nemažėjančia (nedidėjančia), jei kiekvienas tos sekos elementas ne mažesnis (ne didesnis) už pirmesnįjį elementą, t. y., jei su visais numeriais n teisinga nelygybė

$x_{n}\leq x_{n+1}\;\;$ ($x_{n}\geq x_{n+1}$).

Nemažėjančios ir nedidėjančios sekos vadinamos monotoninėmis sekomis. Jei monotoninės sekos $\{x_{n}\}$ elementai su bet kuriuo numeriu n tenkina nelygybę $x_{n}<x_{n+1}\;\;$ ($x_{n}>x_{n+1}$), tai seka $\{x_{n}\}$ vadinama didėjančia (mažėjančia). Didėjančios ir mažėjančios sekos dar vadinamos griežtai monotoninėmis sekomis.

Monotoninės sekos yra aprėžtos arba iš viršaus, arba iš apačios; būtent, nedidėjančios sekos aprėžtos iš viršaus, o nemažėjančios aprėžtos iš apačios pirmaisiais savo elementais. Todėl nedidėjanti seka yra aprėžta iš abiejų pusių, kai ji aprėžta iš apačios, o nemažėjanti seka yra aprėžta iš abiejų pusių, kai ji aprėžta iš viršaus.

Pateiksime tris monotoninių sekų pavyzdžius.

1. Seka $1,\;1,\;{\frac {1}{2}},\;{\frac {1}{2}},\;...,\;{\frac {1}{n}},\;{\frac {1}{n}},\;...\;$ nedidėjanti. Ji aprėžta iš viršaus savo pirmuoju elementu, lygiu vienetui, o iš apačios - skaičiumi 0.

2. Seka $1,\;1,\;2,\;2,\;...,\;n,\;n,\;...\;$ nemažėjanti. Ji aprėžta iš apačios savo pirmuoju elementu, o iš viršaus neaprėžta.

3. Seka ${\frac {1}{2}},\;{\frac {2}{3}},\;{\frac {3}{4}},\;...,\;{\frac {n}{n+1}},\;...\;$ didėjanti. Ji aprėžta iš abiejų pusių: iš apačios pirmuoju elementu ${\frac {1}{2}},$ iš viršaus, pavyzdžiui, skaičiumi 1.

2. Monotoninės sekos konvergavimo požymis.

Įrodysime pagrindinę teoremą.

3.15 teorema. Jei nemažėjanti (nedidėjanti) seka $\{x_{n}\}$ aprėžta iš viršaus (iš apačios), tai ji konverguoja.

Pagal praeitą skirsnį seka $\{x_{n}\},$ atitinkanti 3.15 teoremos sąlygas, yra aprėžta. Todėl 3.15 teoremą galima trumpai formuluoti taip: jei monotoninė seka $\{x_{n}\}$ yra aprėžta iš abiejų pusių, tai ji konverguoja.

Įrodymas. Kadangi seka $\{x_{n}\}$ yra aprėžta, tai jos elementų aibė turi tikslųjį viršutinį rėžį ${\overline {x}}$ ir tikslųjį apatinį rėžį ${\underline {x}}.$ Įrodysime štai ką: jei $\{x_{n}\}$ – nemažėjanti seka, tai jos riba yra minėtasis tikslusis viršutinis rėžis; jei $\{x_{n}\}$ – nedidėjanti seka, tai jos riba yra minėtasis tikslusis apatinis rėžis ${\underline {x}}.$ Išnagrinėsime tik nemažėjančią seką, nes samprotavimai apie nedidėjančią seką būtų analogiški.

Kadangi ${\overline {x}}$ yra sekos $\{x_{n}\}$ elementų aibės tikslusis viršutinis rėžis, tai kiekvieną $\varepsilon >0$ atitinka toks elementas $x_{N},$ kad $x_{N}>{\overline {x}}-\varepsilon$ ir $x_{N}\leq {\overline {x}}$ (bet koks elementas $x_{n}$ ne didesnis už tikslųjį viršutinį rėžį ${\overline {x}},$ t. y. $x_{n}\leq {\overline {x}}$). Iš šitų nelygybių gauname nelygybes $0\leq {\overline {x}}-x_{N}<\varepsilon .$ Kadangi $\{x_{n}\}$ – nemažėjanti seka, tai nelygybės $x_{N}\leq x_{n}\leq {\overline {x}}$ yra teisingos, kai $n\geq N.$ Iš to išplaukia, kad $0\leq {\overline {x}}-x_{n}\leq {\overline {x}}-x_{N},\;$ kai $n\geq N.$ Anksčiau pabrėžėme, kad ${\overline {x}}-x_{N}<\varepsilon ,$ todėl, kai $n\geq N,$ teisingos nelygybės $0\leq {\overline {x}}-x_{n}<\varepsilon ,$ iš kurių gauname nelygybę $|x_{n}-{\overline {x}}|<\varepsilon .$ Vadinasi, įsitikinome, kad ${\overline {x}}$ yra sekos $\{x_{n}\}$ riba. Teorema įrodyta.

1 pastaba. Monotoninės sekos aprėžtumas yra būtinas ir pakankamas jos konvergavimo požymis.

Iš tikrųjų, jei monotoninė seka yra aprėžta, tai pagal 3.15 teoremą ji konverguoja; jei monotoninė seka konverguoja, tai, remiantis 3.8 teorema, ji yra aprėžta.

2 pastaba. Konverguojanti seka gali ir nebūti monotoninė. Pavyzdžiui, seka $\{x_{n}\},$ apibrėžta formule $x_{n}={\frac {(-1)^{n}}{n}},$ konverguoja, o jos riba yra skaičius 0. Tačiau ši seka nėra monotoninė, nes po teigiamo elemento eina neigiamas, o po neigiamo – teigiamas.

3 pastaba. Jei $\{x_{n}\}$ yra nemažėjanti ir aprėžta seka, o ${\overline {x}}$ – tos sekos riba, tai visi jos elementai tenkina nelygybę $x_{n}\leq {\overline {x}}.$ Nedidėjančios ir aprėžtos sekos $\{x_{n}\},$ konverguojančios į ${\underline {x}},$ elementai tenkina nelygybę ${\underline {x}}\leq x_{n}.$ Kad šie teiginiai teisingi, įsitikinome, įrodinėdami 3.15 teoremą.

3.15 teoremos išvada. Sakykime, duota begalinė segmentų sistema $[a_{1};\;b_{1}],\;[a_{2};\;b_{2}],\;...,\;[a_{n};\;b_{n}],\;...\;$ Jei kiekvienas segmentas yra pirmesniajame*, o skirtumas $b_{n}-a_{n}$ (jį vadinsime segmento $[a_{n};\;b_{n}]$ ilgiu) artėja prie nulio, kai $n\to \infty$ (segmentų sistemą, turinčią šią savybę, vadinsime susitraukiančia), tai egzistuoja vienintelis taškas c, priklausantis visiems tos sistemos segmentams.

Įrodymas. Pirmiausia pastebėsime, kad taškas c, priklausantis visiems segmentams, gali būti tik vienas. Iš tikrųjų, jei būtų dar vienas taškas d, priklausantis visiems segmentams, tai segmentas ** [c; d] priklausytų visiems segmentams $[a_{n};\;b_{n}].$ Bet tokiu atveju su bet kokiu numeriu n būtų teisingos nelygybės $b_{n}-a_{n}\geq d-c>0;$ tai neįmanoma, nes $b_{n}-a_{n}\to 0,$ kai $n\to \infty .$ Dabar įsitikinsime, kad taškas c, priklausantis visiems segmentams $[a_{n};\;b_{n}],$ egzistuoja. Kadangi segmentų sistema yra susitraukianti, tai kairiųjų galų seka $\{a_{n}\}$ yra nemažėjanti, o dešiniųjų galų seka $\{b_{n}\},$ – nedidėjanti. Kadangi abi tos sekos yra aprėžtos (visi sekų $\{a_{n}\}$ ir $\{b_{n}\}$ elementai priklauso segmentui $[a_{1};\;b_{1}]$), tai pagal 3.15 teoremą jos konverguoja. Iš to, kad skirtumas $b_{n}-a_{n}$ nyksta, išplaukia, jog minimosios sekos turi bendrą ribą. Tą ribą žymėkime raide c. Remiantis 3 pastaba, su bet kokiu numeriu n teisingos šios nelygybės: $a_{n}\leq c\leq b_{n},$ t. y. taškas c priklauso visiems segmentams $[a_{n};\;b_{n}].$

_______________

* Tai reiškia, kad $a_{n-1}\leq a_{n}<b_{n}\leq b_{n-1}.$

** Siekdami konkretumo, tariame, kad d > c.

3. Keli konverguojančių monotoninių sekų pavyzdžiai.

Išnagrinėsime kelias sekas, kurių ribas apskaičiuosime, remdamiesi 3.15 teorema apie monotoninės sekos ribą. Be to, šiame skirsnyje susipažinsime su vienu bendru sekos ribos ieškojimo metodu, kai seka apibrėžiama rekurentine formule (Rekurentinė formulė (lot. recurrens – grįžtas) – formulė, pagal kurią $(n+1)$-ąjį sekos elementą galima išreikšti jos pirmųjų n elementų reikšmėmis.).

1 pavyzdys. Tirsime seką $\{x_{n}\},$ kurios elementas $x_{n}$ lygus

$x_{n}={\sqrt {a+{\sqrt {a+{\sqrt {a+...+{\sqrt {a}}}}}}}},\;\;a>0$

(n radikalų). Tą pačią seką galima, savaime aišku, nusakyti šitokia rekurentine formule:

$x_{1}={\sqrt {a}},\quad x_{n}={\sqrt {a+x_{n-1}}}.$

Norėdami įsitikinti, kad sekos $\{x_{n}\}$ riba egzistuoja, įrodysime, jog ta seka yra didėjanti ir aprėžta. Kad ta seka didėja, matyti tiesiog. Įsitikinsime, kad seka $\{x_{n}\}$ aprėžta iš viršaus skaičiumi A, kai A – didedsnysis iš dviejų skaičių a ir 2. Jei $x_{n}\leq a,$ tai teiginys jau įrodytas. Jei $x_{n}>a,$ tai, dešinėje nelygybės $x_{n}^{2}=a+x_{n-1}\leq a+x_{n}\;$ pusėje skaičių a pakeitę didesniu už jį skaičiumi $x_{n},$ gauname $x_{n}^{2}<2x_{n},$ o iš čia $x_{n}<2.$ Taigi įrodėme, kad seka $\{x_{n}\}$ aprėžta iš viršaus. Pagal 3.15 teoremą ta seka turi ribą. Pažymėkime tą ribą raide c. Savaime aišku, $c>0.$ Iš rekurentinės formulės turime lygybę

$x_{n}^{2}=a+x_{n-1},$

iš kurios matyti, kad sekos $\{x_{n}^{2}\}$ ir $\{a+x_{n-1}\}$ sutampa. Todėl jos turi bendrą ribą. Kadangi pirmosios sekos riba lygi $c^{2},$ o antrosios $a+c$ (kai n varo į begalybę skaičiams $x_{n}$ ir $x_{n-1},$ tai $x_{n}$ praktiškai nesiskiria nuo $x_{n-1}$ ir abu yra labai arti c reikšmės (abu labai panašūs į c)), tai $c^{2}=a+c.$ Iš čia, turėdami mintyje, kad $c>0,$ randame c:

$c^{2}=a+c,$

$c^{2}-c-a=0,$

$c_{1,2}={\frac {-(-1)\pm {\sqrt {(-1)^{2}-4\cdot (-a)}}}{2}}={\frac {1\pm {\sqrt {1+4a}}}{2}};$

kadangi $c>0,$ tai

$c={\frac {1+{\sqrt {1+4a}}}{2}}.$

2 pavyzdys. Tirsime seką $\{x_{n}\},$ kuria dažniausiai naudojamasi, apskaičiuojant teigiamo skaičiaus a kvadratinę šankį elektronine skaičiavimo mašina. Šita seka apibrėžiama tokia rekurentine formule:

$x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )},\quad n=1,\;2,\;3,\;...$

Pirmuoju elementu $x_{1}$ čia galima laikyti bet kokį teigiamą skaičių.

Įrodysime, kad ta seka konverguoja, o jos riba yra skaičius ${\sqrt {a}}.$ Pirmiausia įsitikinsime, kad seka $\{x_{n}\}$ turi ribą. Tam užtenka įrodyti, kad seka $\{x_{n}\}$ yra aprėžta iš apačios ir kad, pradedant antruoju elementu, ji nedidėja. Iš pradžių įsitikinsime, kad seka $\{x_{n}\}$ aprėžta iš apačios. Sąlygoje pasakyta, jog $x_{1}>0,$ o iš rekurentinės formulės, kai $n=1,$ gauname $x_{2}>0.$ Taip samprotaudami toliau, įsitikiname, kad visi $x_{n}$ yra teigiami.

Dabar įrodysime, kad visi $x_{n},$ kai $n\geq 2,$ tenkina nelygybę $x_{n}\geq {\sqrt {a}}.$ Parašę rekurentinę formulę šitaip: $x_{n+1}={\frac {\sqrt {a}}{2}}{\Big (}{\frac {x_{n}}{\sqrt {a}}}+{\frac {\sqrt {a}}{x_{n}}}{\Big )},$ remsimės beveik savaime aiškia nelygybe $t+{\frac {1}{t}}\geq 2$*, kuri teisinga, kai $t>0$ (imame $t={\frac {x_{n}}{\sqrt {a}}}$). Gausime

$x_{n+1}={\frac {\sqrt {a}}{2}}{\Big (}{\frac {x_{n}}{\sqrt {a}}}+{\frac {\sqrt {a}}{x_{n}}}{\Big )}\geq {\frac {\sqrt {a}}{2}}\cdot 2={\sqrt {a}},$

kai $n\geq 1,$ t. y. $x_{n}\geq {\sqrt {a}},$ pradedant numeriu $n=2.$

Pagaliau įsitikinsime, kad seka $\{x_{n}\},$ kai $n\geq 2,$ nedidėja. Iš rekurentinės formulės matome, kad ${\frac {x_{n+1}}{x_{n}}}={\frac {1}{2}}{\Big (}1+{\frac {a}{x_{n}^{2}}}{\Big )},$ o iš čia, turėdami mintyje, kad $x_{n}\geq {\sqrt {a}},$ gauname ${\frac {x_{n+1}}{x_{n}}}\leq 1,$ arba $x_{n}\geq x_{n+1}$ (kai $n\geq 2$).

Kadangi seka $\{x_{n}\},$ kai $n\geq 2,$ nedidėja ir yra aprėžta iš apačios skaičiumi ${\sqrt {a}}$ (nes $x_{n}\geq {\sqrt {a}}$), tai ji turi ribą, ne mažesnę kaip ${\sqrt {a}}$ (žr. 3.15 ir 3.13 teoremą). Tą ribą pažymėkime raide c ir atsižvelgę į tai, kad $\lim _{n\to \infty }x_{n+1}=c,$ o $\lim _{n\to \infty }{\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )}={\frac {1}{2}}{\Big (}c+{\frac {a}{c}}{\Big )},$ gauname lygybę

$c={\frac {1}{2}}{\Big (}c+{\frac {a}{c}}{\Big )}$**,

$c^{2}={\frac {1}{2}}{\Big (}c^{2}+a{\Big )},$

$2c^{2}=c^{2}+a,$

$2c^{2}-c^{2}=a,$

$c^{2}=a.$

Vadinasi, $c={\sqrt {a}}.$

1 pastaba. Spręsdami pateiktuosius pavyzdžius, naudojomės plačiai taikomu metodu sekos ribai skaičiuoti. Metodo esmė: iš pradžių įsitikinama, kad sekos riba egzistuoja, o paskui iš lygties, kuri gaunama iš rekurentinės formulės, vietoje $x_{n}$ ir $x_{n+1}$ įrašius ieškomąją sekos $\{x_{n}\}$ ribos reikšmę c, apskaičiuojama tos ribos skaitinė reikšmė.

2 pastaba. Rekurentinės formulės dažnai naudojamos skaičiavimo matematikoje, nes jas taikant daug kartų, atliekamos to paties tipo skaičiavimo operacijos, o tai labai patogu skaičiuojant elektroninėmis skaičiavimo mašinomis.

Išnagrinėtoji rekurentinė formulė aprašo, kaip įsitikinome, ${\sqrt {a}}$ skaičiavimo algoritmą (įrodėme, kad $\lim _{n\to \infty }x_{n}={\sqrt {a}}$).

Vėliau bus tiriamas sekos $\{x_{n}\}$ konvergavimo į ${\sqrt {a}}$ greitis. Įrodoma, kad, tinkamai pasirinkus pirmąjį artinį $x_{1},$ kai $a>1,$ jau ketvirtasis artinys $x_{4}$ nuo skaičiaus ${\sqrt {a}}$ skiriasi mažiau kaip $10^{-10}.$

3 pavyzdys. Įsitikinsime, kad sekos $\{c_{n}\},$ kai $c_{n}={\frac {x^{n+1}}{(n+1)!}},$ o x – bet koks fiksuotas skaičius, riba lygi nuliui. Kai n – pakankamai didelis natūralusis skaičius, ${\frac {|x|}{n+1}}<1.$ Todėl, pradedant kuriuo nors numeriu N, turi būti $|c_{n+1}|<|c_{n}|,$ nes $|c_{n+1}|={\frac {|x|^{n}}{n!}}\cdot {\frac {|x|}{n+1}}=|c_{n}|\cdot {\frac {|x|}{n+1}}.$

Vadinasi, pradedant numeriu N, seka $\{|c_{n}|\}$ mažėja. Kadangi, be to, ji aprėžta iš apačios (pavyzdžiui nuliu), tai pagal 3.15 teoremą seka $\{|c_{n}|\}$ konverguoja. Tarkime, kad c yra tos sekos riba. Iš lygybės $|c_{n+1}|=|c_{n}|\cdot {\frac {|x|}{n+1}}$ išplaukia, kad $c=0,$ nes sekos $\{|c_{n+1}|\}$ riba lygi c, o sekos $\{{\frac {|x|}{n+1}}\}$ – nuliui.

_____________________

* Tą nelygybę įrodinėjant, užtenka pastebėti, kad, kai $t>0,$ ji ekvivalenti nelygybei $t^{2}-2t+1\geq 0.$

** Ta lygybė gaunama iš rekurentinės formulės $x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )}.$

Sekos, aproksimuojančios a šaknyje, konvergavimo greitis

Šio skyriaus 3 paragrafo 3 skirsnyje įrodėme, kad sekos $\{x_{n}\}$ apibrėžiamos rekurentine formule

$x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )},\quad n=1,\;2,\;3,\;...,\quad (3.10)$

kai $a>0,$ o $x_{1}$ – bet koks teigiamas skaičius, riba lygi ${\sqrt {a}}.$ Skaičiaus ${\sqrt {a}}$ artiniu galima laikyti bet kurį tos sekos narį $x_{n+1}.$ Tokiu atveju, be abejo, reikia išsiaiškinti, koks turi būti įteracijų* skaičius n, kad ${\sqrt {a}}$ artinys būtų nurodyto tikslumo.

______________

* Iteracija (lot. iteratio – kartojimas) – kokios nors matematinės operacijos pakartojimas. Šiuo atveju viena iteracija yra $x_{n+1}$ apskaičiavimas, kai žinomas $x_{n},$ pagal rekurentinę (3.10) formulę.

______________

Nagrinėsime seką $x_{n},$ apibrėžtą rekurentine (3.10) formule. Tos sekos elementą $x_{n}$ vadinsime skaičiaus $\gamma ={\sqrt {a}}\;$ n-uoju artiniu. Skaičių

$\varepsilon _{n}={\frac {x_{n}-\gamma }{\gamma }}\quad (3.11)$

vadinsime n-ojo artinio santykine paklaida.

Įrodysime teiginį apie santykinės paklaidos $\varepsilon _{n+1}$ įvertinimą, remiantis pirmojo artinio santykine paklaida $\varepsilon _{1}.$

Jei $x_{1}$ taip parinktas, kad $|\varepsilon _{1}|<{\frac {1}{2}},$ tai su bet kokiu $n\geq 1$ teisinga šitokia nelygybė:

$0<\varepsilon _{n+1}<\varepsilon _{1}^{2^{n}}.\quad (3.12)$

[ $0\leq \varepsilon _{n+1}\leq \varepsilon _{1}^{2^{n}},$ kai $\varepsilon _{1}$ gali būti lygus nuliui.]

Įrodymas. Iš (3.11) formulės

$x_{n}=\gamma (1+\varepsilon _{n}).\quad (3.13)$

Remdamiesi (3.10) ir (3.13) formulėmis bei lygybe ${\frac {a}{\gamma }}=\gamma ,$ gauname

$x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )}={\frac {1}{2}}\left[\gamma (1+\varepsilon )+{\frac {a}{\gamma (1+\varepsilon _{n})}}\right]={\frac {1}{2}}\left[{\frac {a(1+\varepsilon _{n})^{2}+a}{\gamma (1+\varepsilon _{n})}}\right]=$

$={\frac {1}{2}}\left[{\frac {a(1+2\varepsilon _{n}+\varepsilon _{n}^{2})+a}{\gamma (1+\varepsilon _{n})}}\right]={\frac {a}{\gamma }}\left[{\frac {(1+2\varepsilon _{n}+\varepsilon _{n}^{2})+1}{2(1+\varepsilon _{n})}}\right]=\gamma \left[{\frac {2+2\varepsilon _{n}+\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}\right]=\gamma \left[1+{\frac {\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}\right].$

Kadangi $x_{n+1}=\gamma (1+\varepsilon _{n+1}),$ tai savaime aišku, kad

$\gamma (1+\varepsilon _{n+1})=\gamma \left[1+{\frac {\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}\right],$

$\varepsilon _{n+1}={\frac {1}{2(1+\varepsilon _{n})}}\cdot \varepsilon _{n}^{2}.\quad (3.14)$

Pagal sąlygą $|\varepsilon _{1}|<{\frac {1}{2}}.$ Iš to išplaukia nelygybės $0<{\frac {1}{2(1+\varepsilon _{n})}}<1.$ Bet tada iš (3.14) lygybės, kai $n=1,$ išplaukia nelygybė $\varepsilon _{2}>0.$ Panašiai naudodamiesi (3.14) lygybe, kai $n=1,\;2,\;...,$ įsitikiname, kad visi $\varepsilon _{n+1}$ yra neneigiami ($n\geq 1$).

Iš (3.14) lygybės, nelygybių $0<{\frac {1}{2(1+\varepsilon _{n})}}<1$ ir iš to, kad visi $\varepsilon _{n},$ kai $n>1,$ yra neneigiami, išplaukia nelygybė $\varepsilon _{n+1}<\varepsilon _{n}^{2},$ kai $n\geq 1.$ Iš čia gauname dešiniąją (3.12) nelygybę. Teiginys įrodytas.

[$|\varepsilon _{1}|<{\frac {1}{2}},\;$ $\varepsilon _{2}<\varepsilon _{1}^{2},\;$ $\varepsilon _{3}<\varepsilon _{2}^{2}<(\varepsilon _{1}^{2})^{2}=\varepsilon _{1}^{4},\;$ $\varepsilon _{4}<\varepsilon _{3}^{2}<(\varepsilon _{1}^{4})^{2}=\varepsilon _{1}^{8}=\varepsilon _{1}^{2^{3}}$ ir taip toliau.]

Iš (3.12) nelygybių matyti, kad ${\sqrt {a}}$ artinio po n iteracijų santykinė paklaida įvertinama, remiantis pirmojo artinio $x_{1}$ santykine paklaida $\varepsilon _{1}$ ir iteracijų skaičiumi n. Vėliau įsitikinsime, jog tuo atveju, kai $a>1$*, pirmąjį artinį galima pasirinkti taip, kad santykinės paklaidos $\varepsilon _{1}$ modulis nebūtų didesnis už 0.05. Savaime aišku, taip pasirinkus $x_{1},$ santykinė paklaida tenkins įrodytojo teiginio sąlygas. Tuomet bus aišku ir koks turi būti iteracijų skaičius n, kad ${\sqrt {a}}$ artinio santykinė paklaida nebūtų didesnė už duotąjį skaičių $\varepsilon$: tą skaičių n galima rasti iš nelygybės* (ta nelygybė gaunama tiesiog iš (3.12) nelygybių)

$(0.05)^{2^{n}}<\varepsilon .\quad (3.15)$

Taigi tarkime, kad $a>1.$ Skaičių a išreikšime šitaip:

$a=2^{2k+i}M;\quad (3.16)$

čia k – sveikasis neneigiamas skaičius, skaičius i lygus arba nuliui, arba vienetui, o skaičius M tenkina sąlygą

$1\leq M<2.\quad (3.17)$

Pabrėšime, kad skaičius a išreiškiamas (3.16) pavidalu vieninteliu budu.

Pirmąjį artinį $x_{1}$ pasirinkime šitaip:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}.\quad (3.18)$

Įsitikinsime, kad tuo atveju, kai M – bet koks skaičius, tenkinantis (3.17) sąlygas, pirmojo ${\sqrt {a}}$ artinio $x_{1},$ rasto pagal (3.18) formulę, santikinės paklaidos $\varepsilon _{1}$ modulis ne didesnis kaip 0.05.

Įrodinėdami šį teiginį, remsimės tikslia santykinės paklaidos išraiška $\varepsilon _{1}={\frac {x_{1}-\gamma }{\gamma }}.$ Iš (3.16) lygybės gauname $\gamma ={\sqrt {a}}=2^{k}{\sqrt {2^{i}M}};$ todėl iš $\varepsilon _{1}$ išraiškos ir iš (3.18) formulės

$\varepsilon _{1}={\frac {x_{1}-\gamma }{\gamma }}={\frac {2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}-2^{k}{\sqrt {2^{i}M}}}{2^{k}{\sqrt {2^{i}M}}}}={\frac {{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}-{\sqrt {2^{i}M}}}{\sqrt {2^{i}M}}}.\quad (3.19)$

3.4 pav.

Alternatyvus 3.4 pav. Oy ašimi grafikas ištemptas tik maždaug 2 kartus daugiau nei Ox ašimi.

Kadangi skaičius i lygus arba nuliui, arba vienetui, o $M\geq 1,$ tai ${\sqrt {2^{i}M}}\geq 1.$ Todėl iš (3.19) lygybės išplaukia nelygybė

$|\varepsilon _{1}|\leq |{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}-{\sqrt {2^{i}M}}|.\quad (3.20)$

Pažymėkime ${\sqrt {2^{i}M}}$ raide X. Kadangi $1\leq M<2,$ o i lygus arba nuliui, arba vienetui, tai visos galimos X reikšmės priklauso intervalui [1; 2):

$1\leq X<2.\quad (3.21)$

Naudodami naujai įvestą žymėjimą, (3.20) nelygybę perrašysime šitaip:

$|\varepsilon _{1}|\leq |{\frac {1}{3}}X^{2}-X+{\frac {17}{24}}|.\quad (3.22)$

Pagal (3.22) nelygybę maksimalioji $|\varepsilon _{1}|$ reikšmė ne didesnė už maksimaliąją reiškinio $|{\frac {1}{3}}X^{2}-X+{\frac {17}{24}}|$ reikšmę, kai X tenkina (3.21) sąlygas. Kad būtų aiškiau, išnagrinėkime funkcijos $f(X)={\frac {1}{3}}X^{2}-X+{\frac {17}{24}}$ grafiką. Iš elementariosios matematikos kurso žinome, kad šios funkcijos grafikas yra parabolė (3.4 pav.)**, kurios viršūnė atitinka $X={\frac {3}{2}}.$

[ $f'(X)={\Big (}{\frac {1}{3}}X^{2}-X+{\frac {17}{24}}{\Big )}'={\frac {2}{3}}X-1,$

${\frac {2}{3}}X-1=0,\;\;\;{\frac {2}{3}}X=1,\;\;\;X={\frac {3}{2}}.$ ]

Kadangi $f(1)=f(2)={\frac {1}{24}},\;$ o $\;f{\Big (}{\frac {3}{2}}{\Big )}=-{\frac {1}{24}},\;$ tai X reikšmės, tenkinančios (3.21) sąlygas, atitinka $f(X)$ reikšmėms, esančioms tarp $-{\frac {1}{24}}$ ir ${\frac {1}{24}}.$ Kitaip sakant,

$|f(X)|=|{\frac {1}{3}}X^{2}-X+{\frac {17}{24}}|\leq {\frac {1}{24}}.$

Iš paskutinės ir (3.22) nelygybių gauname reikiamą $\varepsilon _{1}$ įvertinimą:

$|\varepsilon _{1}|\leq {\frac {1}{24}}<0.05.$

Pastaba. Jei nurodyta, kad santykinė paklaida $\varepsilon$ lygi $10^{-10},$ tai, norint tokiu tikslumu apskaičiuoti bet kokio skaičiaus $a>1$ kvadratinę šaknį, reikia, suradus $x_{1}$ pagal (3.18) formulę, atlikti tik tris iteracijas ($n=3$), nes

$(0.05)^{2^{3}}=0.05^{8}=0.0000000000390625=3.90625\cdot 10^{-11}<10^{-10}.$

[ $(1/24)^{8}=(0.041(6))^{8}\approx 9.0846890375538\cdot 10^{-12}<10^{-11}<10^{-10}.$ Pasirodo, santykinė paklaida $\varepsilon _{4}$ visada bus mažesnė už $10^{-11}$ po 3 iteracijų ($n=3$). ]

_______________

* Jei $a<1,$ tai $a={\frac {1}{b}}$ ir $b>1.$ Tada ${\sqrt {a}}={\frac {1}{\sqrt {b}}}.$

** 3.4 paveiksle ašies Oy mastelis 20 kartų didesnis už ašies Ox mastelį.

Pavyzdžiai

• Apskaičiuosime ${\sqrt {a}}={\sqrt {7}}.$ Iš (3.16) formulės $a=2^{2k+i}M.$ Matome, kad $2^{2k+i}$ turi būti lygus 4, t. y. $2^{2k+i}=4.$ Tada aišku, kad $k=1,\;\;i=0.$ Toliau randame M:

$7=2^{2k+i}M,$

$7=4M,$

$M={\frac {7}{4}}.$

Toliau randame $x_{1}$ pagal (3.18) formulę:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{1}{\Big (}{\frac {1}{3}}\cdot 2^{0}\cdot {\frac {7}{4}}+{\frac {17}{24}}{\Big )}=$

$=2{\Big (}{\frac {1}{3}}\cdot {\frac {7}{4}}+{\frac {17}{24}}{\Big )}={\frac {7}{6}}+{\frac {17}{12}}={\frac {14+17}{12}}={\frac {31}{12}}=$

= 2.58(3).

Skaičių $x_{2}$ randame pagal (3.10) rekurentinę formulę:

$x_{1+1}={\frac {1}{2}}{\Big (}x_{1}+{\frac {a}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {31}{12}}+{\frac {7}{\frac {31}{12}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {31}{12}}+{\frac {84}{31}}{\Big )}={\frac {1}{2}}\cdot {\frac {31^{2}+84\cdot 12}{12\cdot 31}}={\frac {1}{2}}\cdot {\frac {961+1008}{372}}={\frac {1969}{744}}=$

= 2.646505376344086021505376344086.

Toliau rasime $x_{3}$:

$x_{3}={\frac {1}{2}}{\Big (}x_{2}+{\frac {a}{x_{2}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {1969}{744}}+{\frac {7}{\frac {1969}{744}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {1969}{744}}+{\frac {5208}{1969}}{\Big )}={\frac {1}{2}}\cdot {\frac {1969^{2}+5208\cdot 744}{744\cdot 1969}}={\frac {7751713}{2929872}}=$

= 2.6457514184920023809913880196814.

[Atėmę tikslią ${\sqrt {7}}$ reikšmę iš $x_{3}$ reikšmės gauname:

2.6457514184920023809913880196814 - 2.6457513110645905905016157536393 = 1.0742741179048977226604214545141e-7.]

Toliau rasime $x_{4}$:

$x_{4}=x_{3+1}={\frac {1}{2}}{\Big (}x_{3}+{\frac {a}{x_{3}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {7751713}{2929872}}+{\frac {7}{\frac {7751713}{2929872}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {7751713}{2929872}}+{\frac {20509104}{7751713}}{\Big )}=$

$={\frac {1}{2}}\cdot {\frac {7751713^{2}+20509104\cdot 2929872}{2929872\cdot 7751713}}={\frac {1}{2}}\cdot {\frac {60,089,054,434,369+60,089,049,554,688}{2929872\cdot 7751713}}={\frac {120,178,103,989,057}{45,423,053,741,472}}=$

= 120178103989057/45423053741472 = 2.6457513110645927714791488926518.

Atėmę tikslią ${\sqrt {7}}$ reikšmę iš gautos $x_{4}$ reikšmės, gauname:

2.6457513110645927714791488926518 - 2.6457513110645905905016157536393 = 2.1809775331390125063383535564728e-15.

Tai yra gavome paklaidą $\varepsilon \approx 2.180977533139\cdot 10^{-15}<10^{-10}.$

Santykinė paklaida yra tokia:

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

(2.6457513110645927714791488926518 - 2.6457513110645905905016157536393)/2.6457513110645905905016157536393 =

= 8.2433202395785127055094193690687e-16.

Taigi $\varepsilon _{4}\approx 8.2433202395785\cdot 10^{-16}<10^{-10}.$

• Apskaičiuosime ${\sqrt {a}}={\sqrt {39}}.$ Iš (3.16) formulės $a=2^{2k+i}M.$ Pagal (3.17) $1\leq M<2.$ Todėl $2^{2k+i}=32.$ Tada $2k+i=5.$ Iš čia gauname, kad $k=2,\;\;i=1.$

Toliau randame M:

$a=2^{2k+i}M,$

$39=2^{2\cdot 2+1}M,$

$39=2^{5}M,$

$M={\frac {39}{32}}.$

Toliau randame $x_{1}$ pagal (3.18) formulę:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{2}{\Big (}{\frac {1}{3}}\cdot 2^{1}\cdot {\frac {39}{32}}+{\frac {17}{24}}{\Big )}=4{\Big (}{\frac {2}{3}}\cdot {\frac {39}{32}}+{\frac {17}{24}}{\Big )}={\frac {312}{96}}+{\frac {17}{6}}={\frac {312+17\cdot 16}{96}}={\frac {584}{96}}={\frac {73}{12}}=$

= 6.08(3).

[ ${\sqrt {39}}=$ 6.2449979983983982058468931209398. ]

Skaičių $x_{2}$ randame pagal (3.10) rekurentinę formulę:

$x_{2}={\frac {1}{2}}{\Big (}x_{1}+{\frac {a}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {73}{12}}+{\frac {39}{\frac {73}{12}}}{\Big )}=$

= 0.5*(6.0833333333333333333333333333333 + 6.4109589041095890410958904109589) =

= 0.5 * 12.494292237442922374429223744292 = 6.2471461187214611872146118721461.

Toliau rasime $x_{3}$:

$x_{3}={\frac {1}{2}}{\Big (}x_{2}+{\frac {a}{x_{2}}}{\Big )}=$

= 0.5*(6.2471461187214611872146118721461 + 39/6.2471461187214611872146118721461) =

= 6.2449983677207123204232036062421.

[Atėmę tikslią ${\sqrt {39}}$ reikšmę iš $x_{3}$ reikšmės, gauname:

6.2449983677207123204232036062421 - 6.2449979983983982058468931209398 = 3.693223141145763104853022918526e-7.]

Toliau rasime $x_{4}$:

$x_{4}=x_{3+1}={\frac {1}{2}}{\Big (}x_{3}+{\frac {a}{x_{3}}}{\Big )}=$

= 0.5*(6.2449983677207123204232036062421 + 39/6.2449983677207123204232036062421) =

= 6.2449979983984091265040071378015.

Atėmę tikslią ${\sqrt {39}}$ reikšmę iš $x_{4}$ reikšmės, gauname:

6.2449979983984091265040071378015 - 6.2449979983983982058468931209398 = 1.0920657114016861678930991073622e-14.

Taigi, gavome paklaidą $1.0920657114\cdot 10^{-14}.$ O santykinė paklaida $\varepsilon _{4}$ lygi

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

= (6.2449979983984091265040071378015 - 6.2449979983983982058468931209398)/6.2449979983983982058468931209398 =

= 1.7487046620059110065547770263338e-15.

T. y. $\varepsilon _{4}\approx 1.748704662\cdot 10^{-15}<10^{-10}.$

• Apskaičiuosime ${\sqrt {a}}={\sqrt {2}}.$ Iš (3.16) formulės $a=2^{2k+i}M.$ Pagal (3.17) $1\leq M<2.$ Todėl $2^{2k+i}=2.$ Tada $2k+i=1.$ Iš čia gauname, kad $k=0,\;\;i=1.$

Kadangi $2^{2k+i}=2,$ tai $M=1.$

Toliau randame $x_{1}$ pagal (3.18) formulę:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{0}{\Big (}{\frac {1}{3}}\cdot 2^{1}\cdot 1+{\frac {17}{24}}{\Big )}={\frac {2}{3}}+{\frac {17}{24}}={\frac {2\cdot 8+17}{24}}={\frac {33}{24}}=1.375.$

Skaičių $x_{2}$ randame pagal (3.10) rekurentinę formulę:

$x_{2}={\frac {1}{2}}{\Big (}x_{1}+{\frac {a}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}1.375+{\frac {2}{1.375}}{\Big )}=$

= 0.5*(1.375 + 1.4545454545454545454545454545455) = 0.5 * 2.8295454545454545454545454545455 =

= 1.4147727272727272727272727272727.

Toliau rasime $x_{3}$:

$x_{3}={\frac {1}{2}}{\Big (}x_{2}+{\frac {a}{x_{2}}}{\Big )}=$

= 0.5*(1.4147727272727272727272727272727 + 2/1.4147727272727272727272727272727) =

= 1.4142136728733114275282949981745.

[1.4142136728733114275282949981745 - 1.4142135623730950488016887242097 = 0.0000001105002163787266062739648 $\approx 1.105002163787266\cdot 10^{-7}.$]

Toliau rasime $x_{4}$:

$x_{4}={\frac {1}{2}}{\Big (}x_{3}+{\frac {a}{x_{3}}}{\Big )}=$

= 0.5*(1.4142136728733114275282949981745 + 2/1.4142136728733114275282949981745) =

= 1.4142135623730993657935457387426.

Atėmę tikslią ${\sqrt {2}}$ reikšmę iš $x_{4}$ reikšmės, gauname:

1.4142135623730993657935457387426 - 1.4142135623730950488016887242097 = 4.316991857014532913283559342618e-15.

Santykinė paklaida yra tokia:

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

=(1.4142135623730993657935457387426 - 1.4142135623730950488016887242097)/1.4142135623730950488016887242097 =

= 3.0525742164220827729222598916484e-15.

Gavome, kad $\varepsilon _{4}\approx 3.052574216422\cdot 10^{-15}<10^{-10}.$

Šiaip, skaičiuojant su kompiuteriu, santikinę paklaidą $\varepsilon _{4}$ (arba $\varepsilon _{n+1}$) reikia padauginti iš 10, nes procesoriaus x87 FPU be eksponentės dirba su skaičiais nuo 1.00000000000000001 iki 9.99999999999999999. Tada nebus prarastas tikslumas, jei įvertinsime, kad FPU paklaida yra $\delta <10\varepsilon _{n+1}.$

Šaknies traukimas iš didelių skaičių

Bet kokį skaičių a ($a>1$) galima užrašyti taip:

$a=p\cdot 10^{j}\cdot 10^{v}.\quad (1)$

Čia $1<p<10;\;$ skaičius v gali būti nulis arba teigiamas lyginis skaičius (v gali būti 0, 2, 4, 6, 8, 10 ir t. t.); skaičius j lygus arba nuliui arba vienetui.

Jeigu ištrauksime kvadratinę šaknį iš a, tai gausime:

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10^{j}}}\cdot 10^{v \over 2}.\quad (2)$

Jegiu v nelygus nuliui, tai $10^{v/2}$ bus sveikasis skaičius su dvigubai mažiau nulių nei pas $10^{v}.$ Jeigu $j=1,$ tai tereikia apskaičiuoti ${\sqrt {10}},$ kurį bus galima visada taikyti (2) formulėje. Vadinasi, tereikia mokėti greitai ištraukti šaknį iš skaičiaus p.

1) Jeigu $1<p<2,\;$ tai pagal (3.16) formulę

$p=2^{2k+i}M;$

iš čia $k=0,\;\;i=0,\;\;M=p.$

Tada pagal (3.18) formulę randame $x_{1}$:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{0}{\Big (}{\frac {1}{3}}\cdot 2^{0}\cdot p+{\frac {17}{24}}{\Big )}={\frac {1}{3}}p+{\frac {17}{24}}.$

Ir toliau, darant iteracijas, apskaičiuojamas $x_{2},\;x_{3},\;x_{4},\;...$ pagal (3.10) formulę:

$x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {p}{x_{n}}}{\Big )},\quad n=1,\;2,\;3,\;...$

2) Jeigu $2\leq p<4,\;$ tai pagal (3.16) formulę

$p=2^{2k+i}M;$

iš čia $k=0,\;\;i=1.$ Tada $p=2M,\;\;M={\frac {p}{2}}.$

Pagal (3.18) formulę randame $x_{1}$:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{0}{\Big (}{\frac {1}{3}}\cdot 2^{1}\cdot {\frac {p}{2}}+{\frac {17}{24}}{\Big )}={\frac {1}{3}}p+{\frac {17}{24}}.$

Matome, kad šiuo atveju artinio $x_{1}$ formulė tokia pati kaip ir 1) atveju.

Toliau skaičiuojama kaip ir pirmu atveju.

3) Jeigu $4\leq p<8,\;$ tai pagal (3.16) formulę

$p=2^{2k+i}M;$

iš čia $k=1,\;\;i=0.$ Tada $p=4M,\;\;M={\frac {p}{4}}.$

Pagal (3.18) formulę randame $x_{1}$:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{1}{\Big (}{\frac {1}{3}}\cdot 2^{0}\cdot {\frac {p}{4}}+{\frac {17}{24}}{\Big )}={\frac {1}{6}}p+{\frac {17}{12}}.$

Toliau skaičiuojama kaip ir pirmu atveju.

4) Jeigu $8\leq p<10,\;$ tai pagal (3.16) formulę

$p=2^{2k+i}M;$

iš čia $k=1,\;\;i=1.$ Tada $p=8M,\;\;M={\frac {p}{8}}.$

Pagal (3.18) formulę randame $x_{1}$:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{1}{\Big (}{\frac {1}{3}}\cdot 2^{1}\cdot {\frac {p}{8}}+{\frac {17}{24}}{\Big )}={\frac {1}{6}}p+{\frac {17}{12}}.$

Matome, kad šiuo atveju artinio $x_{1}$ formulė tokia pati kaip ir 3) atveju.

Toliau skaičiuojama kaip ir pirmu atveju.

Apibendrinant, skaičiaus ${\sqrt {p}}\,$ ($1<p<10$) apskaičiavimui, taikoma $x_{1}$ reikšmė

$x_{1}={\frac {1}{3}}p+{\frac {17}{24}},\quad (3)$

kai $1<p<4;$

ir taikoma $x_{1}$ reikšmė

$x_{1}={\frac {1}{6}}p+{\frac {17}{12}},\quad (4)$

kai $4\leq p<10.$

• Apskaičiuosime ${\sqrt {a}}={\sqrt {1.5\cdot 10^{8}}}.$ Čia $p=1.5,\;\;j=0,\;\;v=8.$ Tada pagal (1) formulę

$a=p\cdot 10^{j}\cdot 10^{v}=1.5\cdot 10^{0}\cdot 10^{8}.$

Pagal (2) formulę

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10^{j}}}\cdot 10^{v \over 2}={\sqrt {p}}\cdot 10^{8 \over 2}={\sqrt {p}}\cdot 10^{4}.$

Toliau randame $x_{1}$ skaičiui $p=1.5$:

$x_{1}={\frac {1}{3}}p+{\frac {17}{24}}={\frac {1}{3}}\cdot {\frac {3}{2}}+{\frac {17}{24}}={\frac {1}{2}}+{\frac {17}{24}}={\frac {12+17}{24}}={\frac {29}{24}}=$

= 1.2083333333333333333333333333333.

[ ${\sqrt {1.5}}=$ 1.2247448713915890490986420373529. ]

[1.2247448713915890490986420373529 - 1.2083333333333333333333333333333 = 0.0164115380582557157653087040196; tai yra ne santykinė paklaida, o reali paklaida.]

Toliau pagal (3.10) formulę apskaičiuojame $x_{2}$:

$x_{2}=x_{1+1}={\frac {1}{2}}{\Big (}x_{1}+{\frac {p}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {29}{24}}+{\frac {1.5}{\frac {29}{24}}}{\Big )}=$

= 0.5*(1.2083333333333333333333333333333 + 1.5/1.2083333333333333333333333333333) = 1.2248563218390804597701149425287.

Toliau pagal (3.10) formulę apskaičiuojame $x_{3}$:

$x_{3}=x_{2+1}={\frac {1}{2}}{\Big (}x_{2}+{\frac {p}{x_{2}}}{\Big )}=$

= 0.5*(1.2248563218390804597701149425287 + 1.5/1.2248563218390804597701149425287) = 1.2247448764620622240199548319682.

[1.2247448764620622240199548319682 - 1.2247448713915890490986420373529 = 0.0000000050704731749213127946153 = $5.07047\cdot 10^{-9}.$]

Toliau pagal (3.10) formulę apskaičiuojame $x_{4}$:

$x_{4}=x_{3+1}={\frac {1}{2}}{\Big (}x_{3}+{\frac {p}{x_{3}}}{\Big )}=$

= 0.5*(1.2247448764620622240199548319682 + 1.5/1.2247448764620622240199548319682) = 1.2247448713915890595945823395751.

Iš ką tik gautos $x_{4}$ reikšmės atėmę ${\sqrt {1.5}}$ reikšmę, gauname:

1.2247448713915890595945823395751 - 1.2247448713915890490986420373529 = 0.0000000000000000104959403022222 = $1.049594\cdot 10^{-17}.$

Santykinė paklaida yra tokia:

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

(1.2247448713915890595945823395751 - 1.2247448713915890490986420373529)/1.2247448713915890490986420373529 =

= 8.5698993703859497564530032133562e-18.

Taigi $\varepsilon _{4}\approx 8.56989937\cdot 10^{-18}<10^{-10}.$

Iš anksčiau apskaičiuoto

${\sqrt {a}}={\sqrt {p}}\cdot 10^{4}=x_{4}\cdot 10^{4}=$

= 1.2247448713915890595945823395751 * 10000 = 12247.448713915890595945823395751.

Atėmę tikslią ${\sqrt {a}}$ iš 12247.448713915890595945823395751, gauname paklaidą:

12247.448713915890595945823395751 - 12247.448713915890490986420373529 = 1.0495940302222154304017026259951e-13.

O santykinė paklaida skaičiui ${\sqrt {a}}$ yra

(12247.448713915890595945823395751 - 12247.448713915890490986420373529)/12247.448713915890490986420373529 =

= 8.5698993703859497564530032133562e-18.

Tokia pati, kaip skaičiui ${\sqrt {p}}.$

• Apskaičiuosime ${\sqrt {a}}={\sqrt {9.6\cdot 10^{11}}}.$ Čia $p=9.6,\;\;j=1,\;\;v=10.$ Tada pagal (1) formulę

$a=p\cdot 10^{j}\cdot 10^{v}=9.6\cdot 10^{1}\cdot 10^{10}.$

Pagal (2) formulę

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10^{j}}}\cdot 10^{v \over 2}={\sqrt {p}}\cdot {\sqrt {10}}\cdot 10^{10 \over 2}={\sqrt {p}}\cdot {\sqrt {10}}\cdot 10^{5}.$

Šaknį iš 10 lengva apskaičiuoti (nes (4) formulė tinka skaičiams $4\leq a<16$).

${\sqrt {10}}=$ 3.1622776601683793319988935444327.

Toliau randame $x_{1}$ skaičiui $p=9.6$ pagal (4) formulę:

$x_{1}={\frac {1}{6}}p+{\frac {17}{12}}={\frac {1}{6}}\cdot 9.6+{\frac {17}{12}}=$

= 0.16666666666666666666666666666667*9.6 + 1.4166666666666666666666666666667 = 3.0166666666666666666666666666667.

[ ${\sqrt {9.6}}=$ 3.0983866769659335081434123198259. ]

[3.0983866769659335081434123198259 - 3.0166666666666666666666666666667 = 0.0817200102992668414767456531592; tai yra ne santykinė paklaida, o reali paklaida.]

Toliau pagal (3.10) formulę apskaičiuojame $x_{2}$:

$x_{2}=x_{1+1}={\frac {1}{2}}{\Big (}x_{1}+{\frac {p}{x_{1}}}{\Big )}=$

= 0.5*(3.0166666666666666666666666666667 + 9.6/3.0166666666666666666666666666667) = 3.0994935543278084714548802946593.

Toliau pagal (3.10) formulę apskaičiuojame $x_{3}$:

$x_{3}=x_{2+1}={\frac {1}{2}}{\Big (}x_{2}+{\frac {p}{x_{2}}}{\Big )}=$

= 0.5*(3.0994935543278084714548802946593 + 9.6/3.0994935543278084714548802946593) = 3.0983868746074953885623653564026.

[3.0983868746074953885623653564026 - 3.0983866769659335081434123198259 = 0.0000001976415618804189530365767 =

$\approx 1.9764156188\cdot 10^{-7}.$]

Toliau pagal (3.10) formulę apskaičiuojame $x_{4}$:

$x_{4}=x_{3+1}={\frac {1}{2}}{\Big (}x_{3}+{\frac {p}{x_{3}}}{\Big )}=$

= 0.5*(3.0983868746074953885623653564026 + 9.6/3.0983868746074953885623653564026) = 3.0983866769659398117763289613096.

Iš ką tik gautos $x_{4}$ reikšmės atėmę ${\sqrt {9.6}}$ reikšmę, gauname:

3.0983866769659398117763289613096 - 3.0983866769659335081434123198259 = 0.0000000000000063036329166414837 =

$\approx 6.3036329166414837\cdot 10^{-15}.$

Santykinė paklaida yra tokia:

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

(3.0983866769659398117763289613096 - 3.0983866769659335081434123198259)/3.0983866769659335081434123198259 =

= 2.0344887755631127106428103904763e-15.

Taigi $\varepsilon _{4}\approx 2.03448877556\cdot 10^{-15}<10^{-10}.$

Iš ankstesnių skaičiavimų

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10}}\cdot 10^{5}=x_{4}\cdot {\sqrt {10}}\cdot 10^{5}=$

= 3.0983866769659398117763289613096 * 3.1622776601683793319988935444327 * 100000 = 979795.89711327323266266864962311.

Tiksli ${\sqrt {a}}$ reikšmė yra tokia:

${\sqrt {a}}={\sqrt {9.6\cdot 10^{11}}}=$ 979795.89711327123927891362988236.

Atėmę tikslią ${\sqrt {a}}$ reikšmę iš 979795.89711327323266266864962311, gauname paklaidą:

979795.89711327323266266864962311 - 979795.89711327123927891362988236 = 0.00000000199338375501974075 =

$\approx 1.993383755\cdot 10^{-9}.$

Santykinė paklaida skaičiui ${\sqrt {a}}$ yra

(979795.89711327323266266864962311 - 979795.89711327123927891362988236)/979795.89711327123927891362988236 =

= 2.0344887755631126970242116177837e-15,

t. y. tokia pati (su 16 pirmų tokių pačių skaitmenų) kaip skaičiui ${\sqrt {p}}.$

[Čia kažkas negerai su tais 16 skaitmenų, nes iš viso ciferblate yra 32 skaitmenys... Jeigu įdėti (padaryti copy-paste) į kalkuliatorių tokį reiškinį:

(3.0983866769659398117763289613096 - 3.0983866769659335081434123198259)/3.0983866769659335081434123198259 -(979795.89711327323266266864962311 - 979795.89711327123927891362988236)/979795.89711327123927891362988236 =

tai gaunamas toks rezultatas:

= 1.3618598772692596082816173047622e-32.

Kažkokie dvigubi skaičiavimo standartai su "Windows 10" kaluliatorium... Gal po kokių nors dalybos iteracijų labai prasiranda tikslumas, nes ne visi skaitmenys Winows 10 kalkuliatoriuje rodomi (jų yra virš 40, gal apie 47 skaitmenys iš viso, iš kurių rodomi tik 32 skaitmenys). Nors gal ir ne tai, nes ten tik atimtis atliekama... Gal koks bug'as?

Bendrai, viskas teisingai, nes

2.0344887755631127106428103904763e-15 - 2.0344887755631126970242116177837e-15 = 1.36185987726926e-32.

Tai rodo, kad skaičiavimo tikslumas yra 32 skaitmenys, nes jų šituose skaičiavimuose tiek ir yra.]

Jeigu $0<a<1,$ tai pažymėję $b={\frac {1}{a}},$ galime, jei b didelis skaičius, apskaičiuoti ${\sqrt {b}}$ pagal pateiktą aukščiau būdą. Apskaičiavę ${\sqrt {b}},$ galime rasti ${\sqrt {a}}$ taip:

${\sqrt {a}}={\frac {1}{\sqrt {b}}}.$

Tačiau toks ${\sqrt {a}}\;$ ($0<a<1$) skaičiavimo būdas nėra ekonomiškas, nes reikia atlikti dvi [papildomas] dalybos operacijas.

Pateiksime ekonomiškesnį, greitesnį ${\sqrt {a}}$ skaičiavimo būdą. Jei $0<a<1,$ tai a galima užrašyti taip:

$a=p\cdot 10^{j}\cdot 10^{v}.\quad (5)$

Čia $1<p<10;\;$ skaičius v gali būti nulis arba neigiamas lyginis skaičius (v gali būti 0, -2, -4, -6, ...); skaičius j lygus arba nuliui arba $-1.$

Jeigu ištrauksime kvadratinę šaknį iš a, tai gausime:

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10^{j}}}\cdot 10^{v \over 2}.\quad (6)$

Jeigu $j=-1,$ tai tereikia apskaičiuoti $({\sqrt {10}})^{-1},$ kurį bus galima visada taikyti (6) formulėje. Vadinasi, tereikia mokėti greitai ištraukti šaknį iš skaičiaus p. O toliau jau lengva pagal (6) formulę surasti kam lygus ${\sqrt {a}}.$ Kaip rasti ${\sqrt {p}}$ pateikta aukščiau.

Kiek pamenu, skaičiau, kad jeigu norima padauginti skaičių dvejetainėje sistemoje iš 2, tai reikia visus bitus pastumti per vieną bito poziciją į kairę. O jeigu norima skaičių dvejetainėje sistemoje padalinti iš 2, tai reikia visus bitus pastumti viena bito pozicija į dešinę. Tokiu budu, jeigu skaičius j nelygus nuliui, tai reikia skaičių a padauginti iš 10, jei 0<a<1 (arba padauginti iš 0.1, jei a>1). Tada skaičius v bus lyginis skaičius. Toliau norint padalinti v iš 2, tereikia rodiklio v visus bitus pastumti per vieną bitą į dešinę. Taip mes gausime v/2. Ir nereikia procesoriaus FPU daryti sudetingų eksponentės dalybos operacijų (skaičiuojant šaknį iš kokio nors skaičiaus).

Lieka tik paslaptis kodėl procesorių FPU šaknies traukimo operacijai sunaudoja tiek pat arba tik apie 1.5 karto (priklausomai nuo procesoriaus) daugiau ciklų negu dalybos operacijai. Net jeigu artinį $x_{1}$ kažkokiu budu galima gauti su santykine paklaida apie $10^{-6},$ tai vis tiek reikės dviejų iteracijų ir to pasekoje reikės dviejų dalybos operacijų (kiekvienoje iteracijoje yra viena dalybos operacija). Kad užtektų apskaičiuoti ${\sqrt {a}}$ iš vienos iteracijos, artinio $x_{1}$ santykinė paklaida $\varepsilon _{1}$ turi būti apie $10^{-9}$ (kad gauti ${\sqrt {a}}$ su santykine paklaida apie $10^{-18}$). Jeigu kažkokiu budu galima pagal skaičių p pažiūrėt į kažkokią artinių $x_{1}$ lentelę, tai ta lentelė turės apie milijardą artinių ir užims apie $80\;({\text{bit}})\cdot 10^{9}=$ 10 GB, t. y. apie 10 gigabaitų. Į procesorių akivaizdu to neįstatysi, tuo labiau, kad, pavyzdžiui, Pentium II/III naudojo 17...38 ciklus dalybos operacijai ir 27...50 ciklų šaknies traukimo operacijai (tik apie 1.5 karto daugiau ciklų šakniai nei dalybai). Tais laikais ir 1 MB (megabaito) kokios nors atminties (pvz., cache atminties) įdėjimas į procesorių buvo techniškai neįmanomas ar sunkiai pasiekiamas dalykas.

Kadangi artinio $x_{1}$ skaičiavimui [pateiktam aukščiau] nereikia daryti dalybos operacijos, tai galbūt yra, tarkim, apie 100 ar 1000 artinio $x_{1}$ skaičiavimo grupių skaičiui ${\sqrt {p}}.$ Kaip matėme, tik dviejų grupių (1<p<4 ir $4\leq p<10$) užtenka ${\sqrt {p}}$ skaičiavimui, kad gauti pirmo artinio $x_{1}$ santykinę paklaidą $\varepsilon _{1}<1/24<0.05.$ Jeigu sukurti apie 10 tokių artinio $x_{1}$ grupių (1<p<2, 2<p<3, 3<p<4, ..., 9<p<10) arba apie 100 tokių artinio $x_{1}$ grupių (1<p<1.1, 1.1<p<1.2, 1.2<p<1.3, ..., 9.9<p<10), tai peršasi mintis, kad 10-iai grupių tikslumas gali būti apie $10^{-5},$ o šimtui tokių grupių santykinės paklaidos $\varepsilon _{1}$ tikslumas gali būti apie $10^{-10}.$ Tik tai gali užimti darbo suradimui artinio $x_{1}$ formulės kiekvienai grupei (tokia formulė padidinanti tikslumą siauro/trumpo intervalo grupei turėtų egzistuot matematiškai mąstant ir lažinantis...). Tik gal tikslumo padidėjimas gali būti ne toks žymus, kai grupių padaugėja apie 10 kartų (arba 100 kartų), nežinau...

Kitas klausimas yra, kaip kiekvieną skaičių p nusiusti į jį atitinkančią vieną iš 10 (ar vieną iš 100) grupę (kad paskui gauti artinį $x_{1}$). Programuotojas sakytų, kad reikia vidutiniškai praifint 5 kartus, jei yra 10 grupių (arba ~50 kartų, jei yra 100 grupių). Pavyzdžiui, jei yra skaičius 3.5346, tai pirmą kartą paklaust ar skaičius p [jei yra 10 grupių] yra mažiau už 2, antrą kartą paklaust ar skaičius p yra mažiau už 3 ir trečią kartą palyginus 3.5346 su 4 nustatoma, kad skaičius p patenka į trečią grupę, nes yra mažiau už 4. Bet su ifinimais šaknies skaičiavimai bus dar ilgesni. Todėl su tranzistoriais galima sukurt painią schemą, kuri, jeigu skaičiaus p pirmas skaitmuo yra 3, siunčia skaičių p į trečią grupę, kurioje skaičiuoja artinį $x_{1}$ skaičiams nuo 3 iki 4 ir gauna numanomai minėtą tikslumą $10^{-5}.$ O jeigu skaičiaus p pirmas skaitmuo yra 4, tai siunčiasi per tą loginę schemą toks skaičius į ketvirtą grupę (į kurią įeina skaičiai nuo 4 iki 5) ir t. t. Kas ten žino, gal taip ir randamas [skaičiaus p] $x_{1}$ artinys (apie $10^{-10}$ tikslumu) naudojant apie 100 tokių grupių. Ir tada užtenka tik vienos iteracijos (t. y. užtenka vienos dalybos operacijos [kuri yra iteracijoje]), kad gauti apie $10^{-18}$ tikslumą.

Pastebėsime, kad jeigu artinys $x_{1}={\sqrt {a}},$ tai pagal (3.10) formulę

$x_{2}={\frac {1}{2}}{\Big (}x_{1}+{\frac {a}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}{\sqrt {a}}+{\frac {a}{\sqrt {a}}}{\Big )}={\frac {1}{2}}{\Big (}{\sqrt {a}}+{\sqrt {a}}{\Big )}={\sqrt {a}}.$

Čia aiškinama kaip dalybos operacija vyksta procesoriuje:

https://www.quora.com/How-is-integer-division-implemented-in-Computer-Hardware

https://en.wikipedia.org/wiki/Division_algorithm#SRT_division

Pavyzdžiui čia:

https://electronics.stackexchange.com/questions/22410/how-does-division-occur-in-our-computers

skaitant, gali būti neaišku kas yra "reciprocal", "denominator", "numerator". Jie yra atitinkamai 1/x, vardiklis, skaitiklis.

Šis tekstas:

"For floating point numbers, however, it may be quicker to first compute the reciprocal of the denominator and then multiply that times your numerator."

verčiamas per Google translate taip:

"Tačiau slankiojo kablelio skaičiams gali būti greičiau pirmiausia apskaičiuoti vardiklio atvirkštinį skaičių ir tada padauginti jį iš skaitiklio."

Ten pat duotas pavyzdys, jeigu norima apskaičiuoti 1/3.5. Iš pradžių spėjama, koks turi būti g, kad g*3.5 būtų lygus 1. Spėjimas, kad g=0.3. Tada g*3.5=0.3*3.5=1.05. Tada toliau naudojama pagerinta formulė g'=g(2-gd). Gauname 0.3*(2-1.05)=0.285. Jau gana arti 1/3.5=0.28571428571428571428571428571429. Toliau kartojamas procesas: 0.285*(2-0.285*3.5)=0.2857125. Dar toliau: 0.2857125*(2-0.2857125*3.5)=0.285714285703125. Paklaida yra:

0.28571428571428571428571428571429 - 0.285714285703125 = 0.00000000001116071428571428571429 $\approx 1.11607142857\cdot 10^{-11}.$

Čia

https://www.reddit.com/r/explainlikeimfive/comments/yc5qkp/eli5_how_does_a_computer_do_division_efficiently/

dar papildomai rašo apie šitą būdą, kuris ką tik buvo pateiktas skaičiuot 1/D, o suradus kam tai lygu, galima paskui sužinot kam lygu N/D, padauginus N iš apskaičiuoto 1/D (kai D=p, t. y. kai 1<D<10, o jei D>10, tai skaičiaus D eksponentės rodiklis atimamas iš skaičiaus N eksponentės rodiklio...).

Šito metodo išvedimas aprašytas angliškoje Vikipedijoje: https://en.wikipedia.org/wiki/Division_algorithm#Newton–Raphson_division

Jame naudojamas Niutono metodas (liesinių metodas), kuris taikomas polinomo realiųjų šaknų apytikslėms reikšmėms rasti.

Greitas šaknies iš skaičiaus a traukimo būdas

Parodysime praktiškai greičiausią kvadratinės šaknies traukimo būdą. Šitas būdas paremtas Niutono metodu (liestinių metodu). Todėl pirma trumpai apžvelgsime Niutono metodą.

Niutono metodas. Naudojant Niutono metodą galima rasti polinomo (arba kokios nors funkcjos) f(x) šaknį (sprendinį) bet kokiu norimu tikslumu. Nagrinėkime atkarpą [a; b], kurioje yra viena šaknis c (taške c funkcija f(x) kerta ašį Ox). Toliau, priklausomai nuo to ar funkcija f(x) atkarpoje [a; b] įgaubta ar išgaubta bei priklausomai nuo to ar funkcija f(x) didėja ar mažėja atkarpoje [a; b], pasirenkamas tinkamas atkarpos vienas iš galų (a arba b). Pažymėkime tą galą $a_{0}$ ($a_{0}=a$ arba $a_{0}=b$). Toliau taške $a_{0}$ surandame funkcijos f(x) liestinės lygtį, kuri atrodo taip:

$y=f(a_{0})+f'(a_{0})(x-a_{0}).$

Arba taip:

$y-f(a_{0})=f'(a_{0})(x-a_{0}).\quad (1)$

Toliau randamas liestinės y(x) susikirtimo taškas su Ox ašimi. Tą tašką pažymėkime (x; y)=(d; f(d))=(d; 0). Įstatę šio taško koordinates į (1) lygtį gausime:

$0-f(a_{0})=f'(a_{0})(d-a_{0}),$

$-{\frac {f(a_{0})}{f'(a_{0})}}=d-a_{0},$

$d=a_{0}-{\frac {f(a_{0})}{f'(a_{0})}}.\quad (2)$

Dabar turime siauresnį intervalą, kuriame ant Ox ašies yra šaknis c (per kurią pereina funkcjos f(x) grafikas). Ir dabar reikšmė $d=a_{1}$ yra arčiau šaknies c. Pakartoję šį procesą rasime tašką $a_{2},$ esantį dar arčiau prie c taško (ant Ox ašies):

$a_{2}=a_{1}-{\frac {f(a_{1})}{f'(a_{1})}}.$

Taip kartodami (2) formulę galime labai tiksliai rasti funkcijos f(x) susikirtimo su Ox ašimi tašką c. Šios formulės paklaida mažėja kvadratu po kiekvienos iteracijos, kai $a_{n}$ konverguoja į tašką c (po kiekvienos iteracijos, tikslumas didėja daug greičiau, nei jeigu spėliotume šaknį c, dalindami intervalą [a; b] kaskart per pusę ir tikrindami ar skirtinguose intervalų pusėse funkcija f(x) turi priešingus ženklus). Išsamiai apie Niutono metodą galima paskaityti aukštosios algebros arba matematinės analizės vadovėliuose.

Čia pasakojama apie ${\sqrt {a}}$ skaičiavimo metodus: https://www.quora.com/How-do-computers-calculate-square-roots

O čia

https://en.wikipedia.org/wiki/Fast_inverse_square_root#Newton's_method

aprašoma, kaip apskaičiuoti $1/{\sqrt {a}}.$ Kai apytiksliai apskaičiuojama $1/{\sqrt {a}},$ reikia gautą rezultatą padauginti iš a. Tada ir gausime ${\sqrt {a}}.$ Be to $1/{\sqrt {a}}$ rezultatas yra net naudingesnis, nes nereikia dalinti ${\sqrt {a}}$ iš a, kad gauti $1/{\sqrt {a}}.$

Pereikime prie ${\frac {1}{\sqrt {a}}}$ skaičiavimo algoritmo, kuris vadinasi Niutono metodu [šaknies skaičiavimui]. Pažymėkime $y={\frac {1}{\sqrt {a}}},$ o skaičių a pažymėkime per x. Tada

$f(y)={\frac {1}{y^{2}}}-x=0.$

Jeigu $y_{n}$ yra apytikslė skaičiaus y reikšmė, tai $y_{n+1}$ yra dar tikslesnė skaičiaus y reikšmė (pagal liestinių metodą):

$y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}.\quad (3)$

Surandame f(y) išvestinę:

$f'(y)=({\frac {1}{y^{2}}}-x)'={\frac {-2}{y^{3}}}.$

$f'(y_{n})={\frac {-2}{y_{n}^{3}}}.$

Toliau įstatome $f'(y_{n})$ ir $f(y_{n})$ į (3) formulę:

$y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}=y_{n}-{\frac {{\frac {1}{y_{n}^{2}}}-x}{\frac {-2}{y_{n}^{3}}}}=y_{n}+{\frac {{\frac {y_{n}^{3}}{y_{n}^{2}}}-xy_{n}^{3}}{2}}=$

$=y_{n}+{\frac {y_{n}-xy_{n}^{3}}{2}}=y_{n}+{\frac {y_{n}}{2}}-{\frac {xy_{n}^{3}}{2}}={\frac {3y_{n}}{2}}-{\frac {xy_{n}^{3}}{2}}=$

$=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right)={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}.$

Gavome rekurentinę formulę $1/{\sqrt {a}}$ reikšmei skaičiuoti:

$y_{n+1}=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right)={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}.\quad (4)$

Pavyzdžiai

• Apskaičiuosime ${\frac {1}{\sqrt {x}}}={\frac {1}{\sqrt {2.3}}}$ ir rasime kam lygus ${\sqrt {2.3}}={\sqrt {a}}={\frac {x}{\sqrt {x}}}={\frac {2.3}{\sqrt {2.3}}}.$

Pirmiausia turime žinoti $y={\frac {1}{\sqrt {2.3}}}$ apytikslią reikšmę. Tą reikšmę procesoriaus FPU gali gauti iš tam tikros lentelės, kuri galbūt yra procesoriaus L1 Data cache'e arba L2 cache. Tiksli ${\sqrt {2.3}}$ reikšmė yra tokia:

${\sqrt {2.3}}=$ 1.5165750888103101108513650872564.

O $y={\frac {1}{\sqrt {2.3}}}=$ 0.65938047339578700471798482054627.

Mes paimsime, kad $y_{1}=0.6.$ Tada pagal formulę

$y_{n+1}=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right)={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}$

rasime

$y_{2}={\frac {y_{1}\left(3-xy_{1}^{2}\right)}{2}}={\frac {0.6\left(3-2.3\cdot 0.6^{2}\right)}{2}}={\frac {0.6\left(3-0.828\right)}{2}}={\frac {0.6\cdot 2.172}{2}}={\frac {1.3032}{2}}=0.6516.$

Toliau pagal (4) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(3-xy_{2}^{2}\right)}{2}}={\frac {0.6516\left(3-2.3\cdot 0.6516^{2}\right)}{2}}={\frac {0.6516\left(3-2.3\cdot 0.42458256\right)}{2}}={\frac {0.6\cdot 2.023460112}{2}}={\frac {1.3184866089792}{2}}=0.6592433044896.$

[0.65938047339578700471798482054627 - 0.6592433044896 = 0.00013716890618700471798482054627 $\approx 1.37168906187\cdot 10^{-4}.$]

Toliau pagal (4) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(3-xy_{3}^{2}\right)}{2}}=$

= 0.6592433044896*(3 - 2.3*0.6592433044896^2)/2 = 0.65938043059651203399501910136126.

[0.65938047339578700471798482054627 - 0.65938043059651203399501910136126 =

= 0.00000004279927497072296571918501 $\approx 4.279927497\cdot 10^{-8}.$]

Toliau pagal (4) formulę rasime $y_{5}$:

$y_{5}={\frac {y_{4}\left(3-xy_{4}^{2}\right)}{2}}=$

= 0.65938043059651203399501910136126*(3 - 2.3*0.65938043059651203399501910136126^2)/2 = 0.65938047339578283767489142997287.

Taigi, $y_{5}$ paklaida yra tokia:

0.65938047339578700471798482054627 - 0.65938047339578283767489142997287 =

= 0.0000000000000041670430933905734 $=4.1670430933905734\cdot 10^{-15}.$

Po 4 iteracijų gavome paklaidą mažesne nei $10^{-14}.$

Dabar galime rasti a šaknyje:

${\sqrt {2.3}}={\sqrt {a}}={\sqrt {x}}=2.3\cdot y_{5}=$ 2.3*0.65938047339578283767489142997287 = 1.5165750888103005266522502889376.

[1.5165750888103101108513650872564 - 1.5165750888103005266522502889376 = 0.0000000000000095841991147983188 =

$=9.5841991147983188\cdot 10^{-15}.$]

Pažiūrėkime kaip konverguoja $y_{n+1}$ į ${\frac {1}{\sqrt {2.3}}}$ iš dešinės pusės, kai $y_{1}$ parenkamas taip, kad $y_{1}>{\frac {1}{\sqrt {2.3}}}.$

Taigi parinkime $y_{1}=0.66.$ Tada

$y_{2}={\frac {y_{1}\left(3-xy_{1}^{2}\right)}{2}}={\frac {0.66\left(3-2.3\cdot 0.66^{2}\right)}{2}}={\frac {0.66\left(3-2.3\cdot 0.4356\right)}{2}}={\frac {0.66\left(3-1.00188\right)}{2}}={\frac {0.66\cdot 1.99812}{2}}={\frac {1.3187592}{2}}=0.6593796.$

[0.65938047339578700471798482054627 - 0.6593796 = 0.00000087339578700471798482054627 $\approx 8.73395787\cdot 10^{-7}.$]

Toliau pagal (4) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(3-xy_{2}^{2}\right)}{2}}={\frac {0.6593796\left(3-2.3\cdot 0.6593796^{2}\right)}{2}}={\frac {0.6593796\left(3-0.999997350861168\right)}{2}}=$

$={\frac {0.6593796\cdot 2.000002649138832}{2}}=$ 1.3187609467881033886272/2 = 0.6593804733940516943136.

[0.65938047339578700471798482054627 - 0.6593804733940516943136 = 0.00000000000173531040438482054627 $\approx 1.7353104\cdot 10^{-12}.$]

Toliau pagal (4) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(3-xy_{3}^{2}\right)}{2}}=$

0.6593804733940516943136*(3 - 2.3*0.6593804733940516943136^2)/2 = 0.65938047339578700471797797024742.

Paklaida lygi

0.65938047339578700471798482054627 - 0.65938047339578700471797797024742 =

= 0.00000000000000000000000685029885 $=6.85029885\cdot 10^{-24}.$

Po 3 iteracijų gavome paklaidą mažesne nei $10^{-23}.$

Šį kartą, parenkant $y_{1}$ iš dešinės, $y_{1}$ reikšmė buvo daug tikslesnė negu anksčiau. Bet teoriškai Niutono metodas (liestinių metodas) bent jau polinomams gali greitai konverguoti tik iš vienos šaknies y pusės. Čia gi greitas konvergavimas į ${\frac {1}{\sqrt {2.3}}}$ vyksta tiek iš kairės pusės, tiek iš dešinės pusės, kas atrodo kaip prieštaravimas liestinių metodo teorijai. Gal greito konvergavimo tik iš vienos pusės taisyklė taikoma tik polinomams ir kai kurioms funkcjoms, o ši funkcija $f(y)={\frac {1}{y^{2}}}-x$ yra išimtis.

Dabar galime rasti a šaknyje:

${\sqrt {2.3}}={\sqrt {a}}={\sqrt {x}}=x\cdot y_{4}=2.3\cdot y_{4}=$ 2.3*0.65938047339578700471797797024742 = 1.5165750888103101108513493315691.

[1.5165750888103101108513650872564 - 1.5165750888103101108513493315691 = 0.0000000000000000000000157556873 =

$=1.57556873\cdot 10^{-23}.$]

• Apskaičiuosime ${\frac {1}{\sqrt {x}}}={\frac {1}{\sqrt {5.8}}}$ ir rasime kam lygus ${\sqrt {5.8}}={\sqrt {x}}.$

Tiksli ${\sqrt {5.8}}$ reikšmė yra tokia:

${\sqrt {5.8}}=$ 2.4083189157584590960256482060757.

O $y={\frac {1}{\sqrt {5.8}}}=$ 0.41522739926869984414235313897857.

Mes paimsime, kad $y_{1}=0.4.$ Tada pagal (4) formulę

$y_{2}={\frac {y_{1}\left(3-xy_{1}^{2}\right)}{2}}={\frac {0.4\left(3-5.8\cdot 0.4^{2}\right)}{2}}={\frac {0.4\left(3-0.928\right)}{2}}={\frac {0.4\cdot 2.072}{2}}={\frac {0.8288}{2}}=0.4144.$

[0.41522739926869984414235313897857 - 0.4144 = 0.00082739926869984414235313897857 $\approx 8.27399268699844\cdot 10^{-4}.$]

Toliau pagal (4) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(3-xy_{2}^{2}\right)}{2}}=$

= 0.4144*(3 - 5.8*0.4144^2)/2 = 0.4152249278464.

[0.41522739926869984414235313897857 - 0.4152249278464 = 0.00000247142229984414235313897857 $\approx 2.471422299844\cdot 10^{-6}.$]

Toliau pagal (4) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(3-xy_{3}^{2}\right)}{2}}=$

= 0.4152249278464*(3 - 5.8*0.4152249278464^2)/2 = 0.41522739924663512944560082817282.

[0.41522739926869984414235313897857 - 0.41522739924663512944560082817282 = 0.00000000002206471469675231080575 =

$\approx 2.206471469675231\cdot 10^{-11}.$]

Toliau pagal (4) formulę rasime $y_{5}$:

$y_{5}={\frac {y_{4}\left(3-xy_{4}^{2}\right)}{2}}=$

= 0.41522739924663512944560082817282*(3 - 5.8*0.41522739924663512944560082817282^2)/2 = 0.41522739926869984414059439797726.

[0.41522739926869984414235313897857 - 0.41522739926869984414059439797726 = 0.00000000000000000000175874100131 =

$=1.75874100131\cdot 10^{-21}.$]

Tada ${\sqrt {5.8}}={\sqrt {x}}=5.8\cdot y_{5}=$ 5.8*0.41522739926869984414059439797726 = 2.4083189157584590960154475082681.

[2.4083189157584590960256482060757 - 2.4083189157584590960154475082681 = 0.0000000000000000000102006978076 =

$=1.02006978076\cdot 10^{-20}.$]

Po 4 iteracijų gavome paklaidą mažesne nei $10^{-19}$ (beveik $10^{-20}$).

Dabar pažūrėsime kaip konverguos $y_{n}$ į $1/{\sqrt {5.8}}$ iš dešinės pusės (kai $y_{1}>1/{\sqrt {5.8}}$). Parenkame $y_{1}=0.42.$ Tada pagal (4) formulę

$y_{2}={\frac {y_{1}\left(3-xy_{1}^{2}\right)}{2}}={\frac {0.42\left(3-5.8\cdot 0.42^{2}\right)}{2}}={\frac {0.42\left(3-1.02312\right)}{2}}={\frac {0.42\cdot 1.97688}{2}}=0.4151448.$

[0.41522739926869984414235313897857 - 0.4151448 = 0.00008259926869984414235313897857 $\approx 8.2599268699844\cdot 10^{-5}.$]

Toliau pagal (4) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(3-xy_{2}^{2}\right)}{2}}=$

= 0.4151448*(3 - 5.8*0.4151448^2)/2 = 0.4152273746236976005632.

[0.41522739926869984414235313897857 - 0.4152273746236976005632 = 0.00000002464500224357915313897857 $\approx 2.4645002243579\cdot 10^{-8}.$]

Toliau pagal (4) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(3-xy_{3}^{2}\right)}{2}}=$

= 0.4152273746236976005632*(3 - 5.8*0.4152273746236976005632^2)/2 = 0.41522739926869765000924208032344.

[0.41522739926869984414235313897857 - 0.41522739926869765000924208032344 = 0.00000000000000219413311105865513 =

$\approx 2.194133111\cdot 10^{-15}.$]

Toliau pagal (4) formulę rasime $y_{5}$:

$y_{5}={\frac {y_{4}\left(3-xy_{4}^{2}\right)}{2}}=$

= 0.41522739926869765000924208032344*(3 - 5.8*0.41522739926869765000924208032344^2)/2 = 0.41522739926869984414235313896118.

[0.41522739926869984414235313897857 - 0.41522739926869984414235313896118 = 0.00000000000000000000000000001739 =

$=1.739\cdot 10^{-29}.$]

Tada ${\sqrt {5.8}}={\sqrt {x}}=5.8\cdot y_{5}=$ 5.8*0.41522739926869984414235313896118 = 2.4083189157584590960256482059748.

[2.4083189157584590960256482060757 - 2.4083189157584590960256482059748 = 0.0000000000000000000000000001009 =

$=1.009\cdot 10^{-28}.$]

Po 4 iteracijų gavome ${\sqrt {5.8}}$ paklaidą mažesne nei $10^{-27}$ (beveik $10^{-28}$).

[2.4083189157584590960256482059748^2 = 5.7999999999999999999999999995139.]

Pastebėsime, kad abiejuose pavyzdžiuose (skaičiuojant ${\sqrt {2.3}}$ ir ${\sqrt {5.8}}$) skaičiuojant $y_{n}$ iš dešinės (kai $y_{1}>1/{\sqrt {2.3}}$ arba $y_{1}>1/{\sqrt {5.8}}$), jau po pirmos iteracijos $y_{n}<1/{\sqrt {x}}$ (t. y. jau $y_{2}<1/{\sqrt {2.3}}$ arba $y_{2}<1/{\sqrt {5.8}}$). Vadinasi atlikus pirmą iteraciją $y_{n}$ skaičiavimai toliau vyksta iš kairės (pagal Niutono metodą).

Štolco teorema

• Štolco teorema

Trisdešimt antro laipsnio šaknies traukimas taikant Niutono metodą

Skaičiuojant ln(z), kai 0.7<z<1.4, reikalingas greitas šešiolikto ar trisdešimt antro laipsnio šaknies traukimo algoritmas (apie tokius natūrinio logaritmo skaičiavimus parašyta čia https://lt.wikibooks.org/wiki/Matematika/Teiloro_eilutė#Natūrinio_logaritmo_liekamojo_nario_įvertinimas ).

Parodysime kaip apskaičiuoti ${\frac {1}{\sqrt[{32}]{z}}}={\frac {1}{z^{1/32}}}$ taikant Niutono metodą. Pažymėkime $y={\frac {1}{z^{1/32}}}.$ Tada turime funkciją

$f(y)={\frac {1}{y^{32}}}-z=0.$

Pagal liestinių metodą

$y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}.\quad (3)$

Surandame f(y) išvestinę:

$f'(y)=({\frac {1}{y^{32}}}-z)'={\frac {-32}{y^{33}}}.$

$f'(y_{n})={\frac {-32}{y_{n}^{33}}}.$

Toliau įstatome $f'(y_{n})$ ir $f(y_{n})$ į (3) formulę:

$y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}=y_{n}-{\frac {{\frac {1}{y_{n}^{32}}}-z}{\frac {-32}{y_{n}^{33}}}}=y_{n}+{\frac {{\frac {y_{n}^{33}}{y_{n}^{32}}}-zy_{n}^{33}}{32}}=$

$=y_{n}+{\frac {y_{n}-zy_{n}^{33}}{32}}=y_{n}+{\frac {y_{n}}{32}}-{\frac {zy_{n}^{33}}{32}}={\frac {33y_{n}}{32}}-{\frac {zy_{n}^{33}}{32}}=$

$=y_{n}\left({\frac {33}{32}}-{\frac {z}{32}}y_{n}^{32}\right)={\frac {y_{n}\left(33-zy_{n}^{32}\right)}{32}}.$

Gavome rekurentinę formulę ${\frac {1}{\sqrt[{32}]{z}}}$ reikšmei skaičiuoti:

$y_{n+1}={\frac {y_{n}\left(33-zy_{n}^{32}\right)}{32}}.\quad (5)$

Pavyzdžiai

• Apskaičiuosime ${\frac {1}{\sqrt[{32}]{z}}}={\frac {1}{\sqrt[{32}]{1.3}}}$ ir rasime kam lygus ${\sqrt[{32}]{1.3}}={\sqrt[{32}]{z}}.$ Prieš tai dar pastebėsime, kad norint apskaičiuoti ${\sqrt[{32}]{z}}$ reikia 5 kartus ištraukti kvadratinę šaknį iš z. Tiksli ${\sqrt[{32}]{1.3}}$ reikšmė yra tokia:

${\sqrt[{32}]{1.3}}=$ 1.0082325861537202021533989768813.

O $y={\frac {1}{\sqrt[{32}]{1.3}}}=$ 0.99183463590963025992221676151568.

Mes paimsime, kad $y_{1}=0.99.$ Tada pagal (5) formulę

$y_{2}={\frac {y_{1}\left(33-zy_{1}^{32}\right)}{32}}={\frac {0.99\left(33-1.3\cdot 0.99^{32}\right)}{32}}=$

= 0.99*(33 - 1.3*0.72498033595785364231768779197961)/32 = 0.99177969711319507382303549411632.

[0.99183463590963025992221676151568 - 0.99177969711319507382303549411632 = 0.00005493879643518609918126739936 $\approx 5.49\cdot 10^{-5}.$]

Toliau pagal (5) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(33-zy_{2}^{32}\right)}{32}}=$

= 0.99177969711319507382303549411632*(33 - 1.3*0.99177969711319507382303549411632^32)/32 = 0.99183458572688580323137637227142.

[0.99183463590963025992221676151568-0.99183458572688580323137637227142=0.00000005018274445669084038924426 $\approx 5.018\cdot 10^{-8}.$]

Toliau pagal (5) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(33-zy_{3}^{32}\right)}{32}}=$

= 0.99183458572688580323137637227142*(33 - 1.3*0.99183458572688580323137637227142^32)/32 = 0.99183463590958836578366677638205.

(0.99183458572688580323137637227142*(33 - 1.3*0.76922952379474351345468696687987)/32 = 0.99183463590958836578366677638205.)

[0.99183463590963025992221676151568 - 0.99183463590958836578366677638205 = 0.00000000000004189413854998513363 =

$\approx 4.1894138549985\cdot 10^{-14}.$]

Toliau pagal (5) formulę rasime $y_{5}$:

$y_{5}={\frac {y_{5}\left(33-zy_{5}^{32}\right)}{32}}=$

= 0.99183463590958836578366677638205*(33 - 1.3*0.99183463590958836578366677638205^32)/32 = 0.9918346359096302599222167323178.

Atėmę, ką tik gautą $y_{5}$ reikšmę iš tikslios $y={\frac {1}{\sqrt[{32}]{1.3}}}$ reikšmės, gauname:

0.99183463590963025992221676151568 - 0.9918346359096302599222167323178 = 0.00000000000000000000000002919788 =

$\approx 2.919788\cdot 10^{-26}.$

Gavome po keturių iteracijų paklaidą mažesnę už $10^{-25}.$

Padalijus vienetą iš $y_{5}$ reikšmės gauname

${\sqrt[{32}]{1.3}}={\sqrt[{32}]{z}}={\frac {1}{y_{5}}}=$

= 1/0.9918346359096302599222167323178 = 1.0082325861537202021533990065619.

[1.0082325861537202021533989768813 - 1.0082325861537202021533990065619 = -0.0000000000000000000000000296806 =

$\approx -2.96806\cdot 10^{-26}.$]

1.0082325861537202021533990065619^32 = 1.3000000000000000000000012246312.

(kad pakelti 1.0082325861537202021533990065619 32 laipsniu reikia ant kalkuliatoriaus 5 kartus paspasuti mygtuką $x^{2}.$)

Toliau apskaičiuosime ln(1.3). Kadangi

$\ln 1.3=-\ln 1.3^{-1}=-32\ln 1.3^{-1/32}=-32\ln {\frac {1}{z^{1/32}}}\approx -32\ln y_{5},$

tai tereikia apskaičiuoti $\ln y_{5}.$ Gauname

$\ln y_{5}=$ ln(0.9918346359096302599222167323178) = ln(1-0.0081653640903697400777832676822), x=-0.0081653640903697400777832676822.

Pagal Teiloro formulę ( https://lt.wikibooks.org/wiki/Matematika/Teiloro_eilutė#Natūrinio_logaritmo_liekamojo_nario_įvertinimas ) apskaičiuosime -32*ln(1+x) = -32*ln(1-0.0081653640903697400777832676822). Imsime n=9. Tada

$\ln(1+x)=x-{\frac {x^{2}}{2}}+{\frac {x^{3}}{3}}-{\frac {x^{4}}{4}}...+(-1)^{n-1}{\frac {x^{n}}{n}}+r_{n+1}(x);\quad (13)$

čia $|-32r_{n+1}(x)|={\Big |}-32{\frac {x^{n+1}}{(1+c)^{n+1}(n+1)}}{\Big |}\;,\;\;-0.008165364<c<0.$

Vietoje x įstatome -0.0081653640903697400777832676822 į (13) formulę ir su n=9, gauname

$\ln(1+x)=\ln(1-0.0081653640903697400777832676822)=$

= -0.0081653640903697400777832676822 - (-0.0081653640903697400777832676822)^2/2 + (-0.0081653640903697400777832676822)^3/3 - (-0.0081653640903697400777832676822)^4/4 + (-0.0081653640903697400777832676822)^5/5 - (-0.0081653640903697400777832676822)^6/6 + (-0.0081653640903697400777832676822)^7/7 - (-0.0081653640903697400777832676822)^8/8 + (-0.0081653640903697400777832676822)^9/9 =

(galima nuo skaičiaus "-0.0081653640903697400777832676822" iki ženklo "=" nukopijuoti visą eilutę su 9 dėmenim ir įstatyti (padaryti "Paste") į Windows 10 kalkuliatorių ir iškart bus gautas apskaičiuotas atsakymas)

= -0.00819888326460909537597654216045.

[-0.00819888326460909537597654216045 - ln(0.9918346359096302599222167323178) = 1.3273686782585303927157997035195e-22.]

Gavome paklaidą $\approx 1.327368678\cdot 10^{-22}<10^{-21}.$

Įvertinsime liekamąjį narį teoriškai.

$|-32r_{9+1}(x)|<{\Big |}-32\cdot {\frac {(-0.008165364)^{9+1}}{(1-0.008165364)^{9+1}(9+1)}}{\Big |}=$

=| -32*((-0.0081653640903697400777832676822)^10)/(10*(1-0.0081653640903697400777832676822)^10) | =

=| -32*((-0.0081653640903697400777832676822)^10)/(10*0.9918346359096302599222167323178^10) | =

(čia galima padaryti "Copy-Paste", viską nukopijavus, išskyrus modulio ženklus "|" ir paspaudus paskui kalkuliatoriuje ženklą "="; tokiu budu bus gautas žemiau esantis atsakymas)

=| -4.5762824470151633013165981309199e-21 | $<10^{-20}.$

Padauginsime anksčiau per Teiloro eilutę gautą atsakymą -0.00819888326460909537597654216045 iš -32, kad gauti ln(1.3). Taigi,

$\ln(1.3)\approx$ -32*(-0.00819888326460909537597654216045) = 0.2623642644674910520312493491344.

Atėmę apskaičiuotą ln(1.3) [apytikslią] reikšmę iš Windows 10 kalkuliatoriaus tikslios reikšmės, gauname paklaidą:

ln(1.3) - 0.2623642644674910520312493491344 = 0.26236426446749105203549598688095 - 0.2623642644674910520312493491344 =

= 4.2466377465543972041664561314341e-21.

Taigi, gavome paklaidą $4.246637746554\cdot 10^{-21}<|-4.576282447\cdot 10^{-21}|<10^{-20}.$

Ir čia viskas atitinka teoriją.

Palyginsime kiek reikia daugybos ir sudeties operacijų, 5 kartus traukiant kvadratinę šaknį su vieno karto trisdešimt antro laipsnio šaknies traukimo formulę. Priminsime 32-o laipsnio šaknies traukimo formulę (tiksliau 1/z^(1/32))

$y_{n+1}={\frac {y_{n}\left(33-zy_{n}^{32}\right)}{32}}.\quad (5)$

Kad $y_{n}$ pakelti 32 laipsniu, reikia tiesiog $y_{n}$ penkis kartus pakelti kvadratu. Tai jau yra 5 daugybos operacijos. Toliau dar reikia $y_{n}$ padauginti iš z. Jau turime 6 daugybos operacijas. Tada reikia $zy_{n}$ atimti iš 33. Gauname vieną atimties operaciją. Ir galiausiai viską padauginti iš $y_{n}$ ir iš 1/32. Gavosi iš viso 6+2=8 daugybos operacijos ir viena atimties operacija. Truputi tiksliau [nei mes] iš procesoriaus FPU lentelės pasirinkus $y_{1},$ reikšmę užteks 3 iteracijų, kad gauti 18 skaitmenų tikslumą. Todėl kad ištraukti 32-o laipsnio šaknį reikia 8*3=24 daugybos operacijos ir 1*3=3 atimties operacijų. Iš viso 24+3=27 operacijos (nors šiuolaikiniai procesoriai turi FMA (fuse multiply add) instrukcija, kuri leidžia per vieną operaciją sudauginti du skaičius ir pridėti trečią skaičių). Daugybos ir sudeties operacijos x87 FPU skaičiuojamos panašiu greičiu. Čia nereikia vienos dalybos iteracijos, kad atlikti $1/y_{4},$ nes skaičiuojant ln(1+x) be to apseinama kaip galima matyti iš aukščiau.

Išnagrinėsime kiek reikia daugybos ir sudeties operacijų, jei penkis kartus traukti kvadratinę šaknį iš skaičiaus z. Kvadtatinės šaknies (tiksliau 1/x^(1/2)) traukimo vienos iteracijos formulė yra tokia:

$y_{n+1}={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}.\quad (4)$

Kad pakelti $y_{n}$ kvadratu reikia vienos daugybos operacijos. Paskui dar vienos daugybos operacijos reikia, kad padauginti $y_{n}^{2}$ iš x. Jau turime dvi daugybos operacijas. Čia yra dar viena atimties operacija. Dar viena daugybos operacija reiškinio skliaustuose iš $y_{n}$ ir dar viena daugybos operacija visko padauginimas iš 0.5. Taigi, iš viso reikia 2+2=4 daugybos operacijų ir vienos atimties operacijos. Procesoriui gana tiksliai pasirinkus pirmą artinį $y_{1},$ užteks 3 iteracijų. Atlikus 3 iteracijas dar reikia padauginti gautą rezultatą iš x, kad gauti x^(1/2). Šiuolaikiniai procesoriai skaičiuojantis pagal Goldschmidt’o algorithm'ą moka suskaičiuoti vienu metu ir 1/x^(1/2) ir x^(1/2), todėl sutaupoma viena daugybos operacija (palyginus su Niutono metodu), bet nesutaupomi tranzistoriai. Trims iteracijoms reikia 4*3=12 daugybos operacijų ir 1*3=3 atimties operacijos. Taigi, kvadratinę šaknį reikia traukti 5 kartus. Tai iš viso reikia 12*5=60 daugybos operacijų ir 3*5=15 atimties operacijų. Gaunasi, kad iš viso reikia 60+15=75 operacijų, kad penkis kartus ištraukti kvadratinę šaknį (ir tai paskui pritaikyti ln(1+x) skaičiavimuose).

Taigi, kad apskaičiuoti 1/z^(1/32) reikia 27 operacijų. Gauname, kad penkių kartų skaičiavimui 1/x^(1/2) reikia 75/27=2.77778 karto daugiau operacijų negu 1/z^(1/32) skaičiavimui. Ne tiek jau ir daug ir todėl greičiausiai nieks neeikvos papildomų tranzistorių $y_{1}=1/z^{1/32}$ lentelėms bei nekurs naujų instrukcijų, kad greičiau apskaičiuoti ln(x). Tiesa, ln(x), kaip parodyta čia ( https://lt.wikibooks.org/wiki/Matematika/Teiloro_eilutė#e_pakelta_x_liekamojo_nario_įvertinimas ), taikomas skaičiaus kelimui bet kokiu laipsniu. Bet skaičiaus kelimas bet kokiu laipsniu mažai kur reikalingas. Dar yra ir kitokių gana konkurencingų ln(x) skaičiavimo būdų.

Geriau paskaičius, Goldschmidt’o algorithm'as gali būti, kad sutaupo kiekvienoje iteracijoje daugybą iš 0.5. Tai iš viso trijose iteracijose sutaupomos trys daugybos operacijos (arba trys dalybos operacijos iš 2; dvejetainėje sistemoje dalyba iš dviejų yra pastumimas visų bitų į dešinę per vieną bito poziciją; bet kas greičiau ar daugyba, ar pastumimas tai yra dar neaišku).

Va čia gerai paaiškinta apie Goldschmidt’o algorithm'ą: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Goldschmidt’s_algorithm

• Apskaičiuosime ${\frac {1}{\sqrt[{32}]{z}}}={\frac {1}{\sqrt[{32}]{0.75}}}$ ir rasime kam lygus ${\sqrt[{32}]{0.75}}={\sqrt[{32}]{z}}.$ Tiksli ${\sqrt[{32}]{0.75}}$ reikšmė yra tokia:

${\sqrt[{32}]{0.75}}=$ 0.99105022504172762555415281377767.

O $y={\frac {1}{\sqrt[{32}]{0.75}}}=$ 1.0090305967670766121721228310046.

Mes paimsime, kad $y_{1}=1.009.$ Tada pagal (5) formulę

$y_{2}={\frac {y_{1}\left(33-zy_{1}^{32}\right)}{32}}={\frac {1.009\left(33-0.75\cdot 1.009^{32}\right)}{32}}=$

= 1.009*(33 - 0.75*1.3320401627616656376316134288889)/32 = 1.0090305814634409227725711418028.

[1.0090305967670766121721228310046 - 1.0090305814634409227725711418028 = 0.0000000153036356893995516892018 $\approx 1.53\cdot 10^{-8}.$]

Toliau pagal (5) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(33-zy_{2}^{32}\right)}{32}}=$

= 1.0090305814634409227725711418028*(33 - 0.75*1.0090305814634409227725711418028^32)/32 = 1.0090305967670727824366475789071.

[1.0090305967670766121721228310046 - 1.0090305967670727824366475789071 = 0.0000000000000038297354752520975 $\approx 3.8297\cdot 10^{-15}.$]

Toliau pagal (5) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(33-zy_{3}^{32}\right)}{32}}=$

= 1.0090305967670727824366475789071*(33 - 0.75*1.0090305967670727824366475789071^32)/32 = 1.0090305967670766121721228307648.

Atėmę gautą $y_{4}$ reikšmę iš tikslios $y={\frac {1}{\sqrt[{32}]{0.75}}}$ reikšmės, gauname:

1.0090305967670766121721228310046 - 1.0090305967670766121721228307648 = 0.0000000000000000000000000002398 $\approx 2.398\cdot 10^{-28}.$

Tik po 3 iteracijų paklaida yra mažesnė nei $10^{-27}.$ Taip yra todėl, nes parinkome $y_{1}=1.009$ su paklaida

1.0090305967670766121721228310046 - 1.009 = 0.0000305967670766121721228310046 $\approx 3.0596767\cdot 10^{-5}.$

Na o,

${\sqrt[{32}]{z}}={\sqrt[{32}]{0.75}}\approx 1/y_{4}=$ 1/1.0090305967670766121721228307648 = 0.99105022504172762555415281401323.

[0.99105022504172762555415281377767 - 0.99105022504172762555415281401323 = -0.00000000000000000000000000023556 =

$\approx -2.3556\cdot 10^{-28}.$ Čia kažkodėl iš viso, su pirmu nuliu prieš tašką, yra 33 skaitmenys, o ne 32 skaitmenys kaip prieš tai esančiame $y_{4}$ paklaidos atsakyme.]

Toliau apskaičiuosime ln(0.75). Kadangi

$\ln 0.75=-\ln 0.75^{-1}=-32\ln 0.75^{-1/32}=-32\ln {\frac {1}{z^{1/32}}}\approx -32\ln y_{4},$

tai tereikia apskaičiuoti $\ln y_{4}.$ Gauname

$\ln y_{4}=$ ln(1.0090305967670766121721228307648) = ln(1+0.0090305967670766121721228307648), x=0.0090305967670766121721228307648.

Pagal Teiloro formulę apskaičiuosime -32*ln(1+x) = -32*ln(1+0.0090305967670766121721228307648). Imsime n=9. Tada

$\ln(1+x)=x-{\frac {x^{2}}{2}}+{\frac {x^{3}}{3}}-{\frac {x^{4}}{4}}...+(-1)^{n-1}{\frac {x^{n}}{n}}+r_{n+1}(x);\quad (13)$

čia $|-32r_{n+1}(x)|={\Big |}-32{\frac {x^{n+1}}{(1+c)^{n+1}(n+1)}}{\Big |},\;\;0<c<0.009030596767.$

Vietoje x įstatome 0.0090305967670766121721228307648 į (13) formulę ir su n=9, gauname

$\ln(1+x)=\ln(1+0.0090305967670766121721228307648)=$

= 0.0090305967670766121721228307648 - 0.0090305967670766121721228307648^2/2 + 0.0090305967670766121721228307648^3/3 - 0.0090305967670766121721228307648^4/4 + 0.0090305967670766121721228307648^5/5 - 0.0090305967670766121721228307648^6/6 + 0.0090305967670766121721228307648^7/7 - 0.0090305967670766121721228307648^8/8 + 0.0090305967670766121721228307648^9/9 =

= 0.00899006476411815398283337192369.

[0.00899006476411815398283337192369 - ln(1.0090305967670766121721228307648) =

= 0.00899006476411815398283337192369 - 0.00899006476411815398247559369963 =

= 0.00000000000000000000035777822406.]

Gavome paklaidą $3.5777822406\cdot 10^{-22}<10^{-21}.$

Įvertinsime liekamąjį narį teoriškai.

$|-32r_{9+1}(x)|<{\Big |}-32\cdot {\frac {0.009030596767^{9+1}}{(1+0)^{9+1}(9+1)}}{\Big |}=$

=| -32*(0.0090305967670766121721228307648^10)/10 | =

(čia galima padaryti "Copy-Paste", viską nukopijavus į kalkuliatoriaus laukelį kur rodo skaičius, išskyrus modulio ženklus "|" ir paspaudus paskui kalkuliatoriuje ženklą "="; tokiu budu bus gautas žemiau esantis atsakymas)

=| -1.1542888083111938172050915484026e-20 | $<10^{-19}.$

Padauginsime anksčiau per Teiloro eilutę gautą atsakymą ln(1+x)=0.00899006476411815398283337192369 iš -32, kad gauti ln(0.75). Taigi,

$\ln(0.75)\approx$ -32*0.00899006476411815398283337192369 = -0.28768207245178092745066790155808.

Atėmę apskaičiuotą ln(0.75) [apytikslią] reikšmę iš Windows 10 kalkuliatoriaus tikslios reikšmės, gauname paklaidą:

ln(0.75) - (-0.28768207245178092745066790155808) = -0.28768207245178092743921900599383 + 0.28768207245178092745066790155808 =

= 0.00000000000000000001144889556425.

Taigi, gavome paklaidą $1.144889556425\cdot 10^{-20}<|-32r_{n+1}(x)|=|-1.1542888083\cdot 10^{-20}|<10^{-19}.$

Viskas atitinka teoriją.

Goldschmidt dalybos algoritmas

https://en.wikipedia.org/wiki/Division_algorithm#Goldschmidt_division

Goldschmidt'o algoritmas veikia taip, kad skaitiklis ir vardiklis dauginamas iš bendro skaičiaus F_i. Skaičius F_i pasirinktas taip, kad dauginant vardiklį kiekvieną kartą iš F_i, vardiklis konverguoja į vienetą. Taip kelis kartus dauginant skaitiklį ir vardiklį iš F_i skaitiklis konverguoja į Q=N/D; čia N yra dalinys (skaitiklis), D yra daliklis (vardiklis). Formulė Q skaičiuoti yra tokia:

$Q={\frac {N}{D}}{\frac {F_{1}}{F_{1}}}{\frac {F_{2}}{F_{2}}}{\frac {F_{\ldots }}{F_{\ldots }}}.$

Žingsniai Goldschmidt'o dalybai yra:

1. Generuoti F_i .
2. Padauginti dalinį ir daliklį iš F_i .
3. Jeigu daliklis pakankmai arti prie 1, pateikti dalmenį Q, priešintu atveju, kartoti nuo žingsnio 1.

Pirma reikia daliklį (vardiklį) padauginti iš 0.1, kad jisai būtų intervale nuo 0 iki 1 (0<D<1). Kiekvienas F_i randamas iš D_i (D₁ = D):

$F_{i+1}=2-D_{i}.$

Dauginant dalinį $N_{i}$ ir daliklį $D_{i}$ iš faktoriaus F_i duoda $N_{i+1}$ ir $D_{i+1}$:

${\frac {N_{i+1}}{D_{i+1}}}={\frac {N_{i}}{D_{i}}}{\frac {F_{i+1}}{F_{i+1}}}.$

Po pakankmo iteracijų skaičiaus k, $Q=N_{k}$.

Pavyzdžiai

• Apskaičiuosime Q=N/D=3/7. Padauginame 7 iš 0.1: $D_{1}=0.1\cdot 7=0.7.$ Tada

$F_{i+1}=2-D_{i},$

$F_{1+1}=F_{2}=2-D_{1}=2-0.7=1.3.$

Toliau padauginame $N_{1}=3$ ir $D_{1}=0.7$ iš $F_{2}$:

$N_{2}=N_{1}\cdot F_{2}=3\cdot 1.3=3.9,$

$D_{2}=D_{1}\cdot F_{2}=0.7\cdot 1.3=0.91.$

Toliau apskaičiuojame $F_{3}$:

$F_{3}=2-D_{2}=2-0.91=1.09.$

Toliau padauginame $N_{2}$ ir $D_{2}$ iš $F_{3}$:

$N_{3}=N_{2}\cdot F_{3}=3.9\cdot 1.09=4.251,$

$D_{3}=D_{2}\cdot F_{3}=0.91\cdot 1.09=0.9919.$

Toliau apskaičiuojame $F_{4}$:

$F_{4}=2-D_{3}=$ 2 - 0.9919 = 1.0081.

Toliau padauginame $N_{3}$ ir $D_{3}$ iš $F_{4}$:

$N_{4}=N_{3}\cdot F_{4}=$ 4.251*1.0081 = 4.2854331,

$D_{4}=D_{3}\cdot F_{4}=$ 0.9919*1.0081 = 0.99993439.

Toliau apskaičiuojame $F_{5}$:

$F_{5}=2-D_{4}=$ 2 - 0.99993439 = 1.00006561.

Toliau padauginame $N_{3}$ ir $D_{3}$ iš $F_{4}$:

$N_{5}=N_{4}\cdot F_{5}=$ 4.2854331*1.00006561 = 4.285714267265691,

$D_{5}=D_{4}\cdot F_{5}=$ 0.99993439*1.00006561 = 0.9999999956953279.

Taigi $Q=N_{5}\cdot 0.1=$ 4.285714267265691*0.1 = 0.4285714267265691. Tiksli Q=3/7 reikšmė yra tokia:

Q = 3/7 = 0.42857142857142857142857142857143.

[0.42857142857142857142857142857143 - 0.4285714267265691 = 0.00000000184485947142857142857143 $\approx 1.84485947\cdot 10^{-9}.$]

Po 4 iteracijų gavome paklaidą mažesnę nei $10^{-8}.$

Dar padarysime penktą įteraciją, apskaičiuodami $F_{6},\;$ $N_{6}$ ir $D_{6}$:

$F_{6}=2-D_{5}=$ 2 - 0.9999999956953279 = 1.0000000043046721;

$N_{6}=N_{5}\cdot F_{6}=$ 4.285714267265691*1.0000000043046721 = 4.2857142857142856348705633349211,

$D_{6}=D_{5}\cdot F_{6}=$ 0.9999999956953279*1.0000000043046721 = 0.99999999999999998146979811148159.

Taigi $Q=N_{6}\cdot 0.1=$ 4.2857142857142856348705633349211*0.1 = 0.42857142857142856348705633349211.

[0.42857142857142857142857142857143 - 0.42857142857142856348705633349211 = 0.00000000000000000794151509507932 =

$\approx 7.941515\cdot 10^{-18}.$]

Po 5 iteracijų gavome paklaidą mažesnę nei $10^{-17}.$

Goldschmidt'o algoritmas kvadratinės šankies traukimui

https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Goldschmidt's_algorithm

https://lt.wikibooks.org/wiki/Sekos_riba#Greitas_šaknies_iš_skaičiaus_a_traukimo_būdas

https://www.quora.com/How-do-computers-calculate-square-roots

Goldschmidt'o algoritmas galėtų būti taikomas naujesniuose procesoriuose, kai FPU gali atlikti kai kuriuos skaičiavimus vienu metu.

Pirmas būdas užrašyti Goldschmidt'o algoritmą yra toks (skaičius S yra tas skaičius iš kurio reikia ištraukti kvadratinę šaknį):

$b_{0}=S$

$Y_{0}\approx 1/{\sqrt {S}}$ (typically using a table lookup)

$y_{0}=Y_{0}$

$x_{0}=Sy_{0}$

ir daryti iteracijas

$b_{n+1}=b_{n}Y_{n}^{2}$

$Y_{n+1}={\frac {1}{2}}(3-b_{n+1})$

$x_{n+1}=x_{n}Y_{n+1}$

$y_{n+1}=y_{n}Y_{n+1}$

kol $b_{i}$ yra pakankamai arti 1, arba daryti fiksuotą (tam tikrą) kiekį iteracijų. Iteracijos konverguoja prie

$\lim _{n\to \infty }x_{n}={\sqrt {S}}$

ir

$\lim _{n\to \infty }y_{n}=1/{\sqrt {S}}.$

Pastebėsime, kad galima neskačiuoti $x_{n}$ ir $y_{n}$ jeigu abu yra reikalingi. Galima tik paskaičiuoti $x_{n}=Sy_{n}$ pačioje skaičiavimų pabaigoje (tada $y_{n}$ vis tiek reikia skaičiuoti).

Antras būdas, naudojant fused multiply-add operacijas, prasideda

$y_{0}\approx 1/{\sqrt {S}}$ (typically using a table lookup)

$x_{0}=Sy_{0}$

$h_{0}={\tfrac {1}{2}}y_{0}$

ir iteruoti

$r_{n}=0.5-x_{n}h_{n}\;\;$ [ $r_{0}=0.5-Sy_{0}\cdot {\tfrac {1}{2}}y_{0}=0.5(1-Sy_{0}^{2}),\;\;\;$ $r_{1}=0.5-x_{1}h_{1}=0.5-({\tfrac {3}{2}}x_{0}-0.5x_{0}\cdot Sy_{0}^{2})\cdot {\tfrac {1}{2}}({\tfrac {3}{2}}y_{0}-0.5\cdot Sy_{0}^{3})$ ]

$x_{n+1}=x_{n}+x_{n}r_{n}\;\;$ [ $x_{1}=x_{0}+x_{0}r_{0}=x_{0}+x_{0}\cdot 0.5(1-Sy_{0}^{2})={\tfrac {3}{2}}x_{0}-0.5x_{0}\cdot Sy_{0}^{2}$ ]

$h_{n+1}=h_{n}+h_{n}r_{n}\;\;$ [ $h_{1}=h_{0}+h_{0}r_{0}={\tfrac {1}{2}}y_{0}+{\tfrac {1}{2}}y_{0}\cdot 0.5(1-Sy_{0}^{2})={\tfrac {1}{2}}(y_{0}+y_{0}\cdot 0.5(1-Sy_{0}^{2}))={\tfrac {1}{2}}({\tfrac {3}{2}}y_{0}-0.5\cdot Sy_{0}^{3})$ ]

kol $r_{i}$ yra pakankamai arti prie 0, arba daryti fiksuotą skaičių iteracijų. Tai konverguoja prie

$\lim _{n\to \infty }x_{n}={\sqrt {S}},$

ir

$\lim _{n\to \infty }2h_{n}=1/{\sqrt {S}}.$

• Pirmo būdo pavyzdys. Apskaičiuosime ${\sqrt {S}}={\sqrt {3.01}}.$

Sprendimas.

$b_{0}=S=3.01,$

$Y_{0}={\frac {1}{\sqrt {3}}}\approx 1/{\sqrt {S}}\;\;\;$ (parinkome apytikslią S reikšmę lygia 3),

$y_{0}=Y_{0}={\frac {1}{\sqrt {3}}}=$ 0.57735026918962576450914878050196,

$x_{0}=Sy_{0}=3.01\cdot {\frac {1}{\sqrt {3}}}=$ 3.01*0.57735026918962576450914878050196 = 1.7378243102607735511725378293109.

Darome pirmą iteraciją:

$b_{1}=b_{0}Y_{0}^{2}=3.01\cdot ({\frac {1}{\sqrt {3}}})^{2}=3.01\cdot {\frac {1}{3}}=$ 3.01/3 = 1.0033333333333333333333333333333,

$Y_{1}={\frac {1}{2}}(3-b_{1})=$ 0.5*(3 - 1.0033333333333333333333333333333) = 0.99833333333333333333333333333335,

$x_{1}=x_{0}Y_{1}=$ 1.7378243102607735511725378293109*0.99833333333333333333333333333335 = 1.7349279364103389285872502662621,

$y_{1}=y_{0}Y_{1}=$ 0.57735026918962576450914878050196*0.99833333333333333333333333333335 = 0.5763880187409763882349668658678.

Darome antrą iteraciją:

$b_{2}=b_{1}Y_{1}^{2}=$ = 1.0033333333333333333333333333333*0.99833333333333333333333333333335^2 = 0.99999167592592592592592592592593,

$Y_{2}={\frac {1}{2}}(3-b_{2})=$ 0.5*(3 - 0.99999167592592592592592592592593) = 1.000004162037037037037037037037,

$x_{2}=x_{1}Y_{2}=$ 1.7349279364103389285872502662621*1.000004162037037037037037037037 = 1.7349351572446668586552910067157,

$y_{2}=y_{1}Y_{2}=$ 0.5763880187409763882349668658678*1.000004162037037037037037037037 = 0.57639041768925809257650864010487.

Atėmę $x_{2}$ iš ${\sqrt {S}}={\sqrt {3.01}},$ gauname paklaidą:

1.7349351572897472412324994277 - 1.7349351572446668586552910067157 = 4.5080382577208420984281695490435e-11 $\approx 4.508038\cdot 10^{-11}.$

Atėmę $y_{2}$ iš ${\frac {1}{\sqrt {3.01}}},$ gauname paklaidą:

0.57639041770423496386461775006644 - 0.57639041768925809257650864010487 = 1.4976871288109109961569101491839e-11 $\approx 1.497687\cdot 10^{-11}.$

Darome trečią iteraciją:

$b_{3}=b_{2}Y_{2}^{2}=$ = 0.99999167592592592592592592592593*1.000004162037037037037037037037^2 = 0.99999999994803219891278740696279,

$Y_{3}={\frac {1}{2}}(3-b_{3})=$ 0.5*(3 - 0.99999999994803219891278740696279) = 1.0000000000259839005436062965186,

$x_{3}=x_{2}Y_{3}=$ 1.7349351572446668586552910067157*1.0000000000259839005436062965186 = 1.734935157289747241230742381434,

$y_{3}=y_{2}Y_{3}=$ 0.57639041768925809257650864010487*1.0000000000259839005436062965186 = 0.57639041770423496386403401376544.

Atėmę $x_{3}$ iš ${\sqrt {S}}={\sqrt {3.01}},$ gauname paklaidą:

1.7349351572897472412324994277 - 1.734935157289747241230742381434 = 0.000000000000000000001757046266 $\approx 1.757046266\cdot 10^{-21}.$

Atėmę $y_{3}$ iš ${\frac {1}{\sqrt {3.01}}},$ gauname paklaidą:

0.57639041770423496386461775006644 - 0.57639041770423496386403401376544 = 0.000000000000000000000583736301 $\approx 5.83736301\cdot 10^{-22}.$

• Antro būdo pavyzdys. Apskaičiuosime ${\sqrt {S}}={\sqrt {3.01}}.$

Sprendimas.

$y_{0}={\frac {1}{\sqrt {3}}}\approx {\frac {1}{\sqrt {S}}}\;\;\;$ (parinkome apytikslią S reikšmę lygia 3),

$x_{0}=Sy_{0}=3.01\cdot {\frac {1}{\sqrt {3}}}=$ 3.01*0.57735026918962576450914878050196 = 1.7378243102607735511725378293109,

$h_{0}={\tfrac {1}{2}}y_{0}={\tfrac {1}{2}}{\frac {1}{\sqrt {3}}}={\frac {1}{2{\sqrt {3}}}}=$ 0.28867513459481288225457439025098.

Darome pirmą iteraciją:

$r_{0}=0.5-x_{0}h_{0}=0.5-3.01\cdot {\frac {1}{\sqrt {3}}}\cdot {\frac {1}{2{\sqrt {3}}}}=0.5-3.01\cdot {\frac {1}{6}}=$ 0.5 - 3.01/6 = -0.00166666666666666666666666666667,

$x_{1}=x_{0}+x_{0}r_{0}$ = 1.7378243102607735511725378293109 + 1.7378243102607735511725378293109*(-0.00166666666666666666666666666667) = 1.734927936410338928587250266262,

$h_{1}=h_{0}+h_{0}r_{0}=$ 0.28867513459481288225457439025098 + 0.28867513459481288225457439025098*(-0.00166666666666666666666666666667) = 0.28819400937048819411748343293389.

Darome antrą iteraciją:

$r_{1}=0.5-x_{1}h_{1}=$ 0.5 - 1.734927936410338928587250266262*0.28819400937048819411748343293389 = 4.1620370370370370370370370552249e-6,

$x_{2}=x_{1}+x_{1}r_{1}$ = 1.734927936410338928587250266262 + 1.734927936410338928587250266262*4.1620370370370370370370370552249e-6 = 1.7349351572446668586552910067157,

$h_{2}=h_{1}+h_{1}r_{1}=$ 0.28819400937048819411748343293389 + 0.28819400937048819411748343293389*4.1620370370370370370370370552249e-6 = 0.28819520884462904628825432005244.

Atėmę $x_{2}$ iš ${\sqrt {S}}={\sqrt {3.01}},$ gauname paklaidą:

1.7349351572897472412324994277 - 1.7349351572446668586552910067157 = 0.0000000000450803825772084209843 $\approx 4.508038\cdot 10^{-11}.$

Atėmę $h_{2}$ iš ${\frac {1}{2{\sqrt {S}}}}={\frac {1}{2{\sqrt {3.01}}}},$ gauname paklaidą:

0.28819520885211748193230887503322 - 0.28819520884462904628825432005244 = 0.00000000000748843564405455498078 $\approx 7.488435644\cdot 10^{-12}.$

Darome trečią iteraciją:

$r_{2}=0.5-x_{2}h_{2}=$ 0.5 - 1.7349351572446668586552910067157*0.28819520884462904628825432005244 = 2.5983900543606296518569637165611e-11,

$x_{3}=x_{2}+x_{2}r_{2}$ = 1.7349351572446668586552910067157 + 1.7349351572446668586552910067157*2.5983900543606296518569637165611e-11 = 1.7349351572897472412307423814339,

$h_{3}=h_{2}+h_{2}r_{2}=$ 0.28819520884462904628825432005244 + 0.28819520884462904628825432005244*2.5983900543606296518569637165611e-11 = 0.28819520885211748193201700688271.

Atėmę $x_{3}$ iš ${\sqrt {S}}={\sqrt {3.01}},$ gauname paklaidą:

1.7349351572897472412324994277 - 1.7349351572897472412307423814339 = 0.0000000000000000000017570462661 $\approx 1.7570462661\cdot 10^{-21}.$

Atėmę $h_{3}$ iš ${\frac {1}{2{\sqrt {S}}}}={\frac {1}{2{\sqrt {3.01}}}},$ gauname paklaidą:

0.28819520885211748193230887503322 - 0.28819520885211748193201700688271 = 0.00000000000000000000029186815051 $\approx 2.9186815051\cdot 10^{-22}.$

Kiek reikia kilobaitų pastovios atminties, tipo ROM (Read Only Memory), kad išsaugoti lentelę su apytikslėm reikšmėm visų skaičių (be eksponentės) nuo 1 iki 9.99999999999999? Kai tikslumas yra 1, tai reikia 9 arba 10 apytikslių reikšmių, kurios bus naudojamos Goldschmidt'o algoritme. Kadangi FPU naudoja 64 bitų skaičius (be eksponentės; su eksponente naudoja 80 bitų skaičius), tai 10*64 = 640 (bitų) = 80 B.

Jeigu lentelės tikslumas yra 0.1, tai reikia 99 arba 100 apytikslių reikšmių $y_{0},$ kurios bus naudojamos Goldschmidt'o algoritme. Tada reikia ROM atminties, kuri yra tokia pat greita, o gal net greitesnė nei pats greičiausias cache'as, va tiek baitų: 100*64 = 6400 (bitų) = 800 B.

Jeigu lentelės tikslumas yra 0.01 (kaip mūsų atveju), tai reikia 999 arba 1000 apytikslių reikšmių $y_{0},$ kurios bus naudojamos Goldschmidt'o algoritme. Tada reikia ROM atminties, va tiek baitų: 1000*64 = 64000 (bitų) = 8000 B = 8 KB. Reikia 8 kilobaitų ROM (ar panašios) atminties. Tada pagal pirmus 3 skaitmenis [skaičiaus S], kurie gali būti nuo 1.00 iki 9.99 [su žingsniu 0.01], bus nukreiptas signalas į ROM ir duos tas ROM apytikslią $y_{0}\approx 1/{\sqrt {S}}$ reikšmę, kurią toliau panaudos skaičiavimuose. Kadangi 2^10=1024, tai galimas daiktas, kad pirmi 10 bitų nusako skaičius nuo 1.00 iki 9.99. Tie 10 pirmų skaičiaus S bitų per tranzistorių logiką nukreipia, priklausomai nuo skaičiaus S pirmų 3 skaitmenų į beveik 1024 skaičius (į 999 tiksliau), sudarytus iš 64 bitų informacijos kievienas.

Šiaip lentelėje iš 999 (ar 1024) skaičių galima tikriausiai naudoti ne 64 bitų skaičius, o pavyzdžiui, 14 bitų skaičius (2^14 = 16384), kurie nusako $y_{0}\approx 1/{\sqrt {S}}$ pagal skaičiaus S reikšmę. Nes skaičius S vis tiek turi tikslumą iš 3 skaitmenų (nuo 1.00 iki 9.99), pagal kuriuos renkamas skaičius $y_{0}.$ Todėl ir $y_{0}$ tikslumas negali būti daug didesnis negu skaičiaus S tikslumas iš 3 skaitmenų. Todėl suteikus visiems lentelės skaičiams $y_{0}\approx 1/{\sqrt {S}}$ reikšmes iš 4 skaitmenų (nuo 1.000 iki 9.999), to pilnai pakaks. Tada gaunasi, kad tokiai lentelei reikia apie 4 kartus mažiau baitų ROM atminties, nei 64 bitų skaičiams (vietoje 8 KB, reikia apie 2 KB). Pats skaičius S yra 64 bitų tikslumo arba apytiksliai 18 dešimtainių skaitmenų tikslumo, bet pirmi 3 jo skaitmenys per tranzistorių logiką nurodys skaičiaus $y_{0}\approx 1/{\sqrt {S}}$ reikšmę iš 4 skaitmenų. Tada tas 4 skaitmenų $y_{0}$ bus pakrautas į kokį nors registrą, užpildydamas 14 pirmų bitų (gal tik pradžioje palikdamas patį pirmą bitą + ar - ženklo reikšmei) to registro, o likusius 50 ar 49 bitus to registro palikus būti nuliais. Tada tą $y_{0}$ toliau FPU galės naudoti skaičiavimuose pagal Goldschmidt'o algoritmą.

Goldschmidt'o algoritmo įrodymas

Goldschmidt'o algoritmo įrodymas paimtas iš https://www.quora.com/How-do-computers-calculate-square-roots

Iš karto pasakysime, kad Goldschmidt'o algoritmas pagrįstas Newton-Raphson metodu, kuris paremtas Niutono metodu (liestiniu metodu, kai ieškoma funkcijos šaknis).

Kaip pateikta čia [ https://lt.wikibooks.org/wiki/Sekos_riba#Greitas_šaknies_iš_skaičiaus_a_traukimo_būdas ], Newton-Raphson metodas randa $x_{i+1}\approx {\frac {1}{\sqrt {n}}},$ iteruojant pagal formulę:

$x_{i+1}={\tfrac {1}{2}}x_{i}(3-nx_{i}^{2}).$

Skaičiaus n reikšmę mes, žinoma, žinome tiksliai ir norime surasti ${\frac {1}{\sqrt {n}}}.$

Dabar mes galime pažvelgti į Goldschmidt'o algoritmą, kuris skaičiuoja ${\frac {1}{\sqrt {n}}}$ ir ${\sqrt {n}}$ vienu metu. Aišku būtų galima iš pradžių apskaičiuoti ${\frac {1}{\sqrt {n}}}$ ir gautą rezultatą padauginti iš n ir tokiu budu gauti ${\sqrt {n}},$ bet tokiu budu bus sugaišta laiko vienai papildomai daugybai. O Goldschmidt'o algoritmas turėtų naudotis [modernaus] procesoriaus FPU galimybe skaičiuoti dvi užduotis vienu metu, naudojant tipo dvigubą pipline'ą.

Tegu $b_{0}=n.$ Jeigu mes galėsime surasti skaičius $Y_{i},$ tokius, kad $b_{n}=b_{0}Y_{0}^{2}Y_{1}^{2}Y_{2}^{2}\cdots Y_{n-1}^{2}$ artėja prie 1 (kai skaičiaus b indeksas n didėja), tada $y_{n}=Y_{0}Y_{1}Y_{2}\cdots Y_{n-1}$ artės prie ${\frac {1}{\sqrt {b_{0}}}}={\frac {1}{\sqrt {n}}}$ ir $x_{n}=b_{0}Y_{0}Y_{1}Y_{2}\cdots Y_{n-1}$ artės prie ${\sqrt {b_{0}}}={\sqrt {n}}.$

Šitos skaičių $Y_{i}$ serijos yra Newton-Raphson iteracijų žingsniai:

$b_{i}=b_{i-1}Y_{i-1}^{2}$

[ $b_{1}=b_{0}Y_{0}^{2},\;\;\;$ $b_{2}=b_{1}Y_{1}^{2}=b_{0}Y_{0}^{2}Y_{1}^{2},\;\;\;$ $b_{3}=b_{2}Y_{2}^{2}=b_{0}Y_{0}^{2}Y_{1}^{2}Y_{2}^{2},\;\;\;$ ... ]

$Y_{i}={\frac {1}{2}}(3-b_{i})$

[ $Y_{0}$ yra imamas iš lentelės ir lygus apytiksliai ${\frac {1}{\sqrt {n}}}$ reikšmei.

$Y_{1}={\frac {1}{2}}(3-b_{1})={\frac {1}{2}}(3-b_{0}Y_{0}^{2})={\frac {1}{2}}(3-nY_{0}^{2})$ ]

Ir tada mes galime skaičiuoti šaknį iš n taip:

$x_{0}=nY_{0}$

$x_{i}=x_{i-1}Y_{i}$

[ $x_{1}=x_{0}Y_{1}=nY_{0}Y_{1},\;\;\;$ $x_{2}=x_{1}Y_{2}=nY_{0}Y_{1}Y_{2},\;\;\;$ $x_{3}=x_{2}Y_{3}=nY_{0}Y_{1}Y_{2}Y_{3},\;\;\;$ ... (nepamirškim, kad $b_{0}=n$) ]

O šaknis iš skaičiaus 1/n skaičiuojama taip ($y_{i}\approx {\frac {1}{\sqrt {n}}}$):

$y_{0}=Y_{0}$

$y_{i}=y_{i-1}Y_{i}$

[ $y_{1}=y_{0}Y_{1}=Y_{0}Y_{1},\;\;\;$ $y_{2}=y_{1}Y_{2}=Y_{0}Y_{1}Y_{2},\;\;\;$ $y_{3}=y_{2}Y_{3}=Y_{0}Y_{1}Y_{2}Y_{3},\;\;\;$ ...

Be to, iš aukščiau $Y_{1}={\frac {1}{2}}(3-nY_{0}^{2}),$ todėl

$y_{1}=y_{0}Y_{1}=Y_{0}Y_{1}=Y_{0}{\frac {1}{2}}(3-nY_{0}^{2}),$

kas lygu Newton-Raphson iteracijai $x_{i+1}={\tfrac {1}{2}}x_{i}(3-nx_{i}^{2}),$ kurioje $x_{i+1}\approx {\frac {1}{\sqrt {n}}}.$ Tik dabar pas mus raidė x pakeista raide y. ]

[$b_{i}=b_{i-1}Y_{i-1}^{2},\;\;$ $b_{1}=b_{0}Y_{0}^{2},\;$ nes $x_{0}=nY_{0},\;\;$ $y_{0}=Y_{0}\;$ ir t. t., vadinasi, $b_{i}=x_{i-1}y_{i-1}.$]

Kadangi $b_{i}$ artėja prie 1 (kai i artėja prie begalybės), tai iš to, kad $b_{i}=b_{0}Y_{0}^{2}Y_{1}^{2}Y_{2}^{2}\cdots Y_{i-1}^{2}$ ir iš to, kad $y_{i}=Y_{0}Y_{1}Y_{2}\cdots Y_{i}$ bei $x_{i}=b_{0}Y_{0}Y_{1}Y_{2}\cdots Y_{i},$ gauname, kad $b_{i}=x_{i-1}y_{i-1}.$ Todėl mes iš tikrųjų niekada neturime sekti ir skaičiuoti $b_{i}$:

$Y_{i}={\frac {1}{2}}(3-b_{i})=1+{\frac {1}{2}}(1-b_{i})=1+{\frac {1}{2}}(1-x_{i-1}y_{i-1}).$

Ir dabar mums taipogi nereikia sekti ir skaičiuoti $Y_{i}$:

$x_{i}=x_{i-1}Y_{i}=x_{i-1}\left(1+{\frac {1}{2}}(1-x_{i-1}y_{i-1})\right)=x_{i-1}+x_{i-1}{\frac {1}{2}}(1-x_{i-1}y_{i-1}),$

$y_{i}=y_{i-1}Y_{i}=y_{i-1}\left(1+{\frac {1}{2}}(1-x_{i-1}y_{i-1})\right)=y_{i-1}+y_{i-1}{\frac {1}{2}}(1-x_{i-1}y_{i-1}).$

Čia yra kai kurie sudetingi užrašymai, kuriuos galima palengvinti. Skaičiaus ${\frac {1}{\sqrt {n}}}$ apytikslia reikšme iš lentelės pažymėsime raide Y. Tada

$x_{0}=nY,$

$y_{0}=Y.$

Tada toliau iteruojame:

$r_{i}={\frac {1}{2}}(1-x_{i}y_{i}),$

$x_{i+1}=x_{i}+x_{i}r_{i},$

$y_{i+1}=y_{i}+y_{i}r_{i}.$

Nors dalyjimas iš 2 yra pigus, bet mes to galime išvengti užrašę $h_{i}={\tfrac {1}{2}}y_{i}$ vietoje $y_{i}.$ Tai ir yra Goldschmidt'o algoritmas.

Sakykime, Y yra apytiksli reikšmė skaičiaus ${\frac {1}{\sqrt {n}}}.$ Pažymime:

$x_{0}=nY,$

$h_{0}={\frac {Y}{2}}.$

Tada iteruoti:

$r_{i}={\frac {1}{2}}-x_{i}h_{i},$

$x_{i+1}=x_{i}+x_{i}r_{i},$

$h_{i+1}=h_{i}+h_{i}r_{i}.$

Tada $x_{i}$ konverguoja į ${\sqrt {n}},\;$ o $h_{i}$ konverguoja į ${\frac {1}{2{\sqrt {n}}}}.$

[ Iš tikro,

$r_{0}={\frac {1}{2}}-x_{0}h_{0}={\frac {1}{2}}-x_{0}{\tfrac {1}{2}}y_{0}={\frac {1}{2}}(1-x_{0}y_{0}),$

$x_{1}=x_{0}+x_{0}r_{0},$

$h_{1}=h_{0}+h_{0}r_{0}={\tfrac {1}{2}}y_{0}+{\tfrac {1}{2}}y_{0}r_{0}={\tfrac {1}{2}}(y_{0}+y_{0}r_{0})={\tfrac {1}{2}}y_{1}.$

Toliau:

$r_{1}={\frac {1}{2}}-x_{1}h_{1}={\frac {1}{2}}-x_{1}{\tfrac {1}{2}}y_{1}={\frac {1}{2}}(1-x_{1}y_{1}),$

$x_{2}=x_{1}+x_{1}r_{1},$

$h_{2}=h_{1}+h_{1}r_{1}={\tfrac {1}{2}}y_{1}+{\tfrac {1}{2}}y_{1}r_{1}={\tfrac {1}{2}}(y_{1}+y_{1}r_{1})={\tfrac {1}{2}}y_{2}.$

Dar toliau:

$r_{2}={\frac {1}{2}}-x_{2}h_{2}={\frac {1}{2}}-x_{2}{\tfrac {1}{2}}y_{2}={\frac {1}{2}}(1-x_{2}y_{2}),$

$x_{3}=x_{2}+x_{2}r_{2},$

$h_{3}=h_{2}+h_{2}r_{2}={\tfrac {1}{2}}y_{2}+{\tfrac {1}{2}}y_{2}r_{2}={\tfrac {1}{2}}(y_{2}+y_{2}r_{2})={\tfrac {1}{2}}y_{3}.$

Tuo baigtas įrodymas, kad nereikia kai kurioje vietoje dauginti iš 0.5 (skaičiuojant $r_{i}$) ir bus skaičiuojama viskas teisingai. ]

Taigi, kuom Goldschmidt'o algoritmas pranašesnis už Newton-Raphson metodą? Newton-Raphson metodas randa $y_{i+1}\approx {\frac {1}{\sqrt {n}}},$ iteruojant pagal formulę:

$y_{i+1}={\tfrac {1}{2}}y_{i}(3-ny_{i}^{2}).$

Dar jeigu reikia žinoti ${\sqrt {n}},$ tai reikia $y_{i+1}$ padauginti iš n.

Goldschmidt'o gi algoritmas su naujoviška procesoriaus architektūra gali skaičiuoti $x_{i}$ ir $h_{i}$ vienu metu lygiagrečiai, pirma apskaičiavus $r_{i-1}.$ Tiesa, pabaigoje, norint žinoti ${\frac {1}{\sqrt {n}}}$ reikšmę, reikią $h_{i}$ padauginti iš 2. Ir tai galbūt galima atlikti greičiau su šiftinimo visų bitų viena pozicija į kairę operacija negu paprastai padauginus iš n pačioje pabaigoje taikant Newton-Raphson metodą. Be to, taikant Goldschmidt'o algoritmą, mes turėsime ir ${\sqrt {n}}$ reikšmę ir ${\frac {1}{\sqrt {n}}}$ reikšmę (po padauginimo iš 2 ar šiftinimo į kairę). O taikant Newton-Raphson metodą, jei padauginsime $y_{i+1}$ iš n, tai turėsime ${\sqrt {n}},$ bet neturėsime ${\frac {1}{\sqrt {n}}}$ (jeigu kokiuose nors skaičiavimuose reikia jų abiejų). Nebent nukopijuoti iš pat pradžių gautą $y_{i+1}\approx {\frac {1}{\sqrt {n}}}$ reikšmę rezervui, kas gali užimti kažkiek laiko.

Tiesa, Goldschmidt'o algoritmui reikia, greitesniems skaičiavimams, turėti dvi lenteles skaičių. Vienoje lentelėje skaičiai Y, kurie naudojami $x_{0}=nY$ skaičiavimui, o kitoje lentelėje reikia turėti skačius Y/2, kurie naudojami $h_{0}={\frac {Y}{2}}$ gavimui. Bet už tai Newton-Raphson metode kiekvienoje iteracijoje reikia dauginti iš 0.5.

Jeigu neskaičiuoti, kad abiais metodais naudojama FMA operacija, kad išvengti painiavos, tai jau turint $x_{0}$ ir $h_{0},$ iš iteracijų

$r_{i}={\frac {1}{2}}-x_{i}h_{i},$

$x_{i+1}=x_{i}+x_{i}r_{i},$

$h_{i+1}=h_{i}+h_{i}r_{i};$

matome, kad $r_{i}$ skaičiavimui reikia vienos daugybos ir vienos atimties operacijos. $x_{i+1}$ skaičiavimui reikia vienos daugybos ir vienos sudeties operacijos. $h_{i+1}$ skaičiavimui reikia vienos daugybos ir vienos sudeties operacijos. Bet $x_{i+1}$ ir $h_{i+1}$ skaičiuojami lygiagrečiai, todėl iš viso reikia per vieną iteraciją Goldschmidt'o algoritme 2 daugybų ir dviejų sumos operacijų. Iš viso reikia 4 operacijų (nes daugyba ir sudetis užima panašiai ciklų).

Newton-Raphson metode be FMA operacijų, sprendžiant iš iteracijos

$y_{i+1}={\tfrac {1}{2}}y_{i}(3-ny_{i}^{2});$

matome, kad vienoje iteracijoje yra 4 daugybos operacijos ir viena sumos operacija. Taigi, iš viso 5 operacijos. Skirtumas akivaizdus (5 > 4), jeigu po tiek pat iteracijų abiais metodais skaičiuojant gaunamas tas pats tikslumas.

FMA ar FMA3 instrukcijų panaudojimas galėtų Goldschmidt'o algoritmui suteikti didesni pranašumą prieš Newton-Raphson algoritmą (iki 2 kartų; 2 operacijos vs 4 operacijos idealiu atveju...).

Bet Newton-Raphson algoritme irgi galima padaryti viena daugybos operacija mažiau. Tereikia iš pat pradžių padauginti skaičių n iš 0.5,

m=0.5*n.

Tada

$y_{i+1}={\tfrac {1}{2}}y_{i}(3-ny_{i}^{2})=y_{i}({\tfrac {3}{2}}-{\tfrac {1}{2}}ny_{i}^{2})=y_{i}({\tfrac {3}{2}}-my_{i}^{2}).$

• Goldschmidt'o algoritmo pavyzdys. Apskaičiuosime ${\sqrt {n}}={\sqrt {6.309}}$ ir ${\frac {1}{\sqrt {n}}}={\frac {1}{\sqrt {6.309}}}.$

Sprendimas.

$y_{0}={\frac {1}{\sqrt {6.3}}}\approx {\frac {1}{\sqrt {n}}}\;\;\;$ (parinkome apytikslią n reikšmę lygia 6.3; be to, 6.309-6.3=0.009<0.01),

$x_{0}=ny_{0}=6.309\cdot {\frac {1}{\sqrt {6.3}}}=$ 6.309*0.39840953644479787998960572656437 = 2.5135657654302298248544225288946,

$h_{0}={\tfrac {1}{2}}y_{0}={\tfrac {1}{2}}\cdot {\frac {1}{\sqrt {6.3}}}=$ 0.5*0.39840953644479787998960572656437 = 0.19920476822239893999480286328219.

Darome pirmą iteraciją:

$r_{0}=0.5-x_{0}h_{0}=$ 0.5 - 2.5135657654302298248544225288946*0.19920476822239893999480286328219 = -7.1428571428571428571428571428367e-4,

$x_{1}=x_{0}+x_{0}r_{0}$ = 2.5135657654302298248544225288946 + 2.5135657654302298248544225288946*(-7.1428571428571428571428571428367e-4) = 2.5117703613120653749795265128025,

$h_{1}=h_{0}+h_{0}r_{0}=$ 0.19920476822239893999480286328219 + 0.19920476822239893999480286328219*(-7.1428571428571428571428571428367e-4) = 0.19906247910224008360909228980842.

[ ${\sqrt {6.309}}-x_{1}=$ 2.5117722826721374584282937869002 - 2.5117703613120653749795265128025 = 0.0000019213600720834487672740977. ]

[ ${\frac {1}{2{\sqrt {6.309}}}}-h_{1}=$ 0.19906263137360417327851432769854 - 0.19906247910224008360909228980842 = 0.00000015227136408966942203789012 $\approx 1.52271364\cdot 10^{-7}.$ ]

Darome antrą iteraciją:

$r_{1}=0.5-x_{1}h_{1}=$ 0.5 - 2.5117703613120653749795265128025*0.19906247910224008360909228980842 = 7.6494169096209912536443148786058e-7,

$x_{2}=x_{1}+x_{1}r_{1}$ = 2.5117703613120653749795265128025 + 2.5117703613120653749795265128025*7.6494169096209912536443148786058e-7 = 2.5117722826699328655134995747654,

$h_{2}=h_{1}+h_{1}r_{1}=$ 0.19906247910224008360909228980842 + 0.19906247910224008360909228980842*7.6494169096209912536443148786058e-7 = 0.19906263137342945518414166862938.

[${\sqrt {6.309}}-x_{2}=$ 2.5117722826721374584282937869002 - 2.5117722826699328655134995747654 = 0.0000000000022045929147942121348 $\approx 2.204592914794\cdot 10^{-12}.$]

[${\frac {1}{2{\sqrt {6.309}}}}-h_{2}=$ 0.19906263137360417327851432769854 - 0.19906263137342945518414166862938 = 0.00000000000017471809437265906916 $\approx 1.7471809437\cdot 10^{-13}.$]

Darome trečią iteraciją:

$r_{2}=0.5-x_{2}h_{2}=$ 0.5 - 2.5117722826699328655134995747654*0.19906263137342945518414166862938 = 8.7770413345269442665565631474709e-13,

$x_{3}=x_{2}+x_{2}r_{2}$ = 2.5117722826699328655134995747654 + 2.5117722826699328655134995747654*8.7770413345269442665565631474709e-13 = 2.5117722826721374584282908844297,

$h_{3}=h_{2}+h_{2}r_{2}=$ 0.19906263137342945518414166862938 + 0.19906263137342945518414166862938*8.7770413345269442665565631474709e-13 = 0.19906263137360417327851409767236.

[${\sqrt {6.309}}-x_{3}=$ 2.5117722826721374584282937869002 - 2.5117722826721374584282908844297 = 0.0000000000000000000000029024705 $=2.9024705\cdot 10^{-24}.$]

[${\frac {1}{2{\sqrt {6.309}}}}-h_{3}=$ 0.19906263137360417327851432769854 - 0.19906263137360417327851409767236 = 0.00000000000000000000000023002618 $=2.3002618\cdot 10^{-25}.$]

[ $x_{3}^{2}=$ 2.5117722826721374584282908844297^2 = 6.30899999999999999999998541931 .]

• Newton-Raphson algoritmo pavyzdys. Apskaičiuosime ${\sqrt {n}}={\sqrt {6.309}}$ ir ${\frac {1}{\sqrt {n}}}={\frac {1}{\sqrt {6.309}}}.$

Sprendimas.

$y_{0}={\frac {1}{\sqrt {6.3}}}\approx {\frac {1}{\sqrt {n}}}\;\;\;$ (parinkome apytikslią n reikšmę lygia 6.3; be to, 6.309-6.3=0.009<0.01).

$y_{0}={\frac {1}{\sqrt {6.3}}}=$ 0.39840953644479787998960572656437.

m = 0.5*n = 0.5*6.309 = 3.1545.

$y_{i+1}=y_{i}({\tfrac {3}{2}}-my_{i}^{2}).$

Darome pirmą iteraciją:

$y_{1}=y_{0}({\tfrac {3}{2}}-my_{0}^{2})=$ 0.39840953644479787998960572656437*(1.5 - 3.1545*0.39840953644479787998960572656437^2) = 0.39812495820448016721818457961683.

[${\frac {1}{\sqrt {6.309}}}-y_{1}=$ 0.39812526274720834655702865539709 - 0.39812495820448016721818457961683 = 0.00000030454272817933884407578026 $\approx 3.04542728\cdot 10^{-7}.$]

Darome antrą iteraciją:

$y_{2}=y_{1}({\tfrac {3}{2}}-my_{1}^{2})=$ 0.39812495820448016721818457961683*(1.5 - 3.1545*0.39812495820448016721818457961683^2) = 0.39812526274685891036828333725875.

[${\frac {1}{\sqrt {6.309}}}-y_{2}=$ 0.39812526274720834655702865539709 - 0.39812526274685891036828333725875 = 0.00000000000034943618874531813834 $\approx 3.494361887\cdot 10^{-13}.$]

Darome trečią iteraciją:

$y_{3}=y_{2}({\tfrac {3}{2}}-my_{2}^{2})=$ 0.39812526274685891036828333725875*(1.5 - 3.1545*0.39812526274685891036828333725875^2) = 0.39812526274720834655702819534471.

[${\frac {1}{\sqrt {6.309}}}-y_{3}=$ 0.39812526274720834655702865539709 - 0.39812526274720834655702819534471 = 0.00000000000000000000000046005238 $=4.6005238\cdot 10^{-25}.$]

${\sqrt {n}}\approx x_{3}=n\cdot y_{3}=$ 6.309*0.39812526274720834655702819534471 = 2.5117722826721374584282908844298.

$x_{3}^{2}=$ 2.5117722826721374584282908844298^2 = 6.3089999999999999999999854193105.

Norint, pavyzdžiui, greitai apskaičiuoti reiškinį

0.39812526274685891036828333725875*(1.5 - 3.1545*0.39812526274685891036828333725875^2)

reikia tokį reiškinį:

0.39812526274685891036828333725875*(1.5 - 3.1545*0.39812526274685891036828333725875^2) =

nukopijuoti ir įdėti į Windows 10 kalkuliatoriaus langelį, kur rodomi skaičiai (padarius "Paste"). Tada bus gautas rezultatas:

0.39812526274720834655702819534471.

Palyginsime $x_{3}$ reikšmes skaičiuoant dviais algoritmais.

2.5117722826721374584282908844297

2.5117722826721374584282908844298

2.5117722826721374584282937869002

Pirmoji reikšmė yra Goldschmidt'o algoritmo, o antroji reikšmė yra Newton-Raphson algoritmo. Trečioji reikšmė yra tiksli skaičiaus ${\sqrt {6.309}}$ reikšmė.

Vertimas dvejetainės sistemos slankaus kablelio skaičius į dešimtainės sistemos

https://saturncloud.io/blog/how-to-convert-binary-floating-point-to-decimal-the-correct-algorithm/

https://sandbox.mc.edu/~bennet/cs110/flt/ftod.html

https://kyledewey.github.io/comp122-fall17/lecture/week_2/floating_point_interconversions.html

Pavyzdžiui, skaičius 101.001 dvejetainėje sistemoje į dešimtainę sistemą verčiamas taip:

1*2^2 + 0*2^1 + 1*2^0 + 0*2^(-1) + 0*2^(-2) + 1*2^(-3) = 5.125.

Bet double precision ar single precision formate, kuriame ir skaičiuoja procesoriaus FPU (skaičiuoja double extended precision formate, bet verčia į double arba single precision), dvejetainės sistemos skaičiai saugomi taip:

$(s)x.xxxxxxxxxxx\cdot 2^{(s)xxxxx},$

kur (s) reiškia minuso ar pliuso ženklą, o x reiškia vieną bitą. Pirmas bitas prieš kablelį (tašką) visada yra vienetas. Pavyzdžiui, skaičius:

(1)1.001011 * 2^((0)1010)

dvejetainėje sistemoje iškoduojamas taip:

-1.001011 * 2^(1010) = -1.001011 * 2^(1*2^3 + 0*2^2 +1*2^1 + 0*2^0) = -1.001011 * 2^(8+2) = -1.001011 * 2^10 =

= -10010110000 (kablelis paslenkamas 10 pozicijų i dešinę). Toks skaičius reiškia:

- 1*2^10 + 0*2^9 + 0*2^8 + 1*2^7 + 0*2^6 + 1*2^5 + 1*2^4 + 0*2^3 + 0*2^2 + 0*2^1 + 0*2^0 = -848.

O ką daryti, jeigu reikia išversti iš dvejetainės sistemos į dešimtainę sistemą ne sveikąjį skaičių, o slankaus kablelio skaičių (kai po padauginimo iš eksponentės vis tiek gaunamas slankaus kablelio skaičius)?

Pavyzdžiui, čia https://sandbox.mc.edu/~bennet/cs110/flt/ftod.html duotas single precision (32 bit) formate skaičius:

1|01000110|01101011000000000000000

Tai yra skaičius:

-1.01101011000000000000000 * 2^(01000110 - 127) = -1.01101011000000000000000 * 2^(-57) = 1.41796875*2^(-57) =

= 9.8391347153059527386176341678947e-18;

01000110 = 70

-1.01101011000000000000000 = -1.41796875

Vadinasi, jei procesoriaus FPU skaičiuodamas, gavo kažkokį skaičių kaip 9.8391347153059527386176341678947e-18 (po pavertimo į dešimtainę sistemą), tai jis to skaičiaus pirmus tris skaitmenis (9.83) panaudoja žiūrėjimui į lentelę $y_{0}$ reikšmės, kurią toliau taikys Newton-Raphson algoritmo skaičiavimuose.

Vadinasi, būtinai reikia versti iš dvejetainės sistemos į dešimtainę sistemą double precision skaičius, kad pažiūrėti į lentelę kvadratinės šaknies traukimui, nes kitaip nežinosi kokie pirmi to skaičiaus skaitmenys (po pirmo skaičio visada yra kablelis [dešimtainėje sistemoje, jau išvertus]). Neišvertus iš dvejetainės sistemos į dešimtainę, negalima pasakyti kokie yra pirmi to skaičiaus skaitmenys dešimtainėje sistemoje. O dvejetainėje sistemoje, pavyzdžiui dvi vienodos mantisos reikšmės 1.0101101100, bet su skirtinga eksponente reikš visiškai skirtingus skaičius (kurie skirsis ne tik kablelio pozija ar nuliu skaičiumi dešimtainėje sistemoje, bet jų skaitmenys bus visiškai skirtingi)! Todėl nepavertus dvejetainės sistemos slankaus kablelio skaičių į dešimtainę sistemą, negalima sukurti jokios lentelės apytikslės $y_{0}$ reikšmės gavimui, kuri taikoma pirmoje iteracijoje Newton-Raphson algoritme.

Šituo atveju naudinga susidomėti Binary-coded decimal [ https://en.wikipedia.org/wiki/Binary-coded_decimal ]. Binary coded decimal (BCD) sudaro 4 bitai. Normaliai 4 bitai turi 16 reikšmių ir gali saugoti skačius nuo 0 iki 15. Bet BCD atveju 4 bitai saugo tik 10 reikšmių, skaičius nuo 0 iki 9. Visi nešiojami kalkuliatoriai skaičiuoja BCD formate. Taip skaičiai užima daugiau bitų, bet daug lengviau ir greičiau, apskaičiavus, tuos bitus (skaičius ar skaitmenys dešimtainėje sistemoje) parodyti vartotojui.

Taigi, yra dvi galimybės procesoriaus FPU. Arba vardan šaknies traukimo (ir dalybos; beje dalybą galima gauti Newton-Raphson algoritmui baigus skaičiavus ${\frac {1}{\sqrt {n}}},$ rezultatą pakėlus kvadratu, tada gausis 1/n) skaičiuoti viską BCD formate, kas yra sunkiau ir reikalauja daugiau tranzistorių, arba versti dvejetainius skaičius į Binary-coded decimals, pagal kuriuos žiūrėti lentelėje $y_{0}$ reikšmės.

Pastebėsime, kad skaičiuojant 1.41796875*2^(-57),

[1.41796875*2^(-57) = 9.8391347153059527386176341678947e-18]

reikia kažkaip pereiti į BCD sistemą ir per ją skaičiuoti 2^(-57). Tai išeina, vis tiek reikės skaičiuoti BCD sistemoje dvejeto kelimą kažkokiu laipsniu. Tai geriau jau tada apskritai skaičiuoti viską Binary-coded decimal sistemoje.

Binary-coded decimal sistemoje 64 bitai gali išsaugoti 16 dešimtainių skaitmenų (64/4=16). Double extended precision turi 80 bitų. Tai 64 bitai galėtų būti mantisai. Vienas bitas pliuso ar minuso ženklui (mantisos) ir vienas bitas pliuso ar minuso ženklui eksponentės. Tada lieka 80-66=14 bitų. Iš jų for BCD paėmę 12 bitų, gauname 3 dešimtainius skaitmenis eksponentės.

Intel i486 procesorius turi 32 bitų su RAM atmintim data bus. O Pentium ir velesni procesoriai turi 64-bit data bus. Tai Jeigu FPU skaičiuoja su BCD, tai tie skaičiai užima 64+12+2=78 bitus (kiekvienas). Tam neužtenka 64 bitų data bus ir vieno tokio skaičiaus persiuntimas iš atminties RAM ar į atmintį RAM reikalauja dviejų tranzakcijų. i486 procesoriams ir senesniems su 32 bitų data bus reikia trijų tranzakcijų nusiuntimui/paėmimui į/iš RAM (32*3=96 > 78). Gaunas mažesnis pranašumas...

Pavyzdžiui, Free Pascal versijos 2.6 ir 2.4 skaičiavo viską su 15 dešimtainių skaitmenų tikslumu. O naujesnės versijos (kaip 3.2.0) skaičiuoja su 17 dešimtainių skaitmenų tikslumu. Gali būti, kad anksčiau šešioliktas skaitmuo buvo imamas apvalinimui, todėl jo nerodė, o dabar septinioliktas skaitmuo imamas iš lempos (šešioliktas ir septinioliktas skaitmenys dažnai būna neteisingi Free Pascal rezultatuose; oficialiai tai aiškinama dėl vertimo iš dvejetainės į dešimtainę sistemą subtilibių ir čia pačioje teorijoje greičiausiai klaidų nėra). Šešioliktas skaitmuo galėtų būti neteisingas dažnai (jei skaičiuojama su BCD), nes nėra po jo jokio skaitmens kuris garantuotų teisingą tikslumą (ypač po daug visokių operacijų).

Papildymas. Keliant dvejetą teigiamu laipsniu butų galima, duoto pavyzdžiui 1.41796875*2^(-57) atveju, jei eksponentė teigiama, 2^57 neversti iškart į dešimtainę sistemą, o apskaičiavus

2^57 = 144,115,188,075,855,872 = 1.44115188075855872 * 10^17,

t. y. išvertus 2^57 dvejetainėje sistemoje (kas yra 01000110 - 127) į dešimtainę (kas nėra taip jau sunku), jau turėti normalų skaičių dešimtainėje sistemoje

[1.44115188075855872 * 10^17 = 144115188075855872]

1.41796875 * 2^57 = 1.41796875 * 1.44115188075855872*10^17 = 2.04350833091936256 *10^17 = 204350833091936256 (Čia po kablelio yra kaip tik 17 skaitmenų).

Bet jeigu dvejetą reikia kelti neigiamu laipsniu (2^(-57)), tai tada reikia arba 57 kartus šiftinti visus bitus į dešinę arba versti tokį skaičių į dešimtainę sistemą, kas priveda prie dalybos operacijų, kurioms reikia lentelės ir kurios užima daug ciklų ir grįžtame prie tos pačios užduoties nuo kurios pradėjome...

Pašiftinus bitus vargu ar ką išspręs, nes turbut vis tiek tai nepadės konvertavimui į dešimtainę sistemą... Be to, ką tik pateiktam pavyzdyje, reikia atlikti daugybą 1.41796875 * 1.44115188075855872 ir gaunasi, kad vis tiek reikės šituos skaičius sudauginti per Binary-coded decimals.

Papildymas 2. Parodysime kaip pasiekiama lentelės reikšmė ($\approx {\frac {1}{\sqrt {n}}}$) esanti ROM'e, priklausomai nuo skaičiaus n.

Pirmas pavyzdys. Aukščiau gavome tokį skaičių:

9.8391347153059527386176341678947e-18 = $9.839134715305952\cdot 10^{-18}=n.$

Imame pirmus 3 jo skaitmenys: 9.83. Apskaičiuojame ${\sqrt {e}}={\sqrt {10^{-18}}}=10^{-9}.$ Tada ${\frac {1}{\sqrt {e}}}={\frac {1}{10^{-9}}}=10^{9}.$ Lentelėje, kuri yra ROM'e yra jau apskaičiuota $\approx {\frac {1}{\sqrt {9.83}}}$ reikšmė, kuri lygi ${\frac {1}{\sqrt {9.83}}}\approx 0.3189=3.189\cdot 10^{-1}.$ Protingiau eksponentės reikšmės lentelėje nesaugoti, nes ji visada lygi $10^{-1}$ skaičiams n tik su pirmais trim skaitmenim (pagal kuriuos imami skaičiai $y_{0}\approx {\frac {1}{\sqrt {n}}}$), kurie yra nuo 1.00 iki 9.99.

Todėl iš lentelės pagal apytikslę skaičiaus n reikšmę 9.83, gavę 3.189, ją padauginame iš skaičiaus $10^{-1}$ ir iš skaičiaus n eksponentės, iš kurios ištraukta šaknis, tai yra, dar padauginame iš 10^9. Tada gauname:

$y_{0}=3.189\cdot 10^{-1}\cdot 10^{9}=3.189\cdot 10^{8}.$

Ši $y_{0}$ reikšmė ir taikoma pirmoje iteracijoje Newton-Raphson algoritme.

Antras pavyzdys. Pagal skaičių $n=1.4541864575691000\cdot 10^{17},$ rasime $y_{0}$ reikšmę.

${\sqrt {10^{17}}}={\sqrt {10}}\cdot {\sqrt {10^{16}}}=3.162277660168379\cdot 10^{8}.$

Tada

${\frac {1}{\sqrt {10^{17}}}}={\tfrac {1}{3.162277660168379}}\cdot 10^{-8}=0.3162277660168379\cdot 10^{-8}=3.162277660168379\cdot 10^{-9}.$

ROM'o lentelėje yra apytikslė ${\frac {1}{\sqrt {1.45}}}$ reikšmė (be eksponentės $10^{-1}$).

${\frac {1}{\sqrt {1.45}}}\approx 0.8304=8.304\cdot 10^{-1}.$

Tada gauname apytikslę $y_{0}$ reikšmę:

$y_{0}=8.304\cdot 10^{-1}\cdot 3.162277660168379\cdot 10^{-9}=26.25955369003821\cdot 10^{-10}=2.625955369003821\cdot 10^{-9}.$

$y_{0}\approx {\frac {1}{\sqrt {1.4541864575691000\cdot 10^{17}}}}.$

ROM (Read only memory) lentelėje šiame pavyzdyje skaičiui n (su triais pirmais skaitmenim 1.45 ir be eksponentės 10^17) lentelės reikšmė yra 8.304 (kuri yra be eksponentės $10^{-1}$).

Taigi, bet koks skaičius iš lentelės turi 4 dešimtainius skaitmenis. BCD sistemoje vienas toks skaitmuo užima 4 bitus. Vadinasi vienas skaičius iš lentelės užima 4*4=16 bitų. Iš viso yra 999 ar 1000 tokių skaičių lentelėje. Vadinasi visa lentelė ROM'e užima 1000*16=16000 bitų arba 2000 baitų (~2 KB, nes 1 kilobaitas yra 1024*8=8192 bitai, berods).

Kaip pastebėjome, iš eksponentės, kurios laipsnis yra lyginis skaičius, reikia traukti kvadratinę šaknį. Kaip rašyta aukščiau, eksponentei pakaktų trijų skaitmenų (nuo -999 iki +999, t. y. nuo $10^{-999}$ iki $10^{999}$). Dėl to, kad šaknis traukiama tik iš lyginių skaičių, eksponentė šaknyje gali turėti apie 1000/2=500 skirtingų reikšmių. Tam reikia apie vieno kilobaito ROM lentelės. Arba tiksliau 3 BCD skaitmenys užima 3*4=12 bitų. Tada reikia 12*500=6000 bitų ROM lentelės. Pavyzdžiui, jei eksponentė yra $10^{976},$ tai po šaknies traukimo jos reikšmė bus $10^{976/2}=10^{488}.$ Pagal skaičių 976 bus gautas skaičius 488 iš lentelės, kuri yra ROM atmintyje.

Taigi, šaknies traukimas Binary-coded decimal sistemoje yra ilgesnis procesas, nei atrodo iš pirmo žvilgsnio.

Šaknies traukimo vienoje iteracijoje

$y_{i+1}=y_{i}({\tfrac {3}{2}}-my_{i}^{2})\;\;\;$ (m=0.5*n),

yra 4 operacijos. Jeigu daugyba ir sumos operacija užima po 3 ciklus, tai gaunasi 4*3=12 ciklų vienai iteracijai. O iš viso yra 3 iteracijos, tai gaunasi 12*3=36 ciklai.

Čia [ https://lt.wikibooks.org/wiki/Aptarimas:Matematika/Teiloro_eilutė#Skaičiaus_e_skaičiavimas_ir_ne_tik ] gauta, kad šaknies traukimo operacijai reikia apie 29 ciklų. Tada išeitų, kad 4.16 GHz procesoriaus FPU turi gigantiškas ROM lenteles, kurios gali užimti net kelis megabaitus, kad iškart pateikti gana tikslią $y_{0}$ reikšmę ir skaičiuoti tik su dviem ar net viena iteracija (kvadratinę šaknį ar vieneto dalyba iš jos).

Update 1. Po poros dienų pagalvojimų priėjau prie išvados, kad visgi turėtų būti įmanoma sukūrti dvejetainėje sistemoje nedidelę $y_{0}$ lentelę Newton-Raphson algoritmui ROM atmintyje.

Pradžiai, reikėtų žinoti, kad kai dvejetainės sistemos skaičius paverčiamas į dešimtainės sistemos skaičių, tai pirmas [jau dešimtainis] skaitmuo visada yra vienetas. O visi likusieji skaitmenys po kablelio gali būti nuo 0 iki 9. Toks skaičius yra visada mažesnis už 2 (be eksponentės).

Pirmas pavyzdys.

1|01000110|01101011000000000000000

Tai yra skaičius:

-1.01101011000000000000000 * 2^(01000110 - 127) = -1.01101011000000000000000 * 2^(-57) = 1.41796875*2^(-57) =

= 9.8391347153059527386176341678947e-18;

01000110 = 70

-1.01101011000000000000000 = -1.41796875

Čia 127 yra 8 bitų skaičius. Turime

00000001 yra 1,

00000010 yra 2,

00000100 yra 4,

00001000 yra 8,

00010000 yra 16,

00100000 yra 32,

01000000 yra 64,

10000000 yra 128,

11111111 yra 255,

01111111 yra 127.

Tada (pagal 10 - 01 = 01 [2-1=1], ir 10 - 11 =00-01= -01 [2-3=-1]):

01000110

-

01111111

=

01000110 - 01111111 = -00111001. Iš čia atradau, kad dvejetainėje sistemoje atimtis yra paprasčiausiai Exlusive OR arba kitaip XOR. Padarius XOR operaciją su dviais dvejetainiais skaičiais, gaunamas dvejetainis skaičius, kuriame bitai yra vienetai tik ten, kur pirmo ir antro dvejetainio skaičiaus atitinkamų bitų suma lygi vienetui. O jeigu atitinkamų bitų suma tų dviejų skaičių lygi 0 arba 2, tai tada ten gaunama bito reikšmė lygi nuliui.

Dar pridėsiu, kad jeigu pirmas 8 bitų skaičius (iš kurio atimamas antras) didesnis už antrą, tai gaunamas teigiamas 8 bitų skaičius, o jeigu antras 8 bitų skaičius didesnis už pirmą 8 bitų skaičių, tai atlikus XOR dar prirašomas minuso ženklas.

Taigi skaičius -00111001 lygus

-(0*2^7 + 0*2^6 + 1*2^5 + 1*2^4 + 1*2^3 + 0*2^2 + 0*2^1 + 1*2^0) = -(0 + 0 + 32 + 16 + 8 + 0 + 0 + 1) = -57.

Dar pasakysime, kad šitame skaičiuje 1|01000110|01101011000000000000000 (kuris yra single precision arba 32 bitų skaičius) pirmas bitas yra pliuso arba minuso ženklas. 0 reiškia +, o 1 reiškia -. Toliau skaičius 01000110 yra eksponentė, iš kurios visada atimamas skaičius 127; kai eksponentė lygi nuliui, gausime 0-127=-127, t. y. gausime neigiamą eksponentę; o atėme iš maksimalios 8 bitų eksponentės reikšmės, kuri yra 255, skaičių 127, gausime 255-127=128. Vadinasi eksponentė gali būti nuo -127 iki 128.

($2^{-127}=$ 5.8774717541114375398436826861112e-39,

$2^{128}=$ 3.4028236692093846346337460743177e+38).

Skaičius 01101011000000000000000 yra mantisa mūsų 32 bitų skaičiaus. Tiesiog pirmas dvejetainis skaitmuo arba bitas prieš kablelį visada yra vienetas, todėl jis nerašomas. Mantisą sudaro 23 bitai. Su pirmu bitu prieš kablelį (kuris visada yra vienetas), gaunasi 24 bitų dvejetainis skaičius, štai toks: 1.01101011000000000000000.

Skaičius 1.01101011000000000000000, kaip parodyta čia [ https://sandbox.mc.edu/~bennet/cs110/flt/ftod.html ], iškoduojamas taip:

1*2^0 + 0*2^(-1) + 1*2^(-2) + 1*2^(-3) + 0*2^(-4) + 1*2^(-5) + 0*2^(-6) + 1*2^(-7) + 1*2^(-8) =

= 1 + 0 + 0.25 + 0.125 + 0 + 0.03125 + 0 + 0.0078125 + 0.00390625 = 1.41796875.

Kadangi negalima traukti kvadratinės šaknies iš neigiamo skaičiaus, tai laikysime, kad mūsų 32 bitų skaičius yra teigiamas. Tada jis tampa toks:

0|01000110|01101011000000000000000.

Taigi, mes turime single precision (32 bitų skaičių):

01101011000000000000000 * 2^(01000110-01111111) = 01101011000000000000000 * 2^(-00111001) =

[Iš tikro reikia dar kažkur saugoti eksponentės minuso/pliuso ženklą, todėl tikriausiai eksponentė gali būti nuo -126 iki 127, kaip ir parašyta Wikipedijoje. Įprastai, kaip supratau, minuso/pliuso ženklas buna pirmas bitas; todėl pirmas eksponentės bitas tikriausiai yra pliuso-minuso ženklas.]

[Bendrai, angliškoje Vikipedijoje pusiau teisingai. Jeigu eksponentė yra 127=254-127, tai tada gaunamas skaičius 01111111 (o 255-127=128=10000000 ir užimamas minuso-pliuso ženklo bitas). Ir 0-127=-127, ir turint galvoje, kad pirmam eksponentės bitui, 1 reiškia -, o 0 reiškia +, gauname, kad eksponentė gali būti nuo -127 iki 127.]

= 1.41796875 * 2^(-57) = 1.41796875 * 6.9388939039072283776476979255676e-18.

Dabar ištraukime šaknį iš 1.41796875 ir 6.9388939039072283776476979255676e-18.

${\sqrt {1.41796875}}=$ 1.1907849302036031393001193597853.

${\sqrt {6.9388939039072283776476979255676\cdot 10^{-18}}}=$ 2.6341780319308769866855605025025e-9.

Tada mūsų 32 bitų skaičius (be minuso ženklo) lygus:

0|01000110|01101011000000000000000 = 1.1907849302036031393001193597853 * 2.6341780319308769866855605025025e-9 =

= 3.1367395038966740342055127049988e-9.

Skaičiaus 0|01000110|01101011000000000000000 eksponentė yra 01000110-127 = -00111001 = 10111001 = -57.

Dvejetainėje sistemoje šaknis iš eksponentės traukiama paslinkus visus bitus į dešinę viena pozicija (išskyrus pirmą bitą, kuris reiškia eksponentės ženklą). Mūsų atveju nereikia stumti pirmo bito į dešinę, nes jis yra 1 ir reiškia minuso ženklą eksponentės. Paslinkus viena pozicija visus eksponentės bitus, išskyrus pirmą, į dešinę, vienas paskutinis bitas iškris ir bus gautas dvejetainis skaičius:

10011100.

Tas iškritęs bitas buvo 1 ir reiškė 2^1. O pastumus tapo ${\sqrt {2}}.$ Todėl galutiniuose skaičiavimose reikės dar padauginti viską iš $2^{-{\frac {1}{2}}}={\frac {1}{\sqrt {2}}}=$ 0.70710678118654752440084436210485.

Skaičius 10011100 = -00011100 lygus tokiam skaičiui:

-(0*2^0 + 0*2^1 + 1*2^2 + 1*2^3 + 1*2^4 + 0*2^5 + 0*2^6) = -28.

$2^{-28}=$ 2^(-28) = 0.0000000037252902984619140625.

Taigi, pagrindinė įrodymo dalis yra šitas rezultatas:

${\sqrt {1.41796875}}\cdot 2^{-28}=$ 1.1907849302036031393001193597853 * 0.0000000037252902984619140625 =

= 4.436019548042130244151202439462e-9.

Galutinis štrichas yra visko padauginimas iš 0.70710678118654752440084436210485. Tada gauname:

0.70710678118654752440084436210485 * 4.436019548042130244151202439462e-9 = 3.1367395038966740342055127049988e-9.

Taigi, gavome tą patį rezultatą [3.1367395038966740342055127049988e-9] kaip ir aukščiau.

Šis įrodymas nėra visiškai pilnas, nes neparodėme kaip traukti kvadratinę šaknį iš mantisos, kuri yra 1.01101011000000000000000. Bet jau žinome kaip skaičiuoti šaknį Niutono metodu. Todėl galima sudaryti lentelę iš 10 bitų (be pirmo bito prieš kablelį, kuris visada yra 1) skaičiams nuo 1.0000000000 iki 1.1111111111.

10 bitų sudaro dešimtainius skaičius nuo 0 iki 1023 (nes 2^10 = 1024).

Tada pagal pirmus 10 bitų mantisos skaičiaus n, gaunama apytikslė ${\frac {1}{\sqrt {n}}}$ reikšmė (be eksponentės) iš ROM, kuri ir yra $y_{0}$ dvejetainėje sistemoje. Žinant $y_{0}$ galima, taikant Newton-Raphson algoritmą, apskaičiuoti tikslią ${\frac {1}{\sqrt {n}}}$ reikšmę dvejetainėje sistemoje.

Kadangi ROM atmintyje gali būti 1024-ios $y_{0}$ reikšmės, o viena tokia reikšmė užima 10 bitų, tai iš viso reikia 1024*10 = 10240 (bitų) = 1280 baitų ROM atminties.

Pataisymas 1. Pasirodo XOR operacija atimčiai dvejetainėje sistemoje tinka tik, kai vieno iš dvejetainių skaičių [su kuriais atliekama atimtis] visi skaitmenys (bitai) yra vienetai.

Papildymas 1. Single precision (32 bitų) skaičiaus eksponentę butų galima iškart saugoti formate nuo -127 iki 127, kai pirmas bitas lygus 1 reiškia minuso ženklą (o 0 reiškia pliuso ženklą). Neaišku kam dar atiminėti iš eksponentės dvejetainio skaičiaus, kuris gali būti nuo 0 iki 255, skaičių 127, kad gauti eksponentę nuo -127 iki 128 (pagal Wikipedia nuo -126 iki 127). Nebent kartais, kai dirbama su skaičiais didesniais už 1, norima išnaudoti didesnę eksponentės reikšmę (255 vietoje 127). Bet toks nukrypimas nuo tvarkos atrodo duoda mažai naudos. Tuo labiau, kad tokiame eksponentės formate, kai ji gali būti nuo 0 iki 255 negalima sudėti neigiamos ir teigiamos eksponentės, nes visada bus gauta tik didesnė eksponentės reikšmė, nors kartais ji turėtų pasidaryt mažesnė arba neigiama.

Papildymas 2. Čia [ https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding ] lentelėje (ir straipsnyje) pasakyta, kad eksponentė gali būti nuo 1 iki 254 (nuo 00000001 iki 11111110). Tada iš tokio skaičiaus atimamas skaičius 127 (kas lygu dvejetainėje sistemoje 01111111).

254-127=127.

1-127=-126.

Skaičiai 0 (dvejetainėje sistemoje 00000000) ir 255 (dvejetainėje sistemoje 11111111) skirti specialiems atvejams.

Jeigu mantisos visi bitai nuliai ir jeigu eksponentės visi bitai nuliai, tada toks skaičius reiškia +- 0.

O jeigu eksponentės visi bitai nuliai, bet mantisos bitai ne visi nuliai, tai toks skaičius reiškia $(-1)^{\text{sign}}\times 2^{-126}\times 0.{\text{fraction}};$ čia fraction reiškia mantisos bitus.

Jeigu eksponentė yra 255, o mantisos visi bitai nuliai, tada toks skaičius reiškia +- begalybę.

Jeigu eksponentė yra 255, o mantisos ne visi bitai nuliai, tada toks skaičius reiškia NaN. Skaičiumi NaN gali būti 0/0, kvadratinės šaknies traukimas iš neigiamo skaičiaus ir dar kai kurie atvejai.

Išeina, jei pirmas eksponentės bitas yra 1, tai eksponentė yra neigiama, o jei pirmas eksponentės bitas 0, tada eksponentė teigiama. Skaičius 126 dvejetainėje sistemoje yra 127-1 = 01111111-00000001 = 01111110 = 126. Tada -126 = 11111110.

Papildymas 3. Angliškos Vikipedijos single precision straipsnyje yra keli vertingi pavyzdžiai.

Pavyzdys 1. Reikia užrašyti dešimtainės sistemos 1 į dvejetainės sistemos single precision.

$(1)_{10}=(1.0)_{2}\times 2^{0}$

Čia eksponentė yra 0, o pridėjus 127 ji bus 127.

$(127+0)_{10}=(127)_{10}=(0111\ 1111)_{2}$

Po vieneto ir kablelio, mantisoje seka vien tik nuliai. Ir kadangi vienetas [prieš kablelį] nesaugomas single precision formate, tai visas vienetas dešimtainėje sistemoje atrodys taip single precision formate:

$(1)_{10}=(0\ 01111111\ 00000000000000000000000)_{2}=({\text{3F800000}})_{16}$

Pavyzdys 3. Reikia paversti 0.375 į dvejetainę sistemą.

$0.375={(0.011)_{2}}={(1.1)_{2}}\times 2^{-2}$

1.1 = 1*2^0 + 1*2^(-1) = 1 + 1/2 = 1.5.

2^(-2) = 1/4 = 0.25.

1.5 * 0.25 = 0.375.

-2 eksponentė paverčiama į rodos suprantama procesoriui kalbą, t. y. į 127+(-2)=125:

$(127+(-2))_{10}=(125)_{10}=(0111\ 1101)_{2}$

Vadinasi dvejetainėje single precision sistemoje skaičius 0.375, atrodys taip:

$(0.375)_{10}=(0\ 01111101\ 10000000000000000000000)_{2}=({\text{3EC00000}})_{16}$

Pastaba. Šešioliktainė sistema yra hexodecimal'ai. 4 bitai sudaro vieną hexodecimal'ą. Todėl iš viso yra 16 hexodecimal reikšmių:

0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E, F.

Atimtis dvejetainėje sistemoje

https://www.cuemath.com/numbers/binary-subtraction/

https://www.geeksforgeeks.org/binary-subtraction/

Paprasčiausias ir suprantamiausias būdas daryti atimties operaciją yra panašus į sudeties operaciją, tik vietoje carry vyksta borrow (vietoje vieno bito atmintyje, kuri reikia pridėti, yra vienas bitas atmintyje kuri reikia atimti)...

Atimkime iš 8 bitų skaičiaus 01000000, kitą 8 bitų skaičių 00000001, kuris yra tiesiog vienetas.

01000000-

00000001=

00111111.

Atimkime iš 8 bitų skaičiaus 01001000, kitą 8 bitų skaičių 00000001, kuris yra vienetas.

01001000-

00000001=

01000111.

Atimkime iš 8 bitų skaičiaus 01000000, kitą 8 bitų skaičių 00001001.

01000000-

00001001= 00111111-00001000 =

00110111.

Mandresnis būdas atlikti atimtį vadinasi Binary Subtraction Using 1's Complement. Jo esmė yra padaryti XOR operacija tarp skaičiaus, kuris bus atiminėjamas ir skaičiaus sudaryto iš vienetų dvejetainėje sistemoje. Tas rezultatas išėjęs iš šitos XOR operacijos turi būti pridėtas prie skaičiaus iš kurio daroma atimtis. Dar prie visko pridedamas carry bitas (devintas bitas iš kairės, jeigu jis tapo vienetu). Kitaip tariant dar pridedamas skaičius 00000001, kuris reiškia vienetą dvejetainėje sistemoje. Gautas rezultatas ir bus atlikta atimtis tarp dviejų skaičių. Tas rezultatas bus teisingas ir iškart normaliai iškoduojamas, jeigu reiškinyje a - b, teisinga nelygybė $a\geq b.$

Pirmas pavyzdys. Iš 8 bitų skaičiaus 00110110 (kuris yra 2+4+16+32=54) atimsime 8 bitų skaičių 00100010 (kuris yra 2+32 = 34).

Padarome XOR operacija tarp skaičiaus 00100010 ir 8 bitų skaičiaus 11111111.

00100010 XOR

11111111 =

11011101.

Tada toliau sudedame 00110110 (54) su 11011101 (1+4+8+16+64+128 = 221):

_00110110 +

_11011101 =

100010011.

Tada toliau pridedame devintą kairį bitą prie gauto rezultato be devinto bito:

00010011 +

00000001 =

00010100 (4+16 = 20).

Taigi, gavome 8 bitų rezultatą 00010100, kuris reiškia 20 dešimtainėje sistemoje. Ir iš tiesų 54-34=20.

Dar vienas mandras būdas vadinasi Binary Subtraction Using 2’s Complement. Šis budas praktiškai yra tas pats kas Binary Subtraction Using 1's Complement... Skirtumas tas, kad būde Binary Subtraction Using 2’s Complement atlikus XOR operaciją analogiškai kaip Binary Subtraction Using 1's Complement būde, yra pridedamas dar vienetas (jei skaičiuojama su 8 bitais, tai pridedamas 8 bitų skaičius 00000001) nepriklausomai nuo to ar carry [devintas kairis] bitas yra 1 ar 0. Pačio carry bito pridėti nereikia.

Antras pavyzdys. Atimsime iš skaičiaus 01001101 (1+4+8+64 = 77), skaičių 01101010 (2+8+32+64 = 106).

Padarome XOR operacija tarp skaičiaus 01101010 ir 8 bitų skaičiaus 11111111. Tada gauname:

01101010 XOR

11111111 =

10010101 (1+4+16+128 = 149).

Tada pridedame 1 prie 149, t. y. pridedame 00000001 prie 10010101. Gauname:

10010101 +

00000001 =

10010110 (2+4+16+128 = 150).

Tada toliau pridedame 150 prie 77, t. y. pridedame 10010110 prie 01001101. Gauname:

_01001101 +

_10010110 =

011100011 (1+2+32+64+128 = 227).

Matome, kad šiuo atveju carry bitas yra 0, o rezultatas yra 11100011, kas lygu 227 dešimtainėje sistemoje. Na o 77 - 106 = -29. Iš tikro tai rezultatas visada atrodys kaip nesamonė, jeigu iš mažesnio skaičiaus atiminėjamas didesnis.

Bet šitą rezultatą galima atversti į -29 tik be minuso, t. y. galima atversti į 29. Reikia padaryti XOR operaciją tarp mūsų skaičiaus 11100011 (227) ir 8 bitų skaičiaus 11111111. Tada

11100011 XOR

11111111 =

00011100 (4+8+16 =28).

Prie šito rezultato visada dar reikia pridėti vienetą, t. y. 00000001. Tada

00011100 +

00000001 =

00011101 (1+4+8+16 =29).

Taigi gavome 29 be minuso (primename, kad 77 - 106 = -29).

Teoriškai pridėję skaičių 227 (11100011) prie kokio nors 8 bitų skaičiaus ir dar prijėję 1 (00000001), turėtume gauti atimtį iš to 8 bitų skaičiaus. Pavyzdžiui, skaičius 00110110 reiškia 54 (2+4+16+32=54). Tai prie 54 pridėję 227 ir 1, turėtume gauti 54-29=25. Sudedame 54 su 227:

_00110110 +

_11100011 =

100011001.

Pirmą bitą iš kairės, žinoma, ignoruojame. Tada turime atsakymą 00011001, kuris reiškia 1+8+16 = 25. Matome, kad gavome teisingą atsakymą ir net, vadinasi, nereikia dar pridėti vieneto (00000001).

Trečias pavyzdys. Iš Z80 CPU User Manual paimsime pavyzdį NEG operacijos (ar instrukcijos), kuris yra 177 puslapyje ant lapo arba 191 puslapyje pagal Acrobat Reader.

Tas pavyzdys yra NEG instrukcija, kuri 8 bitų skaičių paverčia į neigiamą 8 bitų skaičių. Tas pats kas iš nulio atimti skaičių ir įdėti jį į akumuliatorių.

Apibūdinimas tame PDF yra toks:

The contents of the Accumulator are negated (two’s complement). This method is the same as subtracting the contents of the Accumulator from zero.

Taigi, Zilog Z80 procesorius padaro neigiamą skaičių naudodamas XOR operaciją. Ten yra pavyzdys, kad NEG instrukciją akumuliatoriaus reikšmę paverčia iš 10011000 (8+16+128 = 152) į 01101000 (8+32+64 = 104). Parodysime kaip tai vyksta.

_01101000 XOR

_11111111 =

101100111.

Pirmą kairį bitą atmetame. Atmetę gauname 01100111 (1+2+4+32+64 =103). Toliau dar reikia pridėti vienetą (00000001).

01100111 +

00000001 =

01101000 (8+32+64 = 104).

Taigi gavome žadėtąjį 01101000 (104) po NEG operacijos. Jį galima iškoduoti į -152 analogiškai kaip antrame pavyzdyje. Tai ir padarysime.

01101000 XOR

11111111 =

10010111 (1+2+4+16+128 = 151).

Prie gauto skaičiaus po XOR operacijos, pridedamas vienetas. Tada gauname:

10010111 +

00000001 =

10011000 (8+16+128 = 152).

Gavome 152 be minuso. Atrodo, kam reikalina ta NEG instrukcija, jeigu vis tiek duoda neigiamą skaičių be minuso ir neturintį jokios prasmės? Tas gautas skaičius 01101000 (104) po NEG instrukcijos užkoduotas slypi kaip -152. Todėl jį sudėję, pavyzdžiui, su 183, gausime 183+(-152)=183-152=31. Skaičius 183 = 128+32+16+4+2+1. Dvejetainėje sistemoje 183 yra 10110111.

Taigi, sudėkime 10110111 (183) ir 01101000 (104).

_10110111 +

_01101000 =

100011111.

Pirmas kairis bitas, kuris yra carry, atmetamas. Tada gaunamas skaičius 00011111 (1+2+4+8+16 = 31). Taigi gavome atsakymą 31, kas lygu 183-152.

Toliau bus pateikta truputi istorijos apie Zilog Z80 ir x86 architektūrą.

Pirmą ar beveik pirmą mikroprocesorių intel 4004 sukūrė italas Federico Faggin. Paskui jis sukūrė procesorių 4040. Šitie du procesoriai yra nesusiję su 8008 ar 8080 intel procesoriais ir jų architektūra mirė. Pirmas procesorius, nuo kurio prasidėjo x86 architektūra, buvo 8008. Procesoriaus 8008 architektūrai pagrindą padėjo Datapoint Corporation, originally known as Computer Terminal Corporation (CTC) su savo tam tikrom architektūrom. Paskui Federico Faggin remdamasis tomis architektūromis vienas ar su pagalba sukūrė intel 8080 procesorių (gal dar Federico Faggin prisidėjo kūriant intel 8008 procesorių). Po 8080 procesoriaus sukūrimo Federico Faggin išėjo iš intel ir sukūrė Zilog kompaniją, kurioje sukūrė patobulintą 8080 procesoriaus versiją - procesorių Z80. Tuom jo šlovė ir pasibaigė. Intel patobulino 8080 procesorių, taip sukūrdama intel 8085 procesorių. Dar stipriai patobulinus 8085, gavosi 8086 procesorius, po kurio ir prasidėjo x86 era. Todėl intel 8080 ir Zilog Z80 procesoriai turi daug bendro su x86 architektūra.

Italiano patobulino americano, o paskui americano patobulino italiano, sukurdami 8085 procesorių ir vėliau 8086.

Atrodo dar pačioje pradžioje, kuriant 8008 ar dar prieš jo kūrimą, prisidėjo truputi Texas Instruments, bet kažkas jiems neišėjo.

Neigiamų skaičių sudėtis (2’s Complement)

Jeigu yra 2 neigiami skaičiai 2’s Complement formate, tai juos sudėjus bus gautas taip pat neigiamas skaičius 2’s Complement formate.

Ketvirtas pavyzdys. Turime skaičių 00101101 (1+4+8+32 = 45). Kitas skaičius yra 11001010 (2+8+64+128 = 202). Pirma padarysime juos neigiamais 2’s Complement formate.

00101101 XOR

11111111 =

11010010 (2+16+64+128 = 210).

11010010 + 00000001 = 11010011 (211). Pirmas skaičius 2’s Complement formate yra 211, kas užkoduotai reiškia -45.

11001010 XOR

11111111 =

00110101 (1+4+16+32 = 53).

00110101 + 00000001 = 00110110 (2+4+16+32 = 54). Antras skaičius 2’s Complement formate yra 54, kas užkoduotai reiškia -202.

-45+(-202) = -247. Sudedame:

_11010011 +

_00110110 =

100001001.

Atmetame pirmą bitą (9-tą bitą nuo galo). Tada gauname 00001001 (1+8=9). Toliau atversime šitą skaičių į -247 be minuso, t. y. į 247.

00001001 XOR

11111111 =

11110110.

11110110 + 00000001 = 11110111 (1+2+4+16+32+64+128 = 247).

Pridėkime užkoduotą 2’s Complement formate skaičių -247 (00001001) prie skaičiaus 252 (4+8+16+32+64+128 = 252), Dvejetainėje sistemoje 252 yra 11111100. Sudėjus turėtume gauti 252 + (-247) = 252-247 = 5.

_11111100 +

_00001001 =

100000101.

Atmetame carry bitą. Tada turime skaičių 00000101 (1+4 =5). Taigi, parodėme, kad galima sudėti du neigiamus skaičius 2’s Complement formate ir paskui juos pridėjus prie teigiamo skaičiaus gausime teisingą atsakymą, tarsi būtume atėmę tuos du skaičius (45 ir 202) iš trečio skaičiaus (252).

Kaip galėtų veikti eksponentė single precision formate? Exponentė single precision formate saugoma su reikšmėm nuo 0 iki 254. Bet taip su ja nieko negalima skaičiuoti. Kad suprasti ką ji reiškia, kaip sako Vikipedija ir kitas šaltinis, reikia iš eksponentės 8 bitų reikšmės (kuri buna nuo 0 iki 254) atimti skaičių 127. Iš 128 atėmę 127, gausime 1, t. y. gausime eksponentės reikšmę $2^{1}.$ Jeigu eksponentė yra 130, tai atėmę 127, gausime eksponentės reikšmę $2^{3}\;$ (130 - 127 = 3). Vadinasi, kai pačioje pradžioje duotas skaičius single precision formate ir jeigu jo eksponentė yra nuo 128 iki 254 imtinai, tai pirmas kairis bitas visada lygus vienetui, nes 128 yra 10000000 dvejetainėje sistemoje. O jeigu eksponentė pačioje pradžioje (prieš atimant 127) yra mažiau už 128, tai pirmas kairis bitas visada yra nulis. Todėl procesorius gali suprasti, kada bus operuojama su neigiamais skaičiais ir kada bus operuojama su teigiamais skaičiais.

Jeigu iš eksponentės, kuri yra mažiau už 127 (jeigu eksponentė lygi 127, tai atėmus 127, tai reikš $2^{0}=1$), atimti 127, tai gausime neigiamą skaičių. Skaičius -127 jau yra procesoriaus atmintyje 2’s Complement formate. Tokiame formate jis yra lygus 10000001, nes

01111111 XOR 11111111 = 10000000,

10000000 + 00000001 = 10000001.

Taigi, 2’s Complement formate -127 yra 10000001. Tada tereikia pridėti 10000001 (-127) prie eksponentės (kuri dar yra nuo 1 iki 254; 0 skirtas specialiems atvejams), kad ir kokia ji būtų. Po šitos sudeties bus gautos galimos reikšmės nuo 127 iki -126.

Jeigu pačioje pradžioje (prieš pridedant -127) eksponentė buvo 127, tai 127+(-127)=0 arba

_01111111 +

_10000001 =

100000000.

Atmetus pirmą kairį carry bitą, gausime apvalų 0 (00000000 dvejetainėje sistemoje).

Vadinasi, kai sudauginami du single precision skaičiai, tai jų eksponentės sudedamos ir visada gaunamas arba teigiamas eksponentės skaičius arba užkoduotas (2’s Complement formate) neigiamas eksponentės skaičius. Tai daugybai nieko daugiau daryt ir nereikia. Nes sudaugintų dviejų single precision skaičių, kai prie jų eksponenčių jau pačioje pradžioje pridėtas skaičius -127 (2’s Complement formate), eksponentės bus sudedamos visada teisingai.

Pasinaudoti pirmu kairiu bitu eksponentės (kuri yra 8 bitų single precision formate), prireiks jeigu reikia sudėti du single precision skaičius su skirtingom eksponentėm. Sakėm, kad jeigu [prieš pridedant -127 (2’s Complement formate) prie eksponentės] pirmas kairis bitas yra 1, tai toks eksponentės skaičius yra teigiamas, o jeigu pirmas kairis bitas yra 0, tai toks skaičius yra neigiamas.

Jeigu abu single precision skaičiai, kurie bus sudedami, turi skirtingas, bet teigiamas eksponentes, tai to skaičiaus, kuris turi mažesne eksponentę, mantisos bitai bus paslenkami kažkiek pozicijų į dešinę. Pavyzdžiui, jeigu abu skaičiai turi teigiamas eksponentes ir pirmo skaičiaus eksponentė (po pridėjimo -127 2’s Complement formate) yra 8, o antro skaičiaus eksponentė (po pridėjimo -127 2’s Complement formate) yra 3, tai antro skaičiaus mantisos bitai turi būti pastumi 8-3=5 bitų pozicijom į dešinę. Po šito veiksmo galima sudėti tuos du skaičius ir palikti eksponentę pirmojo skaičiaus. Antrojo skaičiaus išstumti 5 dešiniausi mantisos bitai dingsta.

Neadekvačiai ilgas procesoriaus ar FPU skaičiavimas prasideda, kai reikia sudėti du skaičius su neigiamom eksponentėm, arba kai reikia pridėti skaičių su teigiama eksponente prie skaičiaus su neigiama eksponente. Tokiu atveju skaičiaus su neigiama eksponente, paimamas apdorojimui. Jo eksponentė po XOR su 11111111 ir po pridėjimo 00000001, tampa teigiama (pvz., jeigu buvo 2’s Complement formate -13, tai tapo 13). Tada, kai skaičiaus su neigiama eksponente, eksponentė tapo jau teigiama, ji pridedama prie pirmo skaičiaus eksponentės, kuris yra su teigiama eksponente. Gautas rezultatas reiškia kiek mantisos bito pozicijų reikia pastumti mantisos bitus antro skaičiaus (kuris pradžioje turėjo neigiamą eksponentę).

Pavyzdžiui, jeigu Pirmas skaičius turi teigiamą eksponentę 16, o antras skaičius turi neigiamą eksponentę -4, tai po pavertimo (taikant XOR ir +1) antro skaičiaus eksponentės į +4, šitų dviejų skaičių eksponentės sudedamos (16+4=20) ir antro skaičiaus mantisos visi bitai pastumiami 20 pozicijų į dešinę (20 dešiniausių antro skaičiaus mantisos bitų dingsta, o tuščius bitus mantisos kairėje pakeičia nuliai). Po tokios operacijos galima sudėti tuos du skaičius, palikus eksponentės reikšmę pirmo skaičiaus (kuris turėjo ir turi teigiamą eksponentę lygią 16).

Matome, kad skaičiuojant su skaičiais, po pridėjimo prie eksponenčių -127 2’s Complement formate, tampa neaišku, kurio skaičiaus ekponentė neigiama, o kurio - teigiama. Vadinasi, reikia pridėti 127 prie tokio neaiškaus skaičiaus eksponentės ir tada, jeigu pirmas kairis bitas yra 1, tai eksponentė yra teigiama, o jeigu pirmas kairis bitas yra 0, tai eksponentė šito skaičiaus yra neigiama.

Jeigu abu sudedami skaičiai turi neigiamas eksponentes, tai reikia jas paversti į teigiamus skaičius be minuso ženklo. Pvz., jeigu buvo -9 eksponentė (2’s Complement formate), tai po XOR su 11111111 ir po +1, ši eksponentė taps +9. Bet geriau paprasčiausiai pridėti prie neigiamų eksponenčių (kurios yra 2’s Complement formate) 127. Tada jos taps teigiamos. Žinant, kad, tarkim, pirmo skaičiaus eksponentė po pridėjimo 127 tapo 35, o antro skaičiaus eksponentė po pridėjimo 127 tapo 42, atimame iš didesnio skaičiaus mažesni. Šiuo atveju gauname, 42-35=7. Vadinasi, pirmo skaičiaus su mažesne eksponente visus mantisos bitus reikia pastumti 7 pozicijom į dešinę. Po šitos operacijos tuos du skaičius galima sudėti, išsaugojant antro skaičiaus eksponentę (42). Be to, pridėję 127 prie neigiamų eksponenčių (2’s Complement formate), pagal pirmą kairį bitą (0) žinosime, kad jos neigiamos.

Vadinasi, sudedant du skaičius, reikia dar lyginti, kurio skaičiaus eksponė yra didesnė. Pas modernius procesorius, tas lyginimas gal vyksta momentaliai, lygiagrečiai su kokia nors kita operacija... Skaičius su mažesne eksponente pastumiamas tiek mantisos bito pozicijų į dešinę koks yra [atimties] skirtumas tų eksponenčių. Čia net nereikia nieko daugiau ir daryt.

1. Reikia pridėti 127 prie pirmo ir antro skaičiaus eksponentės.

2. Reikia palyginti, kurio skaičiaus eksponentė yra didensnė (kas gali būti nemažai tranzistorių reikalaujanti schema).

3. Reikia atimti mažesnę eksponetę iš didesnės eksponentės.

4. Reikia pastumti skaičiaus mantisos bitus su mažesne eksponente tiek pozicijų į dešinę, koks yra tų dviejų skaičių eksponenčių skirtumas.

5. Po šitų 4 žingsnių galima atlikti tų dviejų single precision skaičių sudėtį, išsaugant eksponentę to skaičiaus (iš tų dviejų), kuris turėjo ją didesnę.

Pastaba. Single precision formate eksponentės reikšmė 0 (arba 00000000 dvejetainėje sistemoje) skirta specialiems atvejams. Kadangi vykdant skaičiavimus single precision formate, eksponentė yra nuo -126 iki 127 imtinai, tai eksponentė negali būti lygi -127. Kitaip tariant, negalima single precision skaičiams suteikti eksponentės reikšmės -127 ir mažesnės, nes -127+127=0 (skaičius 127 pridedamas, kai, pavyzdžiui, norima išeiti iš skaičiavimų ir papasuoti šitą single precision skaičių kur nors saugojimui ateičiai). O nulis (00000000), kaip sakyta, skirtas specialiems atvejams.

Kita pastaba. Pastumus mantisos bitus kokio nors single precision skaičiaus, reikia dar prirašyti 1 prie kairiausio bito iš buvusios mantisos. Nes single precision formate pirmas bitas prieš kablelį (tašką) visada yra vienetas. Todėl jis nesaugomas single precision 32 bitų gabale. O pastumiant visus mantisos bitus kažkiek pozicijų į dešinę, tas vienetas prieš kablelį neegzistuoja mantisoje, todėl jį kompiuteris turi pats prirašyti mantisoje iš kairės prie pat visų pastumtų mantisos bitų.

Trečia pastaba. Jeigu, pavyzdžiui, FPU daro daug iteracijų ir naudoja ten sudėti su vienu ir tuo pačiu skaičiumi, kai kito skaičiaus eksponentė yra didesnė ir nesikeičia, tai modernus procesoriai galėtų atsiminti, kad to pirmo skaičiaus mantisos bitai visada turi būti pastumti į dešinę tiek pat pozicijų. Tai turėtų pagreitinti sudėtį kai kuriais atvejais. Šiuo atveju reikia lygiagrečiai tikrinti ar antro skaičiaus mantisa nepasikeitė... Kad su neprognozuojamų skaičių sudėtim nebūtų dar ilgesni skaičiavimai, dėl to tai reikia daryti tik lygiagrečiai (du darbus vienu metu).

Ketvirta pastaba. Panašu, kad kai pirdedamas neigiamas skaičius (2’s Complement formate), kuris modulyje yra didesnis už teigiamą skaičių prie kurio pridedamas, tai po sudeties neatsiranda carry vienetas devintame bite skaičiuojant nuo galo, t. y. carry bitas visada buna 0. O jeigu pirdedamas neigiamas skaičius (2’s Complement formate), kuris modulyje yra mažesnis už teigiamą skaičių prie kurio pridedamas, tai visada devintas carry bitas nuo galo (kai sudedami 8 bitų skaičiai) yra 1.

Todėl pagal [devintą] carry bitą galima nustatyti ar po sudeties skaičius yra neigiamas (2’s Complement formate) ar teigiamas. Jeigu carry bitas yra 0, tai gautas skaičius yra neigiamas (2’s Complement formate), o jeigu carry bitas yra 1, tai gautas skaičius (po teigiamo ir neigiamo skaičiaus sudeties) yra teigiamas.

Teoriškai taip galima lyginti kuris sveikasis skaičius yra didesnis. Tik reikia žinoti kokie tai sveikieji skaičiai, ar neigiami ar teigiami... Du neigiamus skaičius (2’s Complement formate) lyginti yra ilgas darbas šituo metodu, nes vieną iš jų reikia paversti į teigiamą skaičių. O jeigu duoti du teigiami skaičiai, tai taipogi vieną iš jų reikia paversti į neigiamą skaičių (2’s Complement formate) ir tada susidarius planą, kas ką reiškia, lyginti juos. Pavyzdžiui, jeigu lyginami du teigiami skaičiai, tai, jeigu pirmas skaičius yra didesnis už antrą, tai antrą skaičių pavertus (po XOR su 11111111 ir po pridėjimo 00000001) į neigiamą skaičių (2’s Complement formate) ir sudėjus su pirmu skaičium, reikia pasižiūrėti į carry bito reikšmę. Šiuo atveju carry bitas bus 1. Tai reikš, kad pirmas skaičius yra didesnis už antrą (kai jie abu buvo teigiami).

Bet greičiausiai yra tranzistorių schemos kurios palygina du teigiamus skaičius momentaliai. Kad lyginti dar ir neigiamus skaičius, t. y., kad momentaliai arba bent greitai lyginti visus sveikuosius skaičius, tranzistorių schemos gali būti gerokai didesnės. Zilog Z80 procesorius net neturi sveikųjų ar teigiamų skaičių palyginimo instrukcijos, kuri įjungtų kokią nors vėliavėlę (1 reiškia įjungta, 0 reiškia išjungta), jei pvz., pirmas skaičius didesnis už antrą. Būna dar instrukcijos, kurios pagal palyginimą (kai atėmus antrą skaičių iš pirmo gaunamas 0) ar kokią nors veliavėlės reikšmę atlieką šuolį iš vieno atminties adreso (Program Counter) į kitą adresą. Z80 procesoriuje tokios pagrindinės veliavėlės yra Zero Flag (pvz., kai akumuliatoriaus reikšmė yra 0 ar atlikus dvieju skaičių atimtį gaunamas 0), Carry Flag (kai carry bitas yra 1, tai Carry Flag is set) ir dar pora kitų veliavėlių pagal kurias galima pakeisti procesoriaus eigą.

Naujesniuose procesoriuose, užmetus akį, tų skaičių lyginimo instrukcijų pilna (ir netgi su ženklais $\geq ,$ ne tik > ar <).

Penkta pastaba. Jeigu sudedami du sveikieji skaičiai, kai pirmas skaičius yra teigiamas, o antro neigiamo (2’s Complement formate) skaičiaus modulis lygus pirmam skaičiui, tai tada juos sudėjus bus gautas 0 ir carry bitas (devintas bitas skaičiuojant nuo galo) bus 1. Todėl iš esmės tokia technika tinka lyginti tik natūraliuosius skaičius, naudojant ženklą $\geq$ (daugiau arba lygu). Skaičius 00000000 po XOR su 11111111 ir po +00000001 bus 00000000 (ir dar priedo įjungs Carry Flag). O bet koks kitas skaičius, pavyzdžiui, 00000001 po XOR 11111111 bus 11111110 ir pridėjus 00000001, bus 11111111 (Carry flag'o neįjungs). Bet sudėjus 00000001 su 00000000 (antras skaičius praėjo vertimą į neigiamą skaičių, bet vis tiek liko tas pats 0) bus gautas 00000001 ir carry bitas bus 0. Ir jeigu sudėsime 00000000 su 11111111 (antras skaičius buvo 00000001 ir praėjo vertimą į neigiamą skaičių, po kurio tapo 11111111), gausime rezultatą 11111111, o carry bitas bus irgi 0. Vadinasi skaičių 00000000 negalima lyginti su kitais teigiamais skaičiais. Todėl šituo budu galima lyginti tik natūraliuosius skaičius, naudojant ženklą $\geq .$

(Palyginus teigiamą skaičių su neigiamu, carry bitas yra 1 (po sudeties teigiamo ir neigiamo skaičiaus 2’s Complement formate), kai teigiamas skaičius yra didesnis už neigiamo skaičiaus modulį. Ir carry bitas yra 0 (po sudeties teigiamo ir neigiamo skaičiaus 2’s Complement formate), kai teigiamas skaičius yra mažesnis už neigiamo skaičiaus modulį).

Minuso ženklas neigiamuose 8 bitų skaičiuose.

Čia [ https://www.zilog.com/docs/z80/um0080.pdf ] 69 puslapyje ant lapo arba 83 puslapyje pagal Acrobat Reader aprašomas Sign Flag:

The Sign Flag (S) stores the state of the most-significant bit of the Accumulator (bit 7). When the Z80 CPU performs arithmetic operations on signed numbers, the binary twoscomplement notation is used to represent and process numeric information. A positive number is identified by a 0 in Bit 7. A negative number is identified by a 1. The binary equivalent of the magnitude of a positive number is stored in bits 0 to 6 for a total range of from 0 to 127. A negative number is represented by the twos complement of the equivalent positive number. The total range for negative numbers is from –1 to –128.

When inputting a byte from an I/O device to a register using an IN r, (C) instruction, the S Flag indicates either positive (S = 0) or negative (S = 1) data.

Akumuliatoriuje ir visur kitur 8 bitų skaičiaus visi bitai žymimi nuo 0 iki 7. Bitas 0 yra pats dešiniausias bitas, o bitas 7 yra pats kairiausias bitas. Minuso ženklas saugomas bite 7 (kairiausiame), kai dirbama su neigiamais skaičiais 2’s Complement formate. Tada teigiami skaičiai gali būti nuo 0 iki 127 imtinai (nors nulis labiau kaip neutralus...), o neigiami skaičiai 2’s Complement formate gali būti nuo -1 iki -128. Jeigu teigiamas skaičius būtų 128, tai jis atrodytų taip 10000000 ir bitas 7 (pats kairiausias) būtų vienetas, kas reikštų neigiamą skaičių (127 yra 01111111, 126 yra 01111110, ...).

Pirmas pavyzdys. Paversime 8 bitų skaičių 128 į neigiamą (2’s Complement formate) skaičių -128 ir pamatysime, kad bitas 7 (pats kairiausias) yra vienetas. Tada pridėsime šį skaičių prie 127, ir gavę -1 (2’s Complement formate) pridėsime prie skaičiaus 4 (00000100). Tada turėtume gauti skaičių 3 (00000011).

10000000 XOR

11111111 =

01111111 (127).

01111111 + 00000001 = 10000000 (vėl 128, tik šikart two's Complement fomate, todėl jis reiškia -128).

Toliau pridėsime -128 (two's Complement fomate) prie 127 ir gausime -1 (two's Complement fomate):

01111111 +

10000000 = 11111111 (255 arba -1 two's Complement fomate).

Toliau pridėsime -1 two's Complement fomate prie 4 ir gausime 3. Taigi,

_00000100 +

_11111111 =

100000011.

Gavome carry bitą vienetą, kas ekvivalentu tam, kad Carry Flag Z80 procesoriuje būtų Set (būtų 1). Be kairiausio carry bito atsakymas yra 00000011 (1+2 = 3). Gavome atsakymą 3, kaip ir sakėme. Be to, 7-tas bitas yra 0, kas reiškia teigiamą skaičių (00000011).

Z80 procesoriuje yra instrukcijų, kurios gali reguoti ar keisti proceso eigą priklausomai nuo to ar Sign Flag (S) yra Set ar Reset (1 ar 0). Tas Sign Flag yra bitas 7 (pats kairiausias iš 8 bitų) procesoriaus akumuliatoriuje, kuriame ir su kuriuo vyksta visokie skaičiavimai (kaip dviejų 8 bitų skaičių sudėtis, kai pirmas skaičius yra akumuliatoriuje, o antras skaičius yra kokiam nors kitame registre ar RAM; gautas rezultatas įrašomas į akumuliatorių).

Nuorodos

• https://en.wikipedia.org/wiki/Methods_of_computing_square_roots

• http://www.techmat.vgtu.lt/konspektai/Dskaiciavimas/t_03.pdf Sekos riba.
• http://www.techmat.vgtu.lt/konspektai/Dskaiciavimas/t_04.pdf Konverguojančios sekos ir jų pagrindinės savybės.