Aptarimas:Matematika/Teiloro eilutė

Teiloro eilutės Benchmark'as

Toks Free Pascal (Compiler Version 2.6.0) kodas:

var a:longint; c:real;
begin
for a:=0 to 1000000000 do
c:=c+a*(1+
sqr(a*1.0)*(-1/6+
sqr(a*1.0)*(1/120+
sqr(a*1.0)*(-1/5040+
sqr(a*1.0)*(1/362880+
sqr(a*1.0)*(-1/39916800+
sqr(a*1.0)*(1/6227020800+
sqr(a*1.0)*(-1/1307674368000+
sqr(a*1.0)*(1/355687428096000+
sqr(a*1.0)*(-1/121645100408832000+
sqr(a*1.0)*(1/51090942171709440000+
sqr(a*1.0)*(-1/25852016738884976640000+
sqr(a*1.0)*(1/15511210043330985984000000+
sqr(a*1.0)*(-1/10888869450418352160768000000+
sqr(a*1.0)/8841761993739701954543616000000))))))))))))));
writeln(c);
Readln;
End.

duoda rezultatą $3.76998768536565\cdot 10^{237}$ po 52 sekundžių ant 2.6 GHz CPU.

Sinuso Teiloro eilutės teisingumo skaičiavimo lygis

Toks Free Pascal kodas:

var a:longint; c:real;
begin
//for a:=0 to 1000000000 do
a:=6;
c:=c+a*(1+
sqr(a*1.0)*(-1/6+
sqr(a*1.0)*(1/120+
sqr(a*1.0)*(-1/5040+
sqr(a*1.0)*(1/362880+
sqr(a*1.0)*(-1/39916800+
sqr(a*1.0)*(1/6227020800+
sqr(a*1.0)*(-1/1307674368000+
sqr(a*1.0)*(1/355687428096000+
sqr(a*1.0)*(-1/121645100408832000+
sqr(a*1.0)*(1/51090942171709440000+
sqr(a*1.0)*(-1/25852016738884976640000+
sqr(a*1.0)*(1/15511210043330985984000000+
sqr(a*1.0)*(-1/10888869450418352160768000000+
sqr(a*1.0)/8841761993739701954543616000000))))))))))))));
writeln(c);
writeln(sin(6));
Readln;
End.

parodo 2 rezultatus (antras rezultatas absoliučiai teisingas):

-0.279415498042951;

-0.27941549819892587.

Dousiu truputi paaiškinimo kodėl antras rezultatas užrašytas didesniu tikslumu. Intel (na ir AMD) procesoriaus x87 FPU (Floating point unit - x87 koprocecosius) skaičiuoja ne double precision (64 bits), o double extended precision (80 bits, 64 bits without mantissa). Tačiau į atminti lengviau krauti 64 bitus, nei 80 bitų, be to rezultatas po daug algoritmo žingsnių ar iteracijų vis tiek prarandamas (dėl divergencijos ar konvergencijos reiškinių), todėl yra numatytas standartas, kad visos programavimo kalbos palaikytų double precision (64 bits) ir single preicision (32 bits), bet nebūtinai double extended precision (80 bits). Greičiau, net reikalaujama, kad galutinis rezultatas būtų pateiktas arba double precision (64 bits) ir single preicision. Tačiau patys skaičiavimai vykstą su double extended precision (80 bits) ir paskui paverčiamas rezultatas sinuso į double precision (64 bits). Originaliai rezultatai pateikti taip:

-2.79415498042951E-001;

-2.7941549819892587E-0001.

Double precision formatas turi bitus nuo 0 bito iki 51 bito, taigi, iš viso 52 bitus užrašyti skaičiui (nuo 52 iki 62 tenka mantisai ir 63 bitas minuso ženklui). Todėl gauname, kad Double precision yra $2^{52}=4503599627370496.$ Mes negalime naudoti 4 be 5, 6, 7, 8, 9. Todėl iš 16 gautų skaitmenų tinka tik 15, t. y. skaičius 999999999999999. Todėl pirmas rezultatas yra double precision ir turi 15 skaitmenų. Kadangi viskas skaičiuojama iš tikro su double extended precision (80 bits), tai suapvalinimui surasti 16 skaitmenį nėra problemų.

