จะคำนวณตัวเลขของจำนวนอตรรกยะทีละตัวได้อย่างไร?

ตามที่ระบุไว้แล้ว คุณต้องเปลี่ยนอัลกอริทึมเป็นตัวเลขต่อหลัก (มีตัวอย่างบางส่วนใน หน้าวิกิพีเดียเกี่ยวกับวิธีการคำนวณรากที่สอง) และใช้ไลบรารีเลขคณิตที่มีความแม่นยำตามอำเภอใจเพื่อทำการคำนวณ (เช่น GMP)

ในตัวอย่างต่อไปนี้ ฉันใช้อัลกอริธึมที่กล่าวถึงก่อนหน้านี้ โดยใช้ GMP (แต่ไม่ใช่ฟังก์ชันสแควร์รูทที่ไลบรารีจัดเตรียมให้) แทนที่จะคำนวณทศนิยมทีละหลัก การใช้งานนี้ใช้ฐานที่ใหญ่กว่า ซึ่งเป็นพหุคูณที่ยิ่งใหญ่ที่สุดของ 10 ที่พอดีกับ unsigned long เพื่อให้สามารถสร้างทศนิยม 9 หรือ 18 หลักในแต่ละครั้งของการวนซ้ำ

นอกจากนี้ยังใช้วิธีนิวตันที่ดัดแปลงเพื่อค้นหา "ตัวเลข" ที่แท้จริง

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gmp.h>

unsigned long max_ul(unsigned long a, unsigned long b)
{
    return a < b ? b : a;   
}

int main(int argc, char *argv[])
{
    // The GMP functions accept 'unsigned long int' values as parameters.
    // The algorithm implemented here can work with bases other than 10,
    // so that it can evaluate more than one decimal digit at a time.
    const unsigned long base = sizeof(unsigned long) > 4
                             ? 1000000000000000000
                             : 1000000000;
    const unsigned long decimals_per_digit = sizeof(unsigned long) > 4 ? 18 : 9;

    // Extract the number to be square rooted and the desired number of decimal
    // digits from the command line arguments. Fallback to 0 in case of errors.
    const unsigned long number = argc > 1 ? atoi(argv[1]) : 0;
    const unsigned long n_digits = argc > 2 ? atoi(argv[2]) : 0;

    // All the variables used by GMP need to be properly initialized before use.
    // 'c' is basically the remainder, initially set to the original number
    mpz_t c;
    mpz_init_set_ui(c, number);

    // At every iteration, the algorithm "move to the left" by two "digits"
    // the reminder, so it multplies it by base^2.
    mpz_t base_squared;
    mpz_init_set_ui(base_squared, base);
    mpz_mul(base_squared, base_squared, base_squared);

    // 'p' stores the digits of the root found so far. The others are helper variables
    mpz_t p;
    mpz_init_set_ui(p, 0UL);    
    mpz_t y;
    mpz_init(y);
    mpz_t yy;
    mpz_init(yy);
    mpz_t dy;
    mpz_init(dy);
    mpz_t dx;
    mpz_init(dx);
    mpz_t pp;    
    mpz_init(pp);

    // Timing, for testing porpuses
    clock_t start = clock(), diff;

    unsigned long x_max = number;
    // Each "digit" correspond to some decimal digits
    for (unsigned long i = 0,
         last = (n_digits + decimals_per_digit) / decimals_per_digit + 1UL;
         i < last; ++i)
    {
        // Find the greatest x such that:  x * (2 * base * p + x) <= c
        // where x is in [0, base), using a specialized Newton method

        // pp = 2 * base * p
        mpz_mul_ui(pp, p, 2UL * base);

        unsigned long x = x_max;
        for (;;)
        {            
            // y = x * (pp + x)
            mpz_add_ui(yy, pp, x);
            mpz_mul_ui(y, yy, x);

            // dy = y - c
            mpz_sub(dy, y, c);

            // If y <= c we have found the correct x
            if ( mpz_sgn(dy) <= 0 )
                break;

            // Newton's step:  dx = dy/y'  where  y' = 2 * x + pp            
            mpz_add_ui(yy, yy, x);
            mpz_tdiv_q(dx, dy, yy);

            // Update x even if dx == 0 (last iteration)
            x -= max_ul(mpz_get_si(dx), 1);
        }        
        x_max = base - 1;

        // The actual format of the printed "digits" is up to you       
        if (i % 4 == 0)
        {
            if (i == 0)
                printf("%lu.", x);
            putchar('\n');
        }
        else
            printf("%018lu", x);

        // p = base * p + x
        mpz_mul_ui(p, p, base);
        mpz_add_ui(p, p, x);

        // c = (c - y) * base^2
        mpz_sub(c, c, y);
        mpz_mul(c, c, base_squared);
    }

    diff = clock() - start;
    long int msec = diff * 1000L / CLOCKS_PER_SEC;
    printf("\n\nTime taken: %ld.%03ld s\n", msec / 1000, msec % 1000);

    // Final cleanup
    mpz_clear(c);
    mpz_clear(base_squared);
    mpz_clear(p);
    mpz_clear(pp);
    mpz_clear(dx);
    mpz_clear(y);
    mpz_clear(dy);
    mpz_clear(yy);
}

คุณสามารถดูตัวเลขที่ส่งออกได้ที่นี่

