12946


以前した質問（[（バイナリ法を用いた）多項式の累乗計算について](https://ja.stackoverflow.com/questions/12910/%E3%83%90%E3%82%A4%E3%83%8A%E3%83%AA%E6%B3%95%E3%82%92%E7%94%A8%E3%81%84%E3%81%9F-%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%81%AE%E7%B4%AF%E4%B9%97%E8%A8%88%E7%AE%97%E3%81%AB%E3%81%A4%E3%81%84%E3%81%A6/12917#12917)）の続きの質問です。

以前の質問で、  
「データ型を指定することで、より大きな累乗計算をできる」  
ことがわかりました。  
それでは、どこまで累乗計算ができるか気になり、次のコードを実行しました。

[code]

    # -*- coding: cp932 -*-
    import numpy
    
    # m次以下のみにする
    def mth_degree_poly(p, m):
        return numpy.poly1d(list(p)[- m - 1:])
    
    def power(f, n, m, type):
        p = numpy.poly1d(numpy.array([1], dtype = type))
        for i in format (n, 'b'):
            p *= p
            p = mth_degree_poly(p, m)
            if i == '1':
                p *= f
                p = mth_degree_poly(p, m)
        return p
    
    # int64型
    type = numpy.int64
    f0 = numpy.poly1d(numpy.array([4, 3, 2, 1], dtype = type))
    print list(reversed((power(f0, 20, 20 * 3, type)).c))
    print list(reversed((power(f0, 21, 21 * 3, type)).c))
    # データ型指定なし
    type = None
    f1 = numpy.poly1d(numpy.array([4, 3, 2, 1]))
    print list(reversed((power(f1, 10, 10 * 3, type)).c))
    print list(reversed((power(f1, 11, 11 * 3, type)).c))
    # object型
    type = numpy.object
    f2 = numpy.poly1d(numpy.array([4, 3, 2, 1], dtype = type))
    print list(reversed((power(f2, 11, 11 * 3, type)).c))
    print list(reversed((power(f2, 12, 12 * 3, type)).c))
    
[/code]

その結果、  
「上記計算では、データ型が  
int64型だと20乗まで、  
指定なしだと（私の環境では）10乗まで、  
object型だと11乗まで  
しか正確に計算できない」  
ことがわかりました。

Numpyでは、  
「N > 20 なる整数に対し、  
(1 + 2x + 3x^2 + 4x^3)^N  
の計算はできない」  
のでしょうか？  
それとも、上記コードでは工夫が足りないということでしょうか？


CC BY-SA 3.0
2015-08-02T05:54:51.057
12943
2015-08-02T07:47:24.600
2017-04-13T12:52:38.920
-1
5363
question
python
Numpyを使ってより大きな累乗計算をさせるには？


質問のコードでは工夫が足りなかったようです。  
argus さんのコメントに従い、以下のように変更。

[code]

    # -*- coding: cp932 -*-
    import numpy
    
    # m次以下のみにする
    def mth_degree_poly(p, m, type):
        return numpy.poly1d(numpy.array(list(p)[- m - 1:], dtype = type))
    
    def power(f, n, m, type):
        p = numpy.poly1d(numpy.array([1L], dtype = type))
        for i in format (n, 'b'):
            p *= p
            p = mth_degree_poly(p, m, type)
            if i == '1':
                p *= f
                p = mth_degree_poly(p, m, type)
        return p
    
    type = numpy.object
    f2 = numpy.poly1d(numpy.array([4L, 3L, 2L, 1L], dtype = type))
    print list(reversed((power(f2, 50, 50 * 3, type)).c))
    
[/code]


CC BY-SA 3.0
2015-08-02T07:38:40.537
12946
2015-08-02T07:47:24.600
2015-08-02T07:47:24.600
5363
5363
12943
answer