محاسبهٔ ثابت‌های ریاضی

هیچ‌کدام از استانداردهای سی‌پلاس‌پلاس و یا سی عدد π را تعریف نکرده‌اند. بنابراین هر برنامه‌نویسی که قصد داشته باشه از توابع مثلثاتی یا آماری استفاده کنه مجبوره خودش π رو تعریف کنه. خوشبختانه فایل‌های سرآیند استاندارد عدد π رو تعریف کرده‌اند، (ثابت M_PI در هدر math.h رو ببینید) با این وجود استفاده از این ثابت و ثابت‌های دیگه بسیار خطرناک هست و برای کاربردهای دقیق باعث بروز خطاهای عددی خواهد شد.

خوشبختانه میشه با صرف هزینهٔ صفر در زمان اجرا π رو به‌صورت استاندارد (با تعاریف مشخص عددی) و البته در زمان کامپال محاسبه کرد. برای این کار از meta-programming به وسیلهٔ template ها در ‪C++‬ استفاده می‌کنیم. از مزایای این روش میشه به قابل حمل بودن و پیروی کامل از استاندارد IEEE 754 اشاره کرد که از لحاظ دقت عددی و یکسان بودن نتایج روی ماشین‌های مختلف بسیار مهم هست.

برای محاسبهٔ عدد π روش‌های مختلفی وجود داره. اکثر این روش‌ها براساس محاسبهٔ یک سری یا دنباله هستند که در بینهایت به عدد پی همگرا میشه. یکی از بهترین روش‌ها فرمول Bailey–Borwein–Plouffe هست. اثبات میشه که:

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

بنابراین با این کد توی ‪C++‬ می‌تونیم به‌سادگی عدد پی رو تا هر تعداد تکرار که میخواهیم حساب کنیم:

 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

#include <cstdint>
#include <limits>
#include <iostream>

template <uint64_t b, uint64_t e>
struct pow {
    static const uint64_t result = b * pow < b, e - 1 >::result;
};

template <uint64_t b>
struct pow<b, 0> {
    static const uint64_t result = 1;
};

template <uint64_t n>
struct bbp {
    constexpr static double pi
            = (1.0 / pow<16, n>::result)
            * (   4.0 / (8 * n + 1.0) - 2.0 / (8 * n + 4.0)
                - 1.0 / (8 * n + 5.0) - 1.0 / (8 * n + 6.0))
            + bbp < n - 1 >::pi;
};

template <>
struct bbp < -1 > {
    constexpr static double pi = 0;
};

int main(int, char* []) {
    std::cout.precision(std::numeric_limits<double>::digits10);
    std::cout << "pi = " << bbp<1>::pi << std::endl;
    std::cout << "pi = " << bbp<5>::pi << std::endl;
    std::cout << "pi = " << bbp<10>::pi << std::endl;
    return 0;
}

نتایج:

pi: 3.14142246642247
pi: 3.14159265322809
pi: 3.14159265358979
pi: 3.14159265358979

چرا؟

روش‌های زیادی برای محاسبهٔ عدد پی وجود داره. چرا دقیقاً این روش بهتره؟ مثلاً روش نیوتون به‌صورت هندسی عدد پی رو به این صورت حساب می‌کنه:

$$\pi=\frac{3}{4}\sqrt{3}+24\int_{0}^{\frac{1}{4}}{\sqrt{x-x^{2}}},dx$$

که خوب از نظر عددی محاسبهٔ پیچیده‌تری داره و زمان کامپایل نمیشه درست حسابش کرد. روش‌های دیگه‌ای هم وجود دارن مثل روش تبدیل همگرایی افزایشی اویلر (OEIS A054387) که اثبات می‌کنه:

$$\pi=\sum_{n=0}^{\infty}{\frac{\left({n!}^{2}2^{n-1}\right)}{\left(2n+1\right)! }}$$

با بررسی تمام این سری‌ها متوجه میشیم که خطای عددی افزایشی در فرمول BBP از تمام فرمول‌های پیوسته و سری‌های گسستهٔ دیگه کمتره.

ثابت‌های دیگر

عدد پی تنها عدد گنگ مورد استفاده در ریاضیات نیست. با این حال در زمینهٔ محاسبات علمی به اعداد دیگه مثل عدد e (پایهٔ لگاریتم طبیعی، عدد اویلر) خیلی کم نیاز پیدا می‌کنیم. (در این مورد خاص دو دلیل داره. یکی این که به‌راحتی قابل تبدیل به محاسبات مختلط هست و دوم این که تابع استاندارد exp در زبان برنامه‌نویسی سی وجود داره (: