Sekos riba

Iš testwiki.
19:38, 9 spalio 2024 versija, sukurta imported>Paraboloid (→Atimtis dvejetainėje sistemoje)
(skirt) ← Ankstesnė versija | Dabartinė versija (skirt) | Vėlesnė versija → (skirt)
Pereiti į navigaciją Jump to search


3. Monotoninės sekos

1. Monotoninių sekų apibrėžimas.

Apibrėžimas. Seka {xn} 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ė
xn≤xn+1 (xn≥xn+1).
Nemažėjančios ir nedidėjančios sekos vadinamos monotoninėmis sekomis. Jei monotoninės sekos {xn} elementai su bet kuriuo numeriu n tenkina nelygybę xn<xn+1 (xn>xn+1), tai seka {xn} 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,12,12,...,1n,1n,... 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 12,23,34,...,nn+1,... didėjanti. Ji aprėžta iš abiejų pusių: iš apačios pirmuoju elementu 12, 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 {xn} aprėžta iš viršaus (iš apačios), tai ji konverguoja.
Pagal praeitą skirsnį seka {xn}, atitinkanti 3.15 teoremos sąlygas, yra aprėžta. Todėl 3.15 teoremą galima trumpai formuluoti taip: jei monotoninė seka {xn} yra aprėžta iš abiejų pusių, tai ji konverguoja.
Įrodymas. Kadangi seka {xn} yra aprėžta, tai jos elementų aibė turi tikslųjį viršutinį rėžį x‾ ir tikslųjį apatinį rėžį x_. Įrodysime štai ką: jei {xn} – nemažėjanti seka, tai jos riba yra minėtasis tikslusis viršutinis rėžis; jei {xn} – nedidėjanti seka, tai jos riba yra minėtasis tikslusis apatinis rėžis x_. Išnagrinėsime tik nemažėjančią seką, nes samprotavimai apie nedidėjančią seką būtų analogiški.
Kadangi x‾ yra sekos {xn} elementų aibės tikslusis viršutinis rėžis, tai kiekvieną ε>0 atitinka toks elementas xN, kad xN>x‾−ε ir xN≤x‾ (bet koks elementas xn ne didesnis už tikslųjį viršutinį rėžį x‾, t. y. xn≤x‾). Iš šitų nelygybių gauname nelygybes 0≤x‾−xN<ε. Kadangi {xn} – nemažėjanti seka, tai nelygybės xN≤xn≤x‾ yra teisingos, kai n≥N. Iš to išplaukia, kad 0≤x‾−xn≤x‾−xN, kai n≥N. Anksčiau pabrėžėme, kad x‾−xN<ε, todėl, kai n≥N, teisingos nelygybės 0≤x‾−xn<ε, iš kurių gauname nelygybę |xn−x‾|<ε. Vadinasi, įsitikinome, kad x‾ yra sekos {xn} 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 {xn}, apibrėžta formule xn=(−1)nn, 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 {xn} yra nemažėjanti ir aprėžta seka, o x‾ – tos sekos riba, tai visi jos elementai tenkina nelygybę xn≤x‾. Nedidėjančios ir aprėžtos sekos {xn}, konverguojančios į x_, elementai tenkina nelygybę x_≤xn. Kad šie teiginiai teisingi, įsitikinome, įrodinėdami 3.15 teoremą.
3.15 teoremos išvada. Sakykime, duota begalinė segmentų sistema [a1;b1],[a2;b2],...,[an;bn],... Jei kiekvienas segmentas yra pirmesniajame*, o skirtumas bn−an (jį vadinsime segmento [an;bn] ilgiu) artėja prie nulio, kai n→∞ (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 [an;bn]. Bet tokiu atveju su bet kokiu numeriu n būtų teisingos nelygybės bn−an≥d−c>0; tai neįmanoma, nes bn−an→0, kai n→∞. Dabar įsitikinsime, kad taškas c, priklausantis visiems segmentams [an;bn], egzistuoja. Kadangi segmentų sistema yra susitraukianti, tai kairiųjų galų seka {an} yra nemažėjanti, o dešiniųjų galų seka {bn}, – nedidėjanti. Kadangi abi tos sekos yra aprėžtos (visi sekų {an} ir {bn} elementai priklauso segmentui [a1;b1]), tai pagal 3.15 teoremą jos konverguoja. Iš to, kad skirtumas bn−an 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: an≤c≤bn, t. y. taškas c priklauso visiems segmentams [an;bn].

_______________

* Tai reiškia, kad an−1≤an<bn≤bn−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ą {xn}, kurios elementas xn lygus
xn=a+a+a+...+a,a>0
(n radikalų). Tą pačią seką galima, savaime aišku, nusakyti šitokia rekurentine formule:
x1=a,xn=a+xn−1.
Norėdami įsitikinti, kad sekos {xn} riba egzistuoja, įrodysime, jog ta seka yra didėjanti ir aprėžta. Kad ta seka didėja, matyti tiesiog. Įsitikinsime, kad seka {xn} aprėžta iš viršaus skaičiumi A, kai A – didedsnysis iš dviejų skaičių a ir 2. Jei xn≤a, tai teiginys jau įrodytas. Jei xn>a, tai, dešinėje nelygybės xn2=a+xn−1≤a+xn pusėje skaičių a pakeitę didesniu už jį skaičiumi xn, gauname xn2<2xn, o iš čia xn<2. Taigi įrodėme, kad seka {xn} 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ę
xn2=a+xn−1,
iš kurios matyti, kad sekos {xn2} ir {a+xn−1} sutampa. Todėl jos turi bendrą ribą. Kadangi pirmosios sekos riba lygi c2, o antrosios a+c (kai n varo į begalybę skaičiams xn ir xn−1, tai xn praktiškai nesiskiria nuo xn−1 ir abu yra labai arti c reikšmės (abu labai panašūs į c)), tai c2=a+c. Iš čia, turėdami mintyje, kad c>0, randame c:
c2=a+c,
c2−c−a=0,
c1,2=−(−1)±(−1)2−4⋅(−a)2=1±1+4a2;
kadangi c>0, tai
c=1+1+4a2.
2 pavyzdys. Tirsime seką {xn}, kuria dažniausiai naudojamasi, apskaičiuojant teigiamo skaičiaus a kvadratinę šankį elektronine skaičiavimo mašina. Šita seka apibrėžiama tokia rekurentine formule:
xn+1=12(xn+axn),n=1,2,3,...
Pirmuoju elementu x1 čia galima laikyti bet kokį teigiamą skaičių.
Įrodysime, kad ta seka konverguoja, o jos riba yra skaičius a. Pirmiausia įsitikinsime, kad seka {xn} turi ribą. Tam užtenka įrodyti, kad seka {xn} yra aprėžta iš apačios ir kad, pradedant antruoju elementu, ji nedidėja. Iš pradžių įsitikinsime, kad seka {xn} aprėžta iš apačios. Sąlygoje pasakyta, jog x1>0, o iš rekurentinės formulės, kai n=1, gauname x2>0. Taip samprotaudami toliau, įsitikiname, kad visi xn yra teigiami.
Dabar įrodysime, kad visi xn, kai n≥2, tenkina nelygybę xn≥a. Parašę rekurentinę formulę šitaip: xn+1=a2(xna+axn), remsimės beveik savaime aiškia nelygybe t+1t≥2*, kuri teisinga, kai t>0 (imame t=xna). Gausime
xn+1=a2(xna+axn)≥a2⋅2=a,
kai n≥1, t. y. xn≥a, pradedant numeriu n=2.
Pagaliau įsitikinsime, kad seka {xn}, kai n≥2, nedidėja. Iš rekurentinės formulės matome, kad xn+1xn=12(1+axn2), o iš čia, turėdami mintyje, kad xn≥a, gauname xn+1xn≤1, arba xn≥xn+1 (kai n≥2).
Kadangi seka {xn}, kai n≥2, nedidėja ir yra aprėžta iš apačios skaičiumi a (nes xn≥a), tai ji turi ribą, ne mažesnę kaip a (žr. 3.15 ir 3.13 teoremą). Tą ribą pažymėkime raide c ir atsižvelgę į tai, kad limn→∞xn+1=c, o limn→∞12(xn+axn)=12(c+ac), gauname lygybę
c=12(c+ac)**,
c2=12(c2+a),
2c2=c2+a,
2c2−c2=a,
c2=a.
Vadinasi, c=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 xn ir xn+1 įrašius ieškomąją sekos {xn} 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, a skaičiavimo algoritmą (įrodėme, kad limn→∞xn=a).
Vėliau bus tiriamas sekos {xn} konvergavimo į a greitis. Įrodoma, kad, tinkamai pasirinkus pirmąjį artinį x1, kai a>1, jau ketvirtasis artinys x4 nuo skaičiaus a skiriasi mažiau kaip 10−10.
3 pavyzdys. Įsitikinsime, kad sekos {cn}, kai cn=xn+1(n+1)!, o x – bet koks fiksuotas skaičius, riba lygi nuliui. Kai n – pakankamai didelis natūralusis skaičius, |x|n+1<1. Todėl, pradedant kuriuo nors numeriu N, turi būti |cn+1|<|cn|, nes |cn+1|=|x|nn!⋅|x|n+1=|cn|⋅|x|n+1.
Vadinasi, pradedant numeriu N, seka {|cn|} mažėja. Kadangi, be to, ji aprėžta iš apačios (pavyzdžiui nuliu), tai pagal 3.15 teoremą seka {|cn|} konverguoja. Tarkime, kad c yra tos sekos riba. Iš lygybės |cn+1|=|cn|⋅|x|n+1 išplaukia, kad c=0, nes sekos {|cn+1|} riba lygi c, o sekos {|x|n+1} – nuliui.

_____________________

* Tą nelygybę įrodinėjant, užtenka pastebėti, kad, kai t>0, ji ekvivalenti nelygybei t2−2t+1≥0.
** Ta lygybė gaunama iš rekurentinės formulės xn+1=12(xn+axn).


Sekos, aproksimuojančios a šaknyje, konvergavimo greitis

Šio skyriaus 3 paragrafo 3 skirsnyje įrodėme, kad sekos {xn} apibrėžiamos rekurentine formule
xn+1=12(xn+axn),n=1,2,3,...,(3.10)
kai a>0, o x1 – bet koks teigiamas skaičius, riba lygi a. Skaičiaus a artiniu galima laikyti bet kurį tos sekos narį xn+1. Tokiu atveju, be abejo, reikia išsiaiškinti, koks turi būti įteracijų* skaičius n, kad a artinys būtų nurodyto tikslumo.

______________

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

______________

Nagrinėsime seką xn, apibrėžtą rekurentine (3.10) formule. Tos sekos elementą xn vadinsime skaičiaus γ=a n-uoju artiniu. Skaičių
εn=xn−γγ(3.11)
vadinsime n-ojo artinio santykine paklaida.
Įrodysime teiginį apie santykinės paklaidos εn+1 įvertinimą, remiantis pirmojo artinio santykine paklaida ε1.
Jei x1 taip parinktas, kad |ε1|<12, tai su bet kokiu n≥1 teisinga šitokia nelygybė:
0<εn+1<ε12n.(3.12)
[ 0≤εn+1≤ε12n, kai ε1 gali būti lygus nuliui.]
Įrodymas. Iš (3.11) formulės
xn=γ(1+εn).(3.13)
Remdamiesi (3.10) ir (3.13) formulėmis bei lygybe aγ=γ, gauname
xn+1=12(xn+axn)=12[γ(1+ε)+aγ(1+εn)]=12[a(1+εn)2+aγ(1+εn)]=
=12[a(1+2εn+εn2)+aγ(1+εn)]=aγ[(1+2εn+εn2)+12(1+εn)]=γ[2+2εn+εn22(1+εn)]=γ[1+εn22(1+εn)].
Kadangi xn+1=γ(1+εn+1), tai savaime aišku, kad
γ(1+εn+1)=γ[1+εn22(1+εn)],
εn+1=12(1+εn)⋅εn2.(3.14)
Pagal sąlygą |ε1|<12. Iš to išplaukia nelygybės 0<12(1+εn)<1. Bet tada iš (3.14) lygybės, kai n=1, išplaukia nelygybė ε2>0. Panašiai naudodamiesi (3.14) lygybe, kai n=1,2,..., įsitikiname, kad visi εn+1 yra neneigiami (n≥1).
Iš (3.14) lygybės, nelygybių 0<12(1+εn)<1 ir iš to, kad visi εn, kai n>1, yra neneigiami, išplaukia nelygybė εn+1<εn2, kai n≥1. Iš čia gauname dešiniąją (3.12) nelygybę. Teiginys įrodytas.
[|ε1|<12, ε2<ε12, ε3<ε22<(ε12)2=ε14, ε4<ε32<(ε14)2=ε18=ε123 ir taip toliau.]
Iš (3.12) nelygybių matyti, kad a artinio po n iteracijų santykinė paklaida įvertinama, remiantis pirmojo artinio x1 santykine paklaida ε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 ε1 modulis nebūtų didesnis už 0.05. Savaime aišku, taip pasirinkus x1, santykinė paklaida tenkins įrodytojo teiginio sąlygas. Tuomet bus aišku ir koks turi būti iteracijų skaičius n, kad a artinio santykinė paklaida nebūtų didesnė už duotąjį skaičių ε: tą skaičių n galima rasti iš nelygybės* (ta nelygybė gaunama tiesiog iš (3.12) nelygybių)
(0.05)2n<ε.(3.15)
Taigi tarkime, kad a>1. Skaičių a išreikšime šitaip:
a=22k+iM;(3.16)
čia k – sveikasis neneigiamas skaičius, skaičius i lygus arba nuliui, arba vienetui, o skaičius M tenkina sąlygą
1≤M<2.(3.17)
Pabrėšime, kad skaičius a išreiškiamas (3.16) pavidalu vieninteliu budu.
Pirmąjį artinį x1 pasirinkime šitaip:
x1=2k(13⋅2iM+1724).(3.18)
Įsitikinsime, kad tuo atveju, kai M – bet koks skaičius, tenkinantis (3.17) sąlygas, pirmojo a artinio x1, rasto pagal (3.18) formulę, santikinės paklaidos ε1 modulis ne didesnis kaip 0.05.
Įrodinėdami šį teiginį, remsimės tikslia santykinės paklaidos išraiška ε1=x1−γγ. Iš (3.16) lygybės gauname γ=a=2k2iM; todėl iš ε1 išraiškos ir iš (3.18) formulės
ε1=x1−γγ=2k(13⋅2iM+1724)−2k2iM2k2iM=13⋅2iM+1724−2iM2iM.(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≥1, tai 2iM≥1. Todėl iš (3.19) lygybės išplaukia nelygybė
|ε1|≤|13⋅2iM+1724−2iM|.(3.20)
Pažymėkime 2iM raide X. Kadangi 1≤M<2, o i lygus arba nuliui, arba vienetui, tai visos galimos X reikšmės priklauso intervalui [1; 2):
1≤X<2.(3.21)
Naudodami naujai įvestą žymėjimą, (3.20) nelygybę perrašysime šitaip:
|ε1|≤|13X2−X+1724|.(3.22)
Pagal (3.22) nelygybę maksimalioji |ε1| reikšmė ne didesnė už maksimaliąją reiškinio |13X2−X+1724| reikšmę, kai X tenkina (3.21) sąlygas. Kad būtų aiškiau, išnagrinėkime funkcijos f(X)=13X2−X+1724 grafiką. Iš elementariosios matematikos kurso žinome, kad šios funkcijos grafikas yra parabolė (3.4 pav.)**, kurios viršūnė atitinka X=32.
[ f′(X)=(13X2−X+1724)′=23X−1,
23X−1=0,23X=1,X=32. ]
Kadangi f(1)=f(2)=124, o f(32)=−124, tai X reikšmės, tenkinančios (3.21) sąlygas, atitinka f(X) reikšmėms, esančioms tarp −124 ir 124. Kitaip sakant,
|f(X)|=|13X2−X+1724|≤124.
Iš paskutinės ir (3.22) nelygybių gauname reikiamą ε1 įvertinimą:
|ε1|≤124<0.05.
Pastaba. Jei nurodyta, kad santykinė paklaida ε lygi 10−10, tai, norint tokiu tikslumu apskaičiuoti bet kokio skaičiaus a>1 kvadratinę šaknį, reikia, suradus x1 pagal (3.18) formulę, atlikti tik tris iteracijas (n=3), nes
(0.05)23=0.058=0.0000000000390625=3.90625⋅10−11<10−10.
[ (1/24)8=(0.041(6))8≈9.0846890375538⋅10−12<10−11<10−10. Pasirodo, santykinė paklaida ε4 visada bus mažesnė už 10−11 po 3 iteracijų (n=3). ]

_______________

* Jei a<1, tai a=1b ir b>1. Tada a=1b.
** 3.4 paveiksle ašies Oy mastelis 20 kartų didesnis už ašies Ox mastelį.

Pavyzdžiai

  • Apskaičiuosime a=7. Iš (3.16) formulės a=22k+iM. Matome, kad 22k+i turi būti lygus 4, t. y. 22k+i=4. Tada aišku, kad k=1,i=0. Toliau randame M:
7=22k+iM,
7=4M,
M=74.
Toliau randame x1 pagal (3.18) formulę:
x1=2k(13⋅2iM+1724)=21(13⋅20⋅74+1724)=
=2(13⋅74+1724)=76+1712=14+1712=3112=
= 2.58(3).
Skaičių x2 randame pagal (3.10) rekurentinę formulę:
x1+1=12(x1+ax1)=12(3112+73112)=12(3112+8431)=12⋅312+84⋅1212⋅31=12⋅961+1008372=1969744=
= 2.646505376344086021505376344086.
Toliau rasime x3:
x3=12(x2+ax2)=12(1969744+71969744)=12(1969744+52081969)=12⋅19692+5208⋅744744⋅1969=77517132929872=
= 2.6457514184920023809913880196814.
[Atėmę tikslią 7 reikšmę iš x3 reikšmės gauname:
2.6457514184920023809913880196814 - 2.6457513110645905905016157536393 = 1.0742741179048977226604214545141e-7.]
Toliau rasime x4:
x4=x3+1=12(x3+ax3)=12(77517132929872+777517132929872)=12(77517132929872+205091047751713)=
=12⋅77517132+20509104⋅29298722929872⋅7751713=12⋅60,089,054,434,369+60,089,049,554,6882929872⋅7751713=120,178,103,989,05745,423,053,741,472=
= 120178103989057/45423053741472 = 2.6457513110645927714791488926518.
Atėmę tikslią 7 reikšmę iš gautos x4 reikšmės, gauname:
2.6457513110645927714791488926518 - 2.6457513110645905905016157536393 = 2.1809775331390125063383535564728e-15.
Tai yra gavome paklaidą ε≈2.180977533139⋅10−15<10−10.
Santykinė paklaida yra tokia:
ε4=x4−γγ=
(2.6457513110645927714791488926518 - 2.6457513110645905905016157536393)/2.6457513110645905905016157536393 =
= 8.2433202395785127055094193690687e-16.
Taigi ε4≈8.2433202395785⋅10−16<10−10.


  • Apskaičiuosime a=39. Iš (3.16) formulės a=22k+iM. Pagal (3.17) 1≤M<2. Todėl 22k+i=32. Tada 2k+i=5. Iš čia gauname, kad k=2,i=1.
Toliau randame M:
a=22k+iM,
39=22⋅2+1M,
39=25M,
M=3932.
Toliau randame x1 pagal (3.18) formulę:
x1=2k(13⋅2iM+1724)=22(13⋅21⋅3932+1724)=4(23⋅3932+1724)=31296+176=312+17⋅1696=58496=7312=
= 6.08(3).
[ 39= 6.2449979983983982058468931209398. ]
Skaičių x2 randame pagal (3.10) rekurentinę formulę:
x2=12(x1+ax1)=12(7312+397312)=
= 0.5*(6.0833333333333333333333333333333 + 6.4109589041095890410958904109589) =
= 0.5 * 12.494292237442922374429223744292 = 6.2471461187214611872146118721461.
Toliau rasime x3:
x3=12(x2+ax2)=
= 0.5*(6.2471461187214611872146118721461 + 39/6.2471461187214611872146118721461) =
= 6.2449983677207123204232036062421.
[Atėmę tikslią 39 reikšmę iš x3 reikšmės, gauname:
6.2449983677207123204232036062421 - 6.2449979983983982058468931209398 = 3.693223141145763104853022918526e-7.]
Toliau rasime x4:
x4=x3+1=12(x3+ax3)=
= 0.5*(6.2449983677207123204232036062421 + 39/6.2449983677207123204232036062421) =
= 6.2449979983984091265040071378015.
Atėmę tikslią 39 reikšmę iš x4 reikšmės, gauname:
6.2449979983984091265040071378015 - 6.2449979983983982058468931209398 = 1.0920657114016861678930991073622e-14.
Taigi, gavome paklaidą 1.0920657114⋅10−14. O santykinė paklaida ε4 lygi
ε4=x4−γγ=
= (6.2449979983984091265040071378015 - 6.2449979983983982058468931209398)/6.2449979983983982058468931209398 =
= 1.7487046620059110065547770263338e-15.
T. y. ε4≈1.748704662⋅10−15<10−10.


  • Apskaičiuosime a=2. Iš (3.16) formulės a=22k+iM. Pagal (3.17) 1≤M<2. Todėl 22k+i=2. Tada 2k+i=1. Iš čia gauname, kad k=0,i=1.
Kadangi 22k+i=2, tai M=1.
Toliau randame x1 pagal (3.18) formulę:
x1=2k(13⋅2iM+1724)=20(13⋅21⋅1+1724)=23+1724=2⋅8+1724=3324=1.375.
Skaičių x2 randame pagal (3.10) rekurentinę formulę:
x2=12(x1+ax1)=12(1.375+21.375)=
= 0.5*(1.375 + 1.4545454545454545454545454545455) = 0.5 * 2.8295454545454545454545454545455 =
= 1.4147727272727272727272727272727.
Toliau rasime x3:
x3=12(x2+ax2)=
= 0.5*(1.4147727272727272727272727272727 + 2/1.4147727272727272727272727272727) =
= 1.4142136728733114275282949981745.
[1.4142136728733114275282949981745 - 1.4142135623730950488016887242097 = 0.0000001105002163787266062739648 ≈1.105002163787266⋅10−7.]
Toliau rasime x4:
x4=12(x3+ax3)=
= 0.5*(1.4142136728733114275282949981745 + 2/1.4142136728733114275282949981745) =
= 1.4142135623730993657935457387426.
Atėmę tikslią 2 reikšmę iš x4 reikšmės, gauname:
1.4142135623730993657935457387426 - 1.4142135623730950488016887242097 = 4.316991857014532913283559342618e-15.
Santykinė paklaida yra tokia:
ε4=x4−γγ=
=(1.4142135623730993657935457387426 - 1.4142135623730950488016887242097)/1.4142135623730950488016887242097 =
= 3.0525742164220827729222598916484e-15.
Gavome, kad ε4≈3.052574216422⋅10−15<10−10.


Šiaip, skaičiuojant su kompiuteriu, santikinę paklaidą ε4 (arba ε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 δ<10εn+1.

Šaknies traukimas iš didelių skaičių

Bet kokį skaičių a (a>1) galima užrašyti taip:
a=p⋅10j⋅10v.(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:
a=p⋅10j⋅10v2.(2)
Jegiu v nelygus nuliui, tai 10v/2 bus sveikasis skaičius su dvigubai mažiau nulių nei pas 10v. Jeigu j=1, tai tereikia apskaičiuoti 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=22k+iM;
iš čia k=0,i=0,M=p.
Tada pagal (3.18) formulę randame x1:
x1=2k(13⋅2iM+1724)=20(13⋅20⋅p+1724)=13p+1724.
Ir toliau, darant iteracijas, apskaičiuojamas x2,x3,x4,... pagal (3.10) formulę:
xn+1=12(xn+pxn),n=1,2,3,...
2) Jeigu 2≤p<4, tai pagal (3.16) formulę
p=22k+iM;
iš čia k=0,i=1. Tada p=2M,M=p2.
Pagal (3.18) formulę randame x1:
x1=2k(13⋅2iM+1724)=20(13⋅21⋅p2+1724)=13p+1724.
Matome, kad šiuo atveju artinio x1 formulė tokia pati kaip ir 1) atveju.
Toliau skaičiuojama kaip ir pirmu atveju.
3) Jeigu 4≤p<8, tai pagal (3.16) formulę
p=22k+iM;
iš čia k=1,i=0. Tada p=4M,M=p4.
Pagal (3.18) formulę randame x1:
x1=2k(13⋅2iM+1724)=21(13⋅20⋅p4+1724)=16p+1712.
Toliau skaičiuojama kaip ir pirmu atveju.
4) Jeigu 8≤p<10, tai pagal (3.16) formulę
p=22k+iM;
iš čia k=1,i=1. Tada p=8M,M=p8.
Pagal (3.18) formulę randame x1:
x1=2k(13⋅2iM+1724)=21(13⋅21⋅p8+1724)=16p+1712.
Matome, kad šiuo atveju artinio x1 formulė tokia pati kaip ir 3) atveju.
Toliau skaičiuojama kaip ir pirmu atveju.
Apibendrinant, skaičiaus p (1<p<10) apskaičiavimui, taikoma x1 reikšmė
x1=13p+1724,(3)
kai 1<p<4;
ir taikoma x1 reikšmė
x1=16p+1712,(4)
kai 4≤p<10.


  • Apskaičiuosime a=1.5⋅108. Čia p=1.5,j=0,v=8. Tada pagal (1) formulę
