3. Αριθμητική κινητής υποδιαστολής#

Η αριθμητική κινητής υποδιαστολής χρησιμοποιείται στα ψηφιακά συστήματα για πράξεις πραγματικών αριθμών. Οι αριθμοί κινητής υποδιαστολής ή αριθμοί μηχανής συνήθως καταλαμβάνουν 32 ή 64 bit και απεικονίζονται με δεκαδική ή επιστημονική γραφή όταν η εξαγωγή γίνεται σε κείμενο.

3.1. Επιστημονική σημειογραφεία#

\[ m\textrm{E}n = \underbrace{m}_\text{σημαντικό μέρος} \cdot {\underbrace{10}_\text{βάση}}^{\overbrace{n}^\text{εκθέτης}} \]

  • \(m\) σημαντικό μέρος (mantissa ή significand): δεκαδικός αριθμός, συνήθως στο διάστημα 1 εως 10.

  • \(n\) εκθέτης (exponent ή radix): ακέραια τιμή

  • \(b\) βάση: 10

Στους υπολογιστές έχει επικρατήσει η μορφή με μικρό e ή κεφαλαίο E μπορεί. Σε ορισμένα εργαλεία (Fortran, Matlab) χρησιμοποιείται επιπλέον το d και D.

Table 3.1 Παραδείγματα επιστημονικής σημειογραφείας#

Δεκαδική μορφή

Επιστημονική μορφή

Επιστημονική μορφή (υπολογιστές)

2

\(2×10^0\)

2e0

300

\(3×10^2\)

3e2

4321.768

\(4.321768×10^3\)

4.321768e3

−53000

\(−5.3×10^4\)

-5.3e4

6720000000

\(6.72×10^9\)

6.72e9

0.2

\(2×10^{−1}\)

2e-1

987

\(9.87×10^2\)

9.87e2

0.00000000751

\(7.51×10^{−9}\)

7.51e-9

3.2. Εσωτερική αναπαράσταση#

Κατά την εξαγωγή αριθμών μηχανής χρησιμοποιείται το δεκαδικό σύστημα και εμφανίζεται ως βάση (base) το 10. Εσωτερικά όμως για την αποθήκευση στην μνήμη χρησιμοποιείται το δυαδικό σύστημα και η βάση του 2. Ο αριθμός κινητής υποδιαστολής διπλής ακρίβειας καταλαμβάνει 64 bit, τα οποία κατανέμονται ως εξής:

  • 1 bit στο πρόσημο (sign),

  • 11 bit στον εκθέτη (exponent),

  • 52 bit στο κλασματικό μέρος (fraction).

https://upload.wikimedia.org/wikipedia/commons/thumb/a/a9/IEEE_754_Double_Floating_Point_Format.svg/618px-IEEE_754_Double_Floating_Point_Format.svg.png

Fig. 3.1 Αναπαράσταση στην μνήμη ενός αριθμού μηχανής διπλής ακρίβειας σύμφωνα με το πρότυπο IEEE 754. (Πηγή: Wikipedia Commons)#

Η εσωτερική αναπαράσταση έχει την μορφή:

\[(-1)^{\overbrace{sign}^\text{πρόσημο}}\underbrace{(1.\overbrace{b_{51}b_{50}\dots b_0}^\text{κλασματικό μέρος})_2}_\text{σημαντικό μέρος} \times {\underbrace{2}_\text{βάση}}^{\overbrace{e}^\text{εκθέτης}-1023} \]

Στο δεκαδικό σύστημα το σημαντικό μέρος παίρνει συνήθως τιμές στο πεδίο \([1,10)\), ενώ κατ’ αναλογία στο δυαδικό σύστημα κυμαίνεται στο πεδίο \([1,2)\). Αυτό συμβαίνει πάντα στις κανονικοποιημένες τιμές στις οποίες το πρώτο ψηφίο είναι υποχρεωτικά 1. Γι’αυτό, από τον τον αριθμό \(1.F\) δηλώνεται ρητά μόνο το κλασματικό μέρος (\(F\)) εξοικονομώντας ένα bit (άρρητη μονάδα). Εξαίρεση αποτελούν οι ειδικές τιμές, οι οποίες όμως διαχωρίζονται με τις ακραίες τιμές του εκθέτη \(e\). Σημειώνεται ότι στην δύναμη περιλαμβάνει ένα όρος μετατόπισης (offset) 1023, έτσι ώστε το \(e\) να είναι θετικός ακέραιος.

