. Aritmetica In Precizie Multipla

Introducere

Acestă lucrare conține următoarele capitole:

Capitolul I – cuprinde prezentarea aritmeticii cu punct mobil, cum ar fi: reprezentarea numerelor, calcule în simplă și dublă precizie, algoritmi folosiți pentru implementarea celor patru operații elementare: adunarea, scăderea, înmulțirea, împărțirea.

Capitolul II – cuprinde calcule în precizie multiplă folosind un pachet de funcții MFLOAT, împreună cu cîteva exemple cum ar calculul funcțiilor Bessel, calculul constantei lui Euler, calculul lui PI.

Capitolul III și IV – cuprind diferite metode de calcul ale funcțiilor transcendente, cum fi: funcția sin, funcția cos, funcția exponențială, funcția arctangentă.

Capitolul V – prezintă descrierea funcțiilor de dază mfloat, precum și informații despre programul mfloat.

Capitolul VI – conține informații despre utilzarea programului precum și listarea unei părți de program.

I. Aritmetica

Scopul principal al acestui capitol este de a efectua un studiu minuțios asupra celor patru procese de bază ale aritmeticii: adunarea, scăderea, înmulțirea, împărțirea. Mulți privesc aritmetica drept un lucru fără însemnătate, pe care îl învață copiii și cu care lucrează calculatoarele, dar aritmetica este un subiect fascinant, avănd multe implicații interesante. Întrucît aritmetica fundamentează atît de multe aplicații ale calculatoarelor, este important să facem un studiu amănunțit asupra metodelor eficiente de calcul cu numere.

Aritmetica este o disciplină vie care a jucat un rol important în istoria omenirii, fiind carecterizată în continuare printr-o rapidă dezvoltare. În acest capitol se analizează algoritmi pentru efectuarea de operații aritmetice asupra mai multor tipuri de cantități, cum sunt numerele cu "punct mobil", numere foarte mari.

I.1. Sisteme poziționale de numerotație

Modul în care studiem aritmetica este strîns legat de modul în care reprezentăm numerele cu care lucrăm. Vom începe studiul asupra acestei discipline cu o discuție privind principalele metode de reprezentare a numerelor.

Notația pozițională folosind baza b este definită prin regula

În cazul sistemului de numerație zecimal b este zece, iar coeficienții a sunt aleși dintre "cifrele zecimale" 0, 1, 2, 3, 4, 5, 6, 7, 8, 9.

Generalizarea sistemului de numerație zecimal se obține cînd b este un întreg mai mare decît 1, iar coeficienții a sunt întregi cuprinși în intervalul 0ak<b. Astfel avem sistem standard de numerație binar (b=2), ternar (b=3), cuaternar (b=4), cinquinar (b=5).

Punctul care apare între a0 și a-1 este denumit punctul bazei, cînd b=10 se numește punct zecimal, iar cînd b=2 poartă denumirea de punct binar etc. Europenii de pe continent folosesc virgula în locul punctului, englezii folosesc în mod frecvent un punct mai ridicat.

Coeficienții a sunt numiți cifre ale reprezentării. Cînd indicele k este mare, cifra ak este „mai semnificativă“ decît cifra ak avînd k mai mic, ne referim la cifra situată la extremitatea din stînga a numărului ca fiind cea mai semnificativă cifră, iar la cifra situată la extremitatea din dreapta ca fiind cea mai puțin semnificativă cifră. În sistemul standard de numerație binar, cifrele binare sînt numite biți, iar în sistemul standard de numerație hexazecimal (cu baza 16) cifrele hexazecimale de la zero la cincisprezece sînt de obicei notate prin simboluri: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E, F.

I.2. Aritmetica numerelor reprezentate cu punct mobil

Secțiunea este destinată studiului principiilor de bază ale aritmeticii numerelor cu "punct mobil", prin analiza mecanismelor interne care stau la baza unor astfel de calcule. Calculatoarele de astăzi au instrucțiuni cablate pentru calcule în sistemul de reprezentare cu punct mobil, fie că posedă subrutinele corespunzătoare.

I.2.1. Calcule cu precizie simplă

Notația cu punct mobil. În aceasta secțiune poziția punctului bazei este în mod

dinamic fie variabila fie "mobilă", iar pe masură ce programul se execută se atașează fiecărui număr o indicație privind poziția curentă a punctului bazei sale. Este necesar în special pentru exprimarea numerelor foarte mari sau a numerelor foarte mici.

În aceasta secțiune vom lucra cu numere reprezentate cu punct mobil în baza b, cu exces de q, pe lungime de p cifre. Numerele vor fi reprezentate prin perechi de valori (e,f), unde e este partea exponențială iar f este partea fracțională. Perechea are semnificația

Aici e este un întreg avînd valoarea specificată, iar f este o fracție cu semn. Vom adopta convenția modul din f fiind mai mic ca unu. Punctul bazei apare la stînga reprezentării poziționale a lui f, iar afirmația că avem numere cu p cifre înseamnă că

B. Calcule normalizate. Un număr cu punct mobil (e,f) este normalizat dacă cifra cea mai semnificativă a reprezentării lui f este diferită de zero, astfel că 1/b | f | < 1, sau dacă f = 0 iar e are cea mai mică valoare posibilă. Pentru a afla care dintre două numere cu virgulă mobilă normalizate este mai mare, comparăm mai întîi părțile exponențiale, iar apoi testăm părțile fracționale numai dacă exponenții sînt egali.

Valoarea (0, .00000001) nu poate fi normalizată fără a produce un exponent negativ dar se obține cîștig de viteză, uniformitate, precum și posibilitatea de a fixa în mod simplu limite pentru erorile relative din calculele noastre.

În continuare vom studia operațiile cu punct mobil normalizat, prezentăm pașii și algoritmul pentru adunarea numerelor cu punct mobil.

Deoarece arirmetica numerelor cu punct mobil este în mod natural aproximativă, nu

exactă, vom utiliza simbolurile „rotunjite”

+, -, ×, ÷

pentru a specifica adunarea, scăderea, înmulțirea șî împărțirea reprezentării cu punct mobil, deoarece aritmetica aceasta este în mod natural aproximată, nu exactă .

Ideea de bază implicată în adunarea numerelor cu punct mobil este destul de simplă: Presupunînd că

se va lua

normalizînd apoi rezultatul .

Algoritmul A ( Adunarea numerelor cu punct mobil )

Avem numerele

reprezentate cu punct mobil in baza b, în exces de q, cu p cifre și normalizate. Algoritmul efectuează suma w = u + v. Acest algoritm se poate folosi și la scădere dacă se înlocuiește v cu (-v).

Pas A1. [ Separarea ] Din reprezentările numerelor u și v se separă părțile fracționale.

Dacă

se schimbă între ele numerele u și v. În multe cazuri, este mai bine să se

combine pasul A2 cu pasul A1, sau cu unul din pașii următori.

diferență mare între exponenți, se introduce

și se trece la Pas A7. În realitate, deoarece presupunem că u este normalizat,

am putea termina algoritmul; dar este uneori util să putem normaliza numerele

posibil nenormalizate.

Pas A5. [ Aliniere la dreapta ] Se deplasează

poziții adică este împărțit prin

Pas A6. [ Adunarea ] Se efectuează

Pas A7. [ Normalizarea ] La acest punct

reprezintă suma lui u cu v, dar fw poate avea mai mult de p cifre și poate fi

mai mare decît unitatea sau mai mic decît 1/b.

Se execută algoritmul N prezentat în continuare, pentru a normaliza și rotunjii

2

| f | ≥ 1

Da

f = 0

Da

1

Nu

Nu

Desenul conține normalizarea lui (e,f). 1- Depășire superioară sau inferioară; 2-Depășire la rotunjire.

Algoritmul N ( Normalizarea )

Un "exponent brut" e și o "fracție brută" f sunt convertite la forma normalizată, rotunjită dacă este necesar la p cifre. Algoritmul presupune că

| f | < b.

Pas N1. [ Testare f ] Dacă | f | 1 ("depășire parte fracționară"), se trece la Pasul N4.

Dacă f = 0, se atribuie lui e cea mai mică valoare pe care o poate avea și se

trece la Pasul N7.

Pas N2. [ Este f normalizat ] Dacă | f | 1/b, se trece la Pasul N5.

Pas N3. [ Aliniere la stînga ] Se deplasează la stînga cu o cifră (adică se înmulțește cu

b), iar e = e – 1. Se revine la Pasul N2.

Pas N4. [ Aliniere la dreapta ] Se deplasează la dreapta cu o cifră (adică se împarte cu

b), iar e = e + 1.

Pas N5. [ Rotunjire ] Se rotunjește f la p poziții. Operația de rotunjire poate conduce la

| f | = 1 ("depășire prin rotunjire"), în astfel de cazuri se revine la Pasul N4.

Pas N6. [ Verificare e ] Dacă e este prea mare, mai mare decît valoarea sa maximă

admisibilă, apare condiția de depășire superioară a exponentului. Dacă e este

prea mic, apare o condiție de depășire inferioară a exponentului.

Pas N7. [ Împachetare ] Se așează e și f împreună în reprezentarea dorită la ieșire.

Operațiile de înmulțirea și împărțirea, sunt mai simple decît adunarea, și se aseamănă între ele.

Algoritmul M (Înmulțirea sau Împărțirea în punct mobil)