สิ่งที่คุณถามคือปัญหายากมาก และเป็นไปได้หรือไม่ที่จะทำ "ทีละคน" (เช่น โดยไม่ต้องใช้พื้นที่ทำงานซึ่งจะปรับขนาดตามระยะทางที่คุณต้องการจะไป) ขึ้นอยู่กับ ทั้งจำนวนอตรรกยะเฉพาะ และฐาน ที่คุณต้องการให้แสดง ตัวอย่างเช่น ในปี 1995 เมื่อ สูตรสำหรับ pi ถูกค้นพบซึ่งช่วยให้สามารถคำนวณเลขฐานสองหลักที่ n ในช่องว่าง O(1) นี่เป็นเรื่องใหญ่จริงๆ มันไม่ใช่สิ่งที่ผู้คนคาดหวังว่าจะเป็นไปได้

หากคุณยินดีที่จะยอมรับช่องว่าง O(n) บางกรณีเช่นกรณีที่คุณพูดถึงนั้นค่อนข้างง่าย ตัวอย่างเช่น หากคุณมี n หลักแรกของรากที่สองของตัวเลขเป็นสตริงทศนิยม คุณสามารถลองเติมเลข 0 ถึง 9 แต่ละหลักต่อท้าย จากนั้นยกกำลังสองสตริงด้วยการคูณยาว (แบบเดียวกับที่คุณเรียนในชั้นประถม) และเลือกอันสุดท้ายที่ไม่หักโหม แน่นอนว่ามันช้ามาก แต่ก็ง่าย วิธีง่ายๆ ในการทำให้มันเร็วขึ้นมาก (แต่ยังคงไม่มีซีมโตติคัลเหมือนเดิม) คือการใช้ไลบรารีคณิตศาสตร์ที่มีความแม่นยำตามอำเภอใจแทนสตริง การทำสิ่งที่ดีกว่าอย่างมีนัยสำคัญต้องใช้แนวทางขั้นสูงกว่า และโดยทั่วไปอาจเป็นไปไม่ได้

ชื่อของคุณพูดว่า:

จะคำนวณตัวเลขของจำนวนอตรรกยะทีละตัวได้อย่างไร?

จำนวนอตรรกยะไม่ได้จำกัดอยู่เพียงรากที่สองส่วนใหญ่ นอกจากนี้ยังรวมถึงตัวเลขในรูปแบบ log(x), exp(z), sin(y) ฯลฯ (ตัวเลขทิพย์) อย่างไรก็ตาม มีปัจจัยสำคัญบางประการที่กำหนดว่าคุณสามารถคำนวณเลขจำนวนอตรรกยะที่กำหนดทีละตัวได้เร็วหรือเร็วแค่ไหน (จากซ้ายไปขวา)

• ไม่สามารถคำนวณจำนวนอตรรกยะทั้งหมดได้ นั่นคือ ไม่มีใครพบวิธีที่จะประมาณความยาวที่ต้องการได้ (ไม่ว่าจะโดยการแสดงออกในรูปแบบปิด อนุกรม หรืออย่างอื่น)
• มีวิธีการแสดงตัวเลขได้หลายวิธี เช่น การขยายเลขฐานสองหรือทศนิยม เป็นเศษส่วนต่อเนื่อง เป็นอนุกรม เป็นต้น และมีอัลกอริธึมที่แตกต่างกันในการคำนวณหลักของตัวเลขที่กำหนด ขึ้นอยู่กับการแทนค่า
• บางสูตรคำนวณตัวเลขที่กำหนดในฐานเฉพาะ (เช่น ฐาน 2) ไม่ใช่ฐานที่กำหนดเอง

ตัวอย่างเช่น นอกจากสูตรแรกที่จะแยกเลข π โดยไม่ต้องคำนวณเลขหลักก่อนหน้าแล้ว ยังมีสูตรอื่นๆ ประเภทนี้ (เรียกว่า สูตรประเภท BBP) ที่แยกตัวเลขของจำนวนอตรรกยะบางจำนวน อย่างไรก็ตาม สูตรเหล่านี้ใช้ได้กับฐานใดฐานหนึ่งเท่านั้น ไม่ใช่ทุกสูตรประเภท BBP ที่มีการพิสูจน์อย่างเป็นทางการ และที่สำคัญที่สุด ไม่ใช่จำนวนอตรรกยะทั้งหมดที่มีสูตรประเภท BBP (โดยพื้นฐานแล้ว เฉพาะค่าคงที่ล็อกและอาร์คแทนบางค่าเท่านั้นที่ใช้ได้ ไม่ใช่ตัวเลขของ แบบฟอร์ม exp(x) หรือ sqrt(x))

ในทางกลับกัน หากคุณสามารถแสดงจำนวนอตรรกยะเป็น เศษส่วนต่อเนื่อง (ซึ่งจำนวนจริงทั้งหมดมี) คุณสามารถแยกตัวเลขจากซ้ายไปขวา และในฐานใดก็ได้ที่ต้องการ โดยใช้ค่าเฉพาะเจาะจง อัลกอริทึม ยิ่งไปกว่านั้น อัลกอริทึมนี้ยังใช้ได้กับค่าคงที่จำนวนจริงใดๆ รวมถึงรากที่สอง เลขชี้กำลัง (e และ exp(x)) ลอการิทึม ฯลฯ ตราบใดที่คุณรู้วิธีแสดงมันเป็นเศษส่วนต่อเนื่อง สำหรับการใช้งาน โปรดดูที่ Digits of pi และ Python เครื่องกำเนิดไฟฟ้า. ดูเพิ่มเติมที่โค้ดสำหรับสร้าง e ทีละหลัก.