Περισσότερες λεπτομέρειες δίνονται στο πρότυπο ΙΕΕΕ 754 και στις αναθεωρήσεις του.

3.3. Ειδικές τιμές#

Οι ακραίες τιμές του εκθέτη χρησιμοποιούνται για τον ορισμό ειδικών τιμών.

Table 3.2 Ειδικές τιμές αριθμητικής κινητής υποδιαστολής#

Ειδική τιμή

Συντομογραφία

Εκθέτης (\(e\))

Κλασματικό μέρος (\(F\))

Σημαντικό μέρος

Signed zero

\(\pm 0\)

\(00000000000_2\)

\(0\)

\(0.F\)

Subnormal numbers

\(-\)

\(00000000000_2\)

\(\neq 0\)

\(0.F\)

Infinite

\(\pm \text{Inf}\)

\(11111111111_2\)

\(0\)

\(0.0\)

Not a Number

\(\text{NaN}\)

\(11111111111_2\)

\(\neq 0\)

\(-\)

Οι μη κανονικοποιημένες τιμές (subnormal) έχουν μηδέν στην θέση της άρρητης μονάδας και χρησιμοποιούνται για την αναπαράσταση πολύ μικρών αριθμών με μειωμένο όμως αριθμό σημαντικών ψηφίων. Αν αυτές οι τιμές κανονικοποιηθούν, δίνουν τιμή μηδέν. Ο χειρισμός των μη κανονικοποιημένων τιμών ελέγχεται από τις παραμέτρους του μεταγλωττιστή (compiler).

Note

  • οι πράξεις με \(\infty\) δίνουν το αναμενόμενο μαθηματικό αποτέλεσμα.

  • Οι πράξεις \(\frac{0}0, \frac{\infty}{\infty},\infty - \infty,0 \times \infty\) έχουν απροσδιόριστο αποτέλεσμα (\(\text{NaN}\)).

  • Οι δυνάμεις \( 0^0, 1^{\infty}, \infty^0\) μαθηματικά δίνουν απροσδιόριστο, αλλά στους υπολογιστές δίνουν 1.

  • Οποιαδήποτε πράξη με \(\text{NaN}\) διαδίδει το \(\text{NaN}\).

Tip

Στις διαιρέσεις πρέπει διαχειρίζεστε κατάλληλα την περίπτωση μηδενικού παρονομαστή.

import numpy as np

print(np.double(1e400))  # Overflow (no error reported)
print(np.double(1e-400))  # Underflow (no error reported)

print(np.double(-0.0))  # Signed zero (negative)
print(1 / np.double(-0.0))  # Minus infinity (warning).
print(np.double(0.0) / np.double(0.0))  # NaN (warning)
print(np.inf * np.double(0.0))  # NaN (warning)

print(np.isnan(np.NaN + 1.))  # NaN propagation

inf
0.0
-0.0
-inf
nan
nan
True

/tmp/ipykernel_321007/947857537.py:7: RuntimeWarning: divide by zero encountered in scalar divide
  print(1 / np.double(-0.0))  # Minus infinity (warning).
/tmp/ipykernel_321007/947857537.py:8: RuntimeWarning: invalid value encountered in scalar divide
  print(np.double(0.0) / np.double(0.0))  # NaN (warning)
/tmp/ipykernel_321007/947857537.py:9: RuntimeWarning: invalid value encountered in scalar multiply
  print(np.inf * np.double(0.0))  # NaN (warning)

3.4. Όρια αριθμών μηχανής#

