Δυστυχώς ο υπολογιστής δεν καταλαβαίνει '0.02' αλλα '3F947AE147AE147B'. Και ο αλγόριθμός σου πάει περίπατο...
Round-off errors, ο εφιάλτης του scientific software.
UPDATE : In-depth (τι σκατά συμβαίνει ;)
Τα πράγματα φαίνονται απλά. Συγκρίνω μια σταθερά η οποία έχει οριστεί μια φορά στην αρχή (f_esu) με μια μεταβλητή που υπολογίζεται ξανά και ξανά (Abse).
Το αποτέλεσμα Abse <= f_esu, όταν και οι δύο όροι είναι ίσοι με 0.02, αντί να είναι μόνιμα true, προκύπτει τυχαίο ! Η προφανής υποψία λοιπόν είναι round-off error στην τιμή της Abse, η οποία επαναυπολογίζεται συνεχώς.
Όπως έδειξα και στην αρχή, η εξέταση των τιμών των μεταβλητών μέσω του Watch List δεν βγάζει κανένα νόημα, οπότε είναι η ώρα να μπούν τα μεγάλα μέσα...
Με το CPU Window μπορεί να δει κανείς πως μεταφράζονται οι high level εντολές του προγράμματος σε assembly (κώδικας μηχανής που λέγαν οι παλιοί), και πλέον δεν μπορεί να ξεφύγει τίποτα. Μπαίνει λοιπόν ένα breakpoint στην εντολή όπου γινεται η σύγκριση Abse <= f_esu (δεύτερο σκέλος μετά το and) και βλέπουμε τα εξής :
Σύντομη επεξήγηση
fld qword ptr [ebp-$10]
Βάλε τo qword (64 bit αριθμός κινητής υποδιαστολής) στην FPU από τη διεύθυνση ebp-$10. Η τιμή αυτή αντιστοιχεί στη μεταβλητή Abse του προγράμματος.
fcomp qword ptr [ebx+$58]
Σύγκρινε την παραπάνω τιμή (Abse) με την qword που υπάρχει στη διεύθυνση ebx+$58, η οποία αντιστοιχεί στη σταθερά f_esu του προγράμματος (αυτό μπορεί κανείς να το διαπιστώσει όταν η τιμή 0.02 γίνει assign στην f_esu, όπου και μπαίνει στην συγκεκριμένη διεύθυνση μνήμης).
fstsw ax
sahf
Στείλε το αποτέλεσμα (flags) από την FPU στον register AX της CPU και αντέγραψε τον register AX στα αντίστοιχα flags της CPU. Με λίγα λόγα μετέφερε το αποτέλεσμα από την FPU στην CPU.
jnbe +$37
Αν το αποτέλεσμα δεν είναι μικρότερο ή ίσο (δηλαδή Abse > f_esu), σήκω και φύγε (πήδα 37 addresses), δηλαδή ουσιαστικά παράκαμψε το begin-end block που ακολουθεί στο πρόγραμμα.
Το πρόβλημα λοιπόν εντοπίζεται στην σύγκριση μεταξύ των δύο τιμών qword ptr [ebp-$10] (Abse) και qword ptr [ebx+$58] (f_esu). Για να προκύπτει λάθος, θα πρέπει η μία εκ των δύο να είναι μεγαλύτερη από 0.02 !
Ας πιάσουμε την ύποπτη μεταβλητή Abse. Το παραπάνω screenshot αντιστοιχεί στην (τυχαία) περίπτωση που το αποτέλεσμα βγαίνει σωστό. Εξετάζοντας την θέση μνήμης ebp-$10 παρατηρούμε (πάνω δεξιά) ότι ο EBP register έχει τιμή $0012FA64 , οπότε :
ebp-$10 = $0012FA54
Στο κάτω δεξιά κουτάκι παρατηρούμε ότι στη θέση $0012FA54 αντιστοιχεί η τιμή $47ΑΕ147Β. Επειδή όμως πρόκειται για qword (64 bit), θα πρέπει να κολλήσουμε από αριστερά και το επόμενο dword (32 bit), δηλαδή αυτό στη θέση μνήμης $0012FA58, που είναι ίσο με $3F947AE1. Έχουμε λοιπόν την τιμή :
$3F947AE147AE147B
Βάζοντάς την στο κατάλληλο calculator προκύπτει ότι $3F947AE147AE147B = 0.02 όπως προβλέπεται !
Τί συμβαίνει όμως με την (τυχαία) περίπτωση που η σύγκριση βγαίνει λάθος ; Το CPU Window σε αυτή τη περίπτωση είναι το εξής :
Βλέπετε την διαφορά ; Η τιμή της (υποτιθέμενης) σταθεράς f_esu είναι πλέον :
$3F947AE147AE147C
Για να τη βάλουμε στο calculator...
Αντί λοιπόν για 0.02, έχουμε 0.020000000000000004 !!! Εδώ λοιπόν φαίνεται καθαρά το round-off error που δεν 'επιασε' το Watch List !
Η πλάκα της όλης υπόθεσης είναι ότι αυτό το φαινομενικά ασήμαντο λάθος, αλλοιώνει τραγικά τα αποτελέσματα ολόκληρου του προγράμματος, και για να αντιμετωπιστεί χρειάζεται αριθμητικές αλχημείες, που σιχαίνομαι.
"To err is human, to really foul things up requires a computer"
ποιό IDE είναι αυτό;
ReplyDeleteΜα τι ρωτας... ;) ;)
ReplyDeleteα, και καλό μου παιδί, αφού κάνεις που κάνεις sci-software και θεσ τρελές ακρίβειες,
ReplyDeleteΓΙΑΤΙ ΔΟΥΛΕΥΕΙΣ ΜΕ ΔΕΚΑΔΙΚΟΥΣ;;; τα κολπάκια με το κάθε μέρος του αριθμού σε long int δεν τα ξέρετε εκεί πάνω στην Τούμπα;
Δεν ζήτησα τρελλές ακρίβειες, απλά θέλω (0.02 <= 0.02) = true
ReplyDeleteΓια πες τωρα τα κολπάκια που ξέρεις ;)
στο γράφω πριν διαβάσω το updated post (μπορεί να το εφαρμόζεις εκεί)
ReplyDeleteεγώ όποτε θέλω μεγάλες ακρίβειες αυτό που κάνω είναι να ορίσω μία custom κλάση μου η οποία θα έχει property 2 long integers (εκεί δεν έχεις σφάλματα, μόνο overflows) έναν για το ακέραιο κι έναν για το δεκαδικό μέρος. Επιπλέον φτιάχνεις μεθόδους για τις αριθμητικές τους πράξεις.
Μα τον Δία !!!
ReplyDeleteMore instructions, slower code !
Εμείς εδώ στην Τούμπα (btw έτσουξε προχθές ε ;) γράφουμε γρήγορο κώδικα ρε, δεν περιμένουμε την Intel και την NVidia να μας σπρώχνουνε κάθε 3 κ λίγο όπως εσάς τα αγγλάκια ;) ;)
Απλή ήταν η λύση :
Abse := Abs(e - 1e-15);
Αριθμητική αλχημεία, αλλά κανένα πρόβλημα μετα.
Καλά τα λες, αλλά δεν μπορούσες να βάλεις καλύτερο κώδικα; να τον προσαρμόσεις λίγο στο πρόβλημα βρε παιδί μου.. Ούτε καν η δήλωση της f_esu δεν φαίνεται..
ReplyDeleteΔεν μπορείς να κάνεις round στο τρίτο-τέταρτο δεκαδικό ψηφίο; εγώ αυτό θα έκανα :)
Τώρα για τα κόμπλεξ της Delphi δεν ξέρω.. Στην Java έχουμε και τα BigDecimal για να είμαστε ακριβής..
Η δήλωση της f_esu είναι
ReplyDeletef_esu := 0.02;
ή αν προτιμάς :
mov eax, [ebp+$08]
mov [ebx+$58], eax
mov eax, [ebp+$0c]
mov [ebx+$5c], eax
Αν κάνεις round στο 3ο-4ο δεκαδικό πάει περίπατο όλος ο αλγόριθμος. Στο 14ο-15ο όμως νομίζω ότι θα ήταν καλό. Έχεις δίκιο :) Είναι πιο όμορφο από το να αφαιρείς ένα 1e-15.
BigDecimal ? Διαβάζω καλά ότι πιάνει 24 bytes ?!?!?!
Ειστε memory και CPU abusers τελικά, καλά τα λέω :)
Φίλε δεν ξέρω τι κάνεις εκεί και τι ακρίβεια θέλεις.. Για αυτό στο είπα το BigDecimal ;)
ReplyDeleteΑυτό που μου πρότεινες είναι πολύ καλύτερο :
ReplyDeletefunction Clean(x : double) : double;
begin
Result := Round(x*1e12)/1e12;
end;
Απλά, όμορφα και γρήγορα :)
Giati den kaneis ena round sta 9 -12 dekadika kai tis dyo metablhtes? Nomizw oti h Delphi exei Round() h mhpws oxi?
ReplyDeleteΜόλις έγινε στο προηγούμενο comment ;)
ReplyDeleteΜπορείς να μου πεις ακριβώς την πράξη που κάνει και βγάζει τέτοιο αριθμό?
ReplyDeleteΗ αμέσως προηγούμενη πράξη είναι
ReplyDeleteAbse := Abs(e);
Η αναζήτηση του χαμένου bit... χαχχα
Τώρα αυτό είναι απάντηση?
ReplyDeleteΟ γενικος κανονας ειναι οτι δεν κανουμε *ποτε* ελεγχο ισοτητας σε floating point μεταβλητες.
ReplyDeleteΚανουμε κατι σαν "if (abs(a-b)<0.0001)" ωστε μετα να εχουμε να πονοκεφαλιαζουμε ΚΑΙ για την ταξη μεγεθους του tolerance, η χρησιμοποιουμε αλχημειες για fixed-point arithmetic.
Ψαξε και βρες το βιβλιο "Numerical Recipes". Ειναι το κλασσικο βιβλιο αναφορας για προγραμματιστριες που πρεπει να κανουν ταυτοχρονως πλεξιμο και αριθμητικη αναλυση. Υπαρχει σιγουρα για Fortran και C, δεν ειμαι σιγουρος για Pascal, αλλα οι αυγοριθμοι ειναι παντου οι ιδιοι και ο κωδικας διαβαζεται ευκολα.
Ναι πλεόν το έμαθα το μάθημά μου και το κάνω έτσι όπως λές. Δηλώνω ZERO = 1e-9 και χρησιμοποιώ αυτή τη σταθερα ως μηδέν.
ReplyDeleteΤο βιβλίο το έχω αγορασμένο από το 2000 :)