14453


<http://worldrecord314.com/pi>  
上記ページには円周率10万桁が載っているのですが、  
9万9901桁から10万桁までを簡単に求めるには  
どうすればよいでしょうか？

参考までにRubyのコードと実行結果を記しておきます。

[code]

    require 'benchmark'
    require 'bigdecimal/math'
    include BigMath
    
    puts Benchmark::CAPTION
    puts Benchmark.measure{
      N = 100000
      M = N * 2
      p BigDecimal(PI(M).to_s).to_s("F")[N - 98..N + 1] # 3.と続くのでN + 1
    }
    
[/code]

実行結果

[code]

          user     system      total        real
    "8575016363411314627530499019135646823804329970695770150789337728658035712790913767420805655493624646"
    118.047000   7.266000 125.313000 (125.424243)
    
[/code]

（追記）  
<https://stackoverflow.com/questions/28284996/python-pi-calculation>  
上記記事のPythonのコードをいじってみたものの、  
Nが10000位で速度がかなり遅くなってしまった。  
なお、N = 1000 ならすぐ求まります。

[code]

    from decimal import Decimal, getcontext
    
    def power(a, n):
        p = 1
        for i in format (n, 'b'):
            p *= p
            if i == '1':
                p *= a
        return p
    
    N = 1000
    getcontext().prec = N + 1
    s = 0
    for k in range(N):
        s += (Decimal(4) / (8 * k + 1) - Decimal(2) / (8 * k + 4) - Decimal(1) / (8 * k + 5) - Decimal(1) / (8 * k + 6)) / power(Decimal(16), k)
    
    print str(s)[N - 98:]
    
[/code]

実行結果

[code]

    5982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
    
[/code]

また、  
π = 2 + 1 / 3( 2 + 2 / 5( 2 + 3 / 7( 2 + ･･･ ( 2 + k / (2k + 1)( 2 +
･･･))･･･)))  
（例えば次を参照してください。<http://xn--w6q13e505b.jp/formula/arctan.html>）  
という式を利用した  
<http://d.hatena.ne.jp/htz/20110613/1307958164>  
に載っているRubyのコード  
（ちなみに、len = 31, 32 のときを比較するとわかるように必ずしも末尾は正しくありません。）  
をいじるとさらに速く求まりました。

[code]

    require 'benchmark'
    
    puts Benchmark::CAPTION
    puts Benchmark.measure{
      N = 10 ** 5
      B = 10 ** N
      B2 = B << 1
      # (i >> 1) / iの部分がk / (2k + 1)にあたる
      p ((N * 8 + 1).step(3, -2).inject(B){|a, i| (i >> 1) * (a + B2) / i} - B).to_s[-100..-1]
    }
    
[/code]

実行結果

[code]

          user     system      total        real
    "8575016363411314627530499019135646823804329970695770150789337728658035712790913767420805655493624646"
     76.547000   7.875000  84.422000 ( 84.476262)
    
[/code]


CC BY-SA 3.0
2015-08-15T18:58:25.173
14452
2015-08-17T16:23:42.967
2017-05-23T12:38:56.467
-1
5363
question
ruby python アルゴリズム
円周率の99901桁から100000桁までを簡単に求めるには？


「99900桁目までを直接計算することなく、99901桁目からの100桁を得たい」という意図と解釈しました。現在、10進数でそのようなことをする方法は知られていないと思います。

16進数でしたら、n桁目までを飛ばしてそれ以降を直接計算する方法が知られています。Bailey–Borwein–Plouffe
formulaを使います。具体的な手順は以下のページの中程に出ています。

<https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula>

2011年のサーベイ論文で他のいくつかの超越数に対しても類似のアルゴリズムが挙げられていますが、10進数には言及がないので、まだ見つかっていないのだと思います。

<https://www.carma.newcastle.edu.au/jon/bbp-bluegene.pdf>


CC BY-SA 3.0
2015-08-15T22:43:01.637
14453
2015-08-15T22:43:01.637
5258
14452
answer


追記のRubyのコードは、  
π ≒ 2 + 1 / 3( 2 + 2 / 5( 2 + 3 / 7( 2 + ･･･ ( 2 + k / (2k + 1) * 2))))  
≒ 2 + 1 / 3( 2 + 2 / 5( 2 + 3 / 7( 2 + ･･･ ( 2 + 1))))  
を利用しています。  
これに対し、  
π ≒ 2 + 1 / 3( 2 + 2 / 5( 2 + 3 / 7( 2 + ･･･ ( 2 + k / (2k + 1) * 2))))  
= 2 + 2 * 1 / 3 + 2 * 1 / 3 * 2 / 5 + ･･･ + 2 * 1 / 3 * ･･･ * k / (2k + 1)  
を利用することを考える。  
π - (2 + 2 * 1 / 3 + 2 * 1 / 3 * 2 / 5 + ･･･ + 2 * 1 / 3 * ･･･ * k / (2k + 1))  
< 2((1 / 2)^(k + 1) + (1 / 2)^(k + 2) + ･･･)  
= (1 / 2)^(k - 1)  
より、よってlog2 10 * d + 1以上のkを選ぶと  
π - (2 + 2 * 1 / 3 + 2 * 1 / 3 * 2 / 5 + ･･･ + 2 * 1 / 3 * ･･･ * k / (2k + 1))  
< 1 / 10^d  
となる。  
よって、除算の切り捨てによる誤差を無視できるなら、  
追記のRubyのコードの

[code]

      p ((N * 8 + 1).step(3, -2)
    
[/code]

の部分は N * 7 + 1 位でよいと考えられる。  
すなわち次のコードになる。

[code]

    require 'benchmark'
    
    puts Benchmark::CAPTION
    puts Benchmark.measure{
      N = 10 ** 5
      # 以下一桁多く計算し、末尾を一つ取り除く
      B = 10 ** (N + 1)
      B2 = B << 1
      # (i >> 1) / iの部分がk / (2k + 1)にあたる
      p (((N + 1) * 7 + (N + 2) % 2).step(3, -2).inject(0){|a, i| (i >> 1) * (a + B2) / i} - B).to_s.chop[-100..-1]
    }
    
[/code]

実行結果

[code]

          user     system      total        real
    "8575016363411314627530499019135646823804329970695770150789337728658035712790913767420805655493624646"
     65.391000   7.563000  72.954000 ( 72.945916)
    
[/code]

なお、「一桁多く計算し、末尾を一つ取り除く」という処理は、  
繰り上がりの問題に対処するために行った。  
しかしこの処理も完全なものではなく、これをすり抜けるものがあることは  
以下のコードを実行することでわかる。

[code]

    n = 2
    str = '1'
    new_str = '14'
    c_new_str = new_str.clone
    n += 1
    while n < 500 && str == c_new_str.chop
      str = new_str
      # 以下一桁多く計算し、末尾を一つ取り除く
      b = 10 ** (n + 1)
      b2 = b << 1
      # (i >> 1) / iの部分がk / (2k + 1)にあたる
      new_str = (((n + 1) * 7 + (n + 2) % 2).step(3, -2).inject(0){|a, i| (i >> 1) * (a + b2) / i} - b).to_s.chop
      p new_str
      c_new_str = new_str.clone
      n += 1
    end
    p n
    
[/code]

実行結果

[code]

    （中略）
    "14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259035"
    "141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360"
    361
    
[/code]


CC BY-SA 3.0
2015-08-17T11:58:04.207
14531
2015-08-17T16:23:42.967
2015-08-17T16:23:42.967
5363
5363
14452
answer