343. CORDIC Method

שבוע מתמטי טוב 2
Math

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

הפעם אני רוצה להסביר על אלגוריתם CORDIC בו משתמשים מחשבונים בד”כ על מנת לממש חישובים של הפונקציות הטריגונומטריות.

נתאר את הבעיה: יש לנו מספר ממשי כלשהו המתאר זווית (לצורך העניין נגביל אותה בגודל מסוים, נגיד לבין 0 ל$ \frac{\pi}{4} $ ($ {45}^{\circ} $)) ואנחנו מעוניינים לחשב את ה$ \sin $/$ \cos $/$ \tan $ שלה. לחלופין יכול להיות שנרצה לחשב את אחת הפונקציות ההפוכות הנ”ל.

ראשית אנחנו מחשבים את $ \arctan x $ עבור x שהוא הופכי של חזקה של 2, כלומר $ \displaystyle{\arctan \frac{1}{2^k}} $ לכל הkים בתחום מסוים.

את זה אנחנו עושים ע”י שימוש בטור טיילור של $ \arctan x $: $ \displaystyle{\arctan{x} = \sum_{n = 0}^{\infty}{\frac{\left( -1 \right)^n}{2n + 1} }x^{2n+1}} $ כאשר $ |x| \le 1 $.

את החישובים האלה אנחנו שומרים בטבלה סטטית בצד (איזשהו HashTable של $ k $ ל$ \displaystyle{\arctan \frac{1}{2^k}} $)

כעת אלגוריתם CORDIC עובד כך:

אנחנו מקבלים זווית $x$ שאנחנו מעוניינים לחשב אותה.

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

  • $ \theta$– מייצג את הזווית שחישבנו את ה $ \cos $וה$ \sin $ שלה עד כה.
  • $ \mathrm{vector} $ – מספר מרוכב בכיוון של $ \cos\theta + i \sin\theta $
  • $ k $ – מספר האיטרציות שעשינו עד כה

המטרה שלנו היא לחשב את $ \mathrm{vector} $ עבור $\theta$ שקרובה מספיק ל$x$ שלנו.

את זה אנחנו עושים באופן הבא: אנחנו מתחילים עם $ \mathrm{vector} = 1 $, $ \theta = 0 $ ו$ k = 0 $.
כעת אנחנו מבצעים את הפעולות הבאות עד שהקירוב מספיק טוב לנו:

  • אם $ x > \theta$:
    • $ \displaystyle{\theta := \theta + \arctan\frac{1}{2^k}} $
    • $ \displaystyle{\mathrm{vector} := \mathrm{vector} \cdot\left(1 + \frac{1}{2^k} \cdot i \right)} $
  • אחרת אם $ x < \theta$:
    • $ \displaystyle{\theta := \theta - \arctan\frac{1}{2^k}} $
    • $ \displaystyle{\mathrm{vector} := \mathrm{vector} \cdot \left(1 - \frac{1}{2^k} \cdot i \right)} $
  • $ k := k + 1 $

הסבר על מה שהולך כאן:

אם $ i $ היחידה המדומה, כלומר המספר המרוכב המקיים $ i^2 = -1 $, אזי מתקיים כי $ (\cos \alpha + i \sin \alpha) \cdot (\cos\beta + i \sin\beta) = \cos(\alpha + \beta) + i \sin(\alpha + \beta) $

כלומר – כפל של שני וקטורים במישור המרוכב הוא וקטור בכיוון סכום הכיוונים שלהם.

נשים לב של $ 1 + \frac{1}{2^k} \cdot i $ יש את אותו כיוון כמו ל$ \cos\arctan\frac{1}{2^k} + i \cdot \sin\arctan\frac{1}{2^k} $. (כי $ \frac{\frac{1}{2^k}}{1} = \frac{1}{2^k} $ ולכן יש להם אותו $ \arctan $)