Τα όρια των αριθμών μηχανής που χρησιμοποιεί ένα ψηφιακό σύστημα μπορούν να γίνουν γνωστά μέσω συναρτήσεων της Python.

import sys
import numpy as np

print(sys.float_info)

a = np.double(1.0)

print(np.finfo(a))

print(np.finfo(a).tiny)
print(np.spacing(a))
print(np.spacing(a * 10))

sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308, min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)
Machine parameters for float64
---------------------------------------------------------------
precision =  15   resolution = 1.0000000000000001e-15
machep =    -52   eps =        2.2204460492503131e-16
negep =     -53   epsneg =     1.1102230246251565e-16
minexp =  -1022   tiny =       2.2250738585072014e-308
maxexp =   1024   max =        1.7976931348623157e+308
nexp =       11   min =        -max
smallest_normal = 2.2250738585072014e-308   smallest_subnormal = 4.9406564584124654e-324
---------------------------------------------------------------

2.2250738585072014e-308
2.220446049250313e-16
1.7763568394002505e-15

Tip

Οι τιμές στο διάστημα [4.9406564584124654e-324,2.2250738585072014e-308) είναι μη κανονικοποιημένες, δηλαδή ξεκινάνε με μηδενικά και όχι την άρρητη μονάδα. Γι’ αυτό το λόγο έχουν λιγότερα σημαντικά ψηφία και μικρότερη ακρίβεια.

Σύμβολο

Παράμετρος

Περιγραφή

Τιμή

max

μέγιστη τιμή στο δεκαδικό σύστημα

1.7976931348623157e+308

max_exp

μέγιστη δύναμη με βάση το 2

1024

max_10_exp

μέγιστη δύναμη με βάση το 10 (\(2^{1024} \approx 2 \cdot 10^{308}\))

308

min

ελάχιστη κανονικοποιημένη τιμή στο δεκαδικό σύστημα

2.2250738585072014e-308

min_exp

ελάχιστη δύναμη με βάση το 2

-1021

min_10_exp

ελάχιστη δύναμη με βάση το 10

-307

\(t\)

dig

σημαντικά ψηφία στο δεκαδικό σύστημα

15

\(t\)

mant_dig

σημαντικά ψηφία στο δυαδικό σύστημα

53

\(ε_M\)

epsilon

έψιλον της μηχανής

2.220446049250313e-16

\(b\)

radix

βάση

2

rounds

σύστημα στρογγυλοποίησης

1

Το έψιλον ή μηδέν της μηχανής (\(ε_M\)) εκφράζει την σχετική απόσταση ανάμεσα σε δύο διαδοχικούς αριθμούς μηχανής. Ένας εναλλακτικός ορισμός του \(ε_M\) είναι το ελάχιστο μέγεθος που πρέπει να προστεθεί στο 1 για να αλλάξει τιμή. Με βάση αυτόν τον ορισμό:

\[ ε_M=b^{1-t}\]

Η συνάρτηση numpy.spacing δίνει την απόλυτη απόσταση ανάμεσα σε έναν αριθμό κινητής υποδιαστολής και τον γειτονικό του.

3.5. Ποσοτικοποίηση των σφαλμάτων#

Για την ποσοτικοποίηση των σφαλμάτων χρησιμοποιούνται δύο προσεγγίσεις:

  • Το απόλυτο σφάλμα

\[Ε=x-x_r\]

  • Το σχετικό σφάλμα

\[ε=\frac{x-x_r}{x_r}\]

όπου:

  • \(x\) η εκτίμηση μίας μεταβλητής και

  • \(x_r\) η πραγματική τιμή της η απλώς μια καλύτερη εκτίμηση.

Κατά τον έλεγχο των σφαλμάτων ορίζεται συνήθως ένα ανώτατο όριο πέρα από το οποίο ο υπολογισμός μας δεν είναι αποδεκτός. Το όριο αυτό ονομάζεται ανοχή \(TOL\) (tolerance) και εφαρμόζεται συνήθως στην σχετική τιμή του σφάλματος.

\[|ε|<TOL_r\]

