15

二百万から三百万桁の平方根の桁を生成したい。

私はNewton-Raphsonを認識していますが、biginteger がサポートされていないため、C または C++ でそれを実装する方法がわかりません。誰かが私を正しい方向に向けることができますか?

また、誰かがPythonでそれを行う方法を知っていれば（私は初心者です）、私も感謝します。

4

9 に答える 9

7

マッピングを使用してみることができます：

a/b -> (a+2b)/(a+b)で始まりa= 1, b= 1ます。これは sqrt(2) に収束します (実際には、連分数表現が得られます)。

重要なポイント: これは、行列の乗算として表すことができます (フィボナッチに似ています)。

a_n と b_n がステップの n 番目の数字の場合

[1 2] [a_n b_n] T = [a_(n+1) b_(n+1)] T
[1 1]

それは今私たちに与えます

[1 2] n [a_1 b_1] T = [a_(n+1) b_(n+1)] T
[1 1]

したがって、2x2 行列が A の場合、A nを計算する必要があります。これは、2 乗を繰り返すことで実行でき、整数演算のみを使用します (したがって、精度の問題について心配する必要はありません)。

また、得られる a/b は常に (gcd(a,b) = gcd(a+2b, a+b) のように) 縮小された形式になることに注意してください。結果、しないでください！

n 番目の分母は (1+sqrt(2))^n のようなものなので、300 万桁を取得するには、3671656番目の項まで計算する必要があります。

約 360 万番目の項を探している場合でも、2 乗を繰り返すと、O(Log n) の乗算と加算で n 番目の項を計算できることに注意してください。

また、これは Newton-Raphson などの反復的なものとは異なり、簡単に並列化できます。

于 2011-03-04T16:08:06.807 に答える
6

EDIT : このバージョンは前のバージョンよりも好きです。これは、整数と小数の両方を受け入れる一般的なソリューションです。n = 2 で精度 = 100000 の場合、約 2 分かかります。Paul McGuire の提案とその他の提案を歓迎します。

def sqrt_list(n, precision):
    ndigits = []        # break n into list of digits
    n_int = int(n)
    n_fraction = n - n_int

    while n_int:                            # generate list of digits of integral part
        ndigits.append(n_int % 10)
        n_int /= 10
    if len(ndigits) % 2: ndigits.append(0)  # ndigits will be processed in groups of 2

    decimal_point_index = len(ndigits) / 2  # remember decimal point position
    while n_fraction:                       # insert digits from fractional part
        n_fraction *= 10
        ndigits.insert(0, int(n_fraction))
        n_fraction -= int(n_fraction)
    if len(ndigits) % 2: ndigits.insert(0, 0)  # ndigits will be processed in groups of 2

    rootlist = []
    root = carry = 0                        # the algorithm
    while root == 0 or (len(rootlist) < precision and (ndigits or carry != 0)):
        carry = carry * 100
        if ndigits: carry += ndigits.pop() * 10 + ndigits.pop()
        x = 9
        while (20 * root + x) * x > carry:
                x -= 1
        carry -= (20 * root + x) * x
        root = root * 10 + x
        rootlist.append(x)
    return rootlist, decimal_point_index
于 2011-03-03T23:33:55.470 に答える
3

[1; 2, 2, ...]おそらく最も良い方法は、2 の平方根の連分数展開を使用することです。

def root_two_cf_expansion():
    yield 1
    while True:
        yield 2

def z(a,b,c,d, contfrac):
    for x in contfrac:
        while a > 0 and b > 0 and c > 0 and d > 0:
            t = a // c
            t2 = b // d
            if not t == t2:
                break
            yield t
            a = (10 * (a - c*t))
            b = (10 * (b - d*t))
            # continue with same fraction, don't pull new x
        a, b = x*a+b, a
        c, d = x*c+d, c
    for digit in rdigits(a, c):
        yield digit

def rdigits(p, q):
    while p > 0:
        if p > q:
           d = p // q
           p = p - q * d
        else:
           d = (10 * p) // q
           p = 10 * p - q * d
        yield d

def decimal(contfrac):
    return z(1,0,0,1,contfrac)

