【问题标题】:How to compute the digits of an irrational number one by one?如何逐个计算无理数的位数?
【发布时间】:2020-02-20 10:42:58
【问题描述】:

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

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

我找到了the following code:

#include<stdio.h>

void main()
{
    int number;

    float temp, sqrt;

    printf("Provide the number: \n");

    scanf("%d", &number);

    // store the half of the given number e.g from 256 => 128
    sqrt = number / 2;
    temp = 0;

    // Iterate until sqrt is different of temp, that is updated on the loop
    while(sqrt != temp){
        // initially 0, is updated with the initial value of 128
        // (on second iteration = 65)
        // and so on
        temp = sqrt;

        // Then, replace values (256 / 128 + 128 ) / 2 = 65
        // (on second iteration 34.46923076923077)
        // and so on
        sqrt = ( number/temp + temp) / 2;
    }

    printf("The square root of '%d' is '%f'", number, sqrt);
}

但是这种方法将结果存储在一个浮点变量中,我不想依赖浮点类型的限制,例如我想提取 10,000 位数字。我还尝试使用本机 sqrt() 函数并使用 this method 将其转换为字符串编号,但我遇到了同样的问题。

【问题讨论】:

  • 你需要一个Arbitrary precision arithmetic
  • 你有两种选择:1)将数字转换为字符串,将字符串中的每个字符一个一个打印出来;或者 2) 使用十进制算术一个一个地得到每个数字(整数截断和与 10 的乘法就是你所需要的)。
  • 在您的问题中,您声明您想逐位“读取”sqrt(5),但在我看来,您真的想编写一个程序来“计算”这个数字,然后按数字打印按数字。请注意,计算机通常可以使用固定数量的有效数字,并且您需要一个特殊的库才能根据需要使用尽可能多的小数位 - 正如@pmg 所指出的那样。 -- 计算机通常固定为一组数字以节省内存,或者至少使其易于管理,以便计算机知道为每个数字分配多少内存......
  • 您的方法近似平方根,直到达到浮点的最大精度,即。直到sqrt 不再改变。要获得更高的精度,您将需要另一种方法。
  • 不可能无限期地使用有限位数来计算 sqrt(5) 的位数进行计算,因为有限位数只有有限数量的状态,所以行为必须开始重复在某些时候,因此必须产生一个重复的十进制数字。重复的十进制数字是有理数,但 sqrt(5) 是无理数。因此,一个“永远”打印 sqrt(5) 数字的程序必须使用任意精度的算法,随着时间的推移使用越来越多的内存(这意味着它也必须在有限的世界中耗尽它的能力)。

标签: c math casting floating-point sqrt


【解决方案1】:

您所问的是一个非常困难的问题,以及是否有可能“一个接一个”地进行(即没有工作空间要求会随着您想走多远而扩展) 取决于特定的无理数和您希望它表示的基数。例如,在 1995 年 formula for pi was discovered that allows computing the nth binary digit in O(1) space 时,这真是一件大事。这不是人们期望的事情。

如果你愿意接受 O(n) 空间,那么像你提到的那种情况是相当容易的。例如,如果您将数字平方根的前 n 位作为十进制字符串,您可以简单地尝试将每个数字附加 0 到 9,然后用长乘法对字符串进行平方(与您在小学时学到的相同),并选择最后一个不会超调的。当然,这很慢,但很简单。使它更快(但仍然渐近同样糟糕)的简单方法是使用任意精度的数学库代替字符串。做得更好需要更先进的方法,一般来说可能是不可能的。

【讨论】:

  • 在“实际复杂度”中可能是 O(1) 空间(假设事物永远不会溢出某些基本字长),但在理论复杂度上不可能是 O(1);您无法在 O(1) 空间中使用任意输入 n 进行计算。
  • @EricPostpischil:在transdichotomous model
  • 嗯,我不确定。在维基百科解释的跨二分模型中,输入值完全独立于输入的数量(n),除了它规定字长至少为 log[2](n),我不确定这是否正确. (当然,至少必须这样才能使列表可以包含要排序的 n 个不同项目,这可能是该规定的唯一目的,但实际上我们需要足够的位来处理实际的输入值。)在使用 Bailey-Borwein-Plouffe 和来计算 π 的数字,k 上的和,其中 k 受 n 影响。
  • @EricPostpischil:此问题中的输入大小是表示 n 的位数,其中 n 是所需的数字位置。跨二分模型只是意味着你的机器的字长足以表示这样的 n。如果是这种情况,那么您可以在 O(1) 空间中计算结果(所需的数字)。
  • 只是一个小调整:'附加每个数字' - 我会以二进制搜索方式进行,不过:从 5 开始,如果平方更大,继续 2,否则7、...
【解决方案2】:

如前所述,您需要将算法更改为逐位数的算法(Wikipedia page about the methods of computing of the square roots 中有一些示例)并使用任意精度的算术库来执行计算(例如,GMP )。

在下面的 sn-p 中,我使用 GMP 实现了前面提到的算法(但不是库提供的平方根函数)。此实现不是一次计算一个十进制数字,而是使用更大的基数,即适合 unsigned long 的 10 的最大倍数,因此它可以在每次迭代中生成 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);
}

可以看到输出的数字here

【讨论】:

  • 这太棒了!我也尝试过 100K 并且工作正常!您能否在代码中添加一些 cmets?谢谢!
  • @harrison4 我会的,当然,但今天晚些时候。
【解决方案3】:

你的标题说:

如何逐个计算无理数的位数?

无理数不限于大多数平方根。它们还包括log(x)exp(z)sin(y) 等形式的数字(超越数字)。但是,有一些重要因素决定了您是否可以或多快地逐一计算给定无理数的位数(即从左到右)。

  • 并非所有无理数都是可计算的;也就是说,没有人找到将它们近似为任何所需长度的方法(无论是通过封闭形式的表达式、系列还是其他方式)。
  • 可以通过多种方式表示数字,例如通过二进制或十进制扩展、连分数、级数等。根据表示形式,有不同的算法可以计算给定数字的位数。
  • 某些公式以特定基数(例如基数 2)而不是任意基数计算给定数字的位数。

例如,除了第一个公式可以在不计算前面的数字的情况下提取 π 的数字,还有其他这种类型的公式(称为BBP-type formulas)可以提取某些无理数的数字。然而,这些公式只适用于特定的基数,并非所有 BBP 类型的公式都有正式证明,最重要的是,并非所有无理数都有 BBP 类型的公式(本质上,只有某些对数和反正切常数,而不是数字形式为exp(x)sqrt(x))。

另一方面,如果你可以将一个无理数表示为一个连分数(所有实数都有),你可以从左到右提取它的数字,并且在任何所需的基础上,使用特定的算法。更重要的是,该算法适用于任何实数常数,包括平方根、指数(eexp(x))、对数等,只要您知道如何将其表示为连分数即可。有关实现,请参见“Digits of pi and Python generators”。另见Code to Generate e one Digit at a Time

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 2011-05-27
    • 2011-01-20
    • 1970-01-01
    • 1970-01-01
    • 2016-08-17
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多