המעשה המופלא בקבוע המסתורי 0x5f3759df (חלק ב' - הקשה)

בואו נמשיך את סיפור המעשה מהפוסט הקודם על קוד מסתורי שמחשב את \( f\left(x\right)=\frac{1}{\sqrt{x}} \) בצורה יעילה עד להפתיע. אנחנו כבר כמעט מסוגלים להבין מה הקוד עושה באופן מלא, רק צריך קודם להבין את העניין הפעוט הזה של איך מספרים מיוצגים במחשב.

פרק רביעי, ובו ביטים עושים דברים

בגדול, אפשר לומר שחלק נכבד מהיקום כולל דברים שמורכבים מדברים. יצירות לגו מפוארות מורכבות מאבני לגו בסיסיות. מולקולות חומר פשוטות ומסובכות בנויות מאטומים (והם בתורם בנויים מ… עזבו, לא פוסט בפיזיקה). המידע הגנטי שלנו שמקודד ב-DNA בנוי מארבע “אותיות בסיסיות” G,C,T,A. התמונה שאתם רואים במסך המחשב מורכבת מפיקסלים - נקודות על המסך שכל אחת מהן היא בעלת צבע אחיד (שבתורו מורכב משלושה צבעים - אדום, ירוק, כחול - בעוצמות משתנות). כאשר מדברים על משהו שמורכב מאבני יסוד בסיסיות לא מספיק לומר מה אבני היסוד - גם צריך להסביר איך הן מתחברות זו לזו כדי ליצור דברים. אצטון ופרופיונאלדהיד הן שתי מולקולות שונות שמורכבות בדיוק מאותם אטומים אבל מחוברים בצורה שונה. כרגע עומד מולי רובוט לגו שבעזרת אותן אבני בניין בדיוק שלו יכלתי להרכיב גם מסוק או טנדר.

כאשר מדובר על לגו, יש אינספור אבנים בסיסיות, אבל אצלנו במדעי המחשב יש בדיוק שני אבני בניין: הספרות 0 ו-1, שבהקשר הזה נקראות ביטים. אנחנו בונים מהן הכל. כל פריט מידע במחשב הוא, בסופו של דבר, ביטים. המספרים השלמים; והמספרים הממשיים; וקומנדר קין והרפתקאותיו והמסמך שאני כותב כרגע וכל המידע שאי פעם נכתב בפייסבוק וכל סרט קולנוע שאי פעם אוכסן במחשב - כולם בסופו של דבר בנויים רק מ-0 ו-1. למה? למה לא לאפשר אבני בניין מורכבות יותר? כי קל, ברמת החומרה של המחשבים, לעבוד רק עם שתי אבני הבניין הללו (בלשון ציורית ולא מדויקת, קל להבדיל ביניהן במערכת אלקטרונית בעזרת “יש זרם חשמלי” ו”אין זרם חשמלי”). גם האופן שבו אנחנו מחברים את 0 ו-1 זה לזה הוא פשוט ביותר - אנחנו פשוט כותבים אותם בשורה. למשל:

011010101

רצף הביטים הזה הוא דוגמא לפריט מידע שמאוחסן במחשב. אבל איזה מידע? ובכן, כאן ההקבלה ללגו או למולקולות קצת משתנה. המחשב יכול לקחת את אותה סדרה של אפסים ואחדות ולחשוב עליה כאילו היא מייצגת דברים שונים ומשונים. היא יכולה לייצג מספר, והיא יכולה באותה מידע בדיוק גם לייצג אות. בשל כך המחשב על פי רוב מבצע איזה שהוא סוג של פירוש כדי להבין איך לחשוב על הסדרה הזו כרגע. זה דומה לאופן שבו מילים נהגות בתור סדרה של הברות בסיסיות, אבל אותו צליל, בשפות שונות, יכול להיות בעל משמעויות שונות. “היא” בעברית ו-he באנגלית נשמעים אותו דבר אבל הם מתפרשים שונה, בהתאם לשומע והשפה שהוא מצפה לשמוע באותו הרגע.

שפות תכנות משתמשות במשתנים. משתנה הוא מקום בזיכרון של המחשב שניתן לו שם קליט בתוך הקוד של התוכניות ובאמצעות השם הזה אפשר לומר לתוכנית לעשות עם המקום הזה דברים - לכתוב שם הרבה פעמים 0, לכתוב שם הרבה פעמים 1, לכתוב 01010101 וכדומה. כדי שלתוכנית יהיה קל להבין מה בדיוק אמור לקרות עם המקום הזה בזכרון, למשתנים בדרך כלל יש טיפוס - משהו שכולל מידע על “מה המשתנה אמור לייצג”. מה זה בדיוק אומר - זה משתנה משפת תכנות לשפת תכנות, ואפילו מסוג אחד של טיפוס לסוג אחר, מבחינת רמת הפירוט שאליה ההגדרה נכנסת. למשל, זה יכול לכלול מידע על כמות הביטים שהמשתנה משתמש בהם (לפעמים בכמה ביטים בדיוק הוא משתמש, ולפעמים בכמה ביטים לכל הפחות הוא אמור להשתמש). פרט לכמות הביטים הטיפוס גם כולל לפעמים מידע על איך אמורים להתייחס אליהם. אותנו מעניינים בהקשר של הקוד שלנו שני טיפוסים שבהם משתמשים בשפת C: הראשון הוא long, שמיועד לתאר ערכים מספריים שלמים, והשני הוא float שנועד לתאר מספרים שיכולים להיות גם שבריים ומיוצגים בייצוג שנקרא נקודה צפה ואסביר בקרוב.

נתחיל בלדבר על long. זו דוגמה לטיפוס שמגדיר את ה”בערך” אבל ההגדרה שלו לא נכנסת לפרטים מדוייקים. אין הגדרה חד משמעית לכמות הביטים שמשתנה מסוג long משתמש בהם, אבל התקן קובע שהוא ישתמש לפחות ב-32 ביטים. בהקשר של הקוד שהופיע ב-Quake אנחנו יודעים שהכוונה הייתה לבדיוק 32 ביטים כי אחרת לא ברור מה הולך שם. לצורך מה שקורה בקוד חשוב שמספר הביטים של ה-long יהיה שווה למספר הביטים ש-float משתמש בו (והמספר הזה הוא חד משמעית 32, כי כך קובע התקן). אין גם הגדרה חד משמעית לאופן שבו הביטים של משתנה מטיפוס long אמורים להתנהג, אבל בפועל מה ש-long תמיד עושה הוא לחשוב על הביטים שלו כמייצגים מספר שלם בבסיס בינארי. יש לי פה הסבר על בסיסי ספירה אבל הנה הרעיון הבסיסי: בבסיס בינארי כל מספר מיוצג על ידי סדרת ביטים שמתארת אותו כסכום של חזקות של 2. למשל, סדרת הביטים 1101 אומרת “זה המספר שמיוצג על ידי הסכום \( 2^{0}+2^{2}+2^{3}=1+4+8=13 \)”. הביט הכי שמאלי מייצג את החזקה הכי גבוהה של 2 שמחברים. כל זה קורה גם בבסיס 10, כמובן: אנחנו רגילים כבר לתרגם אוטומטית משהו כמו 1,089 ל”אלף ועוד שמונים ועוד תשע” בלי אולי לשים לב לכך שאנחנו מחברים חזקות של 10 שנכפלות במקדם כלשהו בבסיס בינארי המקדם הוא רק 0 או 1, אבל הרעיון הוא אותו רעיון.

יש לייצוג מספרים על ידי long רמת סיבוך נוספת שאני חוסך מכם בפוסט הזה כי היא לא רלוונטית - האופן שבו מיוצגים מספרים שליליים. לא ניכנס לזה כרגע. ובמקום זה נעבור לדבר על הייצוג של מספרים על ידי נקודה צפה.

פרק חמישי, ובו נתוודע לפרנקנשטיין של שפות התכנות - הנקודה הצפה

בואו נתחיל מכך שאודה ששיקרתי לכם. קודם הצגתי את העניינים כאילו הרעיון ב-long הוא ייצוג של מספר שלם והרעיון ב-float הוא ייצוג של מספר “ממשי”, בפרט כזה שיכול להיות שבר. ובכן, ראשית כל, float לא יכול לייצג מספר ממשי כללי, למשל את \( \pi \). כל מה שהוא יודע לייצג הוא מספרים רציונליים - שברים שאפשר להציג בתור \( \frac{a}{b} \) כאשר \( a,b \) שניהם שלמים. שנית, אם כל מה שהייתי רוצה הוא לייצג רציונליים הייתי יכול פשוט להשתמש בזוג long, לא לגמרי ברור שצריך טיפוס נתונים חדש בשביל זה. אם כן, “לייצג מספר ממשי” או “לייצג שבר” איננה הסיבה שבגללה אנחנו מתעניינים ב-float. אז מה כן הסיבה?

הסיבה היא שלפעמים לא אכפת לנו אם המספר שלנו לא מיוצג בצורה מדוייקת. לפעמים אפשר לחפף ולעגל קצת, אם זה משתלם לנו. הרעיון ב-float הוא לוותר קצת על הדיוק המושלם ש-long מציע ותחת זאת להרחיב בצורה משמעותית את טווח המספרים שאפשר לייצג באמצעות 32 ביט. אם באמצעות long אפשר לייצג במדויק כל מספר בתחום שבין 0 ל-\( 2^{32}-1 \) (מי שרוצה לנטפק - תזכרו, אמרתי שלא אכנס לשלמים שליליים פה) הרי שבאמצעות float אפשר לייצג מספרים עד בערך \( 2^{127} \), ושברים עד בערך \( 2^{-126} \) ואפילו קטן יותר מכך. המחיר הוא שאי אפשר לייצג את כל המספרים בטווחים הללו; יש לנו מגבלת דיוק. על פי רוב, בשימושים של float שמעניינים אותנו המגבלה הזו לא מפריעה לנו.

ב-float גם כן יכולים להיות מספרים שליליים, והפעם גם אתייחס לאופן שבו מייצגים אותם כי הוא קצת יותר פשוט מאשר ב-long, ולשם שינוי גם מוגדר היטב. בכלל, לנקודה צפה יש יתרון שהיא מוגדרת יחסית טוב בסטנדרט של ה-IEEE ורוב מי שמממש נקודה צפה (בתוכנה/חומרה) יתאים את עצמו לסטנדרט. הקוד של 0x5f3759df מתבסס על זה, כמובן.

מספר בייצוג float מורכב מ-32 ביט, שמחולקים לשלוש קבוצות: הביט הראשון, השמאלי ביותר, הוא הסימן של המספר. אם הוא 0, המספר חיובי; אם הוא 1, המספר שלילי. 8 הביטים הבאים נקראים האקספוננט של המספר, ו-23 הביטים הנותרים נקראים המנטיסה שלו.

כדי להבין את המשמעות של אלו, בואו נראה לרגע על דרכים שונות שבהן אפשר לייצג את המספר \( 314.15 \). אני יכול לכתוב סתם \( 314.15 \), אבל אני גם יכול לכפול בחזקות של 10: למשל, לכתוב \( 31.415\cdot10^{1} \), או \( 3.1415\cdot10^{2} \), או \( 3141.5\cdot10^{-1} \) וכדומה. הבנתם את הרעיון: אני לוקח את המספר הבסיסי, \( 314.15 \), ו”מזיז” את הנקודה העשרונית (“מציף” אותה) כשהמחיר הוא כפל בחזקה מתאימה של 10. הזזתי את הנקודה שמאלה? אני כופל בחזקה חיובית של 10. הזזתי אותה ימינה? אני כופל בחזקה שלילית. באופן הזה אפשר להחליט שכל מספר ייוצג בצורה “נורמלית” שבה יש בדיוק ספרה אחת משמאל לנקודה העשרונית; הייצוג ה”נורמלי” של \( 314.15 \) יהיה, אם כן, \( 3.1415\cdot10^{2} \). האקספוננט של המספר הזה הוא החזקה של 10 בייצוג הנורמלי, והמנטיסה שלו היא המספר שבו מכפילים מצד שמאל.

בואו נראה עוד דוגמה. את המספר \( 1,000 \) קל לייצג עם נקודה צפה: \( 1.0\cdot10^{3} \). מה על מספר ששונה ממנו טיפ-טיפה, נאמר 1,002? אותו אפשר לייצג על ידי \( 1.002\cdot10^{3} \). שימו לב מה קרה - נזקקנו ליותר ספרות במנטיסה כדי לייצג את המספר הזה מאשר את \( 1,000 \) שמיוצג בצורה ישירה באמצעות האקספוננט. באופן דומה, אם אני ארצה לייצג את מיליון זה יהיה קל, אבל אם ארצה לייצג את “מיליון ועוד 2” אצטרך עוד הרבה ספרות במנטיסה. וגם את \( 10^{100} \) קל לייצג, אבל לייצג את \( 10^{100}+2 \) כבר יהיה יותר מדי עבורי - אין לי מספיק מקום במנטיסה בשביל זה כי אצטרך לכתוב \( 1.000\dots2 \) כאשר יש בערך מאה אפסים. הנה כי כן, זו בדיוק מגבלת ה”חוסר דיוק” שדיברתי עליה. את \( 10^{100}+2 \) אני לא יכול לייצג, אבל אני יכול להסתפק ב-\( 10^{100} \) שאותו אני כן יכול לייצג והוא קירוב מצויין ל-\( 10^{100}+2 \). כל עוד אני לא חייב ייצוג מדויק של כל המספרים הללו.

בסדר עד כאן? אז בואו נסבך קצת. הצגתי את מספרי הנקודה הצפה שלי כאילו הם כתובים בבסיס עשרוני, אבל בפועל float מיוצג על ידי ייצוג בינארי (יש גם נקודה צפה של מספרים עשרוניים אבל זה לא רלוונטי לכאן). זה אומר שהאקספוננט והמנטיסה שניהם נכתבים בבסיס בינארי ואנחנו כופלים את המנטיסה בחזקה של 2 ולא 10, אבל זה גם אומר עוד משהו - אין צורך לזכור את הספרה הבודדת שמשמאל לנקודה העשרונית במפורש; אנחנו יודעים שהיא לא 0, כי ככה מוגדרת הצורה הנורמלית של מספר - יש משמאל לנקודה בדיוק ספרה אחת שאינה 0. לכן, עבור מספר נקודה צפה בבסיס בינארי, המנטיסה מתארת רק את מה שקורה מימין לנקודה העשרונית - החלק השברי של המספר.

סיבוך נוסף שטרם דיברתי עליו הוא האופן שבו מאפשרים למספרים שליליים להופיע בתוך האקספוננט. מה שעושים הוא להשתמש במשהו שנקרא bias. יש 8 ביטים של אקספוננט, מה שאומר שאפשר לייצג איתם כל מספר מ-0 עד 255. מכיוון שרוצים חצי חיוביים וחצי שליליים, מגדירים bias של \( B=127 \) ומגדירים שהוא תמיד מחוסר מהאקספוננט. כלומר, אם \( E \) מייצג את \( 200 \) אז האקספוננט של המספר יהיה \( 2^{E-B}=2^{73} \). בצורה הזו האקספוננט הגבוה ביותר הוא לכאורה \( 2^{128} \) והנמוך ביותר הוא \( 2^{-127} \) אבל בפועל הסטנדרט לא מרשה לנו להשתמש באקספוננטים 11111111 ו-0000000 באופן חופשי: את 11111111 שומרים כדי לייצג את אינסוף ואת NaN (ערך שאומר “לא קיבלתי מספר”) ואילו 00000000 שמור כדי לאפשר ייצוג של 0 ושל מספרים נמוכים במיוחד (משהו שנקרא denormalized numbers שאני פשוט לא הולך לדבר עליהם כאן כי לא צריך את זה). לכן טווח האקספוננטים החוקי הוא מ-\( 2^{127} \) ועד \( 2^{-126} \).

בואו נחזור על מה שיש לנו במספר float: יש 32 ביטים בסך הכל. ביט אחד, שנקרא לו \( s \), הוא ביט הסימן. 8 הביטים הבאים, שנקרא להם \( E \), הם הביטים של האקספוננט: אפשר לכתוב \( E=E_{7}E_{6}E_{5}E_{4}E_{3}E_{2}E_{1}E_{0} \) כאשר כל \( E_{i} \) באגף ימין הוא ביט בודד. לסיום, המנטיסה תסומן ב-\( M \) (ברשותכם, לא אכתוב את כל הביטים שלה). עכשיו, בהינתן \( s,E,M \) אפשר לחשב את הערך המפורש של המספר שהם מייצגים ככה:

\( \left(-1\right)^{s}\cdot1.M\cdot2^{E-B} \)

לפעמים במקום \( 1.M \) יותר נוח וקריא לכתוב \( 1+\frac{M}{2^{23}} \) או אפילו לסמן \( m=\frac{M}{2^{23}} \) ואז לכתוב

\( \left(-1\right)^{s}\cdot\left(1+m\right)\cdot2^{E-B} \)

סיימנו עם זה! עכשיו אנחנו מבינים מספרי נקודה צפה ברמה שתספיק להמשך הפוסט. נעבור סוף סוף לשאלת השאלות: מה קורה כשאני לוקח float וחושב על הביטים שלו כמגדירים long?

פרק שישי, ובו שלמים ושברים ולוגריתמים יפים אלו דברים שאותי משמחים

בואו נסתכל על שורת ה”המרה” הידועה לשמצה מהקוד:

i = * ( long * ) &y; // evil floating point bit level hacking

כפי שאמרתי קודם, מה שקורה בשורה הזו איננו המרה - אנחנו לא אומרים לתוכנית לקחת את ה-float שלנו ולעגל אותו עד שיתקבל מספר שלם או משהו. אנחנו עושים משהו ברוטלי ומסוכן באופן כללי: לוקחים את 32 הביטים בזיכרון שמיוצגים על ידי y ואומרים לתוכנית לחשוב עליהם בכוח בתור long. למי שסקרן, זה האופן הטכני שבו זה נעשה: ראשית אנחנו מבקשים “נא לתת לנו את הכתובת בזכרון שבה המידע של y שוכן”. זה מה שעושה האופרטור & כשהוא מוצמד ל-y. אחר כך אנחנו לוקחים את הכתובת הזו, שכרגע התוכנית חושבת עליה בתור “כתובת של float’’, ואנחנו מבצעים עליה פעולה שבשפת c נקראת casting ומתבצעת על ידי ה-long * שבסוגריים (הסוגריים עצמם אומרים לתוכנית שיש כאן פעולת casting). הפעולה הזו אומרת לתוכנית - “נכון שיש לך כתובת שאת חושבת עליה בתור כתובת של float? מעכשיו תחשבי עליה בתור כתובת של long”. הכוכבית האחרונה בתחילת השורה היא דרך לומר “אוקיי, עכשיו נא לתת לי את הערך המספרי שכתוב בתוך כתובת הזיכרון שלך”. התוכנית עושה את הדבר הבא: בשלב הזה, יש לה כתובת זכרון ולידה סימון “הכתובת הזו מכילה long’’. אז התוכנית לוקחת את 32 הביטים מהכתובת, מתייחסת אליהם בתור מספר long ומציבה בתוך i. כל זה ברור למי שמכיר את השפה, ואני מנחש שמי שלא מכיר אותה כבר הלך לאיבוד. לא נורא, לא צריך להבין מה השורה עושה ברמה הטכנית הזו, רק מה האפקט שזה משיג.

ומה זה עושה בפועל, למספר? מיש-מש, זה מה שזה עושה. ביטים שלפני רגע הייתה להם משמעות אחת מקבלים משמעות לא לגמרי קשורה, אבל גם לא לגמרי שונה. אנחנו עדיין יכולים לחשוב על הביטים בתור שלוש קבוצות - \( s,E,M \) - רק שעכשיו כל קבוצה תורמת משהו למספר השלם שמיוצג על ידי ה-long.

המשך הניתוח שאציג מתבסס בעיקר על המאמר הזה שמספק ניתוח יפה מאוד של הסיפור הזה. יש עוד ניתוחים שונים ומשונים שאפשר לעשות וקשה לי לומר מי מהם הוא ה”נכון”, אבל זה שאציג כרגע הוא ללא ספק הפשוט ביותר (והנחמד ביותר, לטעמי) מביניהם.

בגדול, מאוד בגדול, מה שקורה כשמעבירים ככה מספר מ-float אל long הוא שממירים אותו ללוגריתם של עצמו כפול איזה שהוא קבוע. בואו נזכר מה זה לוגריתם בכלל. אם \( x=2^{y} \) אז \( \log x=y \). כלומר, הלוגריתם של \( x \) הוא המספר שאם מעלים את 2 בחזקה שלו, מקבלים את \( x \) (אני מציג פה את מה שנקרא “לוגריתם על בסיס 2” כי זה מה שרלוונטי לנו בפוסט הזה). למשל, \( \log8=3 \) כי \( 2^{3}=8 \). לעומת זאת \( \log7 \) לא הולך לצאת מספר יפה אלא משהו אי-רציונלי שנמצא אי שם בין 2 ל-3.

אנחנו אוהבים לוגריתמים כי הם מבצעים מעין “הורדה בדרגת הקושי” לפעולות חשבוניות מסובכות. כפל וחילוק הופכים להיות חיבור וחיסור, ואילו העלאה בחזקה והוצאת שורש הופכות להיות כפל וחילוק. הנה הכללים המתאימים:

\( \log\left(a\cdot b\right)=\log a+\log b \)

\( \log\left(\frac{a}{b}\right)=\log a-\log b \)

\( \log a^{n}=n\log a \)

\( \log\sqrt[n]{a}=\frac{1}{n}\log a \) (זה בעצם נובע מכך ש-\( \sqrt[n]{a} \) הוא \( a^{\frac{1}{n}} \)).

בימים עברו, לפני המצאת המחשבון, השתמשו בטבלאות לוגריתמים כדי לבצע חישובים מסובכים: טבלת לוגריתמים כללה ערכים של מספרים ולידם את הערך של הלוגריתם שלהם. אם הייתי רוצה לבצע פעולה מסובכת כמו כפל 128 ב-512 מה שהייתי עושה הוא להסתכל בטבלת הלוגריתמים, לראות שהלוגריתמים של שני המספרים הללו הם 7 ו-9 בהתאמה, לחבר את 7 ו-9 לקבלת 16, ואז להסתכל בטבלת הלוגריתמים ולראות שהמספר שהלוגריתם שלו הוא 16 הוא המספר 65536. היו גם כלים מיוחדים בשם סרגלי חישוב שסייעו לעשות את החישוב הזה. בצורה הזו אמנם נדרשה עבודה ראשונית ביצירה של טבלת הלוגריתמים/סרגל החישוב, אבל בעבודה היומיומית הם חסכו הרבה כאב ראש בביצוע פעולות חשבון. מה שאני רוצה לומר כאן הוא שלוגריתמים זה דבר נפלא שכשלומדים אותו בתיכון לפעמים בכלל לא מבינים בשביל מה הוא טוב.

המקרה הנוכחי מושלם עבור לוגריתמים. \( \frac{1}{\sqrt{x}} \) זו דרך אחרת לכתוב \( x^{-\frac{1}{2}} \). נפעיל על זה לוגריתם ונקבל \( \log\left(x^{-\frac{1}{2}}\right)=-\frac{1}{2}\log x \). אז ברמת הלוגריתמים כל פעולת החישוב המסובכת שאנחנו רוצים לבצע היא בסך הכל כפל במינוס חצי. ומה תגידו? כפל כזה מתבצע בפועל!

i = 0x5f3759df - ( i >> 1 ); // what the fuck?

תתעלמו לרגע מהקבוע המסתורי שלנו. מה שיש באגף ימין הוא \( -\frac{i}{2} \). בשביל לראות את זה אני צריך להסביר סוף סוף מה אומר ה-1<< הזה. ובכן,<< זה אופרטור שמופעל על מספר שלם ומבצע הזזה ימינה של הביטים שלו. ה-1 שמצד ימין של האופרטור אומר כמה להזיז ימינה - 1 פירושו להזיז בדיוק פעם אחת. כלומר, המספר \( 0110 \) יהפוך להיות המספר \( 0011 \) וכן הלאה. בפועל הפעולה הזו מבצעת חלוקה ב-2 של המספר השלם (עם עיגול למטה במקרה שמקבלים שבר). אם כן, מה שהשורה המסתורית הזו הוא הוא לכפול את \( i \) במינוס חצי ולהוסיף לו את הקבוע המסתורי בתור… לא לגמרי ברור בתור מה עדיין. אז בואו נמשיך עם הפרטים.

פרק שביעי, שבו התעלומה באה על פתרונה והקוראים מתלוננים על אנטי-קליימקס

כעת, אמרנו שמספר float מיוצג על ידי ביט אחד של \( s \), אחריו ביטים של \( E \) ואחר כך ביטים של \( M \). הביטים של \( M \) הם הראשונים, ולכן הם אכן מייצגים בדיוק את המספר \( M \). הביטים של \( E \), לעומת זאת, מתחילים החל מהמקום ה-24. אם הביט במקום ה-1 מייצג את הספרה שמתאימה ל-\( 2^{0} \), הרי שהביט במקום ה-24 מיצג את הספרה שמתאימה ל-\( 2^{23} \), ולכן \( E \) מייצג את המספר

\( 2^{23}E_{0}+2^{24}E_{1}+\dots+2^{30}E_{7}=2^{23}\left(E_{0}+2^{1}E_{1}+\dots2^{7}E_{7}\right)=2^{23}\cdot E \)

ולבסוף, הביט \( s \) של הסימן מייצג את \( 2^{31} \). כלומר, המספר כולו מתפרש בתור ה-long הבא: \( 2^{31}s+2^{23}E+M \). בפועל, אפשר להתעלם מהביט \( s \) של הסימן: הוצאת שורש היא פעולה שאנחנו מבצעים רק עבור מספרים חיוביים, ולכן הסימן של ה-float הוא חיובי, מה שאומר ש-\( s=0 \) בכל מה שנעשה. על כן המספר מתפרש בתור \( L\cdot E+M \) כאשר \( L=2^{23} \) הוא קבוע שיאפשר לנו לקרוא יותר בקלות מכאן ואילך. זכרו שהמספר המקורי בתור float היה \( x=\left(1+\frac{M}{L}\right)2^{E-B} \). נשאלת כעת השאלה - עד כמה \( L\cdot E+M \) הזה אכן יהיה דומה ל-\( \log x \)? לצורך כך כדאי לקבל הערכה כלשהי לערך של \( \log x \), ולצורך כך אני אשתמש בקירוב מוכר במתמטיקה: זה ידוע שכאשר \( t \) הוא קטן יחסית, אז \( \log\left(1+t\right)\approx t \) (למעוניינים, זה נובע מפיתוח טיילור של \( \log\left(1+t\right) \)). במקרה שלנו, \( \frac{M}{L} \) הוא קטן יחסית (כי הגדול של \( M \) חסום על ידי \( L \)) ולכן אפשר להשתמש בקירוב \( \log\left(1+\frac{M}{L}\right)\approx\frac{M}{L} \). מצד שני, אין סיבה שנשתמש בקירוב הזה באופן עיוור ופשוט נתעלם מכך שאולי כדאי להוסיף “תיקון” כלשהו שיפצה על החלקים שהעפנו מהקירוב. אז נגדיר פרמטר \( \sigma \) שאת הערך שלו נוכל לבחור באופן שרירותי ונשתמש בקירוב הבא: \( \log\left(1+\frac{M}{L}\right)\approx\frac{M}{L}+\sigma \). לא לגמרי ברור בשלב הזה אילו ערכים של \( \sigma \) הם טובים לנו ואיזה לא (אולי \( \sigma=0 \) הוא טוב?) ולכן אנחנו לא מתחייבים על ערך ספציפי עבורו.

כעת, נקבל מהזהויות שקשורות בלוגריתם שראינו למעלה את הדבר הבא:

\( \log\left(x\right)=\log\left(\left(1+\frac{M}{L}\right)2^{E-B}\right)=\log\left(1+\frac{M}{L}\right)+\log2^{E-b} \)

\( \approx\frac{M}{L}+\sigma+E-B \)

עכשיו, אם נסמן \( y=\frac{1}{\sqrt{x}} \), הרי ש-\( y \) הוא המספר שאנחנו מחפשים. בייצוג על ידי נקודה צפה גם הוא ישתמש בפרמטרים \( M,E \), אבל כאלו שיהיו שונים מאלו של \( x \). לכן נשתמש בסימונים כדי להבדיל ביניהם: את \( M,E \) שהשתמשתי בהם עד כה אסמן מעכשיו ב-\( M_{x},E_{x} \), ואילו את האקספוננט והמנטיסה של \( y \), שאותם אני מחפש, אסמן ב-\( E_{y},M_{y} \). אותו חישוב כמו קודם עובד גם עבור \( y \), ולכן יש לנו עכשיו שלוש משוואות:

\( \log\left(x\right)\approx\frac{M_{x}}{L}+\sigma+E_{x}-B \)

\( \log\left(y\right)\approx\frac{M_{y}}{L}+\sigma+E_{y}-B \)

\( \log y=-\frac{1}{2}\log x \)

נשלב את המשוואות הללו יחד:

\( \frac{M_{y}}{L}+\sigma+E_{y}-B\approx-\frac{1}{2}\left(\frac{M_{x}}{L}+\sigma+E_{x}-B\right) \)

נעביר את הקבועים \( \sigma,B \) אגף ונקבל

\( \frac{M_{y}}{L}+E_{y}\approx\left(B+\frac{1}{2}B\right)-\left(\sigma+\frac{1}{2}\sigma\right)-\frac{1}{2}\left(\frac{M_{x}}{L}+E_{x}\right) \)

כלומר

\( \frac{M_{y}}{L}+E_{y}\approx\frac{3}{2}\left(B-\sigma\right)-\frac{1}{2}\left(\frac{M_{x}}{L}+E_{x}\right) \)

לסיום, נכפול את שני האגפים ב-\( L \) ונקבל

\( M_{y}+LE_{y}\approx\frac{3}{2}L\left(B-\sigma\right)-\frac{1}{2}\left(M_{x}+LE_{x}\right) \)

ותראו מה קיבלנו! ה-\( M_{x}+LE_{x} \) באגף ימין הוא בדיוק הערך של המספר שממנו התחלנו, כשמפרשים את הביטים שלו בתור long; והערך באגף שמאל הוא מספר שאם נפרש את הביטים שלו בתור float אז המנטיסה שלו תהיה \( M_{y} \) והאקספוננט שלו יהיה \( E_{y} \). זה גם בדיוק מה שעושים בשורה הבאה:

y = * ( float * ) &i;

ועל כן, המשוואה שלעיל היא בדיוק מה שמנחה את שורת ה-what the fuck? הידועה לשמצה:

i = 0x5f3759df - ( i >> 1 ); // what the fuck?

זה מסביר למה היא נראית ככה וגם מיהו הקבוע המסתורי: הוא פשוט \( \frac{3}{2}L\left(B-\sigma\right) \). זכרו ש-\( L \) הוא פשוט המספר \( 2^{23} \) ו-\( B \) הוא המספר \( 127 \) - אלו פרמטרים שנטועים עמוק בהגדרה של ה-IEEE למהו float, אבל גם אם הערכים שלהם היו שונים היינו עדיין מקבלים משוואה דומה, רק עם “קבוע מסתורי” שונה.

כמובן שעכשיו נשאלת השאלה איזה ערך של פרמטר \( \sigma \) הולך לתת את הקבוע 0x5f3759df מתוך הביטוי \( \frac{3}{2}L\left(B-\sigma\right) \). התשובה היא שזה \( \sigma=0.0450465 \), אבל זה בעצם לא אומר לנו שום דבר. כאן בעצם מגיע החלק המאכזב ביותר בכל הסיפור - מי שיצפה לראות איזה הגיון קוסמי שבזכותו נוצר דווקא המספר 0x5f3759df ולא אחרים לא ימצא אותו - זה ככל הנראה מספר שכותב הקוד הגיע אליו אחרי קצת ניסוי וטעיה - ראה שהוא עובד מספיק טוב, ולא ניסה לשפר יותר. עדיין, אם מישהו רוצה ניתוח קצת יותר מפורט של ערכים אפשריים אחרים, אפשר להסתכל בתזה הזו, שבכלל נמנעת משימוש בלוגריתמים ומסתכלת בצורה מפורשת מאוד על ההבדל בתוך ה-float שגורמות הפעולות שמבצעים עליו. מסבירים שם, למשל, למה קבוע שנותן ערך טוב יותר בתור הקירוב אחרי השורה הזו הוא פחות טוב באופן כללי, כי ניוטון-רפסון מחזיר עליו תוצאה פחות נחמדה, וגם נותנים ערך טוב יותר מ-0x5f3759df בתור קבוע קסם מסתורי עבור הפונקציה. מבחינתי זה מחסל את 0x5f3759df המסכן לגמרי - הוא לא כל כך מעניין אם הבחירה בו הייתה כל כך שרירותית. אולי יום אחד אתבדה ואגלה שהוא נבחר מסיבות מצויינות שאיני מכיר.

פרק שמיני, שבו אנחנו תוהים בשביל מה כל זה היה טוב

הסיפור שלנו מתקרב לסופו, אבל אני רוצה להזכיר למה בכלל נכנסנו אליו מלכתחילה. כזכור, שם המשחק הוא גרפיקה. הגרפיקה הזו:

בשביל לייצר גרפיקה יפה שכזו צריך לדעת לחשב כל מני חישובים. למשל, איך אור משתקף מכל מני משטחים. כשהיינו בימי Wolf3D העליזים כל המשטחים היו פשוטים מאוד - קירות שעמדו בזווית של 90 מעלות ביחס לרצפה וזהו. אבל בעולם תלת-ממדי שנראה טוב, זה לא המצב. יש משטחים באלכסונים, ויש משטחים מעוגלים ועוד ועוד. כשרוצים לחשב איך תתנהג קרן אור שפוגעת במשטח בנקודה כלשהי, אנחנו צריכים לדעת משהו על “הכיוון המקומי” של המשטח באותה נקודה. הכיוון הזה מיוצג באמצעות וקטור יחידה במרחב התלת ממדי. “וקטור יחידה” פירושו שהאורך של הוקטור הוא 1. למה דווקא 1? כי פעולות שמערבות את הוקטור הזה דורשות שהוא יוכפל סקלרית בדברים, ואם האורך שלו הוא לא 1 אז הוא “ינפח” את הדברים הללו באופן מלאכותי. בפועל מה שקורה הוא שקודם כל מוצאים את הכיוון של הוקטור - כלומר, מוצאים וקטור כלשהו שמצביע בכיוון הנכון, ואז מנרמלים את הוקטור - מחלקים אותו באורך של עצמו. אם \( v \) הוא וקטור, אז האורך שלו הוא \( \|v\|\triangleq\sqrt{v\cdot v} \). על כן, הוקטור המנורמל \( \frac{v}{\|v\|} \) שווה ל-\( v\cdot\frac{1}{\sqrt{v\cdot v}} \). הופס! אנחנו צריכים למצוא את ההופכי של שורש של משהו!

אם מסתכלים בקוד ומחפשים שימושים של Q_rsqrt זה בדיוק מה שמוצאים. למשל:

שמופיע בקובץ q_math.c (ואפשר לראות כרגע כאן).

“רגע, זה הכל?” אולי אתם שואלים. ובכן, צריך לזכור שאנחנו מחשבים את הוקטורים הללו עבור אינספור נקודות על כל המשטחים שסביבנו. ככל שיש יותר וקטורים, כך התיאור שלנו של המשטחים נראה יותר ריאליסטי. לכן הפונקציה הזו הולכת להיקרא המון פעמים. ככה בדיוק זה אופטימיזציות: בסופו של דבר צוואר הבקבוק הוא בדיוק בפונקציות הכי קטנות ופשוטות ושם שוברים את הראש על מציאת דרכים טובות יותר לבצע את החישוב.

אז בעצם, מה גרם לחישוב להיות כל כך טוב? שילוב של שני דברים: ראשית, ניוטון-רפסון, שהיא שיטה מגניבה באופן כללי; ושנית, קירוב (קונספטואלי, לכל הפחות) שבוצע ללא שום המרה מפורשת אלא פשוט התייחסות קצת שונה לביטים של הערכים שפעלנו עליהם. אלו הרעיונות המגניבים כאן. ומה עם המספר המסתורי 0x5f3759df? למה בדיוק הוא נבחר? האם יש לו איזו תכונה קסומה שבעטיה הוא נבחר? כנראה שלא, אבל זו תישאר אחת מהתעלומות הקטנות של מדעי המחשב גם לדורות הבאים.