a=p⋅10j⋅10v=1.5⋅100⋅108.
Pagal (2) formulę
a=p⋅10j⋅10v2=p⋅1082=p⋅104.
Toliau randame x1 skaičiui p=1.5:
x1=13p+1724=13⋅32+1724=12+1724=12+1724=2924=
= 1.2083333333333333333333333333333.
[ 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 x2:
x2=x1+1=12(x1+px1)=12(2924+1.52924)=
= 0.5*(1.2083333333333333333333333333333 + 1.5/1.2083333333333333333333333333333) = 1.2248563218390804597701149425287.
Toliau pagal (3.10) formulę apskaičiuojame x3:
x3=x2+1=12(x2+px2)=
= 0.5*(1.2248563218390804597701149425287 + 1.5/1.2248563218390804597701149425287) = 1.2247448764620622240199548319682.
[1.2247448764620622240199548319682 - 1.2247448713915890490986420373529 = 0.0000000050704731749213127946153 = 5.07047⋅10−9.]
Toliau pagal (3.10) formulę apskaičiuojame x4:
x4=x3+1=12(x3+px3)=
= 0.5*(1.2247448764620622240199548319682 + 1.5/1.2247448764620622240199548319682) = 1.2247448713915890595945823395751.
Iš ką tik gautos x4 reikšmės atėmę 1.5 reikšmę, gauname:
1.2247448713915890595945823395751 - 1.2247448713915890490986420373529 = 0.0000000000000000104959403022222 = 1.049594⋅10−17.
Santykinė paklaida yra tokia:
ε4=x4−γγ=
(1.2247448713915890595945823395751 - 1.2247448713915890490986420373529)/1.2247448713915890490986420373529 =
= 8.5698993703859497564530032133562e-18.
Taigi ε4≈8.56989937⋅10−18<10−10.
Iš anksčiau apskaičiuoto
a=p⋅104=x4⋅104=
= 1.2247448713915890595945823395751 * 10000 = 12247.448713915890595945823395751.
Atėmę tikslią a iš 12247.448713915890595945823395751, gauname paklaidą:
12247.448713915890595945823395751 - 12247.448713915890490986420373529 = 1.0495940302222154304017026259951e-13.
O santykinė paklaida skaičiui a yra
(12247.448713915890595945823395751 - 12247.448713915890490986420373529)/12247.448713915890490986420373529 =
= 8.5698993703859497564530032133562e-18.
Tokia pati, kaip skaičiui p.


  • Apskaičiuosime a=9.6⋅1011. Čia p=9.6,j=1,v=10. Tada pagal (1) formulę
a=p⋅10j⋅10v=9.6⋅101⋅1010.
Pagal (2) formulę
a=p⋅10j⋅10v2=p⋅10⋅10102=p⋅10⋅105.
Šaknį iš 10 lengva apskaičiuoti (nes (4) formulė tinka skaičiams 4≤a<16).
10= 3.1622776601683793319988935444327.
Toliau randame x1 skaičiui p=9.6 pagal (4) formulę:
x1=16p+1712=16⋅9.6+1712=
= 0.16666666666666666666666666666667*9.6 + 1.4166666666666666666666666666667 = 3.0166666666666666666666666666667.
[ 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 x2:
x2=x1+1=12(x1+px1)=
= 0.5*(3.0166666666666666666666666666667 + 9.6/3.0166666666666666666666666666667) = 3.0994935543278084714548802946593.
Toliau pagal (3.10) formulę apskaičiuojame x3:
x3=x2+1=12(x2+px2)=
= 0.5*(3.0994935543278084714548802946593 + 9.6/3.0994935543278084714548802946593) = 3.0983868746074953885623653564026.
[3.0983868746074953885623653564026 - 3.0983866769659335081434123198259 = 0.0000001976415618804189530365767 =
≈1.9764156188⋅10−7.]
Toliau pagal (3.10) formulę apskaičiuojame x4:
x4=x3+1=12(x3+px3)=
= 0.5*(3.0983868746074953885623653564026 + 9.6/3.0983868746074953885623653564026) = 3.0983866769659398117763289613096.
Iš ką tik gautos x4 reikšmės atėmę 9.6 reikšmę, gauname:
3.0983866769659398117763289613096 - 3.0983866769659335081434123198259 = 0.0000000000000063036329166414837 =
≈6.3036329166414837⋅10−15.
Santykinė paklaida yra tokia:
ε4=x4−γγ=
(3.0983866769659398117763289613096 - 3.0983866769659335081434123198259)/3.0983866769659335081434123198259 =
= 2.0344887755631127106428103904763e-15.
Taigi ε4≈2.03448877556⋅10−15<10−10.
Iš ankstesnių skaičiavimų
a=p⋅10⋅105=x4⋅10⋅105=
= 3.0983866769659398117763289613096 * 3.1622776601683793319988935444327 * 100000 = 979795.89711327323266266864962311.
Tiksli a reikšmė yra tokia:
a=9.6⋅1011= 979795.89711327123927891362988236.
Atėmę tikslią a reikšmę iš 979795.89711327323266266864962311, gauname paklaidą:
979795.89711327323266266864962311 - 979795.89711327123927891362988236 = 0.00000000199338375501974075 =
≈1.993383755⋅10−9.
Santykinė paklaida skaičiui 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 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=1a, galime, jei b didelis skaičius, apskaičiuoti b pagal pateiktą aukščiau būdą. Apskaičiavę b, galime rasti a taip:
a=1b.
Tačiau toks a (0<a<1) skaičiavimo būdas nėra ekonomiškas, nes reikia atlikti dvi [papildomas] dalybos operacijas.
Pateiksime ekonomiškesnį, greitesnį a skaičiavimo būdą. Jei 0<a<1, tai a galima užrašyti taip:
a=p⋅10j⋅10v.(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:
a=p⋅10j⋅10v2.(6)
Jeigu j=−1, tai tereikia apskaičiuoti (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 a. Kaip rasti 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į x1 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 a iš vienos iteracijos, artinio x1 santykinė paklaida ε1 turi būti apie 10−9 (kad gauti a su santykine paklaida apie 10−18). Jeigu kažkokiu budu galima pagal skaičių p pažiūrėt į kažkokią artinių x1 lentelę, tai ta lentelė turės apie milijardą artinių ir užims apie 80(bit)⋅109= 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 x1 skaičiavimui [pateiktam aukščiau] nereikia daryti dalybos operacijos, tai galbūt yra, tarkim, apie 100 ar 1000 artinio x1 skaičiavimo grupių skaičiui p. Kaip matėme, tik dviejų grupių (1<p<4 ir 4≤p<10) užtenka p skaičiavimui, kad gauti pirmo artinio x1 santykinę paklaidą ε1<1/24<0.05. Jeigu sukurti apie 10 tokių artinio x1 grupių (1<p<2, 2<p<3, 3<p<4, ..., 9<p<10) arba apie 100 tokių artinio x1 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 ε1 tikslumas gali būti apie 10−10. Tik tai gali užimti darbo suradimui artinio x1 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į x1). 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į x1 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] x1 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 x1=a, tai pagal (3.10) formulę
x2=12(x1+ax1)=12(a+aa)=12(a+a)=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 ≈1.11607142857⋅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ą a0 (a0=a arba a0=b). Toliau taške a0 surandame funkcijos f(x) liestinės lygtį, kuri atrodo taip:
y=f(a0)+f′(a0)(x−a0).
Arba taip:
y−f(a0)=f′(a0)(x−a0).(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(a0)=f′(a0)(d−a0),
−f(a0)f′(a0)=d−a0,
d=a0−f(a0)f′(a0).(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=a1 yra arčiau šaknies c. Pakartoję šį procesą rasime tašką a2, esantį dar arčiau prie c taško (ant Ox ašies):
a2=a1−f(a1)f′(a1).
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 an 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 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/a. Kai apytiksliai apskaičiuojama 1/a, reikia gautą rezultatą padauginti iš a. Tada ir gausime a. Be to 1/a rezultatas yra net naudingesnis, nes nereikia dalinti a iš a, kad gauti 1/a.
Pereikime prie 1a skaičiavimo algoritmo, kuris vadinasi Niutono metodu [šaknies skaičiavimui]. Pažymėkime y=1a, o skaičių a pažymėkime per x. Tada
f(y)=1y2−x=0.
Jeigu yn yra apytikslė skaičiaus y reikšmė, tai yn+1 yra dar tikslesnė skaičiaus y reikšmė (pagal liestinių metodą):
yn+1=yn−f(yn)f′(yn).(3)
Surandame f(y) išvestinę:
f′(y)=(1y2−x)′=−2y3.
f′(yn)=−2yn3.
Toliau įstatome f′(yn) ir f(yn) į (3) formulę:
yn+1=yn−f(yn)f′(yn)=yn−1yn2−x−2yn3=yn+yn3yn2−xyn32=
=yn+yn−xyn32=yn+yn2−xyn32=3yn2−xyn32=
=yn(32−x2yn2)=yn(3−xyn2)2.
Gavome rekurentinę formulę 1/a reikšmei skaičiuoti:
yn+1=yn(32−x2yn2)=yn(3−xyn2)2.(4)

Pavyzdžiai

  • Apskaičiuosime 1x=12.3 ir rasime kam lygus 2.3=a=xx=2.32.3.
Pirmiausia turime žinoti y=12.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 2.3 reikšmė yra tokia:
2.3= 1.5165750888103101108513650872564.
O y=12.3= 0.65938047339578700471798482054627.
Mes paimsime, kad y1=0.6. Tada pagal formulę
yn+1=yn(32−x2yn2)=yn(3−xyn2)2
rasime
y2=y1(3−xy12)2=0.6(3−2.3⋅0.62)2=0.6(3−0.828)2=0.6⋅2.1722=1.30322=0.6516.
Toliau pagal (4) formulę rasime y3:
y3=y2(3−xy22)2=0.6516(3−2.3⋅0.65162)2=0.6516(3−2.3⋅0.42458256)2=0.6⋅2.0234601122=1.31848660897922=0.6592433044896.
[0.65938047339578700471798482054627 - 0.6592433044896 = 0.00013716890618700471798482054627 ≈1.37168906187⋅10−4.]
Toliau pagal (4) formulę rasime y4:
y4=y3(3−xy32)2=
= 0.6592433044896*(3 - 2.3*0.6592433044896^2)/2 = 0.65938043059651203399501910136126.
[0.65938047339578700471798482054627 - 0.65938043059651203399501910136126 =
= 0.00000004279927497072296571918501 ≈4.279927497⋅10−8.]
Toliau pagal (4) formulę rasime y5:
y5=y4(3−xy42)2=
= 0.65938043059651203399501910136126*(3 - 2.3*0.65938043059651203399501910136126^2)/2 = 0.65938047339578283767489142997287.
Taigi, y5 paklaida yra tokia:
0.65938047339578700471798482054627 - 0.65938047339578283767489142997287 =
= 0.0000000000000041670430933905734 =4.1670430933905734⋅10−15.
Po 4 iteracijų gavome paklaidą mažesne nei 10−14.
Dabar galime rasti a šaknyje:
2.3=a=x=2.3⋅y5= 2.3*0.65938047339578283767489142997287 = 1.5165750888103005266522502889376.
[1.5165750888103101108513650872564 - 1.5165750888103005266522502889376 = 0.0000000000000095841991147983188 =
=9.5841991147983188⋅10−15.]
Pažiūrėkime kaip konverguoja yn+1 į 12.3 iš dešinės pusės, kai y1 parenkamas taip, kad y1>12.3.
Taigi parinkime y1=0.66. Tada
y2=y1(3−xy12)2=0.66(3−2.3⋅0.662)2=0.66(3−2.3⋅0.4356)2=0.66(3−1.00188)2=0.66⋅1.998122=1.31875922=0.6593796.
[0.65938047339578700471798482054627 - 0.6593796 = 0.00000087339578700471798482054627 ≈8.73395787⋅10−7.]
Toliau pagal (4) formulę rasime y3:
y3=y2(3−xy22)2=0.6593796(3−2.3⋅0.65937962)2=0.6593796(3−0.999997350861168)2=
=0.6593796⋅2.0000026491388322= 1.3187609467881033886272/2 = 0.6593804733940516943136.
[0.65938047339578700471798482054627 - 0.6593804733940516943136 = 0.00000000000173531040438482054627 ≈1.7353104⋅10−12.]
Toliau pagal (4) formulę rasime y4:
y4=y3(3−xy32)2=
0.6593804733940516943136*(3 - 2.3*0.6593804733940516943136^2)/2 = 0.65938047339578700471797797024742.
Paklaida lygi
0.65938047339578700471798482054627 - 0.65938047339578700471797797024742 =
= 0.00000000000000000000000685029885 =6.85029885⋅10−24.
Po 3 iteracijų gavome paklaidą mažesne nei 10−23.
Šį kartą, parenkant y1 iš dešinės, y1 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 į 12.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)=1y2−x yra išimtis.
Dabar galime rasti a šaknyje:
2.3=a=x=x⋅y4=2.3⋅y4= 2.3*0.65938047339578700471797797024742 = 1.5165750888103101108513493315691.
[1.5165750888103101108513650872564 - 1.5165750888103101108513493315691 = 0.0000000000000000000000157556873 =
=1.57556873⋅10−23.]


  • Apskaičiuosime 1x=15.8 ir rasime kam lygus 5.8=x.
Tiksli 5.8 reikšmė yra tokia:
5.8= 2.4083189157584590960256482060757.
O y=15.8= 0.41522739926869984414235313897857.
Mes paimsime, kad y1=0.4. Tada pagal (4) formulę
y2=y1(3−xy12)2=0.4(3−5.8⋅0.42)2=0.4(3−0.928)2=0.4⋅2.0722=0.82882=0.4144.
[0.41522739926869984414235313897857 - 0.4144 = 0.00082739926869984414235313897857 ≈8.27399268699844⋅10−4.]
Toliau pagal (4) formulę rasime y3:
y3=y2(3−xy22)2=
= 0.4144*(3 - 5.8*0.4144^2)/2 = 0.4152249278464.
[0.41522739926869984414235313897857 - 0.4152249278464 = 0.00000247142229984414235313897857 ≈2.471422299844⋅10−6.]
Toliau pagal (4) formulę rasime y4:
y4=y3(3−xy32)2=
= 0.4152249278464*(3 - 5.8*0.4152249278464^2)/2 = 0.41522739924663512944560082817282.
[0.41522739926869984414235313897857 - 0.41522739924663512944560082817282 = 0.00000000002206471469675231080575 =
≈2.206471469675231⋅10−11.]
Toliau pagal (4) formulę rasime y5:
y5=y4(3−xy42)2=
= 0.41522739924663512944560082817282*(3 - 5.8*0.41522739924663512944560082817282^2)/2 = 0.41522739926869984414059439797726.
[0.41522739926869984414235313897857 - 0.41522739926869984414059439797726 = 0.00000000000000000000175874100131 =
=1.75874100131⋅10−21.]
Tada 5.8=x=5.8⋅y5= 5.8*0.41522739926869984414059439797726 = 2.4083189157584590960154475082681.
[2.4083189157584590960256482060757 - 2.4083189157584590960154475082681 = 0.0000000000000000000102006978076 =
=1.02006978076⋅10−20.]
Po 4 iteracijų gavome paklaidą mažesne nei 10−19 (beveik 10−20).
Dabar pažūrėsime kaip konverguos yn į 1/5.8 iš dešinės pusės (kai y1>1/5.8). Parenkame y1=0.42. Tada pagal (4) formulę
y2=y1(3−xy12)2=0.42(3−5.8⋅0.422)2=0.42(3−1.02312)2=0.42⋅1.976882=0.4151448.
[0.41522739926869984414235313897857 - 0.4151448 = 0.00008259926869984414235313897857 ≈8.2599268699844⋅10−5.]
Toliau pagal (4) formulę rasime y3:
y3=y2(3−xy22)2=
= 0.4151448*(3 - 5.8*0.4151448^2)/2 = 0.4152273746236976005632.
[0.41522739926869984414235313897857 - 0.4152273746236976005632 = 0.00000002464500224357915313897857 ≈2.4645002243579⋅10−8.]
Toliau pagal (4) formulę rasime y4:
y4=y3(3−xy32)2=
= 0.4152273746236976005632*(3 - 5.8*0.4152273746236976005632^2)/2 = 0.41522739926869765000924208032344.
[0.41522739926869984414235313897857 - 0.41522739926869765000924208032344 = 0.00000000000000219413311105865513 =
≈2.194133111⋅10−15.]
Toliau pagal (4) formulę rasime y5:
y5=y4(3−xy42)2=
= 0.41522739926869765000924208032344*(3 - 5.8*0.41522739926869765000924208032344^2)/2 = 0.41522739926869984414235313896118.
[0.41522739926869984414235313897857 - 0.41522739926869984414235313896118 = 0.00000000000000000000000000001739 =
=1.739⋅10−29.]
Tada 5.8=x=5.8⋅y5= 5.8*0.41522739926869984414235313896118 = 2.4083189157584590960256482059748.
[2.4083189157584590960256482060757 - 2.4083189157584590960256482059748 = 0.0000000000000000000000000001009 =
=1.009⋅10−28.]
Po 4 iteracijų gavome 5.8 paklaidą mažesne nei 10−27 (beveik 10−28).
[2.4083189157584590960256482059748^2 = 5.7999999999999999999999999995139.]


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

Š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 1z32=1z1/32 taikant Niutono metodą. Pažymėkime y=1z1/32. Tada turime funkciją
f(y)=1y32−z=0.
Pagal liestinių metodą
yn+1=yn−f(yn)f′(yn).(3)
Surandame f(y) išvestinę:
f′(y)=(1y32−z)′=−32y33.
f′(yn)=−32yn33.
Toliau įstatome f′(yn) ir f(yn) į (3) formulę:
yn+1=yn−f(yn)f′(yn)=yn−1yn32−z−32yn33=yn+yn33yn32−zyn3332=
=yn+yn−zyn3332=yn+yn32−zyn3332=33yn32−zyn3332=
=yn(3332−z32yn32)=yn(33−zyn32)32.
Gavome rekurentinę formulę 1z32 reikšmei skaičiuoti:
yn+1=yn(33−zyn32)32.(5)

Pavyzdžiai

  • Apskaičiuosime 1z32=11.332 ir rasime kam lygus 1.332=z32. Prieš tai dar pastebėsime, kad norint apskaičiuoti z32 reikia 5 kartus ištraukti kvadratinę šaknį iš z. Tiksli 1.332 reikšmė yra tokia:
1.332= 1.0082325861537202021533989768813.
O y=11.332= 0.99183463590963025992221676151568.
Mes paimsime, kad y1=0.99. Tada pagal (5) formulę
y2=y1(33−zy132)32=0.99(33−1.3⋅0.9932)32=
= 0.99*(33 - 1.3*0.72498033595785364231768779197961)/32 = 0.99177969711319507382303549411632.
[0.99183463590963025992221676151568 - 0.99177969711319507382303549411632 = 0.00005493879643518609918126739936 ≈5.49⋅10−5.]
Toliau pagal (5) formulę rasime y3:
y3=y2(33−zy232)32=
= 0.99177969711319507382303549411632*(33 - 1.3*0.99177969711319507382303549411632^32)/32 = 0.99183458572688580323137637227142.
[0.99183463590963025992221676151568-0.99183458572688580323137637227142=0.00000005018274445669084038924426 ≈5.018⋅10−8.]
Toliau pagal (5) formulę rasime y4:
y4=y3(33−zy332)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 =
≈4.1894138549985⋅10−14.]
Toliau pagal (5) formulę rasime y5:
y5=y5(33−zy532)32=
= 0.99183463590958836578366677638205*(33 - 1.3*0.99183463590958836578366677638205^32)/32 = 0.9918346359096302599222167323178.
Atėmę, ką tik gautą y5 reikšmę iš tikslios y=11.332 reikšmės, gauname:
0.99183463590963025992221676151568 - 0.9918346359096302599222167323178 = 0.00000000000000000000000002919788 =
≈2.919788⋅10−26.
Gavome po keturių iteracijų paklaidą mažesnę už 10−25.
Padalijus vienetą iš y5 reikšmės gauname
1.332=z32=1y5=
= 1/0.9918346359096302599222167323178 = 1.0082325861537202021533990065619.
[1.0082325861537202021533989768813 - 1.0082325861537202021533990065619 = -0.0000000000000000000000000296806 =
≈−2.96806⋅10−26.]
1.0082325861537202021533990065619^32 = 1.3000000000000000000000012246312.
(kad pakelti 1.0082325861537202021533990065619 32 laipsniu reikia ant kalkuliatoriaus 5 kartus paspasuti mygtuką x2.)
Toliau apskaičiuosime ln(1.3). Kadangi
ln⁡1.3=−ln⁡1.3−1=−32ln⁡1.3−1/32=−32ln⁡1z1/32≈−32ln⁡y5,
tai tereikia apskaičiuoti ln⁡y5. Gauname
ln⁡y5= 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−x22+x33−x44...+(−1)n−1xnn+rn+1(x);(13)
čia |−32rn+1(x)|=|−32xn+1(1+c)n+1(n+1)|,−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ą ≈1.327368678⋅10−22<10−21.
Įvertinsime liekamąjį narį teoriškai.
|−32r9+1(x)|<|−32⋅(−0.008165364)9+1(1−0.008165364)9+1(9+1)|=
=| -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)≈ -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⋅10−21<|−4.576282447⋅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))
yn+1=yn(33−zyn32)32.(5)
Kad yn pakelti 32 laipsniu, reikia tiesiog yn penkis kartus pakelti kvadratu. Tai jau yra 5 daugybos operacijos. Toliau dar reikia yn padauginti iš z. Jau turime 6 daugybos operacijas. Tada reikia zyn atimti iš 33. Gauname vieną atimties operaciją. Ir galiausiai viską padauginti iš yn 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 y1, 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/y4, 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:
yn+1=yn(3−xyn2)2.(4)
Kad pakelti yn kvadratu reikia vienos daugybos operacijos. Paskui dar vienos daugybos operacijos reikia, kad padauginti yn2 iš x. Jau turime dvi daugybos operacijas. Čia yra dar viena atimties operacija. Dar viena daugybos operacija reiškinio skliaustuose iš yn 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į y1, 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ų y1=1/z1/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 1z32=10.7532 ir rasime kam lygus 0.7532=z32. Tiksli 0.7532 reikšmė yra tokia:
0.7532= 0.99105022504172762555415281377767.
O y=10.7532= 1.0090305967670766121721228310046.
Mes paimsime, kad y1=1.009. Tada pagal (5) formulę
y2=y1(33−zy132)32=1.009(33−0.75⋅1.00932)32=
= 1.009*(33 - 0.75*1.3320401627616656376316134288889)/32 = 1.0090305814634409227725711418028.
[1.0090305967670766121721228310046 - 1.0090305814634409227725711418028 = 0.0000000153036356893995516892018 ≈1.53⋅10−8.]
Toliau pagal (5) formulę rasime y3:
y3=y2(33−zy232)32=
= 1.0090305814634409227725711418028*(33 - 0.75*1.0090305814634409227725711418028^32)/32 = 1.0090305967670727824366475789071.
[1.0090305967670766121721228310046 - 1.0090305967670727824366475789071 = 0.0000000000000038297354752520975 ≈3.8297⋅10−15.]
Toliau pagal (5) formulę rasime y4:
y4=y3(33−zy332)32=
= 1.0090305967670727824366475789071*(33 - 0.75*1.0090305967670727824366475789071^32)/32 = 1.0090305967670766121721228307648.
Atėmę gautą y4 reikšmę iš tikslios y=10.7532 reikšmės, gauname:
1.0090305967670766121721228310046 - 1.0090305967670766121721228307648 = 0.0000000000000000000000000002398 ≈2.398⋅10−28.
Tik po 3 iteracijų paklaida yra mažesnė nei 10−27. Taip yra todėl, nes parinkome y1=1.009 su paklaida
1.0090305967670766121721228310046 - 1.009 = 0.0000305967670766121721228310046 ≈3.0596767⋅10−5.
Na o,
z32=0.7532≈1/y4= 1/1.0090305967670766121721228307648 = 0.99105022504172762555415281401323.
[0.99105022504172762555415281377767 - 0.99105022504172762555415281401323 = -0.00000000000000000000000000023556 =
≈−2.3556⋅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 y4 paklaidos atsakyme.]
Toliau apskaičiuosime ln(0.75). Kadangi
ln⁡0.75=−ln⁡0.75−1=−32ln⁡0.75−1/32=−32ln⁡1z1/32≈−32ln⁡y4,
tai tereikia apskaičiuoti ln⁡y4. Gauname
ln⁡y4= 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−x22+x33−x44...+(−1)n−1xnn+rn+1(x);(13)
čia |−32rn+1(x)|=|−32xn+1(1+c)n+1(n+1)|,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⋅10−22<10−21.
Įvertinsime liekamąjį narį teoriškai.
|−32r9+1(x)|<|−32⋅0.0090305967679+1(1+0)9+1(9+1)|=
=| -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)≈ -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⋅10−20<|−32rn+1(x)|=|−1.1542888083⋅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 Fi. Skaičius Fi pasirinktas taip, kad dauginant vardiklį kiekvieną kartą iš Fi, vardiklis konverguoja į vienetą. Taip kelis kartus dauginant skaitiklį ir vardiklį iš Fi skaitiklis konverguoja į Q=N/D; čia N yra dalinys (skaitiklis), D yra daliklis (vardiklis). Formulė Q skaičiuoti yra tokia:
Q=NDF1F1F2F2F…F….

