# 15. Αριθμητική Κινητής Υποδιαστολής: Ζητήματα και Περιορισμοί

Οι αριθμοί κινητής υποδιαστολής αντιπροσωπεύονται στο υλικό του υπολογιστή ως κλάσματα με βάση το 2 (δυαδικά). Για παράδειγμα, το δεκαδικό κλάσμα `0.625` έχει τιμή 6/10 + 2/100 + 5/1000, και με τον ίδιο τρόπο το δυαδικό κλάσμα `0.101` έχει τιμή 1/2 + 0/4 + 1/8. Αυτά τα δύο κλάσματα έχουν πανομοιότυπες τιμές, η μόνη πραγματική διαφορά είναι ότι το πρώτο είναι γραμμένο με κλασματικό συμβολισμό με βάση το 10, και το δεύτερο με βάση το 2.

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

Το πρόβλημα είναι πιο κατανοητό στην αρχή με βάση το 10. Θεωρήστε το κλάσμα 1/3. Μπορεί να το προσεγγίσετε ως κλάσμα βάσης 10:

    0.3

ή, καλύτερα,

    0.33

ή, καλύτερα,

    0.333

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

Με τον ίδιο τρόπο, ανεξάρτητα από το πόσα ψηφία βάσης 2 είστε διατεθειμένοι να χρησιμοποιήσετε, η δεκαδική τιμή 0,1 δεν μπορεί να αναπαρασταθεί ακριβώς ως κλάσμα βάση 2. Στη βάση 2, το 1/10 είναι το κλάσμα που επαναλαμβάνεται

    0.0001100110011001100110011001100110011001100110011...

Σταματήστε σε οποιονδήποτε πεπερασμένο αριθμό bit και λαμβάνετε μια προσέγγιση. Στις περισσότερες μηχανές σήμερα, οι floats προσεγγίζονται χρησιμοποιώντας τα πρώτα 53 bit ξεκινώντας από το πιο σημαντικό bit και με τον παρανομαστή ως δύναμη του δύο. Στην περίπτωση του 1/10, το δυαδικό κλάσμα είναι `3602879701896397 / 2 ** 55` που είναι κοντά αλλά όχι ακριβώς ίσο με την πραγματική τιμή του 1/10.

Πολλοί χρήστες δεν γνωρίζουν την προσέγγιση λόγω του τρόπου με τον οποίο εμφανίζονται οι τιμές. Η Python εκτυπώνει μόνο μια δεκαδική προσέγγιση στην πραγματική δεκαδική τιμή της δυαδικής προσέγγισης που είναι αποθηκευμένη από το μηχάνημα. Στα περισσότερα μηχανήματα, αν η Python τύπωνε την πραγματική δεκαδική τιμή της δυαδικής προσέγγισης που είναι αποθηκευμένη για 0,1, θα πρέπει να εμφανίζει:

    >>> 0.1
    0.1000000000000000055511151231257827021181583404541015625

Αυτά είναι περισσότερα ψηφία από όσα βρίσκουν χρήσιμα οι περισσότεροι, επομένως η Python διατηρεί τον αριθμό των διαχειρίσιμων ψηφίων εμφανίζοντας μια στρογγυλεμένη τιμή αντ” αυτού:

    >>> 1 / 10
    0.1

Απλώς θυμηθείτε, παρόλο που το εκτυπωμένο αποτέλεσμα μοιάζει με την ακριβή τιμή του 1/10, η πραγματική αποθηκευμένη τιμή είναι το πλησιέστερο αναπαραστάσιμο δυαδικό κλάσμα.

Είναι ενδιαφέρον ότι υπάρχουν πολλοί διαφορετικοί δεκαδικοί αριθμοί που μοιράζονται το ίδιο πλησιέστερο κατά προσέγγιση δυαδικό κλάσμα. Για παράδειγμα, οι αριθμοί `0.1` και `0.1000000000000000055511151231257827021181583404541015625` είναι όλα κατά προσέγγιση με `3602879701896397 / 2 ** 55`. Δεδομένου ότι όλες αυτές οι δεκαδικές τιμές μοιράζονται την ίδια προσέγγιση, οποιαδήποτε από αυτές θα μπορούσε να εμφανιστεί διατηρώντας παράλληλα το αμετάβλητο `eval(repr(x)) == x`.

Ιστορικά, το prompt της Python και η ενσωματωμένη συνάρτηση `repr()` θα επέλεγε αυτό με 17 σημαντικά ψηφία, `0.10000000000000001`. Ξεκινώντας με την Python 3.1, η Python (στα περισσότερα συστήματα) είναι πλέον σε θέση να επιλέξει το συντομότερο από αυτά και απλά εμφανίζει το `0.1`.