Jeigu antras rezultatas (-2.7941549819892587E-0001) pateiktas su double extended precision, tada kyla klausimas kodėl yra tik 17 skaitmenų, jeigu žinoma, kad su double extended precision bitai nuo 0 iki 62 yra skaitmenims, 63 bitas yra sveikajam skaičiui ir bitai nuo 64 iki 78 yra mantisai bei 79 bitas skirtas minuso ženklui. Galimi 2 paaiškinimai tokie:

1) $2^{62}=4611686018427387904$ yra 19 skaitmenų, mums tinka tik 18 skaitmenų, vienas skaitmuo nerodomas, nes naudojamas apvalinimui, taigi, gaunasi 17 rodomų skaitmenų (apvalinimui reikia daugiau nei vieno bito, nes apvalinama iš dvejetainės sistemos į dešimtainę);

2) $2^{63}=9223372036854775808$ yra 19 skaitmenų, mums tinka tik 18 skaitmenų (nes reikia, kad būtų 19 devynetų, o ne 9223372036854775808), vienas skaitmuo nerodomas, nes naudojamas apvalinimui, taigi, gaunasi 17 rodomų skaitmenų.

Neteisingo apvalinimo po šaknies traukimo pavyzdis

Va štai toks Free Pascal (Version 1.0.12 2011/12/25; Comiler Version 2.6.0; Debugger GDB 7.2) kodas:

 var
 a:longint;
 c:real;
 begin
 for a:=1 to 1000000000  do
 c:=c+6/(sqr(a*1.0));
 writeln(c);
 writeln(sqrt(c));
 readln;
 end.

duoda 2 rezultatus:

9.86960433860995+000;

3.1415926436458864+0000.