Žingsniai Goldschmidt'o dalybai yra:

  1. Generuoti Fi .
  2. Padauginti dalinį ir daliklį iš Fi .
  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 Fi randamas iš Di (D1 = D):
Fi+1=2−Di.
Dauginant dalinį Ni ir daliklį Di iš faktoriaus Fi duoda Ni+1 ir Di+1:
Ni+1Di+1=NiDiFi+1Fi+1.
Po pakankmo iteracijų skaičiaus k, Q=Nk.

Pavyzdžiai

  • Apskaičiuosime Q=N/D=3/7. Padauginame 7 iš 0.1: D1=0.1⋅7=0.7. Tada
Fi+1=2−Di,
F1+1=F2=2−D1=2−0.7=1.3.
Toliau padauginame N1=3 ir D1=0.7 iš F2:
N2=N1⋅F2=3⋅1.3=3.9,
D2=D1⋅F2=0.7⋅1.3=0.91.
Toliau apskaičiuojame F3:
F3=2−D2=2−0.91=1.09.
Toliau padauginame N2 ir D2 iš F3:
N3=N2⋅F3=3.9⋅1.09=4.251,
D3=D2⋅F3=0.91⋅1.09=0.9919.
Toliau apskaičiuojame F4:
F4=2−D3= 2 - 0.9919 = 1.0081.
Toliau padauginame N3 ir D3 iš F4:
N4=N3⋅F4= 4.251*1.0081 = 4.2854331,
D4=D3⋅F4= 0.9919*1.0081 = 0.99993439.
Toliau apskaičiuojame F5:
F5=2−D4= 2 - 0.99993439 = 1.00006561.
Toliau padauginame N3 ir D3 iš F4:
N5=N4⋅F5= 4.2854331*1.00006561 = 4.285714267265691,
D5=D4⋅F5= 0.99993439*1.00006561 = 0.9999999956953279.
Taigi Q=N5⋅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 ≈1.84485947⋅10−9.]
Po 4 iteracijų gavome paklaidą mažesnę nei 10−8.
Dar padarysime penktą įteraciją, apskaičiuodami F6, N6 ir D6:
F6=2−D5= 2 - 0.9999999956953279 = 1.0000000043046721;
N6=N5⋅F6= 4.285714267265691*1.0000000043046721 = 4.2857142857142856348705633349211,
D6=D5⋅F6= 0.9999999956953279*1.0000000043046721 = 0.99999999999999998146979811148159.
Taigi Q=N6⋅0.1= 4.2857142857142856348705633349211*0.1 = 0.42857142857142856348705633349211.
[0.42857142857142857142857142857143 - 0.42857142857142856348705633349211 = 0.00000000000000000794151509507932 =
≈7.941515⋅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į):
b0=S
Y0≈1/S (typically using a table lookup)
y0=Y0
x0=Sy0
ir daryti iteracijas
bn+1=bnYn2
Yn+1=12(3−bn+1)
xn+1=xnYn+1
yn+1=ynYn+1
kol bi yra pakankamai arti 1, arba daryti fiksuotą (tam tikrą) kiekį iteracijų. Iteracijos konverguoja prie
limn→∞xn=S
ir
limn→∞yn=1/S.
Pastebėsime, kad galima neskačiuoti xn ir yn jeigu abu yra reikalingi. Galima tik paskaičiuoti xn=Syn pačioje skaičiavimų pabaigoje (tada yn vis tiek reikia skaičiuoti).
Antras būdas, naudojant fused multiply-add operacijas, prasideda
y0≈1/S (typically using a table lookup)
x0=Sy0
h0=12y0
ir iteruoti
rn=0.5−xnhn [ r0=0.5−Sy0⋅12y0=0.5(1−Sy02), r1=0.5−x1h1=0.5−(32x0−0.5x0⋅Sy02)⋅12(32y0−0.5⋅Sy03) ]
xn+1=xn+xnrn [ x1=x0+x0r0=x0+x0⋅0.5(1−Sy02)=32x0−0.5x0⋅Sy02 ]
hn+1=hn+hnrn [ h1=h0+h0r0=12y0+12y0⋅0.5(1−Sy02)=12(y0+y0⋅0.5(1−Sy02))=12(32y0−0.5⋅Sy03) ]
kol ri yra pakankamai arti prie 0, arba daryti fiksuotą skaičių iteracijų. Tai konverguoja prie
limn→∞xn=S,
ir
limn→∞2hn=1/S.


  • Pirmo būdo pavyzdys. Apskaičiuosime S=3.01.