decimal((root_two_cf_expansion())すべての 10 進数の反復子を返します。 t1およびt2アルゴリズムでは、次の桁の最小値と最大値です。それらが等しい場合、その数字を出力します。

これは、連分数の負の数などの特定の例外的なケースを処理しないことに注意してください。

(このコードは、あちこちに散らばっている連分数を処理するための Haskell コードの適応です。)

于 2011-03-04T06:05:25.353 に答える
3

作業用？図書館を使おう！

楽しみのために？よかったね ：）

鉛筆と紙で行うことを模倣するプログラムを作成します。最初は 1 桁、次に 2 桁、次に 3、...、...

ニュートンや他の人のことは心配しないでください。あなたのやり方でやってください。

于 2011-03-03T23:04:04.000 に答える
3

任意の大きな数については、The GNU Multiple Precision Arithmetic Library (for C/C++) を参照してください。

于 2011-03-03T22:59:13.893 に答える
2

さて、以下は私が書いたコードです。私にとっては、約 60800 秒で 2 の平方根の小数点以下 100 万桁が生成されましたが、プログラムを実行しているときにラップトップがスリープしていたので、それよりも速いはずです。300 万桁の生成を試みることができますが、それを取得するには数日かかる場合があります。

def sqrt(number,digits_after_decimal=20):
import time
start=time.time()
original_number=number
number=str(number)
list=[]
for a in range(len(number)):
    if number[a]=='.':
        decimal_point_locaiton=a
        break
    if a==len(number)-1:
        number+='.'
        decimal_point_locaiton=a+1
if decimal_point_locaiton/2!=round(decimal_point_locaiton/2):
    number='0'+number
    decimal_point_locaiton+=1
if len(number)/2!=round(len(number)/2):
    number+='0'
number=number[:decimal_point_locaiton]+number[decimal_point_locaiton+1:]
decimal_point_ans=int((decimal_point_locaiton-2)/2)+1
for a in range(0,len(number),2):
    if number[a]!='0':
        list.append(eval(number[a:a+2]))
    else:
        try:
            list.append(eval(number[a+1]))
        except IndexError:
            pass
p=0
c=list[0]
x=0
ans=''
for a in range(len(list)):
    while c>=(20*p+x)*(x):
        x+=1
    y=(20*p+x-1)*(x-1)
    p=p*10+x-1
    ans+=str(x-1)
    c-=y
    try:
        c=c*100+list[a+1]
    except IndexError:
        c=c*100
while c!=0:
    x=0
    while c>=(20*p+x)*(x):
        x+=1
    y=(20*p+x-1)*(x-1)
    p=p*10+x-1
    ans+=str(x-1)
    c-=y
    c=c*100
    if len(ans)-decimal_point_ans>=digits_after_decimal:
            break
ans=ans[:decimal_point_ans]+'.'+ans[decimal_point_ans:]
total=time.time()-start
return ans,total
于 2011-08-09T22:24:20.393 に答える
2

以下は、整数a の平方根を精度の桁数まで計算するための短いバージョンです。これは、2 x桁に 10 を掛けた後のaの整数の平方根を見つけることによって機能します。

def sqroot(a, digits):
    a = a * (10**(2*digits))
    x_prev = 0
    x_next = 1 * (10**digits)
    while x_prev != x_next:
        x_prev = x_next
        x_next = (x_prev + (a // x_prev)) >> 1
    return x_next

いくつかの注意事項があります。

結果を文字列に変換し、正しい位置に小数点を追加する必要があります (小数点を表示したい場合)。

非常に大きな整数を文字列に変換するのは、それほど高速ではありません。

非常に大きな整数の除算も (Python では) あまり高速ではありません。

システムのパフォーマンスによっては、小数点以下 200 万から 300 万の平方根を計算するのに 1 時間以上かかる場合があります。

ループが常に終了することは証明していません。最後の桁が異なる 2 つの値の間で振動する場合があります。またはそうではないかもしれません。

于 2011-03-04T05:00:05.597 に答える
0

Pythonはすでに箱から出して大きな整数をサポートしており、それがC / C ++であなたを阻む唯一のものである場合は、いつでも自分でクイックコンテナクラスを書くことができます。

あなたが言及した唯一の問題は、大きな整数の欠如です。そのためにライブラリを使用したくない場合は、そのようなクラスを作成するためのヘルプを探していますか？

于 2011-03-03T23:16:05.180 に答える
0

これは、すべての場合に終了する、より効率的な整数平方根関数 (Python 3.x) です。平方根にはるかに近い数値から開始するため、ステップ数が少なくなります。int.bit_length には Python 3.1 以降が必要であることに注意してください。簡潔にするために、エラー チェックは省略されています。

def isqrt(n):
    x = (n >> n.bit_length() // 2) + 1
    result = (x + n // x) // 2
    while abs(result - x) > 1:
        x = result
        result = (x + n // x) // 2
    while result * result > n:
        result -= 1
    return result
于 2011-03-04T11:40:36.327 に答える