מה שאנחנו בעצם עושים זה מנסים להתקרב לזווית x כמה שיותר ע”י הוספה או הסרה של$ \displaystyle{ \arctan\frac{1}{2^k} } $ לפי הצורך בכל שלב. בסוף אנחנו מקבלים וקטור עם זווית מאוד קרובה לx. כעת אנחנו יכולים לחשב את $ \sin $ ואת ה$ \cos $ שלו ע”י נירמולו – כלומר ע”י חלוקה בערך המוחלט שלו.

הבעיה בפתרון זה – היא שחילוק יכול להרוס לנו את הדיוק של הפעולות שביצענו עד כה. במקום זאת ניתן להשתמש בטריק הבא: נניח וביצענו אינסוף פעולות, אז קיבלנו ערך מוחלט של: $ \left( 1 + \frac{1}{2^0} \cdot i \right) \cdot \left( 1 + \frac{1}{2^1} \cdot i \right) \cdot \left( 1 + \frac{1}{2^2} \cdot i \right) \cdot \left( 1 + \frac{1}{2^3} \cdot i \right) \cdot \dots $

כלומר הערך המוחלט יהיה $ \sqrt{\left( 1 + \frac{1}{2^{2 \cdot 0}} \right) \cdot \left( 1 + \frac{1}{2^{2 \cdot 1}} \right) \cdot \left( 1 + \frac{1}{2^{2 \cdot 2}} \right) \cdot \left( 1 + \frac{1}{2^{2 \cdot 3}} \right) \cdot \dots} $.

בסוף היינו רוצים לחלק בו, כלומר להכפיל בהופכי שלו: $ \displaystyle{\frac{1}{\sqrt{\left( 1 + \frac{1}{2^{2 \cdot 0}} \right) \cdot \left( 1 + \frac{1}{2^{2 \cdot 1}} \right) \cdot \left( 1 + \frac{1}{2^{2 \cdot 2}} \right) \cdot \left( 1 + \frac{1}{2^{2 \cdot 3}} \right) \cdot \dots}}} $

מסתבר שזה מתכנס לקבוע $ K \approx 0.6072529350088812561694 $, כך שאנחנו יכולים לשמור אותו בצד (עם שאר הקבועים שחישבנו מראש) ובכל פעם להכפיל בו בסוף התהליך.

כאן יש איזשהו מימוש פשוט של CORDIC עם Complex (ראו גם טיפ מספר 128).

שימו לב שלא חישבתי את הקבועים מראש, אבל זה די אפשרי לעשות זאת:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
public static double Sin(double x)
{
    Complex vector = 1;
    double theta = 0;
    Complex oldVector = 0;
    
    for (int k = 0; oldVector != vector; k++)
    {
        oldVector = vector;
        double power = 1/Math.Pow(2, k); // Get it faster..
        double powerAngle = Math.Atan(power); // Get it from a hash..
        if (x > theta)
        {
            theta += powerAngle;
            vector *= (1 + power * Complex.ImaginaryOne);
        }
        else if (x < theta)
        {
            theta -= powerAngle;
            vector *= (1 - power * Complex.ImaginaryOne);
        }
    }
    vector /= Complex.Abs(vector); // Multiply by a constant instead.
    return vector.Imaginary;
}

היתרון של CORDIC על אלגוריתמים אחרים הוא שהוא מאוד מהיר בעולם בו אין כפל - אם תסתכלו על האלגוריתם, הוא למעשה משתמש רק בחיבור/חיסור, כפל בחזקה של 2 (שהוא מאוד זול – הוא Shift) וכפל חד פעמי בתוצאה. (אמנם נראה כאן שיש כפל, אבל אפשר לשים לב ש$ \left(1+ai\right)\left(1+bi\right)=1-ab+i\left(a+b\right) $ ואם a או b חזקות של 2 אז המכפלה היחידה שיש פה היא Shift)

תוכלו לקרוא על זה עוד כאן.

סופ"ש נומרי טוב!

שתף