Η επιλογή του σφάλματος και της ανοχής μπορεί να διαφέρουν κατά περίπτωση. Για παράδειγμα, όταν οι τιμές μας προσεγγίζουν το μηδέν, ενδείκνυται η χρήση του απόλυτου σφάλματος.

Exercise 3.1

Υπολογίστε το απόλυτο και σχετικό σφάλμα της τιμής 3.14 σε σχέση με τον αριθμό π της βιβλιοθήκης numpy. Ελέγξτε αν το σχετικό σφάλμα ικανοποιεί μια ανοχή 0.1%.

import numpy as np

E = 3.14 - np.pi
e = (3.14 - np.pi) / np.pi
TOL=0.1/100

print(E)
print(e)
print(f"{e:.2%}")
print(abs(e)<TOL)

-0.0015926535897929917
-0.0005069573828972128
-0.05%
True

3.6. Σφάλμα στρογγυλοποίησης#

Οι διαδοχικοί αριθμοί κινητής υποδιαστολής δεν έχουν πάντα σταθερή απόσταση μεταξύ τους όπως οι ακέραιοι, καθώς εκτός από το τελευταίο σημαντικό ψηφίο, μπορεί να αλλάξει κι ο εκθέτης. Όπως φαίνεται στο παρακάτω γράφημα, οι αριθμοί μηχανής είναι πιο πυκνοί σε μικρές απόλυτες τιμές και αραιώνουν σε μεγαλύτερες.

_images/4e0deea82f089e619fa08c3005f5330ab0aba473ab04eb45282a1fee627b22be.png

Εκτός ιδανικών καταστάσεων, ένας πραγματικός αριθμός βρίσκεται στο ενδιάμεσο διάστημα στο διάστημα μεταξύ δύο διαδοχικών αριθμών μηχανής. Κατά την αποθήκευση του πραγματικού αριθμού ως αριθμού μηχανής ο υπολογιστής καλείται να επιλέξει μία από τις δύο τιμές εφαρμόζοντας:

  • Τεμαχισμό, δηλαδή να τεμαχίσει τον πραγματικό αριθμό απορρίπτοντας όσα ψηφία δεν χωράνε στην εσωτερική αναπαράσταση με σχετικό σφάλμα:

    \[ |ε|\leq b^{1-t} \]

  • Στρογγυλοποίηση, δηλαδή να στρογγυλοποιήσει τον πραγματικό αριθμό στον κοντινότερο αριθμό μηχανής με σχετικό σφάλμα:

    \[ |ε|\leq\frac{1}{2}b^{1-t} \]

    Το μέγεθος \(u=\frac{1}{2}b^{1-t}\) καλείται και μοναδιαίο σφάλμα στρογγυλοποίησης.

Στην ακραία περίπτωση της ισότητας το σφάλμα της στρογγυλοποίησης είναι το μισό του σφάλματος τεμαχισμού και γι’ αυτό εφαρμόζεται σχεδόν πάντα στους υπολογιστές. Στο υπόλοιπο μέρος του συγγράμματος θα αναφερθούμε αποκλειστικά στο σφάλμα στρογγυλοποίησης (round-off error). Σημειώνεται ότι στρογγυλοποίηση εφαρμόζεται και μετά από οποιαδήποτε μαθηματική πράξη στον υπολογιστή.

print(10 - 9.99 == 0.01)
print(10 - 9.99)

False
0.009999999999999787

Το σφάλμα στρογγυλοποίησης μεγενθύνεται μετά από διαδοχικούς υπολογισμούς.

x = 1
dx = 1 / 3
n = 100000
for i in range(n):
    x = x - dx
for i in range(n):
    x = x + dx
print(x)

0.9999999999980784

Λόγω του σφάλματος στρογγυλοποίησης συνήθως δεν έχει νόημα να ελέγχουμε την ισότητα δύο αριθμών μηχανής, αλλά το σφάλμα, σχετικό ή/και απόλυτο. Για αυτή την περίπτωση η numpy παρέχει την συνάρτηση isclose με προεπιλεγμένα όρια rtol=1e-05, atol=1e-08.