Sprendimas.
b0=S=3.01,
Y0=13≈1/S (parinkome apytikslią S reikšmę lygia 3),
y0=Y0=13= 0.57735026918962576450914878050196,
x0=Sy0=3.01⋅13= 3.01*0.57735026918962576450914878050196 = 1.7378243102607735511725378293109.
Darome pirmą iteraciją:
b1=b0Y02=3.01⋅(13)2=3.01⋅13= 3.01/3 = 1.0033333333333333333333333333333,
Y1=12(3−b1)= 0.5*(3 - 1.0033333333333333333333333333333) = 0.99833333333333333333333333333335,
x1=x0Y1= 1.7378243102607735511725378293109*0.99833333333333333333333333333335 = 1.7349279364103389285872502662621,
y1=y0Y1= 0.57735026918962576450914878050196*0.99833333333333333333333333333335 = 0.5763880187409763882349668658678.
Darome antrą iteraciją:
b2=b1Y12= = 1.0033333333333333333333333333333*0.99833333333333333333333333333335^2 = 0.99999167592592592592592592592593,
Y2=12(3−b2)= 0.5*(3 - 0.99999167592592592592592592592593) = 1.000004162037037037037037037037,
x2=x1Y2= 1.7349279364103389285872502662621*1.000004162037037037037037037037 = 1.7349351572446668586552910067157,
y2=y1Y2= 0.5763880187409763882349668658678*1.000004162037037037037037037037 = 0.57639041768925809257650864010487.
Atėmę x2 iš S=3.01, gauname paklaidą:
1.7349351572897472412324994277 - 1.7349351572446668586552910067157 = 4.5080382577208420984281695490435e-11 ≈4.508038⋅10−11.
Atėmę y2 iš 13.01, gauname paklaidą:
0.57639041770423496386461775006644 - 0.57639041768925809257650864010487 = 1.4976871288109109961569101491839e-11 ≈1.497687⋅10−11.
Darome trečią iteraciją:
b3=b2Y22= = 0.99999167592592592592592592592593*1.000004162037037037037037037037^2 = 0.99999999994803219891278740696279,
Y3=12(3−b3)= 0.5*(3 - 0.99999999994803219891278740696279) = 1.0000000000259839005436062965186,
x3=x2Y3= 1.7349351572446668586552910067157*1.0000000000259839005436062965186 = 1.734935157289747241230742381434,
y3=y2Y3= 0.57639041768925809257650864010487*1.0000000000259839005436062965186 = 0.57639041770423496386403401376544.
Atėmę x3 iš S=3.01, gauname paklaidą:
1.7349351572897472412324994277 - 1.734935157289747241230742381434 = 0.000000000000000000001757046266 ≈1.757046266⋅10−21.
Atėmę y3 iš 13.01, gauname paklaidą:
0.57639041770423496386461775006644 - 0.57639041770423496386403401376544 = 0.000000000000000000000583736301 ≈5.83736301⋅10−22.


  • Antro būdo pavyzdys. Apskaičiuosime S=3.01.
