3

以前した質問((バイナリ法を用いた)多項式の累乗計算について)の続きの質問です。

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

# -*- 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))

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

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

3
  • 1
    type = numpy.object;f2 = numpy.poly1d(numpy.array([4, 3, 2, 1], dtype = type)) の部分の array を numpy.array([4L, 3L, 2L, 1L], dtype=typ) としてみて下さい(L(long) suffix を追加)。これで多倍長計算ができる様になります。
    – user9156
    2015年8月2日 6:37
  • N = 21 のときはうまくいくのですが、N = 22 とすると、負の数が出てきておかしくなるのはなぜでしょう?
    – TOM
    2015年8月2日 6:58
  • 1
    def mth_degree_poly(p, m): return numpy.poly1d(list(p)[- m - 1:]) => numpy.poly1d(numpy.array(list(p)[- m - 1:], dtype = numpy.object)), p = numpy.poly1d(numpy.array([1], dtype = type)) => p = numpy.poly1d(numpy.array([1L], dtype = type)) にする必要があります。
    – user9156
    2015年8月2日 7:14

1 件の回答 1

1

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

# -*- 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))

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。