12917


以前した質問 "[冪乗計算を高速に行うには？](https://ja.stackoverflow.com/q/10006)" の多項式版についての質問です。

次のように「（バイナリ法を用いた）多項式の累乗計算」を行うコードを  
Ruby 2.2 および Python 2.7 で作成し、それぞれ (1 + x)^n を計算させたのですが、  
Python 2.7版の方が n が34以上で結果がおかしくなってしまいました。  
どのように直せばよいか教えていただけないでしょうか？

(Ruby 2.2)

[code]

    # (f_ary×b_ary)のm次以下を取り出す
    def mul(f_ary, b_ary, m)
      s1, s2 = f_ary.size, b_ary.size
      ary = Array.new(s1 + s2 - 1, 0)
      s10 = [s1 - 1, m].min
      (0..s10).each{|i|
        s20 = [s2 - 1, m - i].min
        (0..s20).each{|j|
          ary[i + j] += f_ary[i] * b_ary[j]
        }
      }
      ary
    end
    
    # ary^nのm次以下を取り出す
    def power(ary, n, m)
      return [1] if n == 0
      k = power(ary, n >> 1, m)
      k = mul(k, k, m)
      return k if n & 1 == 0
      return mul(k, ary, m)
    end
    
    f0 = [1, 1]
    p power(f0, 33, 33)
    p power(f0, 34, 34)
    
[/code]

出力結果:

[code]

    [1, 33, 528, 5456, 40920, 237336, 1107568, 4272048, 13884156, 38567100, 92561040
    , 193536720, 354817320, 573166440, 818809200, 1037158320, 1166803110, 1166803110
    , 1037158320, 818809200, 573166440, 354817320, 193536720, 92561040, 38567100, 13
    884156, 4272048, 1107568, 237336, 40920, 5456, 528, 33, 1]
    [1, 34, 561, 5984, 46376, 278256, 1344904, 5379616, 18156204, 52451256, 13112814
    0, 286097760, 548354040, 927983760, 1391975640, 1855967520, 2203961430, 23336062
    20, 2203961430, 1855967520, 1391975640, 927983760, 548354040, 286097760, 1311281
    40, 52451256, 18156204, 5379616, 1344904, 278256, 46376, 5984, 561, 34, 1]
    
[/code]

(Python 2.7)  
NumPyを使っています。

[code]

    import numpy
    
    def power(f, n):
        p = numpy.poly1d([1])
        for i in list(format (n, 'b')):
            p *= p
            if i == '1':
                p *= f
        return p
    
    f0 = numpy.poly1d([1, 1])
    print list(reversed((power(f0, 33)).c))
    print list(reversed((power(f0, 34)).c))
    
[/code]

出力結果:

[code]

    [1, 33, 528, 5456, 40920, 237336, 1107568, 4272048, 13884156, 38567100, 92561040
    , 193536720, 354817320, 573166440, 818809200, 1037158320, 1166803110, 1166803110
    , 1037158320, 818809200, 573166440, 354817320, 193536720, 92561040, 38567100, 13
    884156, 4272048, 1107568, 237336, 40920, 5456, 528, 33, 1]
    [1, 34, 561, 5984, 46376, 278256, 1344904, 5379616, 18156204, 52451256, 13112814
    0, 286097760, 548354040, 927983760, 1391975640, 1855967520, -2091005866, -196136
    1076, -2091005866, 1855967520, 1391975640, 927983760, 548354040, 286097760, 1311
    28140, 52451256, 18156204, 5379616, 1344904, 278256, 46376, 5984, 561, 34, 1]
    
[/code]

（追記）  
次のようにバイナリ法を用いずとも全く同じ問題が発生している。

[code]

    import numpy
    
    def power(f, n):
        i = 1
        p = f
        while i < n:
            p *= f
            i += 1
        return p
    
    f0 = numpy.poly1d([1, 1])
    print list(reversed((power(f0, 33)).c))
    print list(reversed((power(f0, 34)).c))
    
[/code]

出力結果:  
質問のPython 2.7版の出力結果と同じ


CC BY-SA 4.0
2015-08-01T07:07:09.320
12910
2021-03-10T05:57:34.323
2021-03-10T05:57:34.323
3060
5363
question
ruby python アルゴリズム
バイナリ法を用いた多項式の累乗計算について


numpy では、整数は C
などと同じフォーマットで格納されているので、桁あふれを起こしているのでしょう。20億を越えたあたりでおかしくなっているので、32bitの符号付き整数だとすると辻褄が合います。ただ、手元の
Linux Mint で試したところ、64bitのためか、問題が再現しませんでした。

試しに

[code]

    import numpy
    
    def power(f, n):
        p = numpy.poly1d(numpy.array([1], dtype=numpy.int64))
        for i in format (n, 'b'):
            p *= p
            if i == '1':
                p *= f
        return p
    
    f0 = numpy.poly1d(numpy.array([1, 1], dtype=numpy.int64))
    print list(reversed((power(f0, 33)).c))
    print list(reversed((power(f0, 34)).c))
    
[/code]

などと、データ型を指定してみたら、どうでしょうか。


CC BY-SA 3.0
2015-08-01T08:25:23.663
12917
2015-08-01T08:25:23.663
3605
12910
answer