Sprendimas.
y0=13≈1S (parinkome apytikslią S reikšmę lygia 3),
x0=Sy0=3.01⋅13= 3.01*0.57735026918962576450914878050196 = 1.7378243102607735511725378293109,
h0=12y0=1213=123= 0.28867513459481288225457439025098.
Darome pirmą iteraciją:
r0=0.5−x0h0=0.5−3.01⋅13⋅123=0.5−3.01⋅16= 0.5 - 3.01/6 = -0.00166666666666666666666666666667,
x1=x0+x0r0 = 1.7378243102607735511725378293109 + 1.7378243102607735511725378293109*(-0.00166666666666666666666666666667) = 1.734927936410338928587250266262,
h1=h0+h0r0= 0.28867513459481288225457439025098 + 0.28867513459481288225457439025098*(-0.00166666666666666666666666666667) = 0.28819400937048819411748343293389.
Darome antrą iteraciją:
r1=0.5−x1h1= 0.5 - 1.734927936410338928587250266262*0.28819400937048819411748343293389 = 4.1620370370370370370370370552249e-6,
x2=x1+x1r1 = 1.734927936410338928587250266262 + 1.734927936410338928587250266262*4.1620370370370370370370370552249e-6 = 1.7349351572446668586552910067157,
h2=h1+h1r1= 0.28819400937048819411748343293389 + 0.28819400937048819411748343293389*4.1620370370370370370370370552249e-6 = 0.28819520884462904628825432005244.
Atėmę x2 iš S=3.01, gauname paklaidą:
1.7349351572897472412324994277 - 1.7349351572446668586552910067157 = 0.0000000000450803825772084209843 ≈4.508038⋅10−11.
Atėmę h2 iš 12S=123.01, gauname paklaidą:
0.28819520885211748193230887503322 - 0.28819520884462904628825432005244 = 0.00000000000748843564405455498078 ≈7.488435644⋅10−12.
Darome trečią iteraciją:
r2=0.5−x2h2= 0.5 - 1.7349351572446668586552910067157*0.28819520884462904628825432005244 = 2.5983900543606296518569637165611e-11,
x3=x2+x2r2 = 1.7349351572446668586552910067157 + 1.7349351572446668586552910067157*2.5983900543606296518569637165611e-11 = 1.7349351572897472412307423814339,
h3=h2+h2r2= 0.28819520884462904628825432005244 + 0.28819520884462904628825432005244*2.5983900543606296518569637165611e-11 = 0.28819520885211748193201700688271.
Atėmę x3 iš S=3.01, gauname paklaidą:
1.7349351572897472412324994277 - 1.7349351572897472412307423814339 = 0.0000000000000000000017570462661 ≈1.7570462661⋅10−21.
Atėmę h3 iš 12S=123.01, gauname paklaidą:
0.28819520885211748193230887503322 - 0.28819520885211748193201700688271 = 0.00000000000000000000029186815051 ≈2.9186815051⋅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ų y0, 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ų y0, 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ą y0≈1/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 y0≈1/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 y0. Todėl ir y0 tikslumas negali būti daug didesnis negu skaičiaus S tikslumas iš 3 skaitmenų. Todėl suteikus visiems lentelės skaičiams y0≈1/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 y0≈1/S reikšmę iš 4 skaitmenų. Tada tas 4 skaitmenų y0 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ą y0 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 xi+1≈1n, iteruojant pagal formulę:
xi+1=12xi(3−nxi2).
Skaičiaus n reikšmę mes, žinoma, žinome tiksliai ir norime surasti 1n.
Dabar mes galime pažvelgti į Goldschmidt'o algoritmą, kuris skaičiuoja 1n ir n vienu metu. Aišku būtų galima iš pradžių apskaičiuoti 1n ir gautą rezultatą padauginti iš n ir tokiu budu gauti 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 b0=n. Jeigu mes galėsime surasti skaičius Yi, tokius, kad bn=b0Y02Y12Y22⋯Yn−12 artėja prie 1 (kai skaičiaus b indeksas n didėja), tada yn=Y0Y1Y2⋯Yn−1 artės prie 1b0=1n ir xn=b0Y0Y1Y2⋯Yn−1 artės prie b0=n.
Šitos skaičių Yi serijos yra Newton-Raphson iteracijų žingsniai:
bi=bi−1Yi−12
[ b1=b0Y02, b2=b1Y12=b0Y02Y12, b3=b2Y22=b0Y02Y12Y22, ... ]
Yi=12(3−bi)
[ Y0 yra imamas iš lentelės ir lygus apytiksliai 1n reikšmei.
Y1=12(3−b1)=12(3−b0Y02)=12(3−nY02) ]
Ir tada mes galime skaičiuoti šaknį iš n taip:
x0=nY0
xi=xi−1Yi
[ x1=x0Y1=nY0Y1, x2=x1Y2=nY0Y1Y2, x3=x2Y3=nY0Y1Y2Y3, ... (nepamirškim, kad b0=n) ]
O šaknis iš skaičiaus 1/n skaičiuojama taip (yi≈1n):
y0=Y0
yi=yi−1Yi
[ y1=y0Y1=Y0Y1, y2=y1Y2=Y0Y1Y2, y3=y2Y3=Y0Y1Y2Y3, ...
Be to, iš aukščiau Y1=12(3−nY02), todėl
y1=y0Y1=Y0Y1=Y012(3−nY02),
kas lygu Newton-Raphson iteracijai xi+1=12xi(3−nxi2), kurioje xi+1≈1n. Tik dabar pas mus raidė x pakeista raide y. ]
[bi=bi−1Yi−12, b1=b0Y02, nes x0=nY0, y0=Y0 ir t. t., vadinasi, bi=xi−1yi−1.]
Kadangi bi artėja prie 1 (kai i artėja prie begalybės), tai iš to, kad bi=b0Y02Y12Y22⋯Yi−12 ir iš to, kad yi=Y0Y1Y2⋯Yi bei xi=b0Y0Y1Y2⋯Yi, gauname, kad bi=xi−1yi−1. Todėl mes iš tikrųjų niekada neturime sekti ir skaičiuoti bi:
Yi=12(3−bi)=1+12(1−bi)=1+12(1−xi−1yi−1).
Ir dabar mums taipogi nereikia sekti ir skaičiuoti Yi:
xi=xi−1Yi=xi−1(1+12(1−xi−1yi−1))=xi−1+xi−112(1−xi−1yi−1),
yi=yi−1Yi=yi−1(1+12(1−xi−1yi−1))=yi−1+yi−112(1−xi−1yi−1).
Čia yra kai kurie sudetingi užrašymai, kuriuos galima palengvinti. Skaičiaus 1n apytikslia reikšme iš lentelės pažymėsime raide Y. Tada
x0=nY,
y0=Y.
Tada toliau iteruojame:
ri=12(1−xiyi),
xi+1=xi+xiri,
yi+1=yi+yiri.
Nors dalyjimas iš 2 yra pigus, bet mes to galime išvengti užrašę hi=12yi vietoje yi. Tai ir yra Goldschmidt'o algoritmas.
Sakykime, Y yra apytiksli reikšmė skaičiaus 1n. Pažymime:
x0=nY,
h0=Y2.
Tada iteruoti:
ri=12−xihi,
xi+1=xi+xiri,
hi+1=hi+hiri.
Tada xi konverguoja į n, o hi konverguoja į 12n.
[ Iš tikro,
r0=12−x0h0=12−x012y0=12(1−x0y0),
x1=x0+x0r0,
h1=h0+h0r0=12y0+12y0r0=12(y0+y0r0)=12y1.
Toliau:
r1=12−x1h1=12−x112y1=12(1−x1y1),
x2=x1+x1r1,
h2=h1+h1r1=12y1+12y1r1=12(y1+y1r1)=12y2.
Dar toliau:
r2=12−x2h2=12−x212y2=12(1−x2y2),
x3=x2+x2r2,
h3=h2+h2r2=12y2+12y2r2=12(y2+y2r2)=12y3.
Tuo baigtas įrodymas, kad nereikia kai kurioje vietoje dauginti iš 0.5 (skaičiuojant ri) ir bus skaičiuojama viskas teisingai. ]