Fie numerele

reprezentate cu punct mobil în baza b, în exces de q, cu p cifre și normalizate, acest algoritm efectuează produsul w = u × v sau cîtul w = u ÷ v.

Pas M1. [ Separarea ] Din reprezentările numerelor u și v se separă părțile fracționale.

Uneori este convenabil să se testeze la acest pas dacă operanzii sunt nuli.

Pas M2. [ Operarea ] Se efectuează

(1) – pentru înmulțire

(2) – pentru împărțire

Dacă este necesar, reprezentarea lui fw poate fi redusă la p+2 sau p+3 cifre

la acest pas.

Pas M3. [ Normalizarea ] Se execută "Algoritmul N" asupra lui

pentru a normaliza, rotunji și împacheta rezultatul.

Normalizarea este mai simplă în acest caz, deoarece alinierea la stînga apare cel mult o dată, iar depășirea de rotunjire nu poate să apară după împărțire.

I.2.2. Acuratețea aritmeticii cu punct mobil

Calculul în notație cu punct mobil este inexact și nu este deloc greu ca acesta să fie folosit greșit. O problemă principală a analizei matematice este determinare gradului de acuratețe pe care îl au rezultatele diferitelor metode numerice. Avem problema "intervalului de credibilitate" nu știm care din răspunsurile calculatorului sunt corecte.

Matematicienii care au încercat să analizeze operațiile cu punct mobil, au observat că sarcina este destul de dificilă, încît s-au mulțumit cu argumente plauzibile.

În acestă secțiune studiem cîteva dintre caracteristicile erorilor în aritmetica cu punct mobil, avînd scopul de a descoperii modalități de efectuare a calculelor cu punct mobil, pentru a ajuta la o analiză rezonabilă a propagării erorilor.

Prin reprezentarea numărului real exact y, folosim o aproximare ŷ = y(1+ε), cantitatea ε = (ŷ – y)/y se numește eroare relativă de aproximare. Operațiile de înmulțire și împărțire cu punct mobil nu măresc cu mult eroarea relativă, dar scăderea unor cantități aproape egale poate să creasor cantități aproape egale poate să crească cu mult eroarea relativă.

Unele legi ca de exemplu legea asociativități în adunare nu mai este valabilă :

( u + v ) + w u + ( v + w ) pentru mulțimile u, v, w.

Exemplu :

(11111113. + -11111111.) + 7.511111 = 2.0000000 + 7.5111111 = 9.5111111

11111113. + (-11111111. + 7.511111) = 11111113. + -11111103 = 10.000000.

Notațiile matematice "a1 + a2 + a3" sau " ak", 1kn sunt în mod inerent bazate pe presupunerea unei asociativități, astfel încît programatorul să fie atent în mod special asupra faptului că el nu presupune implicit valabilitatea legii asociativității.

A.O abordare axiomatică. Deși legea asociativității nu mai este valabilă, legea comutativității u + v = v + u este valabilă. Rutinele cu punct mobil trebuie să fie concepute astfel încît să prezerve, pe cît posibil, cît mai multe dintre legile obișnuite ale matematicii.

Considerăm cîteva dintre legile de bază care sunt valabile pentru operațiile cu punct mobil normalizate, avem astfel:

u – v = u + (-v) ;

-(u + v) = (-u) + (-v) ;

u + v = 0 dacă și numai dacă v = -u ;

u + 0 = u.

Se pot deduce și alte identități : u – v = -( v – u ).

Dacă u v atunci u + w v + w.

Ideea este ca operațiile cu punct mobil să satisfacă relațiile:

u + v = round(u + v),

u – v = round(u – v),

u v = round(u v),

u ÷ v = round(u / v).

round(x) – este cea mai bună aproximație cu punct mobil a lui x.

round(-x) = – round(x)

x ≤ y implică round(x) ≤ round(y).

Putem scrie de asemenea alte identități :

u × v = v × u, (-u ) × v = -(u × v), 1 × v = v ;

u × v = 0 dacă și numai dacă u = 0 sau v = 0 ;

(-u) ÷ v = u ÷ (-v) = -( u ÷ v ) ;

0 ÷ v = 0, u ÷ 1 = u, u ÷ u = 1 .

Dacă u v și w > 0, atunci uw v w și u ÷w v/w și w÷u w÷v .

Dacă u+v = u+v, atunci (u + v) – v = u; și dacă uv = uv ≠ 0, atunci ( u v) ÷ v = u .

Legea distributivă dintre × și + poate produce greșeli destul de mari.

Fie u = 20000.000, v = 6.0000000 și w = 6.0000003 atunci :

( u × v ) + ( u × w ) = -120000.00 + 120000.01 = 0.01

u × ( v + w ) = 20000.00 × .00000030000000 = 0.006

deci

u × ( v + w ) ( u × v ) + ( u × w ).

Pe de altă parte avem b × ( v + w ) = ( b × v ) + ( b × w ), cînd b este baza notației cu punct mobil, atunci round(bx) = b round(x).

Funcția round(x) nu este definită cînd | x | este prea mic sau prea mare.

O cale pentru a calcula medii și deviații standard folosind aritmetica cu punct mobil este de a utiliza formulele recurente:

Cînd

în care eroarea relativă este limitată independent de x:

Un exemplu de proceduri tipice pentru estimarea erorilor putem considera legea asociativă pentru înmulțire.

Numărului

i s-a acordat denumirea specială de ulp, avînd semnificația unei "unități pe ultima poziție" a părții fracționale. Operațiile cu punct mobil sînt corecte pînă la o jumătate de ulp. Legea asociativă pentru înmulțire este valabilă pînă la aproximativ doi ulp de eroare relativă.

S-a arătat că ( u × v ) × w este apoximativ egal cu u × ( v × w ), exceptînd cazul cînd depășirea superioară sau inferioară a exponentului reprezintă o problemă.

Compararea numerelor cu punct mobil ne ajută la stabilirea valorilor relative ale două mărimi reprezentate cu punct mobil. Pentru numerele cu punct mobil în baza b, exces q, ca de exemplu:

avem următoarele definiții :

uv () dacă și numai dacă

u v () dacă și numai dacă

uv () dacă și numai dacă

u v () dacă și numai dacă

Definițiile se aplică atît valorilor nenormalizate cît și celor normalizate.

Condițiile u v (strict mai mic decît), u v (aproximativ egal cu), u v (strict mai mare decît) trebuie întotdeauna să aibă loc pentru orice pereche de valori u și v. Relația

u v este mai puternică decît u v și ar trebui interpretată ca "u este în mod esențial egal cu v". Toate relațiile sunt date în termenii unui număr real pozitiv care măsoară gradul de aproximare avut în vedere.

Iată cîteva consecințe simple ale definiției de mai sus:

dacă u v () atunci u v () ;

dacă u v () atunci u v () ;

u u () ;

dacă u v () atunci u < v ;

dacă u v (1) și 1 2 atunci u v (2) ;

dacă u v (1) și 1 2 atunci u v (2) ;

dacă u v (1) și 1 2 atunci u v (2) ;

dacă u v (1) și v w (2) atunci u w (1+2) ;

dacă u v (1) și v w (2) atunci u w (1+2) .

putem arăta că:

| u – v | | u | și | u – v | | v | implică u v () ;

| u – v | | u | sau | u – v | | v | implică u v () ;

pentru numerele u și v normalizate cu punct mobil, cînd <1

u v () implică | u – v | b | u | și | u – v | b | v | ;

u v () implică | u – v | b | u | sau | u – v | b | v | .

Relațiile , , și sînt utile în cadrul algoritmilor numerici .

I.2.3. Calcule cu dublă precizie

În aritmetica cu punct mobil folosim un cuvînt pentru memorarea valorilor cu punct mobil. Precizia poate fi mărită folosind tehnici convenabile de programare care folosesc două, sau mai multe cuvinte de memorie pentru reprezentarea fiecărui număr. Dubla precizie este aplicabilă la mai multe probleme care nu necesită precizie extrem de mare. Dubla precizie este folosită nu numai pentru a extinde precizia părții fracționale a numărului cu punct mobil, dar și pentru a extinde rangul părții exponențiale.

Ne vom opri asupra versiunii cu punct mobil ale aritmeticii cu dublă precizie. Rutinele în dublă precizie nu vor rotunji rezultatele, uneori se admite strecurarea unei erori mici.

Rotunjirea corectă în dublă precizie necesită programare suplimentară, rutinele cu acuratețe integrală au nevoie de două ori mai mult spațiu, și cu 50% mai mult timp.

Aritmetica în precizie simplă este un "aliment de bază" pentru cei ce utilizează aritmetica cu punct mobil, dubla precizie se utilizează și în situații cînd rezultatele nu sunt importante. Diferența dintre acuratețea pe șapte și cea pe opt poziții (cifre) este considerabilă, dar rareori ne pasă de diferența între acuratețea cu 15 și cea cu 16 poziții. Dubla precizie este folosită pentru pașii intermediari din timpul calculării rezultatelor în precizie simplă.

Aceste proceduri sunt lipsite de acuratețe, deoarece ele sunt tipice pentru felul de simplificări făcute în general în rutinele incorecte.

Adunarea este efectuată adunînd împreună jumătățile cele mai puțin semnificative și separat jumătățile cele mai semnificative, propagînd în mod corespunzător transportul.

