如何一一计算非理性数字? [英] How to compute the digits of an irrational number one by one?

问题描述

我想逐位读取C中5的sqrt的小数. 5的平方根是2,23606797749979 ...，所以这是预期的输出:

I want to read digit by digit the decimals of the sqrt of 5 in C. The square root of 5 is 2,23606797749979..., so this'd be the expected output:

2
3
6
0
6
7
9
7
7
...

我发现

但是这种方法将结果存储在float变量中，并且我不想依赖float类型的限制，例如，我想提取10,000个数字.我还尝试使用本机sqrt()函数，并使用此方法将其强制转换为字符串号，但是我遇到了相同的问题问题.

But this approach stores the result in a float variable, and I don't want to depend on the limits of the float types, as I would like to extract like 10,000 digits, for instance. I also tried to use the native sqrt() function and casting it to string number using this method, but I faced the same issue.

推荐答案

如前所述，您需要将算法更改为一个数字一位( GMP ).

As already noted, you need to change the algorithm into a digit-by-digit one (there are some examples in the Wikipedia page about the methods of computing of the square roots) and use an arbitrary precision arithmetic library to perform the calculations (for instance, GMP).

在以下代码段中，我使用GMP(但未提供库提供的平方根函数)实现上述算法.此实现无需一次计算一个十进制数字，而是使用了一个更大的基数，该最大基数可以容纳unsigned long内的10的最大倍数，因此每次迭代可以生成9或18个十进制数字.

In the following snippet I implemented the before mentioned algorithm, using GMP (but not the square root function that the library provides). Instead of calculating one decimal digit at a time, this implementation uses a larger base, the greatest multiple of 10 that fits inside an unsigned long, so that it can produce 9 or 18 decimal digits at every iteration.

它还使用改编的牛顿方法来查找实际的数字".

It also uses an adapted Newton method to find the actual "digit".

#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);
}

您可以在此处看到输出的数字.

You can see the outputted digits here.

这篇关于如何一一计算非理性数字?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！