Λάβετε υπόψη ότι αυτό είναι στην ίδια τη φύση του δυαδικού κινητής υποδιαστολής: αυτό δεν είναι σφάλμα στην Python, ούτε είναι σφάλμα στον κώδικα σας. Θα δείτε το ίδιο πράγμα σε όλες τις γλώσσες που υποστηρίζουν την αριθμητική κινητής υποδιαστολής του υλικού σας (αν και ορισμένες γλώσσες μπορεί να μην εμφανίζουν τη διαφορά από προεπιλογή ή σε όλες τις λειτουργίες εξόδου).

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

    >>> format(math.pi, '.12g')  # give 12 significant digits
    '3.14159265359'

    >>> format(math.pi, '.2f')   # give 2 digits after the point
    '3.14'

    >>> repr(math.pi)
    '3.141592653589793'

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

Μια ψευδαίσθηση μπορεί να γεννήσει μια άλλη. Για παράδειγμα, καθώς το 0,1 δεν είναι ακριβώς το 1/10, το άθροισμα τριών τιμών του 0,1 μπορεί να μην αποφέρει ακριβώς 0,3:

    >>> 0.1 + 0.1 + 0.1 == 0.3
    False

Επίσης, εφόσον το 0.1 δεν μπορεί να πλησιάσει την ακριβή τιμή του 1/10 και το 0.3 δεν μπορεί να πλησιάσει την ακριβή τιμή του 3/10, τότε η προ-στρογγυλοποίηση με τη συνάρτηση `round()` δεν μπορεί να βοηθήσει:

    >>> round(0.1, 1) + round(0.1, 1) + round(0.1, 1) == round(0.3, 1)
    False

Αν και οι αριθμοί δεν μπορούν να γίνουν πιο κοντά στις προβλεπόμενες ακριβείς τιμές τους, η συνάρτηση `math.isclose()` μπορεί να είναι χρήσιμη για τη σύγκριση ανακριβών τιμών:

    >>> math.isclose(0.1 + 0.1 + 0.1, 0.3)
    True

Εναλλακτικά, η συνάρτηση `round()` μπορεί να χρησιμοποιηθεί για τη σύγκριση χονδρικών προσεγγίσεων:

    >>> round(math.pi, ndigits=2) == round(22 / 7, ndigits=2)
    True

Η αριθμητική δυαδικής κινητής υποδιαστολής επιφυλάσσει πολλές εκπλήξεις όπως αυτή. Το πρόβλημα με το «0.1» εξηγείται λεπτομερώς παρακάτω, στην ενότητα «Σφάλμα Αναπαράστασης». Δείτε Παραδείγματα προβλημάτων κινητής υποδιαστολής για μια ευχάριστη περίληψη του τρόπου λειτουργίας της δυαδικής κινητής υποδιαστολής και των ειδών προβλημάτων που αντιμετωπίζονται συνήθως στην πράξη. Δείτε επίσης Οι κίνδυνοι του Floating Point <http://www.indowsway.com/floatingpoint.htm>`_ για μια πιο ολοκληρωμένη περιγραφή άλλων κοινών εκπλήξεων.

Όπως λέει αυτό κοντά στο τέλος, «δεν υπάρχουν εύκολες απαντήσεις.» Ωστόσο, μην είστε υπερβολικά επιφυλακτικοί με την κινητή υποδιαστολή! Τα σφάλματα στις λειτουργίες κινητής υποδιαστολής Python κληρονομούνται από το υλικό κινητής υποδιαστολής και στα περισσότερα μηχανήματα δεν είναι της τάξης του 1 σε 2**53 ανά πράξη. Αυτό είναι παραπάνω από επαρκές για τις περισσότερες εργασίες, αλλά πρέπει να έχετε κατά νου ότι δεν είναι δεκαδική αριθμητική και ότι κάθε λειτουργία float μπορεί να υποστεί νέο σφάλμα στρογγυλοποίησης.

Ενώ υπάρχουν παθολογικές περιπτώσεις, για την πιο περιστασιακή χρήση της αριθμητικής κινητής υποδιαστολής, θα δείτε στο τέλος το αποτέλεσμα που περιμένετε εάν απλώς στρογγυλοποιήσετε την εμφάνιση των τελικών αποτελεσμάτων σας στον αριθμό των δεκαδικών ψηφίων που περιμένετε. Το `str()` συνήθως αρκεί, και για καλύτερο έλεγχο δείτε τους προσδιοριστές μορφής της μεθόδου `str.format()` σε Format string syntax.

Για περιπτώσεις χρήσης που απαιτούν ακριβή δεκαδική αναπαράσταση, δοκιμάστε να χρησιμοποιήσετε το module `decimal` που εφαρμόζει δεκαδική αριθμητική κατάλληλη για λογιστικές εφαρμογές και εφαρμογές υψηλής ακριβείας.

