342. Pi with the BBP formula

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

לכבוד יום הפאי – נדבר היום על נוסחת BBP לחישוב פאי.

קיימות לא מעט נוסחאות המחשבות את פאי המבוססות על אנליזה מתמטית.

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

ב1995 גילה בחור בשם Simon Plouffe באמצעות תוכנת מחשב שכתב, נוסחה מעט מסובכת לחישוב פאי:

$ \displaystyle{\sum_{k = 0}^{\infty}\left[ \frac{1}{16^k} \left( \frac{4}{8k + 1} - \frac{2}{8k + 4} - \frac{1}{8k + 5} - \frac{1}{8k + 6} \right) \right]} $

היתרון הגדול של הנוסחה הוא שהיא מאפשרת לחשב את הספרה הnית של פאי אחרי הנקודה ללא חישוב קודמותיה.

אבל יש פה קאץ’. הקאץ’ הוא שנקבל את הספרה הnית של הייצוג של פאי בבסיס 16. (הקסה-דצימלי)

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

לרגל יום הפאי, יכתב כאן קוד שמחשב את הספרה הnית של פאי אחרי הנקודה בבסיס 16 לפי נוסחת BBP. שלבי פיתוח הנוסחה (לא כל כך מסובכים) נמצאים כאן.

הקוד:

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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
public static char PiDigit(int digitIndex)
{
    digitIndex--;
    double approximation =
        4 * PiSeries(1, digitIndex) -
        2 * PiSeries(4, digitIndex) -
        PiSeries(5, digitIndex) -
        PiSeries(6, digitIndex);
    approximation = approximation % 1.0;
    approximation = (1.0 + approximation) % 1.0;
    long digit = (long)(16 * approximation) % 16;
    return digit.ToString("X").First();
}
public static double PiSeries(int remainder, int digitIndex)
{
    return PiFiniteSeries(remainder, digitIndex) + PiInfiniteSeries(remainder, digitIndex);
}
private static double PiFiniteSeries(int remainder, int digitIndex)
{
    double sum = 0;
    for (int i = 0; i <= digitIndex; i++)
    {
        long currentNum = 8 * i + remainder;
        sum += PowerModulo(16, digitIndex - i, currentNum) / (double)currentNum;
        sum = sum % 1.0;
    }
    return sum;
}
private static double PiInfiniteSeries(int remainder, int digitIndex)
{
    double sum = 0;
    double previousSum = 1;
    for (long i = digitIndex + 1; sum != previousSum; i++)
    {
        long currentNum = 8 * i + remainder;
        previousSum = sum;
        sum += Math.Pow(16, digitIndex - i) / currentNum;
    }
    return sum;
}
public static long PowerModulo(long a, long n, long m)
{
    long result = 1;
    long currentPower = a;
    while (n > 0)
    {
        if (n % 2 == 1)
        {
            result *= currentPower;
            result = result % m;
        }
        currentPower *= currentPower;
        currentPower = currentPower % m;
        n /= 2;
    }
    return result;
}

בגדול: מנסים לחשב את הספרה הראשונה אחרי הנקודה של $ \pi \cdot 16^n $. את זה עושים ע"י חישוב שני טורים: אחד סופי ואחד אינסופי. בדרך מתעלמים מהחלק השלם של כל המספרים (שהרי מעניין אותנו רק הספרה הראשונה אחרי הנקודה של המספר)

שימו לב שגם יש פה אלגוריתם חישוב חזקה מודולו מספר שדומה לאלגוריתם שראינו פעם (ע"י פירוק המעריך לבינארית)

היתרון העיקרי של הנוסחה היא שאפשר לחשב ספרה בלי להיסחב עם כל הספרות האחרות.

המשך יום פאי שמח!

שתף