Как вычислить цифры иррационального числа одну за другой?

Как уже отмечалось, алгоритм нужно изменить на поразрядный (примеры есть в страница Википедии о методах вычисления квадратных корней) и использовать арифметическую библиотеку произвольной точности для выполнения вычислений (например, 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 году, когда была обнаружена формула для числа пи, позволяющая вычислять 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)), логарифмы и т. д., если вы знаете, как выразить их в виде непрерывной дроби. Для реализации см. Цифры числа пи и Python генераторы. См. также Код для генерации одной цифры за раз.