Este posibil să adunăm jumătățile cele mai puțin semnificative și să obținem semnul greșit. Cea mai simplă soluție este să se anticipeze corect semnul, în Pasul A2 (comform Algoritmului aduarea numerelor cu punct mobil) nu vom presupune numai Eu>Ev, vom mai presupune că | u | | v |, astfel știim că semnul final va fi semnul lui u. Adunarea în dublă precizie se aseamănă cu cea în simplă precizie, exceptînd faptul că totul este făcut de două ori.

Scăderea este convertită în adunare prin schimbarea semnului celui de-al doilea operand.

Ne sunt necesari doar opt baiți situați la extremitatea din stînga, este convenabil să se lucreze numai cu cifre situate la stînga barei verticale din schemă, înseamnă că nu mai este necesar să calculăm produsul celor două jumătați mai puțin semnificative.

Împărțirea cu punct mobil în dublă precizie este mai dificilă. Avem numerele care trebuie împărțite sub forma (um + εul)/(vm + εvl), unde ε este inversa mărimii maxime reprezentabile, iar vm se presupune că este normalizată.

Fracția se poate descompune astfel:

Întrucît 0 | vl | < 1 iar 1 / b | vm | < 1, avem | vl / vm | < b , iar eroarea datorată eliminării termenilor care implică ε² poate fi neglijată. Metoda constă din a calcula:

apoi scădem din rezultat produsul dintre:

Tabela următoare conține timpul mediu aproximativ de calcul pentru subrutinele în dublă precizie, și subrutinele în precizie simplă care apar în această capitol:

I.3. Aritmetica în precizie multiplă

Să considerăm acum operații pe numere, care au o precizie arbitrar de înaltă. Pentru simplitate presupunem că lucrăm cu numere întregi, în loc de numere care au în interior un punct al bazei.

I.3.1. Algoritmii clasici

În această secțiune vom discuta algoritmii pentru:

adunarea sau scăderea întregilor avînd n- poziții , dînd un răspuns avînd n- poziții și un transport ;

înmulțirea unui întreg avînd n- poziții cu un întreg avînd m- poziții, dînd un răspuns avînd ( n+m )- poziții ;

împărțirea unui întreg avînd (n+m)- poziții la un întreg avînd n- poziții, dînd un cît avînd (m+1)- poziții și un rest avînd n- poziții.

Aceștia pot fi numiți "algoritmi clasici", deoarece cuvîntul "algoritm" a fost folosit multe secole numai în legătură cu aceste procese. Termenul "întreg avînd n- poziții" înseamnă orice întreg mai mic decît bˆn, unde b este baza notației poziționale obișnuite în care sunt exprimate numerele, care pot fi scrise utilizînd cel mult n "poziții" în această notație.

În continuatre vom studia algoritmi care efectuează operațiile menționate înainte, pentru întregi exprimați în notația cu baza b, unde b 2 număr întreg dat.

Numerele cu precizie extinsă pot fi privite ca numere scrise în notație cu baza w, în care w este mărimea cuvîntului calculatorului. De exemplu, un întreg care se extinde pe 10 cuvinte ale unui calculator avînd mărimea cuvîntului w = 10¹º are 100 cuvinte, dar noi îl vom considera un număr cu 10 poziții în baza 10¹º. Acest punct de vedere este justificat prin aceleași motive ca și conversia de la notația binară la cea octală în care am grupat biții împreună.

În acești termeni, avem următoarele operații primitive:

aº) adunarea sau scăderea întregilor avînd 1- poziție, dînd un răspuns avînd

1- poziție și un transport;

bº) înmulțirea unui întreg avînd 1- poziție cu un întreg avînd 1- poziție dînd un

răspuns avînd 2- poziții;

cº) împărțirea unui întreg avînd 2- poziții la un întreg avînd 1- poziție, cu condiția

ca rezultatul să fie un cît întreg avînd 1- poziție și un rest avînd 1- poziție.

Prin ajustarea mărimii cuvîntului aproape toate calculatoarele vor avea la dispoziție aceste trei operații. Astfel vom construi algoritmii (a), (b) și (c) menționați anterior, în termenii operațiilor primitive (aº), (bº) și (cº).

Cînd b=10 și ne imaginăm că efectuăm operațiile manual, atunci, operația (aº) este analoagă cu memorizarea tabelei de adunare, (bº) este analoagă cu memorizarea tabelei înmulțirii iar (cº) este de fapt memorizarea tabelei înmulțirii inversată. Operațiile (a), (b) și (c) asupra numerelor cu precizie înaltă pot fi efectuate folosind simple adunări, scăderi, înmulțiri și împărțiri în format lung proceduri care le-am învațat în școala elementară.

Pentru a evita munca suplimentară pentru calcularea semnelor vom presupune că toate numerele cu care vom lucra în această secțiune sînt nenegative.

Vom trata pentru început adunarea, deoarece ideile care apar aici le vom folosi și la ceilalți algoritmi.

Algoritm A ( Adunarea întregilor nenegativi )

Avem întregii nenegativi cu n-poziții

acest algoritm alcătuiește suma lor în baza b,

Unde

este transportul, care va fi întotdeauna egal fie cu 0, fie cu 1.

Pas A1. [Inițializarea] Se efectuează j n, k 0. Variabila j va trece prin variantele

poziții ale cifrelor, iar variabila k păstrează urma depășirilor la fiecare pas.

Pas A2. [Adunarea cifrelor] Se efectuează:

k ia valoarea 0 sau 1, depinzînd dacă s-a produs sau nu un „transport“, adică

dacă:

sau nu. În timpul adunării celor două numere, cel mult un transport este

posibil, deoarece avem întotdeauna:

prin inducție asupra calculului.

Pas A3. [Ciclu pe j] j = j – 1.

Dacă j > 0 , mă întorc la Pasul A2

Altfel efectuez

și termin algoritmul .

Algoritmul scăderii este prezentat în continuare, el este similar cu cel al adunării.

Algoritmul S ( Scăderea întregilor nenegativi )

Avem întregii nenegativi cu n-poziții

acest algoritm produce diferența lor nenegativă în baza b,

Pas S1. [Inițializarea] Se efectuează j n , k 0.

Pas S2. [Scăderea cifrelor] Se efectuează:

k ia valoarea -1 sau 0, depinzînd dacă apare sau nu un "împrumut", adică:

sau nu. La calculul lui wj trebuie notat că avem:

Pas S3. [Ciclu pe j] j = j – 1.

Dacă j > 0 , mă întoarc la Pasul S2

Altfel termin algoritmul.

La înmulțire vom duce mai departe ideile folosite în algoritmul A.

Algoritmul M ( Înmulțirea întregilor nenegativi )

Avem întregii nenegativi

acest algoritm alcătuiește produsul lor în baza b,

Pas M1. [Inițializare] Toate variabilele

se aduc la valoarea zero. Se efectuează j m.

Pas M2. [Multiplicator zero ?]

Dacă

efectuez

și trec la pasul M6.

Pas M3. [Inițializare i ] Se efectuează i n, k 0.

Pas M4. [Multiplicare și adunare] Se afectuează:

apoi se face:

Aici transportul k va fi întotdeauna în domeniul 0 k < b;

Pas M5. [Ciclu pe i ] i = i -1.

Dacă i > 0 , revin la pasul M4;

Altfel efectuez

Pas M6. [Ciclu pe j ] j = j -1.

Dacă j > 0 , revin la pasul M2;

Altfel algoritmul se termină.

Algoritmul M este ilustrat în Tabelul 1, presupunînd că b=10, prin prezantarea stărilor calculului la începutul pașilor M5 și M6.

Tabelul 1

Înmulțirea lui 914 cu 84

Ultimul algoritm prezentat aici este împărțirea în format lung, în cadrul căreia dorim să împărțim întregi avînd (n+m)- poziții, la întregi avînd n- poziții. Metodele obișnuite cu creionul și hîrtia implică o oarecare intuiție și ingeniozitate din partea celui care efectuează împărțirea.

O analiză a împărțiri în format lung ne arată că problema generală se descompune în pași mai simpli, fiecare dintre acești pași fiind împărțirea unui număr u avînd (n+1)- poziții la un divizor v avînd n- poziții, unde 0 u / v < b; restul r după fiecare pas este mai mic decît v, astfel că vom putea folosi cantitatea rb + (următoarea poziție a dividentului) ca fiind noul u în urmatorul pas. De exemplu, dacă împărțim 3142 la 47, vom împărți mai întîi pe 314 la 47 obținînd 6 și un rest de 32; apoi împărțim pe 322 la 47 obținînd 6 și un rest de 40; astfel avem cîtul de 66 iar restul de 40. Această idee este valabilă în general însă dorim să găsim un algoritm de împătțire mai rapid și mai eficient, astfel căutările noastre se reduc la următoarea problemă:

Fie:

întregi nenegativi în notație cu baza b, astfel încît u/v<b. Să se găsească un algoritm pentru a determina [ u / v ].

q

v1v2 … vn)u0u1u2 … un

qv

r

Se caută o modalitate mai rapidă pentru determinarea lui q .

Observăm echivalența condiției

cu condiția:

adică:

aceasta este, pur și simplu, condiția ca:

Dacă notăm r = u – qv, atunci q este unicul întreg astfel ca 0 r < v.