print(10 - 9.99 == 0.01)
print(np.isclose(10 - 9.99, 0.01, rtol=1e-07))

False
True

Αντίστοιχα για τον έλεγχο ισότητας δύο πινάκων (στοιχείο προς στοιχείο) χρησιμοποιείται η συνάρτηση allclose.

3.7. Πηγές σφαλμάτων#

Τόσο στους αναλυτικούς υπολογισμούς, όσο και στους υπολογισμούς με αριθμούς μηχανής υπεισέρχονται σφάλματα. Αν και γενικώς προσπαθούμε να τα περιορίσουμε, αυτό δεν είναι απόλυτα εφικτό. Oι πηγές σφαλμάτων στους αριθμητικούς υπολογισμούς είναι:

  • Η υπερχείλιση (overflow).

    • Αποτελέσματα μεγαλύτερα από το άνω όριο παίρνουν την ειδική τιμή Inf.

  • Η υποχείλιση (underflow).

    • Αποτελέσματα μικρότερα από το κάτω όριο παίρνουν την τιμή 0.

  • Η στρογγυλοποίηση

    • Όλα τα αποτελέσματα πράξεων στρογγυλοποιούνται στους αντίστοιχους αριθμούς κινητής υποδιαστολής.

  • Η αποκοπή

    • Σφάλματα που εισάγει η αριθμητική μέθοδος λόγω της χρήσης πεπερασμένων όρων ή αντί άπειρων.

Tip

Σε πολλά συγγράμματα ο όρος σφάλμα αποκοπής (truncation error) χρησιμοποιείται ως συνώνυμο του σφάλματος τεμαχισμού (chopping error). Το παρόν σύγγραμμα υιοθετεί την ορολογία των Chapra και Canale [2], που διαχωρίζει τις δύο έννοιες.

Στα πλαίσια της αριθμητικής ανάλυσης είναι απαραίτητη η ποσοτικοποίηση των σφαλμάτων και η κατανόηση των παραγόντων που τα επηρεάζουν.

Οι πρώτες 3 πηγές σφαλμάτων μπορούν να περιοριστούν χρησιμοποιώντας αριθμούς κινητής υποδιαστολής υψηλότερης ακρίβειας, π.χ. διπλής ή αντί για μονής, καθώς επίσης και με την αναδιατύπωση ορισμένων αναλυτικών σχέσεων.

Η τελευταία πηγή σφαλμάτων μπορεί να περιοριστεί με την κατάλληλη επιλογή της αριθμητικής μεθόδου και των παραμέτρων της, όπως θα δούμε στα παρακάτω κεφάλαια.

Exercise 3.2

Υπολογίστε τις ρίζες της συνάρτησης \(f(x)=x^2-100000.1x+10000\).

import numpy as np

a = 1
b = -100000.1
c = 10000

x1 = (-b+np.sqrt(b**2-4*a*c))/(2*a)
x2 = (-b-np.sqrt(b**2-4*a*c))/(2*a)
# x1=(2*c)/(-b-np.sqrt(b**2-4*a*c))
# x2=(2*c)/(-b+np.sqrt(b**2-4*a*c))

print(x1, x2)

100000.0 0.10000000000582077

Ακόμη κι όταν υπάρχουν αναλυτικές λύσεις χρειάζεται προσοχή στην κωδικοποίηση του υπολογισμού. Στην παραπάνω άσκηση οι ρίζες του τριωνύμου είναι γνωστές 100000 και 0.1. Όπως φαίνεται από τα αποτελέσματα, ο υπολογισμός της δεύτερης ρίζας αποκλίνει ελαφρώς και αυτό οφείλεται στο σφάλμα στρογγυλοποίησης κατά την αφαίρεση στον αριθμητή. Ένας τρόπος να περιοριστεί το σφάλμα είναι να χρησιμοποιηθεί η εναλλακτική σχέση:

\[ x_{1,2}=\frac{2c}{-b\mp \sqrt{b^2-4ac}} \]

