精确的n次方根

我正在寻找Python第n个根函数/算法，但在你发布之前:NO INTEGER root, HELL!
我在哪里可以获得至少一个指南，如何编程 n次根函数，产生精确的float/Decimal ?
对于root(125, 1756482845) ， 不返回1或0的函数(第一个参数是数字，第二个参数是根深度(或其他东西))

编辑:所以，你给了我这个解决方案:n ** (1.0 / exp)，我知道当我问这个问题，但它只是不工作，例如，exp = 3。1/3不能用有理数表示，所以125 ** (1/3)给出了不正确的结果4.999999...。我要求一些"智能"算法，它可以为如此漂亮的数字提供正确的结果，并且至少为有理exp提供4个小数点精确的结果。如果没有这样的函数或算法，我将使用这个(n ** (1/exp))

我会尝试gmpy2库。

>>> import gmpy2
>>> gmpy2.root(125,3)
mpfr('5.0')
>>>

gmpy2使用MPFR库来执行正确的舍入浮点运算。默认精度为53位，但可以增加。

>>> gmpy2.root(1234567890123456789**11, 11)
mpfr('1.2345678901234568e+18')  # Last digits are incorrect.
>>> gmpy2.get_context().precision=200
>>> gmpy2.root(1234567890123456789**11, 11)
mpfr('1234567890123456789.0',200)
>>>

免责声明:我维护gmpy2 .

您可以对答案进行二分搜索。如果你想找到等于N的k次方根的X，你可以在X上做一个二分搜索，在二分搜索的每一步测试X^k是否等于N +-一个小常数，以避免精度问题。

代码如下:

import math
N,K = map(float,raw_input().split()) # We want Kth root of N
lo = 0.0
hi = N
while 1:
    mid = (lo+hi)/2
    if math.fabs(mid**K-N) < 1e-9: # mid^K is really close to N, consider mid^K == N
        print mid
        break
    elif mid**K < N: lo = mid
    else: hi = mid

对于(N,K) =(125,3)，打印正确答案5.0。您可以通过更改1e-9常量使其更精确，但Python中

有与浮点变量精度限制相关的精度限制

这是math模块的pow函数。

import math
math.pow(4, 0.5)

将返回4的平方根，即2.0。

对于root(125, 1756482845)，需要做的是

math.pow(125, 1.0 / 1756482845)

你的意思是:

>>> 125**(1/9.0)
1.7099759466766968

你可能会感兴趣的是bigfloat模块(个人没有使用过，只是知道它存在:)-实际上在过去安装它有问题-可能是OS X的错误)

在Squeak Smalltalk中，如果整数接收器恰好是某个整数的n次幂，则有一个nthRoot:消息回答确切的Integer结果。然而，如果解是一个代数根，那么实现不会回退到一个朴素的n**(1/exp);该方法通过对残差的适当处理，舍入到最接近的浮点数。

此处复制相关代码(MIT许可)。基本算法是用Newton-Raphson来搜索整数的截断n次根:

Integer>>nthRootTruncated: aPositiveInteger
    "Answer the integer part of the nth root of the receiver."
    | guess guessToTheNthMinusOne nextGuess |
    self = 0 ifTrue: [^0].
    self negative
        ifTrue:
            [aPositiveInteger even ifTrue: [ ArithmeticError signal: 'Negative numbers don''t have even roots.' ].
            ^(self negated nthRootTruncated: aPositiveInteger) negated].
    guess := 1 bitShift: self highBitOfMagnitude + aPositiveInteger - 1 // aPositiveInteger.
    [
        guessToTheNthMinusOne := guess raisedTo: aPositiveInteger - 1.
        nextGuess := (aPositiveInteger - 1 * guess * guessToTheNthMinusOne + self) // (guessToTheNthMinusOne * aPositiveInteger).
        nextGuess >= guess ] whileFalse:
            [ guess := nextGuess ].
    ( guess raisedTo: aPositiveInteger) > self  ifTrue:
            [ guess := guess - 1 ].
    ^guess

这不是特别聪明，因为在指数很大的情况下收敛会很慢，但是，它是有效的。然后，相同的根从零舍入:

Integer>>nthRootRounded: aPositiveInteger
    "Answer the integer nearest the nth root of the receiver."
    | guess |
    self = 0 ifTrue: [^0].
    self negative
        ifTrue:
            [aPositiveInteger even ifTrue: [ ArithmeticError signal: 'Negative numbers don''t have even roots.' ].
            ^(self negated nthRootRounded: aPositiveInteger) negated].
    guess := self nthRootTruncated: aPositiveInteger.
    ^self * 2 > ((guess + 1 raisedTo: aPositiveInteger) + (guess raisedTo: aPositiveInteger))
        ifTrue: [guess + 1]
        ifFalse: [guess]

然后在nroot中测试准确性:

Integer>>nthRoot: aPositiveInteger
    "Answer the nth root of the receiver.
    Answer an Integer if root is exactly this Integer, else answer the Float nearest the exact root."
    | guess excess scaled nBits |
    guess := self nthRootRounded: aPositiveInteger.
    excess := (guess raisedTo: aPositiveInteger) - self.
    excess = 0 ifTrue: [ ^ guess ].
    nBits := Float precision - guess highBitOfMagnitude.
    nBits <= 0 ifTrue: [ ^(Fraction numerator: guess * 4 - excess sign denominator: 4) asFloat].
    scaled := self << (nBits * aPositiveInteger).
    guess := scaled nthRootRounded: aPositiveInteger.
    excess := (guess raisedTo: aPositiveInteger) - scaled.
    ^(Fraction numerator: guess * 4 - excess sign denominator: 1 << (nBits + 2)) asFloat

这也可以应用于Fraction，但最接近的float有点复杂，Squeak的实现目前还很幼稚。

适用于大整数，如:

(10 raisedTo: 600) nthRoot: 300 -> 100 "exact"
(10 raisedTo: 600) + 1 nthRoot: 300 -> 100.0 "inexact"

如果你没有这样的期望，最初的猜测可以使用不精确的朴素n**(1/exp)。

代码应该易于在Python中移植，并留下大量优化空间。

我没有检查Python中可用的内容，但也许您需要正确四舍五入的LargeInteger -> Float和Fraction -> Float，就像这里解释的那样(Smalltalk也是，抱歉，但语言并不重要)。

http://smallissimo.blogspot.fr/2011/09/clarifying-and-optimizing.html
http://smallissimo.blogspot.fr/2011/09/reviewing-fraction-asfloat.html

相关内容

最新更新

热门标签：