SPOJ `CUBERT` 中的错误答案
Wrong answer in SPOJ `CUBERT`
我在 SPOJ 上 this problem 的解决方案得到 错误答案。
该题要求计算一个整数（最多150位）的立方根，并输出截断至小数点后10位的答案。
它还要求计算答案模 10 中所有数字的总和作为 'checksum' 值。
这里是确切的问题陈述：
Your task is to calculate the cube root of a given positive integer.
We can not remember why exactly we need this, but it has something in
common with a princess, a young peasant, kissing and half of a kingdom
(a huge one, we can assure you).
Write a program to solve this crucial task.
Input
The input starts with a line containing a single integer t <= 20, the
number of test cases. t test cases follow.
The next lines consist of large positive integers of up to 150 decimal
digits. Each number is on its own separate line of the input file. The
input file may contain empty lines. Numbers can be preceded or
followed by whitespaces but no line exceeds 255 characters.
Output
For each number in the input file your program should output a line
consisting of two values separated by single space. The second value
is the cube root of the given number, truncated (not rounded!) after
the 10th decimal place. First value is a checksum of all printed
digits of the cube root, calculated as the sum of the printed digits
modulo 10.
Example
Input:
5
1
8
1000
2 33076161
Output:
1 1.0000000000
2 2.0000000000
1 10.0000000000
0 1.2599210498
6 321.0000000000
这是我的解决方案：
from math import pow
def foo(num):
num_cube_root = pow(num, 1.0 / 3)
# First round upto 11 decimal places
num_cube_root = "%.11f" % (num_cube_root)
# Then remove the last decimal digit
# to achieve a truncation of 10 decimal places
num_cube_root = str(num_cube_root)[0:-1]
num_cube_root_sum = 0
for digit in num_cube_root:
if digit != '.':
num_cube_root_sum += int(digit)
num_cube_root_sum %= 10
return (num_cube_root_sum, num_cube_root)
def main():
# Number of test cases
t = int(input())
while t:
t -= 1
num = input().strip()
# If line empty, ignore
if not num:
t += 1
continue
num = int(num)
ans = foo(num)
print(str(ans[0]) + " " + ans[1])
if __name__ == '__main__':
main()
它非常适合示例案例：Live demo。
谁能告诉我这个解决方案有什么问题？
您可以尝试使用精度值足够大的decimal模块。
编辑：感谢@DSM，我意识到 decimal 模块不会产生非常精确的立方根。我建议您检查是否所有数字都是 9，如果是，则将其四舍五入为整数。
此外，我现在也用 Decimals 执行 1/3 除法，因为将 1/3 的结果传递给 Decimal 构造函数会导致精度降低。
import decimal
def cbrt(n):
nd = decimal.Decimal(n)
with decimal.localcontext() as ctx:
ctx.prec = 50
i = nd ** (decimal.Decimal(1) / decimal.Decimal(3))
return i
ret = str(cbrt(1233412412430519230351035712112421123121111))
print(ret)
left, right = ret.split('.')
print(left + '.' + ''.join(right[:10]))
输出：
107243119477324.80328931501744819161741924145124146
107243119477324.8032893150
cbrt(10) 的输出是：
9.9999999999999999999999999999999999999999999999998
您的解决方案有两个问题，都与浮点运算的使用有关。第一个问题是 Python floats 只携带大约 16 位有效的小数位精度，所以只要你的答案需要超过 16 位左右的有效位（所以在点之前超过 6 位，和后面的 10 位数字），您几乎没有希望获得正确的尾随数字。第二个问题更微妙，甚至会影响 n 的小值。那就是您四舍五入到 11 位小数然后删除最后一位的方法会因 double rounding 而遭受潜在错误。例如，以 n = 33 为例。 n 的立方根，小数点后 20 位左右，是：
3.20753432999582648755...
四舍五入到该点后的 11 位时，您最终得到
3.20753433000
现在删除最后一位数字得到 3.2075343300，这不是您想要的。问题是舍入到小数点后 11 位最终会影响第 11 位数字左侧的数字。
那么你能做些什么来解决这个问题呢？好吧，您可以完全避免浮点数并将其简化为纯整数问题。我们需要一些整数 n 的立方根到小数点后 10 位（将最后一位向下舍入）。这相当于将 10**30 * n 的立方根计算为最接近的整数，再次向下舍入，然后将结果除以 10**10。所以这里的基本任务是计算任何给定整数 n 的立方根的底数。我无法找到任何关于计算整数立方根的现有 Stack Overflow 答案（在 Python 中更少），所以我认为值得详细说明如何这样做。
事实证明，计算整数的立方根非常容易（借助一点点数学知识）。有多种可能的方法，但一种既高效又易于实现的方法是使用 Newton-Raphson 方法的纯整数版本。在实数上，牛顿求解方程的方法 x**3 = n 对 n 的立方根进行近似 x，并迭代到 return 改进的近似。所需的迭代是：
x_next = (2*x + n/x**2)/3
在实际情况下，您会重复迭代，直到达到某个所需的公差。事实证明，在整数上，本质上相同的迭代有效，并且在正确的退出条件下，它会给我们 完全 正确的答案（不需要容差）。整数情况下的迭代是：
a_next = (2*a + n//a**2)//3
（注意使用底除运算符 // 代替上面通常的真正除法运算符 /。）从数学上讲，a_next 正好是 [=36 的底=].
下面是一些基于这次迭代的代码：
def icbrt_v1(n, initial_guess=None):
"""
Given a positive integer n, find the floor of the cube root of n.
Args:
n : positive integer
initial_guess : positive integer, optional. If given, this is an
initial guess for the floor of the cube root. It must be greater
than or equal to floor(cube_root(n)).
Returns:
The floor of the cube root of n, as an integer.
"""
a = initial_guess if initial_guess is not None else n
while True:
d = n//a**2
if a <= d:
return a
a = (2*a + d)//3
一些示例使用：
>>> icbrt_v1(100)
4
>>> icbrt_v1(1000000000)
1000
>>> large_int = 31415926535897932384626433
>>> icbrt_v1(large_int**3)
31415926535897932384626433
>>> icbrt_v1(large_int**3-1)
31415926535897932384626432
icbrt_v1 中存在一些问题和效率低下的问题，我们很快就会解决。但首先，简要解释一下为什么上面的代码有效。请注意，我们从假定大于或等于立方根底的初始猜测开始。我们将证明这个 属性 是一个循环不变量：每次 我们到达 while 循环的顶部时，a 至少是 floor(cbrt(n)) .此外，每次迭代产生的值 a 严格小于旧值，因此我们的迭代保证最终收敛到 floor(cbrt(n))。为了证明这些事实，请注意当我们进入 while 循环时，有两种可能性：
情况 1。a 严格大于 n 的立方根。然后a > n//a**2，代码继续下一次迭代。写a_next = (2*a + n//a**2)//3，则我们有：
a_next >= floor(cbrt(n))。这是因为 (2*a + n/a**2)/3 至少是 n 的立方根，而 n 又是 AM-GM inequality 应用于 a、a 和 n/a**2：这三个量的几何平均数恰好是n的立方根，所以算术平均数必须是至少[=22=的立方根].所以我们的循环不变量被保留用于下一次迭代。
a_next < a：由于我们假设 a 大于立方根 n/a**2 < a，因此 (2a + n/a**2) / 3 是比 a 小 ，因此 floor((2a + n/a**2) / 3) < a。这保证了我们在每次迭代中都朝着解决方案取得进展。
情况2。a小于或等于n的立方根。那么a <= floor(cbrt(n))，但是从上面建立的循环不变量我们也知道a >= floor(cbrt(n))。所以我们完成了：a 是我们追求的价值。 while 循环此时退出，因为 a <= n // a**2.
上面的代码有几个问题。首先，从 n 的初始猜测开始是低效的：代码将花费其前几次迭代（粗略地）将 a 的当前值除以 3 每次直到它进入解的邻域。初始猜测的更好选择（并且在 Python 中很容易计算）是使用超过 n.
的立方根的 2 的第一个幂
initial_guess = 1 << -(-n.bit_length() // 3)
更好的是，如果 n 足够小以避免溢出，则使用浮点运算来提供初始猜测，例如：
initial_guess = int(round(n ** (1/3.)))
但这给我们带来了第二个问题：我们算法的正确性要求初始猜测不小于实际整数立方根，并且随着 n 变大，我们不能保证对于上面基于浮点数的 initial_guess（尽管对于足够小的 n，我们可以）。幸运的是，有一个非常简单的修复方法：对于 any 正整数 a，如果我们执行一次迭代，我们总是会得到一个至少为 floor(cbrt(a)) 的值（使用我们上面使用的相同 AM-GM 参数）。所以我们所要做的就是在开始测试收敛性之前至少执行一次迭代。
考虑到这一点，下面是上述代码的更高效版本：
def icbrt(n):
"""
Given a positive integer n, find the floor of the cube root of n.
Args:
n : positive integer
Returns:
The floor of the cube root of n, as an integer.
"""
if n.bit_length() < 1024: # float(n) safe from overflow
a = int(round(n**(1/3.)))
a = (2*a + n//a**2)//3 # Ensure a >= floor(cbrt(n)).
else:
a = 1 << -(-n.bit_length()//3)
while True:
d = n//a**2
if a <= d:
return a
a = (2*a + d)//3
手头有 icbrt，很容易将所有内容放在一起计算立方根，精确到小数点后十位。在这里，为简单起见，我将结果输出为字符串，但您可以轻松构造一个 Decimal 实例。
def cbrt_to_ten_places(n):
"""
Compute the cube root of `n`, truncated to ten decimal places.
Returns the answer as a string.
"""
a = icbrt(n * 10**30)
q, r = divmod(a, 10**10)
return "{}.{:010d}".format(q, r)
示例输出：
>>> cbrt_to_ten_places(2)
'1.2599210498'
>>> cbrt_to_ten_places(8)
'2.0000000000'
>>> cbrt_to_ten_places(31415926535897932384626433)
'315536756.9301821867'
>>> cbrt_to_ten_places(31415926535897932384626433**3)
'31415926535897932384626433.0000000000'
我在 SPOJ 上 this problem 的解决方案得到 错误答案。
该题要求计算一个整数（最多150位）的立方根，并输出截断至小数点后10位的答案。
它还要求计算答案模 10 中所有数字的总和作为 'checksum' 值。
这里是确切的问题陈述：
Your task is to calculate the cube root of a given positive integer. We can not remember why exactly we need this, but it has something in common with a princess, a young peasant, kissing and half of a kingdom (a huge one, we can assure you).
Write a program to solve this crucial task.
Input
The input starts with a line containing a single integer t <= 20, the number of test cases. t test cases follow.
The next lines consist of large positive integers of up to 150 decimal digits. Each number is on its own separate line of the input file. The input file may contain empty lines. Numbers can be preceded or followed by whitespaces but no line exceeds 255 characters.
Output
For each number in the input file your program should output a line consisting of two values separated by single space. The second value is the cube root of the given number, truncated (not rounded!) after the 10th decimal place. First value is a checksum of all printed digits of the cube root, calculated as the sum of the printed digits modulo 10.
Example
Input:
5
18
1000
2 33076161
Output:
1 1.0000000000
2 2.0000000000
1 10.0000000000
0 1.2599210498
6 321.0000000000
这是我的解决方案：
from math import pow
def foo(num):
num_cube_root = pow(num, 1.0 / 3)
# First round upto 11 decimal places
num_cube_root = "%.11f" % (num_cube_root)
# Then remove the last decimal digit
# to achieve a truncation of 10 decimal places
num_cube_root = str(num_cube_root)[0:-1]
num_cube_root_sum = 0
for digit in num_cube_root:
if digit != '.':
num_cube_root_sum += int(digit)
num_cube_root_sum %= 10
return (num_cube_root_sum, num_cube_root)
def main():
# Number of test cases
t = int(input())
while t:
t -= 1
num = input().strip()
# If line empty, ignore
if not num:
t += 1
continue
num = int(num)
ans = foo(num)
print(str(ans[0]) + " " + ans[1])
if __name__ == '__main__':
main()
它非常适合示例案例：Live demo。
谁能告诉我这个解决方案有什么问题？
您可以尝试使用精度值足够大的decimal模块。
编辑：感谢@DSM，我意识到 decimal 模块不会产生非常精确的立方根。我建议您检查是否所有数字都是 9，如果是，则将其四舍五入为整数。
此外，我现在也用 Decimals 执行 1/3 除法，因为将 1/3 的结果传递给 Decimal 构造函数会导致精度降低。
import decimal
def cbrt(n):
nd = decimal.Decimal(n)
with decimal.localcontext() as ctx:
ctx.prec = 50
i = nd ** (decimal.Decimal(1) / decimal.Decimal(3))
return i
ret = str(cbrt(1233412412430519230351035712112421123121111))
print(ret)
left, right = ret.split('.')
print(left + '.' + ''.join(right[:10]))
输出：
107243119477324.80328931501744819161741924145124146
107243119477324.8032893150
cbrt(10) 的输出是：
9.9999999999999999999999999999999999999999999999998
您的解决方案有两个问题，都与浮点运算的使用有关。第一个问题是 Python floats 只携带大约 16 位有效的小数位精度，所以只要你的答案需要超过 16 位左右的有效位（所以在点之前超过 6 位，和后面的 10 位数字），您几乎没有希望获得正确的尾随数字。第二个问题更微妙，甚至会影响 n 的小值。那就是您四舍五入到 11 位小数然后删除最后一位的方法会因 double rounding 而遭受潜在错误。例如，以 n = 33 为例。 n 的立方根，小数点后 20 位左右，是：
3.20753432999582648755...
四舍五入到该点后的 11 位时，您最终得到
3.20753433000
现在删除最后一位数字得到 3.2075343300，这不是您想要的。问题是舍入到小数点后 11 位最终会影响第 11 位数字左侧的数字。
那么你能做些什么来解决这个问题呢？好吧，您可以完全避免浮点数并将其简化为纯整数问题。我们需要一些整数 n 的立方根到小数点后 10 位（将最后一位向下舍入）。这相当于将 10**30 * n 的立方根计算为最接近的整数，再次向下舍入，然后将结果除以 10**10。所以这里的基本任务是计算任何给定整数 n 的立方根的底数。我无法找到任何关于计算整数立方根的现有 Stack Overflow 答案（在 Python 中更少），所以我认为值得详细说明如何这样做。
事实证明，计算整数的立方根非常容易（借助一点点数学知识）。有多种可能的方法，但一种既高效又易于实现的方法是使用 Newton-Raphson 方法的纯整数版本。在实数上，牛顿求解方程的方法 x**3 = n 对 n 的立方根进行近似 x，并迭代到 return 改进的近似。所需的迭代是：
x_next = (2*x + n/x**2)/3
在实际情况下，您会重复迭代，直到达到某个所需的公差。事实证明，在整数上，本质上相同的迭代有效，并且在正确的退出条件下，它会给我们 完全 正确的答案（不需要容差）。整数情况下的迭代是：
a_next = (2*a + n//a**2)//3
（注意使用底除运算符 // 代替上面通常的真正除法运算符 /。）从数学上讲，a_next 正好是 [=36 的底=].
下面是一些基于这次迭代的代码：
def icbrt_v1(n, initial_guess=None):
"""
Given a positive integer n, find the floor of the cube root of n.
Args:
n : positive integer
initial_guess : positive integer, optional. If given, this is an
initial guess for the floor of the cube root. It must be greater
than or equal to floor(cube_root(n)).
Returns:
The floor of the cube root of n, as an integer.
"""
a = initial_guess if initial_guess is not None else n
while True:
d = n//a**2
if a <= d:
return a
a = (2*a + d)//3
一些示例使用：
>>> icbrt_v1(100)
4
>>> icbrt_v1(1000000000)
1000
>>> large_int = 31415926535897932384626433
>>> icbrt_v1(large_int**3)
31415926535897932384626433
>>> icbrt_v1(large_int**3-1)
31415926535897932384626432
icbrt_v1 中存在一些问题和效率低下的问题，我们很快就会解决。但首先，简要解释一下为什么上面的代码有效。请注意，我们从假定大于或等于立方根底的初始猜测开始。我们将证明这个 属性 是一个循环不变量：每次 我们到达 while 循环的顶部时，a 至少是 floor(cbrt(n)) .此外，每次迭代产生的值 a 严格小于旧值，因此我们的迭代保证最终收敛到 floor(cbrt(n))。为了证明这些事实，请注意当我们进入 while 循环时，有两种可能性：
情况 1。a 严格大于 n 的立方根。然后a > n//a**2，代码继续下一次迭代。写a_next = (2*a + n//a**2)//3，则我们有：
a_next >= floor(cbrt(n))。这是因为(2*a + n/a**2)/3至少是n的立方根，而n又是 AM-GM inequality 应用于a、a和n/a**2：这三个量的几何平均数恰好是n的立方根，所以算术平均数必须是至少[=22=的立方根].所以我们的循环不变量被保留用于下一次迭代。a_next < a：由于我们假设a大于立方根n/a**2 < a，因此(2a + n/a**2) / 3是比a小 ，因此floor((2a + n/a**2) / 3) < a。这保证了我们在每次迭代中都朝着解决方案取得进展。
情况2。a小于或等于n的立方根。那么a <= floor(cbrt(n))，但是从上面建立的循环不变量我们也知道a >= floor(cbrt(n))。所以我们完成了：a 是我们追求的价值。 while 循环此时退出，因为 a <= n // a**2.
上面的代码有几个问题。首先，从 n 的初始猜测开始是低效的：代码将花费其前几次迭代（粗略地）将 a 的当前值除以 3 每次直到它进入解的邻域。初始猜测的更好选择（并且在 Python 中很容易计算）是使用超过 n.
initial_guess = 1 << -(-n.bit_length() // 3)
更好的是，如果 n 足够小以避免溢出，则使用浮点运算来提供初始猜测，例如：
initial_guess = int(round(n ** (1/3.)))
但这给我们带来了第二个问题：我们算法的正确性要求初始猜测不小于实际整数立方根，并且随着 n 变大，我们不能保证对于上面基于浮点数的 initial_guess（尽管对于足够小的 n，我们可以）。幸运的是，有一个非常简单的修复方法：对于 any 正整数 a，如果我们执行一次迭代，我们总是会得到一个至少为 floor(cbrt(a)) 的值（使用我们上面使用的相同 AM-GM 参数）。所以我们所要做的就是在开始测试收敛性之前至少执行一次迭代。
考虑到这一点，下面是上述代码的更高效版本：
def icbrt(n):
"""
Given a positive integer n, find the floor of the cube root of n.
Args:
n : positive integer
Returns:
The floor of the cube root of n, as an integer.
"""
if n.bit_length() < 1024: # float(n) safe from overflow
a = int(round(n**(1/3.)))
a = (2*a + n//a**2)//3 # Ensure a >= floor(cbrt(n)).
else:
a = 1 << -(-n.bit_length()//3)
while True:
d = n//a**2
if a <= d:
return a
a = (2*a + d)//3
手头有 icbrt，很容易将所有内容放在一起计算立方根，精确到小数点后十位。在这里，为简单起见，我将结果输出为字符串，但您可以轻松构造一个 Decimal 实例。
def cbrt_to_ten_places(n):
"""
Compute the cube root of `n`, truncated to ten decimal places.
Returns the answer as a string.
"""
a = icbrt(n * 10**30)
q, r = divmod(a, 10**10)
return "{}.{:010d}".format(q, r)
示例输出：
>>> cbrt_to_ten_places(2)
'1.2599210498'
>>> cbrt_to_ten_places(8)
'2.0000000000'
>>> cbrt_to_ten_places(31415926535897932384626433)
'315536756.9301821867'
>>> cbrt_to_ten_places(31415926535897932384626433**3)
'31415926535897932384626433.0000000000'