Η πρώτη ρίζα υπολογίζεται με \(-\) και η δεύτερη με \(+\).

3.8. Ταχύτητα υπολογισμού#

Exercise 3.3

Χρησιμοποιώντας την μονάδα κώδικα timeit, μετρήστε τον χρόνο που χρειάζονται οι 4 βασικοί αριθμητικοί τελεστές, η δύναμη και η συνάρτηση exp της numpy και αδιαστατοποιήστε τα αποτελέσματα σε σχέση με τον χρόνο της πρόσθεσης. Προσέξτε ότι τα αποτελέσματα διαφοροποιούνται ελαφρώς αναλόγως τις αριθμητικές τιμές που θα επιλεγούν για τις πράξεις.

import timeit
import numpy as np

a = np.arange(0, 20, 0.1)
b = 2.5
iterations = 100000

dt_addition=timeit.timeit(stmt="c=a+b", number=iterations, globals=globals())
dt_addition=timeit.timeit(stmt="c=a+b", number=iterations, globals=globals())   #first call is always longer. Second call is more accurate

print("Πρόσθεση\t",1.0)
print("Αφαίρεση\t",timeit.timeit(stmt="c=a-b", number=iterations, globals=globals())/dt_addition)
print("Πολλαπλασιασμός\t",timeit.timeit(stmt="c=a*b", number=iterations, globals=globals())/dt_addition)
print("Διαίρεση\t",timeit.timeit(stmt="c=a/b", number=iterations, globals=globals())/dt_addition)
print("Δύναμη\t\t",timeit.timeit(stmt="c=a**b", number=iterations, globals=globals())/dt_addition)
print("Εκθετική συν.\t",timeit.timeit(stmt="c=np.exp(a)", number=iterations, globals=globals())/dt_addition)

Πρόσθεση	 1.0

Αφαίρεση	 0.9779591186230746
Πολλαπλασιασμός	 0.9736175814420275

Διαίρεση	 1.0407038819739813

Δύναμη		 4.858316690286215
Εκθετική συν.	 1.6874739281396394

Ο αριθμός flop/s ή FLOPS ορίζεται ως το πλήθος των πράξεων κινητής υποδιαστολής που εκτελεί ένας υπολογιστής στην μονάδα του χρόνου, συνήθως σε ένα δευτερόλεπτο. Η μέτρηση των FLOPS μπορεί να διαφοροποιηθεί σημαντικά ανάλογα το τρόπο μέτρησης για αυτό θέλει προσοχή όταν αναφέρονται τέτοιες τιμές. Οι υπολογιστές υψηλών επιδόσεων έχουν ξεπεράσει το \(1Εflop/s=1\times 10^{18} flop/s\) σύμφωνα με την λίστα που διατηρεί ο δικτυακός τόπος TOP500.

Exercise 3.4

Ένας απλός τρόπος εκτίμησης του αριθμού FLOPS συνίσταται στην μέτρηση του χρόνου που χρειάζεται ένας υπολογιστής για να εκτελέσει έναν πολλαπλασιασμό και μία πρόσθεση μαζί (multiply-and-add instruction), δηλαδή \(y←a+x*y\), μαζί με τον χρόνο που απαιτείται για την ανάκτηση από την μνήμη των δεδομένων που εμπλέκονται στις δύο αυτές πράξεις. Εκτιμήστε τον αριθμό FLOPS του υπολογιστή σας με βάση τον παραπάνω ορισμό.

import timeit
import numpy as np

def FLOPfunction(iterations:int):
    x = 1+1e-10
    y = 1.
    a = 0.5
    for i in range(iterations):
        y=a+x*y 
    #print(y)   #Make sure y does not reach infinity

iterations = 10000000
dt = timeit.timeit(stmt="FLOPfunction(iterations)", number=1, globals=globals())
FLOPS = 2*iterations / dt   #MAD instruction counts as 2 iterations
print(f"{FLOPS*1e-6:.2f} MFLOPS")

46.91 MFLOPS