Pentru a aborda această problemă facem o presupunere asupra lui q, bazată pe cele mai semnificative cifre ale lui u și v.

Astfel presupunem că:

Pe ˆq îl obținem prin împărțirea primelor două cifre cele mai semnificative ale lui u la cea mai semnificativă cifră a lui v, iar dacă rezultatul este b sau mai mare decît b, îl vom înlocui prin (b – 1).

Valoarea ˆq este întotdeauna, o foarte bună aproximare a răspunsului q, atîta timp cît v1 este suficient de mare. Pentru a arăta cît de mult se apropie ˆq de q vom enunța două teoreme importante.

Teorema A. În notația de mai sus, ˆq q.

Teorema B. Dacă:

În final concluzia este independentă de b, nu mai are importanță cît este b de mare, cîtul de încercare ˆq nu va fi niciodată eronat cu mai mult de 2.

Condiția

este condiția normalizării binare cu punct mobil într-un calculator binar. O modalitate simplă de a asigura ca v1 să fie suficient de mare, este multiplicarea atît a lui u cît și a lui v cu

acest lucru nu schimbă valoarera lui u / v, și nici nu mărește numărul de poziții în v.

Algoritmul folosește o alegere puțin mai îmbunătățită a lui ˆq în pasul D3, care garantează că q = ˆq sau ˆq-1; de fapt alegerea îmbunătățită a lui ˆq făcută aici se bucură aproape întotdeauna de acuratețe.

Algoritmul D ( Împărțirea întregilor nenegativi )

Avem întregii nenegativi

vom forma în baza b cîtul

Pas D1. [ Normalizare ] Se efectuează:

Pas D2. [ Inițializare j ] Se efectuează j 0. Ciclul pe j, pașii D2 la D7. vor fi de fapt

o împărțire a lui

pentru a obține o singură cifră a cîtului; comform figurii următoare:

Împărțirea în format lung

1 – q ≠ ˆq 2 – 1 ≤ j ≤ m

2

1

Pas D3. [ Calculul ˆq ]

Dacă

se pune ˆq b -1; dacă nu, se efectuează:

Se testează apoi dacă:

dacă este adevărat, se diminuează ˆq cu 1 și se repetă acest test.

Pas D4. [ Înmulțire și scădere ] Se înlocuiește

Acest pas constă dintr-o simplă multiplicare cu un număr avînd 1- poziție,

combinată cu o scădere. Cifrele

trebuie să fie menținute pozitive; dacă rezultatul acestui pas este însă negativ,

trebuie să fie lăsat ca valoare adevărată plus

adică ca și cînd ar fi complementul lui b al valorii reale și trebuie să

memorăm un "împrumut" la stînga.

Pas D5. [ Testarea restului ] Se introduce qj ˆq. Dacă rezultatul pasului D4 a fost

negativ, se trece la pasul D6; dacă nu, se trece la pasul D7

Pas D6. [ Adunare înapoi ] Se diminuează qj cu 1 și se adună

Va apărea o depășire la stînga lui uj, dar aceasta trebuie ignorată întrucît se

anulează reciproc cu "împrumutul" de la pasul D4.

Pas D7. [ Ciclu pe j ] Se majorează j cu 1. Acum, dacă j ≤ m, se revine la pasul D3.

Pas D8. [ Denormalizarea ]

Acum

este cîtul dorit, iar restul operației poate fi obținut împărțind

II. Calculul preciziei multiple folosind Mfloat

Precizia pentru standardul IEEE-754 de reprezentare a numerelor în virgulă mobilă este destul de bună pentru majoritatea aplicațiilor științifice și tehnice. Dar uneori un om de știință sau un matematician are probleme cu calculele deoarece algoritmul sau problema de matematică este prost determinată (de exemplu găsirea rădăcinii unui polinom). Deoarece numerele în virgulă mobilă sunt numai o reprezentare a numerelor reale și raționale, există cîte o eroare adițională la fiecare operație cu virgulă mobilă.

Nu există nici o metodă care să reprezinte fiecare număr real exact într-o formă binară și de aceea o aproximare este necesară. O mertodă practică este aceea de a extinde mantisa numerelor în virgulă mobilă, astfel rezulatatul devine mai exact cu fiecare operație. Dar există o limită pentru că fiecare număr în virgulă mobilă este necesar o memorie mai extinsă iar timpul necesar rezolvării este mai mic. Din această cauză un om de știința care are probleme numerice de acest gen caută niște subrutine rapide pentru numere în virgulă mobilă.

Acest pachet reprezintă un asemenea pachet de subrutine. Informația este folosită în formă binară pentru că există instrucțiuni optine pentru manipularea numerelor binare în setul de instructiuni 8086 și 80386. Un format binar are avantajul, că folosește memoria minimizată (mai restrînsă). Singurul dezavantaj este efortul mai mare care trebuie depus pentru intrările și ieșirile numerelor în formatul ASCII.

Scopul este un timp mai scurt pentru numere foarte precise. Așa că nu se pierde timpul pentru rotunjiri (fiecare rezultat este trunchiat) și gestionarea valorilor speciale (plus, minus infinit sau rezultat nedeterminat). Singurul număr special este zero.

În eventualitatea unei depășiri sau altor erori de calcul se afișează un indicator de eroare. Bine înțeles că în asemenea cazuri rezultatul numeric este greșit. Pentru dezvoltarea unui program care folosește numerele în virgulă mobilă, programul poate fii întrerupt la o eroare, dacă folosiți un alt indicator ("useInterruptAtError=1"). Pe ecran apare un mesaj de eroare indicînd eroarea (de exemplu: "Eroare MFLOAT: împărțire cu zero"). În cazul unei subdepășiri rezultatul este pus la zero fără nici un mesaj.

Există o gamă largă pentru acest exponent. Cel mai mare număr care poate fi reprezentat este +/-7.07E+9863, iar cel mai mic este aproximativ +/-7.07E-9865.

Pentru a evita timpul de a căuta erorile de supra încărcare a stivei este inclus un test de stivă în procedura de asamblare.

II.1. Calculularea lui PI

Există multe metode de a calcula numărul pi care reprezintă relația dintre perimetrul și diametrul unui cerc.

Funcția arctan(x) se calculează folosind seria:

Avantajul primelor două ecuații este acela că converg foarte rapid și fiecare element al seriei poate fi calculat recursiv din vechiul element numai prin două diviziuni ale unui întreg. O divizare print-un întreg este foarte rapidă pentru numerele în virgulă mobilă lungi. A treia ecuație este alcătuită din cunoscuta valoare

Această metodă este mai înceată decît cea precedentă ei. A patra ecuație are doar importanță teoretică, în folosirea numerică dar ecuația este foarte greșită. Programul exemplu compară diferite metode de a calcula pi.

Alte formule pentru pi:

Aceste formule sunt mai lente decît prima a lui John Machin. Toarte rezultatele programului test sunt corecte pînă la 74 sau 75 cifre zecimale după virgulă.

O bună metodă de convergență pentru a calcula pi este însemnătatea aritmetico-geometrică. Această metodă a fost populară printre cei care au folosit pi cu 10^9 cifre pe super calculatoare. Metoda folosește multiplicatori FFT.

Se consideră secvența următoare:

Apoi secvența a(k) converge la 1/pi mai repede decît pătratic, adică, obținem mai mult decît dublu numărului cifrelor a fiecărei iterație.

Se obține:

y(0) = 7.07E-1

y(1) = 1.72E-1

y(2) = 7.47E-3

y(3) = 1.40E-5

y(4) = 4.86E-11

y(5) = 5.92E-22

y(6) = 8.75E-44

y(7) = 1.91E-87

y(8) = 9.16E-175

y(9) = 2.08E-349

y(10) = 1.10E-698

.

.

Dacă seria converge foarte bine, schimbarea lui a(k) este următoarea :

Această metodă este demonstrată în programul "PIGM.Pas". Rezultatul pi este corect cu 74 cifre zecimale pentru acest exemplu.

II.2. Calculul Funcției Bessel

Evaluarea funcției Bessel pentru argumente extinse: Există unele serii unde unele elemente sunt foarte extinse iar rezultatul este foarte mic. În acest caz erorile numerice sunt mari. Un exemplu în funcția Bessel este funcția JO(x) pentru un argument x=100,0 .

Seria converge foarte bine dar există probleme privind acuratețea calculului. Programul test arată influența exactității calculului față de rezultat.

II.3. Calculul constantei lui Euler

Constanta este definită de următoarele limite:

Evaluarea directă a limitei prin punerea lui n foarte sus nu are succes deoarece converge în mod jalnic. Formulele Euler-Maclaurin ajută la calculul aceastor limite.

Constantele B(k) sunt numere Bernoulli, calculate prin:

Dacă se folosește

Ceea ce rămîne poate fi estimat:

Folosind această formulă constanta gamma poate fi calculată pîna la 100 zecimale foarte ușor. Eroarea se datorează estimării a ceea ce rămîne și descrește precum "n" crește la o valoare optimă, dar nu există convergență. Pentru valori mai mari a lui "n" eroarea crește și nu există limită sus. Seria este una asimtotică. Pentru o acuratețe mai mare numărul "m" poate fi mărit, dar efortul de calcul crește și el. 1000 de zecimale nu sunt o problemă pentru numerele în virgulă mobilă dar nu și pentru această versiune (limitată la 77 zecimale). Rezultatul programului test este corect pînă la 74 zecimale.

