如何一位一位地计算无理数的数字？
How to compute the digits of an irrational number one by one?
我想在 C 中逐位读取 5 的开方的小数。
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 将其转换为字符串数字，但我遇到了同样的问题。
你问的是一个非常难的问题，是否有可能做到"one by one"（即不工作space要求与你想走多远成比例）取决于特定的无理数 和你希望它表示的基数 。例如，在 1995 年 formula for pi was discovered that allows computing the nth binary digit in O(1) space，这个真的很重要。这不是人们期望的事情。
如果您愿意接受 O(n) space，那么您提到的某些情况就相当容易了。例如，如果你有一个数字的平方根的前 n 位数字作为十进制字符串，你可以简单地尝试将每个数字附加 0 到 9，然后用长乘法对字符串进行平方（与你在小学学到的一样），并选择最后一个不会超调的。当然这很慢，但是很简单。使它更快（但仍然渐近地同样糟糕）的简单方法是使用任意精度的数学库代替字符串。做得更好需要更先进的方法，一般来说可能是不可能的。
如前所述，您需要将算法更改为逐位算法（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 个小数位。
它还使用了一种改进的牛顿法来找到实际的 "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);
}
您可以看到输出的数字here。
您的标题是：
How to compute the digits of an irrational number one by one?
无理数不限于大多数平方根。它们还包括 log(x)、exp(z)、sin(y) 等形式的数字（超越数）。但是，有一些重要因素决定了您能否或多快地逐一（即从左到右）计算给定的无理数的数字。
- 并非所有无理数都是可计算的；也就是说，没有人找到一种方法将它们近似为任何所需的长度（无论是通过封闭形式表达式、级数还是其他方式）。
- 数字的表示方式有很多种，例如通过二进制或十进制展开式、连续分数式、级数等。根据表示形式，可以使用不同的算法来计算给定数字的位数。
- 一些公式计算给定数字在特定基数（例如基数 2）中的数字，而不是在任意基数中。
例如，除了第一个公式不计算前面的数字来提取π的数字外，还有其他此类公式（称为BBP-type formulas）提取某些无理数的数字。然而，这些公式只适用于特定的基数，并非所有 BBP-type 公式都有正式证明，最重要的是，并非所有无理数都有 BBP-type 公式（本质上，只有某些 log 和 arctan 常数做，而不是 exp(x) 或 sqrt(x)).
形式的数字
另一方面，如果您可以将无理数表示为 连分数（所有实数都有），您可以从左到右提取它的数字，并且在任何需要的基础上，使用特定的算法。此外，该算法适用于任何实数常数，包括平方根、指数（e 和 exp(x)）、对数等，只要您知道如何将其表示为连分数即可.有关实现，请参阅“Digits of pi and Python generators". See also Code to Generate e one Digit at a Time.
我想在 C 中逐位读取 5 的开方的小数。 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 将其转换为字符串数字，但我遇到了同样的问题。
你问的是一个非常难的问题，是否有可能做到"one by one"（即不工作space要求与你想走多远成比例）取决于特定的无理数 和你希望它表示的基数 。例如，在 1995 年 formula for pi was discovered that allows computing the nth binary digit in O(1) space，这个真的很重要。这不是人们期望的事情。
如果您愿意接受 O(n) space，那么您提到的某些情况就相当容易了。例如，如果你有一个数字的平方根的前 n 位数字作为十进制字符串，你可以简单地尝试将每个数字附加 0 到 9，然后用长乘法对字符串进行平方（与你在小学学到的一样），并选择最后一个不会超调的。当然这很慢，但是很简单。使它更快（但仍然渐近地同样糟糕）的简单方法是使用任意精度的数学库代替字符串。做得更好需要更先进的方法，一般来说可能是不可能的。
如前所述，您需要将算法更改为逐位算法（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 个小数位。
它还使用了一种改进的牛顿法来找到实际的 "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);
}
您可以看到输出的数字here。
您的标题是：
How to compute the digits of an irrational number one by one?
无理数不限于大多数平方根。它们还包括 log(x)、exp(z)、sin(y) 等形式的数字（超越数）。但是，有一些重要因素决定了您能否或多快地逐一（即从左到右）计算给定的无理数的数字。
- 并非所有无理数都是可计算的；也就是说，没有人找到一种方法将它们近似为任何所需的长度（无论是通过封闭形式表达式、级数还是其他方式）。
- 数字的表示方式有很多种，例如通过二进制或十进制展开式、连续分数式、级数等。根据表示形式，可以使用不同的算法来计算给定数字的位数。
- 一些公式计算给定数字在特定基数（例如基数 2）中的数字，而不是在任意基数中。
例如，除了第一个公式不计算前面的数字来提取π的数字外，还有其他此类公式（称为BBP-type formulas）提取某些无理数的数字。然而，这些公式只适用于特定的基数，并非所有 BBP-type 公式都有正式证明，最重要的是，并非所有无理数都有 BBP-type 公式（本质上，只有某些 log 和 arctan 常数做，而不是 exp(x) 或 sqrt(x)).
另一方面，如果您可以将无理数表示为 连分数（所有实数都有），您可以从左到右提取它的数字，并且在任何需要的基础上，使用特定的算法。此外，该算法适用于任何实数常数，包括平方根、指数（e 和 exp(x)）、对数等，只要您知道如何将其表示为连分数即可.有关实现，请参阅“Digits of pi and Python generators". See also Code to Generate e one Digit at a Time.