Taigi, kuom Goldschmidt'o algoritmas pranašesnis už Newton-Raphson metodą? Newton-Raphson metodas randa yi+1≈1n, iteruojant pagal formulę:
yi+1=12yi(3−nyi2).
Dar jeigu reikia žinoti n, tai reikia yi+1 padauginti iš n.
Goldschmidt'o gi algoritmas su naujoviška procesoriaus architektūra gali skaičiuoti xi ir hi vienu metu lygiagrečiai, pirma apskaičiavus ri−1. Tiesa, pabaigoje, norint žinoti 1n reikšmę, reikią hi 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 n reikšmę ir 1n reikšmę (po padauginimo iš 2 ar šiftinimo į kairę). O taikant Newton-Raphson metodą, jei padauginsime yi+1 iš n, tai turėsime n, bet neturėsime 1n (jeigu kokiuose nors skaičiavimuose reikia jų abiejų). Nebent nukopijuoti iš pat pradžių gautą yi+1≈1n 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 x0=nY skaičiavimui, o kitoje lentelėje reikia turėti skačius Y/2, kurie naudojami h0=Y2 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 x0 ir h0, iš iteracijų
ri=12−xihi,
xi+1=xi+xiri,
hi+1=hi+hiri;
matome, kad ri skaičiavimui reikia vienos daugybos ir vienos atimties operacijos. xi+1 skaičiavimui reikia vienos daugybos ir vienos sudeties operacijos. hi+1 skaičiavimui reikia vienos daugybos ir vienos sudeties operacijos. Bet xi+1 ir hi+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
yi+1=12yi(3−nyi2);
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
yi+1=12yi(3−nyi2)=yi(32−12nyi2)=yi(32−myi2).


  • Goldschmidt'o algoritmo pavyzdys. Apskaičiuosime n=6.309 ir 1n=16.309.