O altă metodă pentru calcularea constantei Euler e sugerată de Sweeney și o formulă este găsită de Brent.

Această formulă este probabil foarte precisă, mai rapidă și nu necesită multă memorie, dar acuratețea rezultatelor intermediare trebuie să fie mai mare decît cea a rezultatelor. De aceea rezultatul ce folosește numerle în virgulă mobilă are 51 de zecimale corecte (aproximativ 2/3 din precizia maximă).

II.4. Calculul Funcției Gamma pentru valori complexe

Funcția Gamma este definită prin:

Există și alte definiți. Această integrală nu se folosește pentru calculele numerice. Se folosește formula:

Mărginirea restului :

Aceste formule sunt foarte bune pentru valori mari pentru Re{z}. Pentru a avea o valoare înaltă, pentru partea reală a lui "z", formula recurentă este folosită în mod repetat.

Aceasta se face, pănă cînd o parte reală a lui "z" este realizată.

Metoda este folosită de programul "Gamma.PAS" utilizat în PASCAL pentru argumente reale.

III. Metode MFLOAT

Un număr în precizie multiplă este alcătuit dintr-un tablou de cuvinte (1 cuvînt = 1 byte = 16 biți). Primul cuvînt reprezintă exponentul bazei doi. Și este reprezentat printr-un deplasament de 8000 Hex (=32768) cu domeniul de la 1 la FFFF Hex (=65535). Valoarea specială 0 (zero) este utilizată pentru a reprezenta numărul în virgulă mobilă zero. În acest caz restul tabloului este arbitrar. Al doilea cuvînt din tablou reprezintă primul cuvînt al mantisei. Mantisa se normalizează pentru fiecare număr în virgulă mobilă. De aceea bitul cel mai semnificativ al primului cuvînt de mantisă este întotdeauna setat. Bitul de semn este unu, dacă numărul în virgulă mobilă este negativ, și zero, dacă numărul în virgulă mobilă este pozitiv. Următoarele cuvinte ale numărului în precizie multiplă sunt cuvinte ale mantisei. Lungimea mantisei poate fi între 1 și 16 cuvinte și poate fi aleasă de utilizator.

Exemplu: 8002 490F Hex reprezintă 3,141540528…

Toate cuvintele pentru mantisă reprezintă în medie 4,8165 (=16*log 10(2)) cifre ale numărului zecimal. Deci avem nevoie de 16 cuvinte pentru mantisă pentru o acuratețe de 77 cifre zecimale iar tot numărul necesită 17 cuvinte sau 34 bytes.

Și acum despre metoda disponibilă în acest format:

Fundamentul MFLOAT se compune din elemente de algebră combinate cu asamblarea. Scopul principal este de a reduce efortul de calcul, deoarece este o problemă importantă pentru aritmetica în precizie înaltă, dacă nici o componentă hardware ca coprocesorul nu este folosit. Doar unele probleme matematice pot fi rezolvate, doar dacă avem disponibile subrutine putenice (de viteză) pentru operații de bază.

Programele de nivel înalt ca Pascalul sunt relativ mai lente în comparație cu asamblorul, deoarece un program scris în nivel înalt nu are acces la bitul de transport al procesorului. De aceea bitul de transport trebuie calculat folosind o altă metodă consumatoare de timp. Este tipic, un factor de timp cîștigat de la 6 la 20, dacă folosim asamblarea în locul programelor de nivel înalt. Mai multă abilitate și mai mult timp este necesar în programarea 8086, dacă viteza este maximă. Dimensiunea codului utilizat nu este o problemă deoarece în ziua de astăzi există destulă memorie. De aceea nu este o problemă să scriem o serie sinus într-o singură subrutină în limbaj de asamblare.

O altă idee este să nu pierdem timpul cu rotunjirea rezultatului. Astfel sunt folosite numai trunchierile.

Este considerată mai întîi adunarea a două numere. Primul pas este de a deplasa mantisa numărului, care are valoarea absolută mai mică, iar exponentul numărului în virgulă mobilă este același. Majoritatea timpului este pierdut cu această operație, doar dacă nu este făcută corect. Ideea de bază a acestui algoritm este de a evita alocarea masivă de memorie, deoarece această operație este mare consumatoare de timp.

Multiplicarea a două numere este împărțită în subprobleme de multiplicare a unui cuvînt pentru mantisă a numărului "a" cu un cuvînt pentru mantisă a numărului "b" și adăugarea acestuia la rezultat. Prin adăugarea acestor produse doar un singur cuvînt mantisă este calculat pentru formatul final (cuvîntul suplimentar) iar toate celelalte produse care au efecte minore asupra rezultatului sunt neglijate. Pentru numere în precizie înaltă este obținut un factor de viteză de doi și doar puțină acuratețe este pierdută.

Exemplu (pentru numere zecimale):

7,34 * 2,14

–––––

29??

0734?

1468

–––––

1570??

? – rerprezintă cifre care nu se iau în calcul (sunt neglijate) la rezultat

rezultatul este: 15,7 folosind 3 cifre zecimale.

Pentru un număr egal de cuvinte mantisă și folosind instrucțiunile 80386 este posibilă impărțirea întregi probleme în multiplicări de cîte două cuvinte (32 biti). Folosind această metodă se obține un factor de viteză de la 1,5 pînă la 3.

Împărțirea a două numere "mfloat" este făcută ca și în calcul manual. Este estimată prima cifră iar restul este calculat. Dacă restul este mai mare decît divizorul sau dacă este negativ estimarea este corectată. Pentru o cifră se folosesc numere de la 0 la 65535 (1 cuvînt = 16 biți). Pentru estimarea, primelor două cuvinte mantisă de număr "a" sunt împărțite cu primul cuvînt mantisă al numărului "b". De vreme ce numerele sunt într-un format în virgulă mobilă normalizat, primul cuvînt mantisă pentru "b" este cuprins între [32768, 65535]. Împărțirea poate fi făcută de instrucțiunea DIV al 8086 în cîteva cicluri. Poate fi arătat că această eliminare este mai mare sau egală cu rezultatul corect și poate fi maxim de două ori mai mare decît rezultatul exact. Încă o dată un cuvînt suplimentar previne o acumulare de erori de trunchiare. Pentru un număr egal de cuvinte mantisă și folosind instrucțiunile 80386, este obținut un factor de viteză de la 1,5 pînă la 3.

Cea mai complicată subrutină este rădăcina pătratică. Algoritmul este același cu cel utilizat în calculul manual. Primul cuvînt mantisă este calculat bit cu bit de către metoda estimării și a corectării. Celelalte cuvinte mantisă sunt estimate ca cuvinte întregi și corecte ca și algoritmul împărțirii.

Prima dată cuvîntul mantisă al rezultatului este estimat prin împărțirea primelor două cuvinte cu primul cuvînt al rezultatului cu ajutorul instrucțiuni DIV din 8086, iar restul este calculat prin multiplicări repetate ale cuvintelor mantisă. Dacă rezultatul este negativ estimarea este corectată pînă cînd devine corectă. Un singur cuvînt suplimentar este folosit pentru a prevenii mari acumulări de erori de trunchiat. Interesant este că pentru numerele în precizie înaltă calcularea rădăcinii pătratice este mai rapidă decît calculare împărțirii. Pentru un număr egal de cuvinte mantisă și folosind instrucțiunile 80386, este obținut un factor de viteză de la 1,5 pînă la 3.

Și acum funcțiile transcendente:

Dar putem folosi de asemenea:

O serie mult mai simplă trebuie calculată deoarece ea constă în jumătate de termeni. Majoritatea timpului este folosit pentru calcularea seriei. Astfel că ea este implementată într-o subrutină în limbaj de asamblare. Prima dată se calculează x*x, ("*" – reprezintă înmulțirea) apoi se folosește o formulă pentru temenii seriei. Doar produsele esențiale sunt calculate.

Seria sinh(x) trebuie calculată avînd argumentul în domeniul [-0,346, 0,346] (= 0.5*ln(2) ), deoarece este folosită următoarea formulă:

Calcularea lui 2 la puterea (k) este o simplă manipulare a exponentului a două numere în virgulă mobilă. Funcțiile sinh(x), cosh(x), tanh(x), coth(x) și exp(x) pot fi derivate din această serie de bază și ecuația:

Este important să notăm că viteza rădăcinii pătratice este un mare avantaj pentru această subrutină.

Algoritmi similari sunt folosiți pentru funcția sin(x). este necesar să calculăm sin(x) în domeniul [-0,3926, 0,3926] (pi/8).

Formulele folosite sunt:

Logaritmii se pot calcula după formula:

Este necesar să calculăm seria în intrvalul [-0,1715, 0,1716] (=3-2*sqrt(2)), deoarece următoarea ecuație este folosită:

și "x" poate fi ales din intervalul [0,7071, 1,4142].

Seria de bază pentru inversa funcției trigonometrice este seria arctan:

Din nou este important să reducem argumentul seriei, deaorece comportamentul de convergență nu este mai bun decît cel al seriei artanh(x).

Formulele folosite sunt:

Dacă această formulă se folosește în mod repetat argumentul seriei este redus considerabil. Această formula se folosește pînă cînd valoarera absolută a argumentului este mai mică decît 0,0625. Este posibil să găsim un optim al numărului de reducere, dacă este investigat timpul de lucru, dar depinde de numărul de cuvinte pentru mantisă folosite și de argumentul însuși. Din această serie de bază derivă următoarele funcții: arcsin(x), arccos(x), arccot(x), atan2(x).

IV. Precizia MFLOAT

MFLOAT este aritmetica în virgulă mobilă și de aceea rezultatul normal nu este exact. Pentru a cîștiga viteză, nu trebuie să se piardă timp cu rotunjirea rezultatului, deoarece se cîștigă aproximativ un bit suplimentar sau 0,301 cifre zecimale. Este mai eficient să folosim cu un cuvînt de mantisă mai mult și se va cîștiga 16 biți sau 4,816 cifre zecimale. MFLOAT folosește metoda de rotunjire a lui zero (cu excepția unor cazuri speciale de adunare și de scădere), dar rezultatele nu sunt rotunjite corect întotdeauna. Erorile operațiilor algebrice sunt mai mici decît +/-1,032 LSB (least signifiant bit = bitul cel mai puțin semnificativ) pentru o acuratețe de 10000 cifre zecimale. Eroarea maximă este calculată de formula:

Pentru funcții transcendente ( sqrt(x), sin(x), exp(x), … ) erorile sunt ceva mai mari și depind de argumentul funcției.

V. Versiunea MFLOAT pentru PASCAL

Se include următoare linie ca prima linie a codul sursă, în fața oricărei declarații sau funcții : {$M 20000, 0, 650000}.

Această comandă instalează o stivă mare (20000 bytes) pentru program. Aceasta se folosește dacă avem mai multe variabile locale. Mărimea maximă permisă a stivei este aproximativ 65000.

Se definește lungimea maximă a numerelor în virgulă mobilă, dacă vrei mai puțină precizie dar mai multe variabile.

Schimbă linia a 12 din "PFLOAT.PAS", de exemplu: CONST MfloatWords =10

Acum există 9 cuvinte mantisă și un cuvînt de exponent.

Domeniul constantei MfloatWords este: [2 .. max+1].

"max"este definit în MFLOATA.OBJ și nu poate fi schimbat.

77 cifre zecimale … max = 16 (cuvinte mantisă)

2408 cifre zecimale … max = 500 (cuvinte mantisă)

Dacă folosiți instrucțiunile 80386, este bine să se folosească o valoare impară pentru "mfloatwords" (rezultă numere pare de cuvinte mantisă –> calculare mai rapidă). Nu are sens să se definească precizie mai mare decît "max+1", deoarece toate procedurile pot lucra numai cu numere de "max" cuvinte mantisă.

Se include instrucțiunea "USES" la începutul codului sursă. Toate funcțiile pentru tipul mfloat sunt incluse dacă este folosită sub forma:

USES pfloat;

Definește tratarea erorii, dacă nu vrei să se întrerupă programul la o eroare. Se include la începutul codului sursei:

BEGIN

.

useInterruptAtError=false; { definit în "PFLOAT.PAS" }

.

END.

Valoarea implicită este adevărată, aceasta înseamnă că programul se termină la o eroare. Desigur se poate schimba acest lucru în fișierul "PFLOAT.PAS".

Definește precizia numerelor de tip MFLOAT:

Se include la începutul codului sursă :

BEGIN

.

SetMantissawords(5); { definit în "PFLOAT.PAS" }

.

END.

Valoarea implicită este 1, înseamnă că, MFLOAT este calculat cu precizie minimă. Valoarea variabilei "mantissawords" nu se modifică direct, deoarece se poate depăși domeniul permis care este de 16 cuvinte mantisă.

Se alege opțiunile de compilare:

Menu: Options/Compiler/Runtime errors: Stack checking.

Această opțiune este utilă în faza de elaborare a programelor. Verificarea stivei din secvența scrisă în asamblare este independentă față de această verificare.

Iată un exemplu pentru începutul codului sursă:

{$M 20000, 0, 650000}

PROGRAM myprog(input,output);

.

.

.

BEGIN

{ program principal }

useInterruptAtError := false;

Setmantissawords(10);

.

.

.

END.

V.1. Descrierea funcțiilor de bază mfloat

Tipurile următoare sunt disponibile în "PFLOAT.PAS":

TYPE mfloat = ARRAY[0..MfloatWords-1] OF integer;

Acest tip de declarație pentru numere în precizie multiplă setează o precizie maximă de 16 cuvinte mantisă sau aproximativ 77 cifre zecimale. Precizia maximă este o constantă fixată în momentul asamblării. În această bibliotecă de funcții precizia maximă este fixată la 16 cuvinte mantisă. Este posibil să se schimbe declarația de tip pentru a reduce numărul de cuvinte mantisă. Aceasta duce la o economie de memoriei și la o viteză de execuție mai mare cu costul reducerii preciziei. Este posibil să se ajusteze precizia dinamic folosind "PROCEDURE Setmantissawords".

V.1.1 Funcțiile de bază mfloat

PROCEDURE Setmantissawords(number : integer);

Această procedură stabilește lungimea mantisei numărului în virgulă flotantă. Numărul de cuvinte trebuie să fie cuprins între [1…min(MfloatWords-1,max)]. altfel se limitează la aceste valori. max = 16 pentru 77 cifre zecimale versiune max = 500 pentru 2408 cifre zecimale versiune. Dacă numărul de cuvinte mantisă este mărit și numărul calculelor precedente sunt folosite mai departe, cuvintele adiționale nu sunt inițializate cu zero, valorile lor nu sunt cunoscute. Eroarea numărului datorată următoarelor cuvinte mantisă este foarte mic și are o acuratețe foare apropiată de cea obținută dacă ar fi inițializate cu zero. Dacă e important ca cuvintele mantisă să fie zero, programatorul trebuie să le inițializeze. Dacă numărul de cuvinte de mantisă e redus, fiecare număr poate fi văzut ca trunchiat, deși, nimic nu se schimbă la aceste numere.

FUNCTION GetMantissawords : integer;

Această funcție returnează numărul actual de cuvinte mantisă.

PROCEDURE ResetError;

Procedură de resetare a indicator de eroare.

FUNCTION GetError : boolean;

Se returnează indicatorul de eroare. Dacă există o eroare de execuție în timpul calculului, indicatorul de eroare este setat pe true.

PROCEDURE mfloatError(Nr : integer);

Indicatorul de eroare este setat. Dacă variabila booleană "useInterruptAtError" este adevărată, un mesaj de eroare va fi afișat pe ecran și programul se va termina.

PROCEDURE equm(VAR a, b : mfloat); { *** a <– b *** }

Valoarea lui "b" este copiată în variabila "a" cu precizia definită. Valoarea lui "b" nu se schimbă.

PROCEDURE addm(VAR a, b : mfloat); { *** a <– a + b *** }

Valoarea lui "b" este adunată la valoarea lui "a" cu precizia definită. Rezultatul este depus în variabila "a". Valoarea lui "b" nu se schimbă.

PROCEDURE subm(VAR a, b : mfloat); { *** a <– a – b *** }

Valoarea lui "b" este scăzută din valoarea lui "a" cu precizia definită. Rezultatul este depus în variabila "a". Valoarea lui "b" nu se schimbă.

PROCEDURE multm(VAR a, b : mfloat); { *** a <– a * b *** }

Valoarea lui "a" este înmulțită cu valoarea lui "b" cu precizia definită. Rezultatul este depus în variabila "a". Valoarea lui "b" nu se schimbă.

PROCEDURE divm(VAR a, b : mfloat); { *** a <– a / b *** }

Valoarea lui "a" este împărțită cu valoarea lui "b" cu precizia definită. Rezultatul este depus în variabila "a". Valoarea lui "b" nu se schimbă.

PROCEDURE multi(VAR a : mfloat; b : integer); { *** a <– a * b *** }

Valoarea lui "a" este înmulțită cu valoarea întreagă a lui "b" cu precizia definită. Rezultatul este depus în variabila "a".

PROCEDURE divi(VAR a : mfloat; b : integer); { *** a <– a / b *** }

Valoarea lui "a" este împărțită cu valoarea lui "b" cu precizia definită. Rezultatul este depus în variabila "a".

PROCEDURE inversm(VAR a : mfloat); { *** a <– 1 / a *** }

Numărul unu este împărțit la numărul "a" cu precizia definită. Rezultatul este depus în variabila "a".

PROCEDURE negm(VAR a : mfloat); { *** a <– – a *** }

Această procedură schimbă semnul numărului variabilei "a".

FUNCTION eqZero(VAR a : mfloat) : boolean; { *** eqZero <– a = 0 *** }

Această funcție verifică, dacă numărul în precizie multiplă "a" este zero. Valoarea lui "a" nu este schimbată.

FUNCTION gtZero(VAR a : mfloat) : boolean; { *** gtZero <– a > 0 *** }

Această funcție verifică, dacă numărul în precizie multiplă "a" este mai mare decît zero. Valoarea lui "a" nu este schimbată.

FUNCTION geZero(VAR a : mfloat) : boolean; { *** geZero <– a >= 0 *** }

Această funcție verifică, dacă numărul în precizie multiplă "a" este mai mare sau egal cu zero. Valoarea lui "a" nu este schimbată.

FUNCTION gtm(VAR a, b : mfloat) : boolean; { *** gtm <– a > b *** }