Μια άλλη μορφή ακριβούς αριθμητικής υποστηρίζεται από το module `fractions`, η οποία υλοποιεί την αριθμητική με βάση τους ορθολογικούς αριθμούς (έτσι οι αριθμοί όπως το 1/3 μπορούν να αναπαρασταθούν ακριβώς).

Εάν είστε ένα εντατικός χρήστης πράξεων κινητής υποδιαστολής, θα πρέπει να ρίξετε μια ματιά στο πακέτο NumPy και πολλά άλλα πακέτα για μαθηματικές και στατιστικές πράξεις που παρέχονται από το project SciPy. Δείτε <https://scipy.org>.

Η Python παρέχει εργαλεία που μπορεί να βοηθήσουν σε εκείνες τις σπάνιες περιπτώσεις που πραγματικά θέλετε να μάθετε την ακριβή τιμή ενός float. Η μέθοδος `float.as_integer_ratio()` εκφράζει την τιμή ενός float ως κλάσμα:

    >>> x = 3.14159
    >>> x.as_integer_ratio()
    (3537115888337719, 1125899906842624)

Δεδομένου ότι η αναλογία είναι ακριβής, μπορεί να χρησιμοποιηθεί για την αναδημιουργία χωρίς απώλειες της αρχικής τιμής:

    >>> x == 3537115888337719 / 1125899906842624
    True

Η μέθοδος `float.hex()` εκφράζει ένα float σε δεκαεξαδικό (βάση 16), δίνοντας πάλι την ακριβή τιμή που έχει αποθηκευτεί στον υπολογιστή σας:

    >>> x.hex()
    '0x1.921f9f01b866ep+1'

Αυτή η ακριβής δεκαεξαδική αναπαράσταση μπορεί να χρησιμοποιηθεί για την ανακατασκευή της τιμής float ακριβώς:

    >>> x == float.fromhex('0x1.921f9f01b866ep+1')
    True

Δεδομένου ότι η αναπαράσταση είναι ακριβής, είναι χρήσιμη για την αξιόπιστη μεταφορά τιμών σε διαφορετικές εκδόσεις της Python (ανεξαρτησία πλατφόρμας) και την ανταλλαγή δεδομένων με άλλες γλώσσες που υποστηρίζουν την ίδια μορφή (όπως Java και C99).

Ένα άλλο χρήσιμο εργαλείο είναι η συνάρτηση `sum()` που βοηθά στον μετριασμό της απώλειας ακρίβειας κατά την άθροιση. Χρησιμοποιεί εκτεταμένη ακρίβεια για ενδιάμεσα βήματα στρογγυλοποίησης καθώς οι τιμές προστίθενται σε ένα τρέχον σύνολο. Αυτό μπορεί να κάνει τη διαφορά στη συνολική ακρίβεια, ώστε τα σφάλματα να μην συσσωρεύονται στο σημείο που επηρεάζουν το τελικό σύνολο:

    >>> 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 == 1.0
    False
    >>> sum([0.1] * 10) == 1.0
    True

Το `math.fsum()` πηγαίνει πιο μακριά και παρακολουθεί όλα τα «χαμένα ψηφία» καθώς οι τιμές προστίθενται σε ένα τρέχον σύνολο, έτσι ώστε το αποτέλεσμα να έχει μόνο μία στρογγυλοποίηση. Αυτό είναι πιο αργό από το `sum()` αλλά θα είναι πιο ακριβές σε ασυνήθιστες περιπτώσεις όπου οι είσοδοι μεγάλου μεγέθους ακυρώνουν η μία την άλλη αφήνοντας ένα τελικό άθροισμα κοντά στο μηδέν:

    >>> arr = [-0.10430216751806065, -266310978.67179024, 143401161448607.16,
    ...        -143401161400469.7, 266262841.31058735, -0.003244936839808227]
    >>> float(sum(map(Fraction, arr)))   # Exact summation with single rounding
    8.042173697819788e-13
    >>> math.fsum(arr)                   # Single rounding
    8.042173697819788e-13
    >>> sum(arr)                         # Multiple roundings in extended precision
    8.042178034628478e-13
    >>> total = 0.0
    >>> for x in arr:
    ...     total += x                   # Multiple roundings in standard precision
    ...
    >>> total                            # Straight addition has no correct digits!
    -0.0051575902860057365

## 15.1. Σφάλμα Αναπαράστασης

Αυτή η ενότητα εξηγεί το παράδειγμα «0.1» λεπτομερώς και δείχνει πώς μπορείτε να εκτελέσετε μια ακριβή ανάλυση περιπτώσεων όπως αυτή μόνοι σας. Υποτίθεται ότι έχετε βασική εξοικείωση με την αναπαράσταση δυαδικής κινητής υποδιαστολής.