Sprendimas.
y0=16.3≈1n (parinkome apytikslią n reikšmę lygia 6.3; be to, 6.309-6.3=0.009<0.01),
x0=ny0=6.309⋅16.3= 6.309*0.39840953644479787998960572656437 = 2.5135657654302298248544225288946,
h0=12y0=12⋅16.3= 0.5*0.39840953644479787998960572656437 = 0.19920476822239893999480286328219.
Darome pirmą iteraciją:
r0=0.5−x0h0= 0.5 - 2.5135657654302298248544225288946*0.19920476822239893999480286328219 = -7.1428571428571428571428571428367e-4,
x1=x0+x0r0 = 2.5135657654302298248544225288946 + 2.5135657654302298248544225288946*(-7.1428571428571428571428571428367e-4) = 2.5117703613120653749795265128025,
h1=h0+h0r0= 0.19920476822239893999480286328219 + 0.19920476822239893999480286328219*(-7.1428571428571428571428571428367e-4) = 0.19906247910224008360909228980842.
[ 6.309−x1= 2.5117722826721374584282937869002 - 2.5117703613120653749795265128025 = 0.0000019213600720834487672740977. ]
[ 126.309−h1= 0.19906263137360417327851432769854 - 0.19906247910224008360909228980842 = 0.00000015227136408966942203789012 ≈1.52271364⋅10−7. ]
Darome antrą iteraciją:
r1=0.5−x1h1= 0.5 - 2.5117703613120653749795265128025*0.19906247910224008360909228980842 = 7.6494169096209912536443148786058e-7,
x2=x1+x1r1 = 2.5117703613120653749795265128025 + 2.5117703613120653749795265128025*7.6494169096209912536443148786058e-7 = 2.5117722826699328655134995747654,
h2=h1+h1r1= 0.19906247910224008360909228980842 + 0.19906247910224008360909228980842*7.6494169096209912536443148786058e-7 = 0.19906263137342945518414166862938.
[6.309−x2= 2.5117722826721374584282937869002 - 2.5117722826699328655134995747654 = 0.0000000000022045929147942121348 ≈2.204592914794⋅10−12.]
[126.309−h2= 0.19906263137360417327851432769854 - 0.19906263137342945518414166862938 = 0.00000000000017471809437265906916 ≈1.7471809437⋅10−13.]
Darome trečią iteraciją:
r2=0.5−x2h2= 0.5 - 2.5117722826699328655134995747654*0.19906263137342945518414166862938 = 8.7770413345269442665565631474709e-13,
x3=x2+x2r2 = 2.5117722826699328655134995747654 + 2.5117722826699328655134995747654*8.7770413345269442665565631474709e-13 = 2.5117722826721374584282908844297,
h3=h2+h2r2= 0.19906263137342945518414166862938 + 0.19906263137342945518414166862938*8.7770413345269442665565631474709e-13 = 0.19906263137360417327851409767236.
[6.309−x3= 2.5117722826721374584282937869002 - 2.5117722826721374584282908844297 = 0.0000000000000000000000029024705 =2.9024705⋅10−24.]
[126.309−h3= 0.19906263137360417327851432769854 - 0.19906263137360417327851409767236 = 0.00000000000000000000000023002618 =2.3002618⋅10−25.]
[ x32= 2.5117722826721374584282908844297^2 = 6.30899999999999999999998541931 .]


  • Newton-Raphson algoritmo pavyzdys. Apskaičiuosime n=6.309 ir 1n=16.309.