Această funcție verifică, dacă numărul în precizie multiplă "a" este mai mare decît numărul în precizie multiplă "b". Valoarea lui "a" și "b" nu este schimbată.

FUNCTION eqm(VAR a, b : mfloat) : boolean; { *** eqm <– a = b *** }

Această funcție verifică, dacă numărul în precizie multiplă "a" este egal cu numărul în precizie multiplă "b". Valoarea lui "a" și "b" nu este schimbată.

PROCEDURE GetZerom(VAR a : mfloat); { *** a <- 0 *** }

Această procedură copiază valoarea zero = 0,0000.. în variabila "a".

PROCEDURE GetOnem(VAR a : mfloat); { *** a <- 1 *** }

Această procedură copiază valoarea unu = 1,0000.. în variabila "a".

PROCEDURE GetPim(VAR a : mfloat); { *** a <- pi *** }

Această procedură copiază valoarea pi = 3,14159.. în variabila "a".

PROCEDURE GetLn2m(VAR a : mfloat); { *** a <- ln(2) *** }

Această procedură copiază valoarea ln(2) = 0,69314.. în variabila "a".

FUNCTION strtomf(VAR a : mfloat; { *** a <– string *** }

b : string) : integer;

Această procedură convertește șirul "b" într-un număr în precizie multiplă "a". Dacă conversia se realizează corect, se returnează valoarea zero. Altfel este returnată poziția caracterului unde apare eroarea. Pentru conversie, precizia calculării este stabilită cu un mantissaword mai mult față de cum era specificat și apoi este repus la valoarea originală. Această funcție nu este optimizată pentru viteză, doar pentru precizie, deoarece, în mod normal nu este folosită foarte des.

FUNCTION mftoa(VAR a : mfloat; { *** string <– a *** }

len : integer) : string;

Această procedură convertește un număr în precizie multiplă "a" la un șir "str" folosind numărul de caractere ASCII, care sunt definite în "len". Dacă numărul caracterelor este prea mic, șirul este setat cu '*'. Formatul folosit pentru această conversie este ',32767F' (vezi FUNCTION mftostr). Pentru conversie, precizia calculului se stabilește la o valoare, astfel ca eroarea să nu aibă efect asupra cifrelor zecimale calculate iar timpul conversiei să nu fie prea mare. Apoi, se setează înapoi la valoarea originală. Această funcție nu este optimizată pentru viteză, doar pentru precizie, deoarece, în mod normal nu este folosită foarte des. Lungimea maximă a șirului este 255 caractere. Dacă specificați "len" mai mare decît 255, lungimea este trunchiată la 255 caracter.

FUNCTION mftostr(VAR a : mfloat; { *** string <– a *** }

len : integer;

format : string) : string;

Această procedură convertește un număr în precizie multiplă "a" la un șir "str" folosind numărul maxim de caracte ASCII, care sunt definite în "len". Dacă numărul caracterelor este prea mic, șirul este setat cu '*'. Ordinea caracterelor în șirul format este arbitrară. Este de asemenea permisă folosirea literelor mici pentru format, dar în acest caz litera pentru exponentul zece este de asemenea litera 'e'. Pentru conversie, precizia calculului este stabilită la o valoare, astfel ca eroarea să aibă puțin efect pentru zecimale arătate și timpul conversiei să nu fie prea mare. Apoi, se setează înapoi la valoarea originală. Această funcție nu este optimizată pentru viteză, doar pentru precizie, deoarece, în mod normal nu este folosită foarte des. Lungimea maximă a șirului este 255 caractere. Dacă specificați "len" mai mare decît 255, lungimea este trunchiată la 255 caracter.

FUNCTION MfToD(VAR a : mfloat) : double; { *** MfToD <- a *** }

Această funcție convertește un număr în precizie multiplă într-un număr în dublă precizie. Valoarea lui "a" nu este schimbată.

FUNCTION MfToLd(VAR a : mfloat) : extended; { *** MfToLd <- a *** }

Această funcție convertește un număr în precizie multiplă într-un număr de tip extended (număr în virgulă mobilă reprezentat pe zece octeți). Valoarea lui "a" nu este schimbată.

PROCEDURE DToMf(VAR a : mfloat; b : double); { *** a <- b *** }

Această procedură convertește un număr în dublă precizie într-un număr în precizie multiplă.

PROCEDURE LdToMf(VAR a : mfloat; b : extended); { *** a <- b *** }

Această procedură convertește un număr de tip extended într-un număr în precizie multiplă.

PROCEDURE IntToMf(VAR a : mfloat; b : integer); { *** a <- b *** }

Această procedură convertește un număr întreg într-un număr în precizie multiplă.

V.1.2. Funcții standard de bibliatecă (Borland C: MATH.H)

PROCEDURE atanm(VAR a : mfloat); { *** a <- arctan(a) *** }

Această procedură calculează inversul funcției tan(x) avînd argumentul "a". Rezultatul este depus în variabila "a".

-pi/2 a pi/2 domeniul rezultatului

PROCEDURE atan2m(VAR a, b : mfloat); { *** a <- arg(b + ia) *** }

Acestă funcție calculează funcția extended arctan phi = arctan(a/b). Dacă "a" este zero sau mai puțin decît zero, valoarea corectă "phi" este returnată, care este cuprinsă în intervalul (-pi, pi].

a = r * sin(phi)

b = r * cos(phi)

r = sqrt(sqr(a) + sqr(b)) ≥ 0

Dacă a=0 și b=0, valoarea "phi" nu este unică. Indicatorul de eroare este setat și numărul zero este returnat. Variabla "b" nu se schimbă.

PROCEDURE ceilm(VAR a : mfloat); { *** a <– ceil(a) *** }

Ultimul numărul intreg cel mai mic, care este mai mare sau egal cu numărul mfloat "a", este calculat și returnat variabilei "b".

PROCEDURE cosm(VAR a : mfloat); { *** a <- cos(a) *** }

Procedura calculează funcția cos(x) avînd argumentul "a". Rezultatul este depus în variabila "a".

PROCEDURE expm(VAR a : mfloat); { *** a <- exp(a) *** }

Această procedură calculează funcția exponențială exp(x) avînd argumentul "a".

Rezultatul este depus în variabila "a".

a < 22712 domeniul argumentului fără depășire

a < -22712 avem o subdepășire rezultatul = 0

a 0 domeniul rezultatului

PROCEDURE fabsm(VAR a : mfloat); { *** a <– fabs(a) *** }

Această procedură returnează valoarea absolută a numărului "a". Rezultatul este depus în variabila "a".

PROCEDURE floorm(VAR a : mfloat); { *** a <– floor(a) *** }

Cel mai mare număr intreg care e mai mic sau egal cu numărul în virgulă mobilă "a", este calculat și returnat variabilei "a".

PROCEDURE fmodm(VAR a, b : mfloat); { *** a <– fmod(a,b) *** }

Aceasta este o funcție modul a unui punct mobil. Numărul "a" este împărțit la "b", așa că rezultatul este un întreg. Rezultatul împărțirii este returnat variabilei "a". Variabila "b" nu se schimbă.

a = k * b + c;

"k" este un întreg;

0 c / b < 1; "c" este atribuit lui "a"

Dacă variabila "b" este zero, indicatorul de eroare este setat. Pentru "b" negativ rezultatul este negativ sau zero. Pentru un "k" foarte mare care nu poate fi reprezentat de un număr în virgulă mobilă , rezultatul este zero.

PROCEDURE frexpm(VAR a : mfloat; { *** a <– frexp(a,b) *** }

VAR b : integer);

Numărul mfloat "a" este împărțit în exponenții săi și mantisă:

a = c * (2^b);

0.5 |c| < 1 dacă a ≠ 0;

c este atribuit lui a;

Această procedură returnează două rezultate. Dacă numărul mfloat este zero, rămîne zero și "b" este setat pe zero.

PROCEDURE hypotm(VAR a, b : mfloat); { *** a <– hypot(a,b) *** }

Această funcție este folosită pentru a calcula valoarea absolută a numerelor complexe:

c = sqrt(sqr(a) + sqr(b));

"c" este atribuit lui "a",

"b" rămîne neschimbat.

|a + i*b| = hypot(a,b);

Calculează rezultatul corect chear dacă a pătrat nu poate fi reprezentat, dacă este prea mare sau prea mic.

PROCEDURE ldexpm(VAR a : mfloat; b : integer); { *** a <– a * 2^b *** }

Numărul mfloat "a" este multiplicat de puterile lui doi: 2^b.

a a * (2^b);

Acest calcul este foarte rapid, deoarece doar exponentul lui doi este manipulat.

PROCEDURE logm(VAR a : mfloat); { *** a <- log(a) *** }

Această procedură calculează logaritmul natural ln(x) avînd argumentul "a". Rezultatul este depus în variabila "a".

a > 0 domeniul argumentului

-22712 a 22712 domeniul rezultatului

PROCEDURE modfm(VAR a, b : mfloat); { *** a, b <– modf(a) *** }

Această procedură desparte numărul "a" într-o parte întreagă și o parte fracțională:

a = k + b;

"k" este un intreg;

0 b < 1;

"k" este atribuit lui "a"

PROCEDURE powm(VAR a, b : mfloat); { *** a <– a^b *** }