Το Σφάλμα αναπαράστασης (Representation error) αναφέρεται στο γεγονός ότι ορισμένα (τα περισσότερα, στην πραγματικότητα) δεκαδικά κλάσματα δεν μπορούν να αναπαρασταθούν ακριβώς ως δυαδικά (βάση 2) κλάσματα. Αυτός είναι ο κύριος λόγος για τον οποίο η Python (ή Perl, C, C++, Java, Fortran, και πολλές άλλες) συχνά δεν εμφανίζουν τον ακριβή δεκαδικό αριθμό που περιμένετε.

Γιατί συμβαίνει αυτό; Το 1/10 δεν μπορεί να αναπαρασταθεί ακριβώς ως δυαδικό κλάσμα. Από το 2000 τουλάχιστον, σχεδόν όλες οι μηχανές χρησιμοποιούν δυαδική αριθμητική κινητής υποδιαστολής IEEE 754 και σχεδόν όλες οι πλατφόρμες αντιστοιχίζουν τα Python floats σε IEEE binary64 «διπλής ακρίβειας» τιμές. Οι τιμές του IEEE 754 binary64 περιέχουν 53 bits ακρίβειας, επομένως κατά την είσοδο ο υπολογιστής προσπαθεί να μετατρέψει το 0,1 στο πλησιέστερο κλάσμα που μπορεί να έχει τη μορφή J/2**N όπου J είναι ένας ακέραιος που περιέχει ακριβώς 53 bits. Ξαναγράφεται

    1 / 10 ~= J / (2**N)

ως

    J ~= 2**N / 10

και υπενθυμίζοντας ότι το J έχει ακριβώς 53 bits (είναι `>= 2**52` αλλά `< 2**53`), η καλύτερη τιμή για το N είναι 56:

    >>> 2**52 <=  2**56 // 10  < 2**53
    True

Δηλαδή, το 56 είναι η μόνη τιμή για το N που αφήνει το J με ακριβώς 53 bits. Η καλύτερη δυνατή τιμή για το J είναι τότε αυτό το πηλίκο στρογγυλοποιημένο:

    >>> q, r = divmod(2**56, 10)
    >>> r
    6

Δεδομένου ότι το υπόλοιπο είναι περισσότερο από το μισό του 10, η καλύτερη προσέγγιση επιτυγχάνεται με στρογγυλοποίηση προς τα επάνω:

    >>> q+1
    7205759403792794

Ως εκ τούτου, η καλύτερη δυνατή προσέγγιση στο 1/10 στο IEEE 754 διπλής ακρίβειας είναι:

    7205759403792794 / 2 ** 56

Η διαίρεση του αριθμητή και του παρονομαστή με δύο μειώνει το κλάσμα σε:

    3602879701896397 / 2 ** 55

Λάβετε υπόψη ότι από τη στιγμή που κάναμε στρογγυλοποίηση, αυτό είναι στην πραγματικότητα λίγο μεγαλύτερο από το 1/10· αν δεν είχαμε στρογγυλοποιήσει προς τα πάνω, το πηλίκο θα ήταν λίγο μικρότερο από το 1/10. Αλλά σε καμία περίπτωση δεν μπορεί να είναι ακριβώς 1/10!

Έτσι ο υπολογιστής δεν «βλέπει» ποτέ 1/10: αυτό που βλέπει είναι το ακριβές κλάσμα που δίνεται παραπάνω, η καλύτερη διπλή προσέγγιση IEEE 754 που μπορεί να πάρει:

    >>> 0.1 * 2 ** 55
    3602879701896397.0

Αν πολλαπλασιάσουμε αυτό το κλάσμα με 10**55, μπορούμε να δούμε την τιμή με 55 δεκαδικά ψηφία:

    >>> 3602879701896397 * 10 ** 55 // 2 ** 55
    1000000000000000055511151231257827021181583404541015625

που σημαίνει ότι ο ακριβής αριθμός που είναι αποθηκευμένος στον υπολογιστή είναι ίσος με την δεκαδική τιμή 0.1000000000000000055511151231257827021181583404541015625. Αντί να εμφανιστεί η πλήρης τιμή, πολλές γλώσσες (συμπεριλαμβανομένων των παλαιότερων εκδόσεων της Python), στρογγυλοποιούν το αποτέλεσμα σε 17 σημαντικά ψηφία:

    >>> format(0.1, '.17f')
    '0.10000000000000001'

Τα modules `fractions` και `decimal` κάνουν αυτούς τους υπολογισμούς εύκολους:

    >>> from decimal import Decimal
    >>> from fractions import Fraction

    >>> Fraction.from_float(0.1)
    Fraction(3602879701896397, 36028797018963968)

    >>> (0.1).as_integer_ratio()
    (3602879701896397, 36028797018963968)

    >>> Decimal.from_float(0.1)
    Decimal('0.1000000000000000055511151231257827021181583404541015625')

    >>> format(Decimal.from_float(0.1), '.17')
    '0.10000000000000001'