Sprendimas.
y0=16.3≈1n (parinkome apytikslią n reikšmę lygia 6.3; be to, 6.309-6.3=0.009<0.01).
y0=16.3= 0.39840953644479787998960572656437.
m = 0.5*n = 0.5*6.309 = 3.1545.
yi+1=yi(32−myi2).
Darome pirmą iteraciją:
y1=y0(32−my02)= 0.39840953644479787998960572656437*(1.5 - 3.1545*0.39840953644479787998960572656437^2) = 0.39812495820448016721818457961683.
[16.309−y1= 0.39812526274720834655702865539709 - 0.39812495820448016721818457961683 = 0.00000030454272817933884407578026 ≈3.04542728⋅10−7.]
Darome antrą iteraciją:
y2=y1(32−my12)= 0.39812495820448016721818457961683*(1.5 - 3.1545*0.39812495820448016721818457961683^2) = 0.39812526274685891036828333725875.
[16.309−y2= 0.39812526274720834655702865539709 - 0.39812526274685891036828333725875 = 0.00000000000034943618874531813834 ≈3.494361887⋅10−13.]
Darome trečią iteraciją:
y3=y2(32−my22)= 0.39812526274685891036828333725875*(1.5 - 3.1545*0.39812526274685891036828333725875^2) = 0.39812526274720834655702819534471.
[16.309−y3= 0.39812526274720834655702865539709 - 0.39812526274720834655702819534471 = 0.00000000000000000000000046005238 =4.6005238⋅10−25.]
n≈x3=n⋅y3= 6.309*0.39812526274720834655702819534471 = 2.5117722826721374584282908844298.
x32= 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 x3 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 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⋅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ę y0 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 y0 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 1n, 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 y0 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ė (≈1n) esanti ROM'e, priklausomai nuo skaičiaus n.
Pirmas pavyzdys. Aukščiau gavome tokį skaičių:
9.8391347153059527386176341678947e-18 = 9.839134715305952⋅10−18=n.
Imame pirmus 3 jo skaitmenys: 9.83. Apskaičiuojame e=10−18=10−9. Tada 1e=110−9=109. Lentelėje, kuri yra ROM'e yra jau apskaičiuota ≈19.83 reikšmė, kuri lygi 19.83≈0.3189=3.189⋅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 y0≈1n), 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:
y0=3.189⋅10−1⋅109=3.189⋅108.
Ši y0 reikšmė ir taikoma pirmoje iteracijoje Newton-Raphson algoritme.
Antras pavyzdys. Pagal skaičių n=1.4541864575691000⋅1017, rasime y0 reikšmę.
1017=10⋅1016=3.162277660168379⋅108.
Tada
11017=13.162277660168379⋅10−8=0.3162277660168379⋅10−8=3.162277660168379⋅10−9.
ROM'o lentelėje yra apytikslė 11.45 reikšmė (be eksponentės 10−1).
11.45≈0.8304=8.304⋅10−1.
Tada gauname apytikslę y0 reikšmę:
y0=8.304⋅10−1⋅3.162277660168379⋅10−9=26.25955369003821⋅10−10=2.625955369003821⋅10−9.
y0≈11.4541864575691000⋅1017.
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 10999). 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 10976, tai po šaknies traukimo jos reikšmė bus 10976/2=10488. 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
yi+1=yi(32−myi2) (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ą y0 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ę y0 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,
2128= 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.
1.41796875= 1.1907849302036031393001193597853.
6.9388939039072283776476979255676⋅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 2. Todėl galutiniuose skaičiavimose reikės dar padauginti viską iš 2−12=12= 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:
1.41796875⋅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ė 1n reikšmė (be eksponentės) iš ROM, kuri ir yra y0 dvejetainėje sistemoje. Žinant y0 galima, taikant Newton-Raphson algoritmą, apskaičiuoti tikslią 1n reikšmę dvejetainėje sistemoje.
Kadangi ROM atmintyje gali būti 1024-ios y0 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)sign×2−126×0.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×20
Č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=(3F800000)16
Pavyzdys 3. Reikia paversti 0.375 į dvejetainę sistemą.
0.375=(0.011)2=(1.1)2×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=(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≥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ę 21. Jeigu eksponentė yra 130, tai atėmę 127, gausime eksponentės reikšmę 23 (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š 20=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 ≥, 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ą ≥ (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ą ≥.
(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

Gauta iš „https://lt.wikibooks.beta.math.wmflabs.org/w/index.php?title=Sekos_riba&oldid=99“