Această procedură calculează a la puterea b. Rezultatul este depus în "a", variabila "b" nu se schimbă. Dacă exponentul este un întreg, calculul este făcut prin multiplicări repetate. Numărul multiplicărilor este proporțional cu log(|b|). Dacă "b" este fracțional, calculul este făcut cu ajutorul logaritmilor. În acest caz, baza trebuie să fie pozitivă, altfel indicatorul de eroare este setat, pentru că rezultatul nu este real. Dacă a=0 și b0 rezultatul nu este definit și indicatorul de eroare este setat.

PROCEDURE powi(VAR a : mfloat; b : integer); { *** a <– a^b *** }

Această procedură calculează "a" la puterea "b" unde "b" este un exponent întreg. Rezultatul este depus în "a", variabila "b" nu se schimbă. Calculul este făcut prin multiplicări repetate. Numărul multiplicărilor este proporțional cu log(|b|). Dacă a=0 și b0 rezultatul nu este definit și indicatorul de eroare este setat.

PROCEDURE sinm(VAR a : mfloat); { *** a <- sin(a) *** }

Această procedură calculează funcția sin(x) avînd argumentul "a". Rezultatul este depus în variabila "a".

PROCEDURE sqrtm(VAR a : mfloat); { *** a <- sqrt(a) *** }

Această procedură calculează rădăcina pătrată a numărului "a". Rezultatul este depus în variabila "a". Dacă "a" este negativ, indicatorul de eroare este setat.

V.1.3. Funcții standard extinse

PROCEDURE cossinm(VAR a, b : mfloat); { *** a <– cos(a), b <– sin(a) *** }

Această procedură calculează funcția sin(x) și cos(x) prin doar o evaluare a seriilor, folosite pentru funcția sinus. Această procedură este mai rapidă decît două apeluri individuale. Este recomandată de exemplu pentru calcularea funcției exponențiale cu argumentele complexe.

PROCEDURE sqrm(VAR a : mfloat); { *** a <- sqr(a) *** }

Această procedură calculează radical din "a". Rezultatul este depus în variabila "a".

sqr(a) = a*a.

PROCEDURE truncm(VAR a : mfloat); { *** a <– trunc(a) *** }

Această procedură returnează partea întreagă a numărului "a" (trunchiat). Rezultatul este depus în variabila "a". Semnul nu influențează cu nimic și rămîne neschimbat.

VI. Analiza Aplicației

VI.1. Listarea programului unit Bessel.pas

{$N+} {$M 20000, 0, 650000}

unit bessel;

interface

FUNCTION J0(x : extended) : extended;

FUNCTION J0x(x : extended) : extended;

procedure bes_jo;

{Acest program test ne arata avantajele numerelor mfloat}

{Functia bessel JO(x) este calculata cu ajutorul seriilor}

implementation

USES pfloat, crt,graph,grafica;

{––––––––––––––––––-}

FUNCTION J0(x : extended) : extended;

CONST epsi = 1e-20;

VAR sum, sqrx, prod, mepsi : mfloat;

i : integer;

BEGIN

ldtomf(mepsi, epsi);

ldtomf(sqrx, x);

ldexpm(sqrx, -1);

sqrm(sqrx);

Getonem(sum);

Getonem(prod);

i := 0;

REPEAT

i := i + 1;

multm(prod, sqrx);

divi(prod, i);

divi(prod, i);

IF odd(i) THEN

subm(sum, prod)

ELSE

addm(sum, prod);

UNTIL gtm(mepsi, prod);

J0 := mftold(sum);

END;

{––––––––––––––––––-}

FUNCTION J0x(x : extended) : extended;

CONST epsi = 1e-20;

VAR sum, sqrx, prod : extended;

i : integer;

BEGIN

sqrx := sqr(x / 2);

sum := 1;

prod := 1;

i := 0;

REPEAT

i := i + 1;

prod := prod * sqrx / i / i;

IF odd(i) THEN

sum := sum – prod

ELSE

sum := sum + prod

UNTIL epsi > prod;

J0x := sum;

END;

{––––––––––––––––––}

VAR accuracy : integer;

x : extended;

y : integer;

d : string;

procedure bes_jo;

begin

setcolor(black);

outtextxy(xx-153,yy-180,'Calculul Functiei Bessel JO(x)

‘);

setcolor(white);

x := 100;

outtextxy(xx-290,yy-140,'Cuvant Mantisa Rezultat

(x=100.0)');

FOR accuracy := 1 TO 15 DO

BEGIN

Setmantissawords(accuracy);

str(accuracy,d);

outtextxy(xx-250,yy-130+accuracy*10,d);

outtextxy(xx-250,yy-130+accuracy*10,'JO(100.0)=');

str(J0(x),d);

outtextxy(xx,yy-130+accuracy*10,d);

END;

outtextxy(xx-250,yy+30,'Aritmetica IEEE JO(100.0)=');

{este cea implementata de coprocesor}

str(J0x(X),d); {conversie extended la string}

outtextxy(xx,yy+30,d);

outtextxy(xx-250,yy+70,'Este nevoie de o precizie de

12 cuvinte mantisa pentru a obtine');

outtextxy(xx-300,yy+80,' un rezultat corect.');

END;

begin

end.

VI.2. Listarea programului unit PI.pas

{$N+} {$M 20000, 0, 650000}

unit pi;

{** Se prezinta citeva metode pentru calculul lui PI **}

interface

procedure cal_pi;

implementation

USES pfloat, crt, graph, grafica;

VAR a,b,c, bound : mfloat;

i : integer;

procedure cal_pi;

begin

setcolor(white);

Setmantissawords(16);

Getonem(bound);

ldexpm(bound, -16 * Getmantissawords – 16);

setcolor(black);

outtextxy(xx-55,yy-180,'Calculul lui PI');

setcolor(white);

outtextxy(xx-298,yy-140,'PI= 6*arctan(1/5)-

– 4*arctan(1/239) :');

Getonem(a);

multi(a,16);

divi(a,5);

equm(b,a);

i := 0;

REPEAT

i := i + 1;

divi(b, 25);

equm(c, b);

divi(c, 2 * i + 1);

IF odd(i) THEN subm(a, c)

ELSE addm(a, c);

UNTIL gtm(bound, c);

getonem(b);

multi(b, 4);

divi(b, 239);

subm(a, b);

i := 0;

REPEAT

i := i + 1;

divi(b, 239);

divi(b, 239);

equm(c, b);

divi(c, 2 * i + 1);

IF odd(i) THEN addm(a, c)

ELSE subm(a, c);

UNTIL gtm(bound, c);

outtextxy(xx-298,yy-120,mftostr(a,75, '.F77'));

outtextxy(xx-298,yy-80,'PI = 48*arctan(1/18)+

+ 32*arctan(1/57)-20*arctan(1/239):');

getonem(a);

multi(a, 8);

divi(a, 3);

equm(b, a);

i := 0;

REPEAT

i := i + 1;

divi(b, 324);

equm(c, b);

divi(c, 2 * i + 1);

IF odd(i) THEN subm(a, c)

ELSE addm(a, c);

UNTIL gtm(bound, c);

getonem(b);

multi(b, 32);

divi(b, 57);

addm(a, b);

i := 0;

REPEAT

i := i + 1;

divi(b, 3249);

equm(c, b);

divi(c, 2 * i + 1);

IF odd(i) THEN subm(a, c)

ELSE addm(a, c);

UNTIL gtm(bound, c);

getonem(b);

multi(b, 20);

divi(b, 239);

subm(a, b);

i := 0;

REPEAT

i := i + 1;

divi(b, 239);

divi(b, 239);

equm(c, b);

divi(c, 2 * i + 1);

IF odd(i) THEN addm(a, c)

ELSE subm(a, c);

UNTIL gtm(bound, c);

outtextxy(xx-298,yy-60,mftostr(a,75, '.F77'));

outtextxy(xx-298,yy-10,'PI = 6*arctan(1/sqrt(3)):');

Getonem(a);

multi(a, 3);

sqrtm(a);

multi(a, 2);

equm(b, a);

i := 0;

REPEAT

i := i + 1;

divi(b, 3);

equm(c, b);

divi(c, 2 * i + 1);

IF odd(i) THEN subm(a, c)

ELSE addm(a, c);

UNTIL gtm(bound, c);

outtextxy(xx-298,yy+10,mftostr(a,75, '.F77'));

outtextxy(xx-298,yy+50,'Valoarea interna memorata PI:');

Getpim(a);

outtextxy(xx-298,yy+70,mftostr(a,75, '.F77'));

END;

begin

end.

Bibliografie

Donald E. Knuth – Tratat de programare a calculatoarelor. Algoritmi seminumerici, vol.2, Editura Tehnică, București 1983, pp.203-300.

Brent R.P. – Fast Multiple-Precision Evaluațion of Elementary Funcțion.

Beyer W.A. and Waterman M.S. – Error Analiysis of Computațion of Euler's Constant, Mathematics of Computațion, vol.28, Aprile 1974, pp.599-604.

Salamin, E., – Computațion of Pi Using Arithmetic-Geometric Mean, Mathematics of Computațion, vol.30, Iulie 1976, pp.565-570.

Russel Rector and Gorge Alexy – The 8086 Book, McGraw-Hill Inc.1980.

Turbo Pascal 6.0, Ghid de utilizare, Cluj-Napoca 1992.