Šie rezultatai gaunami po 11 sekundžių su 2.6 GHz CPU ir DDR-800 (400 MHz) RAM (procesorius yra dual-core, bet vis tiek Free Pascal naudoja tik vieną core; procesorius turi 512 KB L2 per core, bet x87 FPU turi 8 80 bitų (double-extended precision) registrus, todėl galimas daiktas, kad L2 cashe, kaip ir L1 cashe net neegzistuoja; procesorius taip pat turi 8 16/32/64 bitų General Purpose registrus ir 16 128 bitų SSE registrų, taigi, smulkioms programoms ir net GPU (kurie turi iš tikro šiuo metu apie 8-32 core'ius, o ne 640-2048 shader-core'ius) blur'o effektams užtenka 8 arba 16 registrų (iš kvadrato 4*4 pikselių judama į dešinį šona trinant 4 pikseliu iš kairės ir vietoje jų įdėjus iš RAM 4 naujus pikselius iš dešinės ir kiekvieną kart sudedant 16 pixelių), net jei RAM veikia iki 8 ar 16 kartų lėčiau nei CPU).

Tikroji $\pi$ reikšmė yra 3,1415926535897932384626433832795. O įrašius į Microsoft Windows kalkuliatorių reikšmę 9,86960433860995 ir ištraukius su Microsoft Windows kalkuliatoriumi šaknį gaunamas rezultas 3,1415926436458864995793203208207. Palyginus matome, kad suapvalinta blogai:

3.1415926436458864+0000;

3,14159264364588649957932.

Gal čia Free Pascal (kuris greičiausiai naudojasi CPU koprocesoriaus x87 FPU šaknies traukimo instrukcija FSQRT, 151 puslapis arba 5.2.2; http://www.intel.com/content/www/us/en/processors/architectures-software-developer-manuals.html) apvalina nevisada teisingai tik traukiant šaknį, o ne, pavyzdžiui, dauginant.

Toks Free Pascal (Version 1.0.12 2011/12/25; Comiler Version 2.6.0; Debugger GDB 7.2) kodas:

 var
 a:longint;
 c:real;
 begin
 for a:=3 to 1000000000  do
 c:=c+6/(sqr(a*1.0-2));
 writeln(c);
 writeln(sqrt(c));
 readln;
 end.

duoda rezultatus

9.86960433860995+000 ir

3.1415926436458864+0000 po 11 sekundžių su 2.6 GHz CPU ir DDR-800 (400 MHz) RAM. Atimties instrukcija neprailgina skaičiavimo laiko ir tai greičiausiai ne dėl savaiminio Free Pascal kodo optimizavimosi, o dėl procesoriaus architektūros ypatumų. Profesionaliai kalbant toks ypatumas būtų pavadintas CISC architekturos pranašumu palyginus su RISC architektūra, taučiau tiesa vis tiek turbūt yra tokia, kad visi procesoriai yra RISC tik, kadangi Intel ir AMD procesoriai yra populiariausi, tai daugiausiai visokių palengvinančių instrukcijų prikūrta, kaip MMX, SSE, AVX, kurios greičiausiai tik sutrumpina programos rašymą, bet nepagreitina niekada patį programos veikimo procesą. Teoriškai MMX, SSE arba AVX instrukcijos turėtų pagreitintį, pavyzdžiui, grafikoje vektorių t. y. geometrijos ir pikselių skaičiavimą 3-4 kartus, bet kadangi spalvų yra 3, o ne keturios ir trimatė erdvė yra iš 3 koordinačių, tai kai kurios užduotys iš tiesų galbūt gali būti padarytos truputi greičiau be visokių SSE, nes reikia 3 operacijų vietoje 4. Cisc architektūra teoriškai sukuria daug ilgų pipline'ų, kada per viena operaciją ir atimama ir sudauginama. Tačiau daug didesnė tikimybė, kad loop'inant kažkaip savaime atėmimo operacija gaunasi nemokamai ir nėra ten jokių piplainų; visi CPU veikia tuo pačiu principu ir RISC ir CISC ir GPU. Taip pat labai galimas variantas, kad konstantos, kuri yra CPU registre, atėmimui reikia nuo 1 iki 4 ciklų todėl atėmimo operacija užima nuo 1/2.6=0.38 sekundės iki 4/2.6=1.54 sekundės. Apytiksliai vienai daugybos arba sudeties operacijai reikia 3,5 ciklo, tačiau sekunės skirtumas nejaučiamas, gal dėl to, kad sudėtis (ir atimtis) daroma truputi greičiau nei daugyba (be abejonės yra sudeties, daugybos ir dalybos pipline'ai 80 bitų (arba 64 bitų sveikiems skaičiams), o jų kombinacijų greičiasiai nėra).

Kaip [in software] pagreitinti sinuso Teiloro eilutę

Kad apskaičiuoti sinusą yra skitos Intel CPU x87 FPU sinuso instrukcijos (galima iš viso skaičiuoti sinusą arba kosinusą arba tangenta arba dar su specifinę instrukcija ir sinusą ir kosinusą vienu metu, kas greičiau nei atskirai, bet lėčiau nei kiekvieną atskirai). Šios x87 FPU (x87 FPU yra nuo Intel 8087, 287, 387; koprocesorius 487 nededamas kartu su CPU tik SX versijai, o toliau nuo intel 486 CPU visi procesoriai turi integruotus x87 FPU į tą patį čipą) koprocesoriaus instrukcijos cos, sin, tan, sincos yra in hardware, t. y., sinuso skaičiavimą Free Pascal kalba naudoja iš x87 FPU intrukcijų per asamblerį, o ne rašo asambleryje Teiloro eilutes. Bet galbūt sinusą in software galima apskaičiuoti greičiau nei in hardware, ypač jei skaičiai yra nuo 0 iki $2\pi ,$ nes tada nereikia skaičiuoti frac() (arba ceil()) funkcijos, kuriai reikia apie 65 ciklų, tuo tarpu, kai sinusui in hardware reikia 122 ciklų.

Kad pagreitinti sinuso skaičiavimą reikia skaičiuoti ne nuo 0 iki $2\pi ,$ o nuo $-\pi$ iki $-\pi$ (šitas triukas netinka kosinusui, nes kosinusas kitaip simetriškas koordinačių ašims nei sinusas). Taigi, normalai bet kokį didelį skaičių a nuo $-\infty$ iki $\infty$ mes paverstume į skaičių 0 iki $2\pi$ taip:

$b=6.283185307179586477*frac(a/6.283185307179586477).$

Skaičių b mes jau galime dėti į Teiloro eilutę ir gausime tą patį rezultatą, kokį gautume įdėjus į begalo ilgą Teiloro eilutę skaičių a. Skaičiuojant su skaičiu b vietoje a Teiloro eilutė labai sutrumpėja, kad gauti tą patį tikslumą.

Greitesnis būdas apskaičiuoti yra toks:

$b_{2}=3.141592653589793238-6.283185307179586477*frac(a/6.283185307179586477).$

Įdėjus $b_{2}$ į sinuso Teiloro eilutę atsakymo tikslumas turėtų padidėti dvigubai palyginus su tuo, ką gautume įdėjus b (tarkim vietoje 5 teisingų skaimenų po kablelio, gautume 10 teisingų skaimenų po kablelio).

Toks Free Pascal kodas:

var a:longint; b,c:real;
begin
//for a:=0 to 1000000000 do
a:=6;
b:=3.141592653589793238-6.283185307179586477*frac(a/6.283185307179586477);
c:=c+b*(1+
sqr(b*1.0)*(-1/6+
sqr(b*1.0)*(1/120+
sqr(b*1.0)*(-1/5040+
sqr(b*1.0)*(1/362880+
sqr(b*1.0)*(-1/39916800+
sqr(b*1.0)*(1/6227020800+
sqr(b*1.0)*(-1/1307674368000+
sqr(b*1.0)*(1/355687428096000+
sqr(b*1.0)*(-1/121645100408832000+
sqr(b*1.0)*(1/51090942171709440000+
sqr(b*1.0)*(-1/25852016738884976640000+
sqr(b*1.0)*(1/15511210043330985984000000+
sqr(b*1.0)*(-1/10888869450418352160768000000+
sqr(b*1.0)/8841761993739701954543616000000))))))))))))));
writeln(c);
writeln(sin(6));
Readln;
End.

duoda rezultatus:

-2.79415498198926E-001;

-2.7941549819892587E-0001.

Tikslumas išaugo dvigubai (o gal ir dar daugiau?), kaip ir buvo tikėtasi.

Update 1 [2024 04 05]. Su naujasne Free Pascal versija (Free Pascal IDE Version 1.0.1 [2020/06/04]; Compiler Version 3.2.0) paskutinis kodas duoda tokius rezultatus:

-2.7941549819892575E-001

-2.79415498198925872810E-0001

Tiksli sin(6) reikšmė iš kalkuliatoriaus yra tokia (čia yra 6 radianai):

-0.27941549819892587281155544661189.

Kiek suprantu, kodo pirmas 17 skaitmenų rezultatas yra Double precision (64 bit), o antras 21 skaitmens rezultatas yra Double Extended precision (80 bit). Double precision (64 bit) garantuoja 15 pirmų teisingų skaitmenų iš 17 (likę du skaitmenys kartais gali būti teisingi dėl konvertavimo iš dvejetainės į dešimtainę sistemą). O Double Extended precision (80 bit) berods garantuoja 18 pirmų teisingų skaitmenų. Matyt, kad kai kuriais atvejais gali būti teisingas ir 21 skaitmuo Double Extended precision (80 bit) formate.

Toks Free Pascal kodas (kuris skiriasi tik tuo, kad vietoje 6 yra 38):

var a:longint; b,c:real;
begin
a:=38;
b:=3.141592653589793238-6.283185307179586477*frac(a/6.283185307179586477);
c:=c+b*(1+
sqr(b*1.0)*(-1/6+
sqr(b*1.0)*(1/120+
sqr(b*1.0)*(-1/5040+
sqr(b*1.0)*(1/362880+
sqr(b*1.0)*(-1/39916800+
sqr(b*1.0)*(1/6227020800+
sqr(b*1.0)*(-1/1307674368000+
sqr(b*1.0)*(1/355687428096000+
sqr(b*1.0)*(-1/121645100408832000+
sqr(b*1.0)*(1/51090942171709440000+
sqr(b*1.0)*(-1/25852016738884976640000+
sqr(b*1.0)*(1/15511210043330985984000000+
sqr(b*1.0)*(-1/10888869450418352160768000000+
sqr(b*1.0)/8841761993739701954543616000000))))))))))))));
writeln(c);
writeln(sin(38));
Readln;
End. 

duoda rezultatus:

2.9636857870938516E-001

2.96368578709385317445E-0001

Vėl gavome 17 skaitemnų ir 21 skaitmens rezultatus. Kalkuliatoriaus sin(38) reikšmė yra tokia (38 radianai):

0.29636857870938531739229664984902.

Pirmas rezultatas duoda pirmus 15 teisingų skaitmenų, o antras rezultatas duoda 18 pirmų teisingų skaitmenų.

Kode įrašius 136 vietoje 38, gaunami tokie rezultatai:

-7.9043320672288875E-001

-7.90433206722888758095E-0001

Kalkuliatoriaus sin(136) reikšmė yra tokia (136 radianai):

-0.79043320672288875799022474878835.

Pirmas rezultatas duoda visus 17 teisingus skaitmenis. O antras rezultatas duoda 17 pirmų teisingų skaitmenų, bet galima sakyti 18 pirmų teisingų skaitmenų, nes ...22888758 beveik tas pats kas ...2288875799.

Taigi, konvertuojant iš dvejetainės į dešimtainę sistemą gaunami dešimtainėje sistemoje paskutiniai pora/keli skaitmenys kartais teisingi, o kartais neteisingi. Double precision (64 bit) garantuoja 15 teisingų pirmų skaitmenų, o Double Extended precision (80 bit) garantuoja, ko gero, 18 pirmų teisingų skaitmenų.

Kode funkcija frac() duoda skaičių be sveikos dalies, pavyzdžiui, frac(56423.465684681)=0.465684681.

Kode pirmą eilutę galima pakeisti tokia:

var a,b,c:real;

tada skaičius a galės būti ne tik sveikas skaičius (galės būti bet koks realusis skaičius).

Čia pasinaudota redukcijos formule (Redukcijos formulės)

$\sin(\pi -\alpha )=\sin \alpha ,\;\;{\text{kai}}\;\;0\leq \alpha <2\pi .$

Va čia ( https://en.wikipedia.org/wiki/Extended_precision#Working_range ) iš angliškos Vikipedijos apie Double Extended precision (80 bit):

"The 80-bit floating-point format has a range (including subnormals) from approximately 3.65×10⁻⁴⁹⁵¹ to 1.18×10⁴⁹³². Although log₁₀(2⁶⁴) ≅ 19.266, this format is usually described as giving approximately eighteen significant digits of precision (the floor of log₁₀(2⁶³), the minimum guaranteed precision)."

Be to, $\log _{10}2^{63}\approx 18.9648897=$

=log₁₀(9,223,372,036,854,775,808) = 18.964889726830815298465550367643,

t. y. beveik 19 skaitmenų dauguma atveju, bet kartais garantuoti (teisingi) buna tik 18 pirmų skaitmenų. Matyt, taip.

Kodėl yra 21 skaitmuo Double Extended precision (80 bit) formate, galima paaiškinti taip (analogiškai, kaip tai aiškinama Double precision (64 bit) Vikipedijos straipsnyje: https://en.wikipedia.org/wiki/Double-precision_floating-point_format#IEEE_754_double-precision_binary_floating-point_format:_binary64 ):

$2^{-64}\approx 5.42101\cdot 10^{-20}=$ 5.4210108624275221700372640043497e-20.

Šis skaitmuo užims 21 skaitmenį ir atrodys maždaug taip 0.00000000000000000005, nes pavyzdžiui $10^{-3}=0.001,$ t. y. užima 4 skaitmenis. Todėl ir yra 21 skaitmuo ir kai kuriais atvejais visi 21 skaitmenys gali būti teisingi, konvertuojant iš dvejetainės į dešimtainę sistemą.

Beje, $2^{-63}\approx 1.084202\cdot 10^{-19}=$ 1.0842021724855044340074528008699e-19. Tada reikia 20 skaitmenų tokiam skaičiui. Tiksliai negaliu pasakyti kodėl aukščiau skaičiuojama su 63 bitais, o ne 64 bitais ($\log _{10}2^{63}\approx 18.9648897$), nes Vikipedijos Extended precision (80 bit) straipsnyje ( https://en.wikipedia.org/wiki/Extended_precision#x86_extended_precision_format ) paveikslėlyje rodoma, kad yra 63 bitai ir vienas bitas integer (sveikasis skaičius), kuris, kaip manau, stovi prieš kablelį pats pirmas (tai gaunas kaip ir 1+63=64 bitai).

Bet panašu (iš to paties straipsnio, lentelės ir paaiškinimo po lentele), kad bitai 0-62 veikia su visais koprocesoriais įskaitant 8087 ir 80287 ir vėlesnius kaip 80387, o 63-čias bitas veikia tik su 8087 ir 80287, jei jo reikšmė yra 0. Kai 63-čio bito reikšmė yra 1, tai koprocesoriaus skaičiavimai veikia normaliai su visais koprocesoriais nuo 8087 iki bet kokio velesnio, kaip 80387. Iš 63-čio bito naudos gali išgauti tik 8087 ir 80287 koprocesoriai. If bit 63 is zero then "Unnormal. Only generated on the 8087 and 80287. The 80387 and later treat this as an invalid operand.". If bit 63 is one then "Normalized value.".

Iš čia ( https://en.wikipedia.org/wiki/Extended_precision#Working_range ):

"Bounds on conversion between decimal and binary for the 80-bit format can be given as follows: if a decimal string with at most 18 significant digits is correctly rounded to an 80-bit IEEE 754 binary floating-point value (as on input) then converted back to the same number of significant decimal digits (as for output), then the final string will exactly match the original; while, conversely, if an 80-bit IEEE 754 binary floating-point value is correctly converted and (nearest) rounded to a decimal string with at least 21 significant decimal digits then converted back to binary format it will exactly match the original."

Ten rašo, kad jeigu iš 18 dešimtainių skaitmenų paversti į 80 bitų (IEEE 754) dvejetainės sistemos skaičių, o paskui tą dvejetainį skaičių paversti vėl į dešimtainį skaičių, tai gausime tą patį dešimtainį skaičių su 18 tokių pačių skaitmenų. O Jeigu 80 bitų (IEEE 754 - 80-bit floating-point format) [Double Extended precision] formato dvejetainės sistemos skaičių paversti į 21 skaitmens dešimtainį skaičių ir paskui šitą 21 skaitmens dešimtainės sistemos skaičių paversti vėl į dvejetainį skaičių (80-bit floating-point formate), tai tas dvejetainis skaičius bus lygiai toks pats koks buvo pradžioje. Tai gal dėl to ir yra 21 [dešimtainio] skaitmens atsakymas, kai FPU skaičiuoja sin(x) iš FPU instrukcijų.

Skaičiaus e skaičiavimas ir ne tik

Toks Free Pascal kodas (skaičiuojantis e):

var a:longint; c:real;
begin
c:=1;
for a:=1 to 1000000000 do
c:=c*1.000000001;
writeln(c);
Readln;
End.

duoda atsakymą "2.7182818271394744E+000" po 6 sekundžių su 4.16 GHz veikiančių procesorium (pirmą kartą duoda šitą atsakymą po maždaug 85 sekundžių). Čia yra milijardas iteracijų ir milijardas daugybos operacijų. Todėl gaunasi, kad vienai daugybos operacijai reikia 6*4.16=24.96=~25 ciklų. Bet negali būti tiek daug ciklų vienai daugybos operacijai, todėl tikriausiai daug ciklų suvalgo pati iteracijos operacija.

Kalkuliatoriaus reikšmė yra tokia:

$(1.000000001)^{10^{9}}=$ 2.7182818270999043223766440238603. Gavome su FP kodu 10 (beveik 11) pirmų teisingų skaitmenų. Vadinasi, po daug dauginimų skaičiavimo tikslumas mažėja.

Tiksli skaičiaus e reikšmė yra tokia: e=2.7182818284590452353602874713527.

Toks Free Pascal kodas:

var a:longint; c:real;
begin
for a:=1 to 1000000000 do
c:=c+1/a;
writeln(c);
Readln;
End.

duoda atsakymą "2.1300481502506980E+001" po 6 sekundžių su 4.16 GHz veikiančių procesorium (pirmą kartą duoda šitą atsakymą po maždaug 103 sekundžių). Kode yra milijardas iteracijų, 1 milijardas sudeties operacijų ir 1 milijardas dalybos operacijų. Sudeties operacija užima apie pora ciklų, o dalybos operacija užima daug daugiau negu sudeties ar daugybos operacija. Čia mes gavome, kad vienai dalybos operacijai reikia 6*4.16=24.96=~25 ciklų. Gali būti, kad kol planuojama nauja iteracija, FPU (floating point unit) nesnaudžia ir daro dalybos operaciją ar kažkas tokio.

Su 2.6 GHz procesorium buvo skaičiuotas toks pat kodas čia: https://lt.wikibooks.org/wiki/Matematika/Kreiviniai_integralai#Kreivės_masė

ten buvo gautas atsakymas 21,3004815025070 po 8 sekundžių su 2,6 GHz procesoriumi. Išeina, kad tenai 2.6 GHz procesoriui vienai dalybos operacijai reikėjo 8*2.6=20.8=~21 ciklo. Arba viena sekunde ten netiksliai paskaičiuota (tada 9*2.6=23.4 ciklo), arba 2.6 GHz procesoriui reikėjo truputi mažiau ciklų šitame kode dalybai (ir viskam kitam) negu 4.16 GHz dažniu veikiančiam procesoriui.

Toks Free Pascal kodas:

var a:longint; c:real;
begin
for a:=1 to 1000000000 do
c:=c+sqrt(a);
writeln(c);
Readln;
End.

duoda atsakymą "2.1081851083598383E+013" po 7 sekundžių su 4.16 GHz veikiančių procesorium (pirmą kartą [paleidus] duoda šitą atsakymą po maždaug 90 sekundžių). Kode yra milijardas iteracijų, 1 milijardas sudeties operacijų ir 1 milijardas kvadratinės šaknies traukimo operacijų. Oficialiai [kvadratinės] šaknies traukimo operacija užima tik truputi daugiau ciklų nei dalybos operacija. Čia mes gavome, kad vienai šaknies traukimo operacijai reikia 7*4.16=29.12=~29 ciklų.

Kodas skaičiuoja integralą

$\int _{0}^{10^{9}}{\sqrt {x}}\;dx=\int _{0}^{10^{9}}x^{1 \over 2}\;dx={\frac {x^{3/2}}{\frac {3}{2}}}|_{0}^{10^{9}}={\frac {2}{3}}(10^{9})^{3/2}-0={\frac {2}{3}}{\sqrt {10^{27}}}\approx 2.1081851067789\cdot 10^{13}=$

=21,081,851,067,789.195546659290296218.

Su Free Pascal kodo skaičiavimu gavome 9 pirmus teisingus skaitmenis.

Toks Free Pascal kodas:

var a:longint; c:real;
begin
for a:=1 to 1000000000 do
c:=c+a*1.0*a;
writeln(c);
Readln;
End.

duoda atsakymą "3.3333333383324648E+026" po 5 sekundžių su 4.16 GHz veikiančių procesorium (pirmą kartą [paleidus] duoda šitą atsakymą po maždaug 101 sekundės). Kode yra milijardas iteracijų ir 1 milijardą kartų skaičius a keliamas kvadratu. Galima eilutę "c:=c+a*1.0*a;" pakeisti eilute "c:=c+sqr(a*1.0);". Tuomet skaičiavimas užtruks irgi 5 sekundes. Jei a nebus padaugintas iš 1.0 tada bus klaidos... Gaunasi, kad vienai sudečiai ir vienai [naudingai] daugybos operacijai reikia 5*4.16=20.8=~21 ciklo. Žinoma, čia greičiausiai kažkaip iteracijos paima daug ciklų, o pati daugybos operacija kokius 4 ciklus naudoja.

Kodas skaičiuoja integralą

$\int _{0}^{10^{9}}x^{2}\;dx={\frac {x^{3}}{3}}|_{0}^{10^{9}}={\frac {(10^{9})^{3}}{3}}-0={\frac {10^{27}}{3}}\approx 3.3333333333\cdot 10^{26}=$

=333,333,333,333,333,333,333,333,333.33333.

Su FP kodu gavome 9 pirmus teisingus skaitmenis.

Sinuso ir kosinuso greitesnis skaičiavimas

$\sin 85^{\circ }-\cos 5^{\circ }=$ -2.7725921425635885468885319116808e-48 $\approx -2.772592\cdot 10^{-48}.$ Toks atsakymas gaunamas su Windows 10 kalkuliatorium. Toks greičiausiai ir yra didžiausias Windows 10 kalkuliatoriaus tikslumas ($\approx 10^{-48}$).

$\sin 73^{\circ }-\cos 17^{\circ }=$ -9.3009073833400679282048050079536e-48 $\approx -9.30090738334\cdot 10^{-48}.$

Taigi, galima apskaičiuoti, pavyzdžiui, sinusą su Makloreno eilute nuo 0 iki 45 laipsnių, o sin(x) reikšmes nuo 45 laipsnių iki 90 laipsnių (kai 45<x<90) galima gauti taip:

$\sin x={\sqrt {1-\cos ^{2}(90-x)}}.$

O kosinuso reikšmes, kai 45<x<90 (laipsnių), galima gauti, apskaičiavus sinuso Makloreno eilutę su $0<x<{\frac {\pi }{4}}$ ir paskui atlikti tokias operacijas (kelimo kvadratu ir šaknies traukimo operacijas):

$\cos x={\sqrt {1-\sin ^{2}({\frac {\pi }{2}}-x)}}.$

Tik truputi neaišku, kas procesoriaus FPU yra greičiau ar ištraukti kvadratinę šaknį, ar skaičiuot maždaug turbut dvigubai ilgesnę Makloreno eilutę.

Sinuso liekamojo nario formulė yra tokia:

$|r(x)|=|{\frac {x^{2k+3}}{(2k+3)!}}\cos c|<|{\frac {x^{2k+3}}{(2k+3)!}}|.$

Tarkime, kad k=10, o $x={\frac {\pi }{2}}\approx 1.5707963.$ Tada paklaida bus

$|r(x)|<|{\frac {x^{2k+3}}{(2k+3)!}}|=|{\frac {x^{2\cdot 10+3}}{(2\cdot 10+3)!}}|=|{\frac {x^{23}}{23!}}|=$

= 1.5707963^23/25852016738884976640000 = 32415.819192852196702666149390257/25852016738884976640000 =

= 1.2538990485834848423867495105035e-18 $\approx 1.253899\cdot 10^{-18}.$

O dabar tarkime, kad k=5, o $x={\frac {\pi }{4}}\approx 0.785398.$ Tada paklaida bus

$|r(x)|<|{\frac {x^{2k+3}}{(2k+3)!}}|=|{\frac {x^{2\cdot 5+3}}{(2\cdot 5+3)!}}|=|{\frac {x^{13}}{13!}}|=$

= 0.785398^13/6227020800 = 0.04326804604255175660997065487462/6227020800 =

= 6.9484344813095463901406359321326e-12 $\approx 6.94843448\cdot 10^{-12}.$

Matome, kad dvigubai ilgesnė Makloreno eilutė (su $x={\frac {\pi }{2}}\approx 1.5707963$) sinusui skaičiuoti duoda geresnį tikslumą ($\approx 1.253899\cdot 10^{-18}$) negu trumpesnė Makloreno eilutė (su $x={\frac {\pi }{4}}\approx 0.785398$ ir kurios tikslumas $\approx 6.94843448\cdot 10^{-12}$).

Gal maždaug pusantro karto pailgėja Makloreno eilutė, kai $x={\frac {\pi }{4}}$ pakeičiamas $x={\frac {\pi }{2}},$ kad gauti tą patį tikslumą. Sakyčiau, ilgesnę Makloreno eilutę paskaičiuot turėtų būti greičiau procesoriaus FPU, negu skaičiuoti trumpesnę Makloreno eilutę ir traukuti šaknį. Bet jeigu reikia vienu metu apskaičiuoti ir sinusą ir kosinusą, tai tikriausiai greičiau apskaičiuoti sinusą pagal Makloreno eilutę, o kosinusą paskui apskaičiuoti pagal formulę (naudojant šaknies traukimą):

$\cos x={\sqrt {1-\sin ^{2}x}}.$

Taip kompiuteriui bus greičiau nei skaičiuoti ir sinuso, ir kosinuso Makloreno eilutes.

Pradžioje ten klaidos (nereikia traukti jokios šaknies!) ir gaunasi, kad galima labai greitai apskaičiuoti sin(x), kai 45<x<90 laipsnių. Štai taip:

$\sin x=\cos(90^{\circ }-x),\quad 45^{\circ }<x<90^{\circ }.$

Pavyzdžiui, jeigu x=90 laipsnių, tai cos(90-x)=cos(90-90)=cos(0)=1 ir sin(x)=sin(90)=1.

Kosinusas gi, apskaičiuojamas taip:

$\cos x=\sin(90^{\circ }-x),\quad 45^{\circ }<x<90^{\circ }.$