前置き

python3系では標準で整数は無限桁扱うことができ，mpmathパッケージを利用すれば実数(複素数)の多倍長数値演算が可能である．個人的に，numpy.arrayと組み合わせて利用しているのだが，いくつか注意が必要なので個人的な備忘録として注意すべき点を記録する．

多倍長評価のdemo

>>> import python 
>>> math.factorial(100)                                                                                     
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
>>> import mpmath
>>> mpmath.mp.dps = 100
>>> print(mpmath.pi)                                                                                            
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

精度の確認

python 標準のfloatに関する精度は次のコマンドで確認可能

>>> import sys 
>>> sys.float_info                                                          
sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308,
			   min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, 
			   dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)

epsilonは計算機イプシロンで1の次に大きい最小の数とのこと．numpyのarray(ndarray)はC言語で実装されていると何かで読んだので，C言語の型の制約を受ける．実際numpy.arrayを作る際にdtypeが指定される．

numpyの場合finfo(浮動小数点), iinfo(整数)で確認することができnumpy:overflow-error

>>> import numpy as np                                                          
>>> np.finfo(np.float)                                                             
finfo(resolution=1e-15, min=-1.7976931348623157e+308, max=1.7976931348623157e+308, dtype=float64)
>>> np.finfo(np.float128)                                                                                                                                                             
finfo(resolution=1e-18, min=-1.189731495357231765e+4932, max=1.189731495357231765e+4932, dtype=float128)
>>> np.iinfo(np.int)                                                                                                
iinfo(min=-9223372036854775808, max=9223372036854775807, dtype=int64)

から確認できる．float128はCで言うlong doubleであってfloat128(4倍浮動小数点)ではない．Overflowの挙動は

>>> np.array([np.finfo(np.float).max], dtype=np.float) * 2  
array([inf])
>>> np.array([np.iinfo(np.int).max], dtype=np.int) + 1                                                              
array([-9223372036854775808])

となる．pythonも浮動小数点は倍精度の壁があるので，floatに関してoverflowするのは仕方ないとおもう．しかし，python3はintを無限精度で扱えるのでうっかりoverflowや意図せぬ精度不足に陥ることがある．numpy.arrayでnp.intの精度限界を超えた整数を扱いたければ，dtypeを無指定にするか，dtypeにobjectを指定すれば良い．

つまり

import math
>>> np.array([math.factorial(25)], dtype=np.int)                                                                    
Traceback (most recent call last):
  File "<ipython-input-46-0451d8765c1b>", line 1, in <module>
    np.array([math.factorial(25)], dtype=np.int)
OverflowError: Python int too large to convert to C long

だが

>>> np.array([math.factorial(25)])                                                                                  
array([15511210043330985984000000], dtype=object)

となる．dtypeをobjectに指定(無指定)するとmpmathのような多倍長浮動小数点も扱える．

多倍長浮動小数(mpmath)でnumpy.arrayをつかう

intの場合と同じ様にdtypeを無指定(or dtype=object)にすれば良い．

>>> import mpmath                                                                                                   
>>> mpmath.mp.dps = 50                                                                                             
>>> np.array([mpmath.mpf("1")])/3                                                                                   
array([mpf('0.33333333333333333333333333333333333333333333333333311')], dtype=object)

昔は演算の順序を変えると未定義だと怒られることがあったが，今はdtypeがobjectなのかnp.floatやnp.intなのか意識せずほぼ同等四則演算等の操作が可能(arrayの要素が四則演算可能なら)．当然であるが，dtypeがnp.float(やnp.int)のarrayと演算を施すと精度が失われる．

>>> np.array(mpmath.linspace(0, 1, 10)) + np.linspace(0, 10, 10)                                                    
array([mpf('0.0'),
       mpf('1.2222222222222222715654677611180684632725185818142358'),
       mpf('2.4444444444444445431309355222361369265450371636284716'),
       mpf('3.6666666666666668146964032833542053898175557454427101'),
       mpf('4.8888888888888890862618710444722738530900743272569433'),
       mpf('6.1111111111111109137381289555277261469099256727430567'),
       mpf('7.3333333333333336293928065667084107796351114908854202'),
       mpf('8.555555555555556345047484177889095412360297309027773'),
       mpf('9.7777777777777781725237420889445477061801486545138865'),
       mpf('11.0')], dtype=object)

計算速度の比較

numpy.arrayで dtype=np.intやnp.floatは計算速度のbest perfomanceがでるだろう．dtype=objectにするとどの程度遅くなるのか確認しておく．

import numpy as np
import mpmath
import time

def addsub(x, imax):
    x0 = np.copy(x)
    time0 = time.time()
    for n in range(imax):
        xx = x + x
        x = xx - x
    print(time.time()-time0, "[sec]", sum(x-x0), x.dtype)

def powpow(x, imax):
    x0 = np.copy(x)
    time0 = time.time()
    for n in range(imax):
        xx = x**(3)
        x = xx**(1/3)
    print(time.time()-time0, "[sec]", float(sum(x-x0)), x.dtype)

imax = 10000
sample = 10000

x = np.arange(0, sample, 1, dtype=np.int)
addsub(x, imax)

x = np.arange(0, sample, 1, dtype=object)
addsub(x, imax)


x = np.arange(0, sample, 1, dtype=np.float)
powpow(x, imax)

x = np.arange(0, sample, 1, dtype=object)
powpow(x, imax)

次が実行結果

0.06765484809875488 [sec] 0 int64
3.387320041656494 [sec] 0 object
8.931994915008545 [sec] -0.00024138918300131706 float64
16.408015489578247 [sec] -0.00024138918300131706 object

arrayの要素がintの場合は，objectに指定すると圧倒的に遅くなる傾向にあるが，一方でfloatに関しては2倍程度と許容の範囲だろう．さてmpmathの場合を考えよう．mpmathはデフォルトの仮数部は15桁なので,15桁と仮数部100桁の場合で計算速度を比較しよう．ただ，上の設定だと数時間ほど掛かりそうなのでimaxを2桁落とす．そのため，上の結果と比較するには計算時間を100倍すれば，おおよそ比較可能．

imax = 100

mpmath.mp.dps=100
x = np.array(mpmath.arange(0, sample, 1), dtype=object)
powpow(x, imax)

次が実行結果

15.098075151443481 [sec] -2.4139151442170714e-06 object

mpmath，仮数部15桁で計算量を1/100にしているのでだいたい100倍程度かかる．

ちなみにpowpowの定義では精度は改善されない．なぜなら$x^{1/3}$とした際に指数部がpythonのfloatで評価されてしまうため15桁以降は意味の無い数字となる．そこでpowpowを少し修正して

def papawpowpow(x, imax):
    x0 = np.copy(x)
    time0 = time.time()
    for n in range(imax):
        xx = x**(3)
        x = xx**(mpmath.mpf("1/3"))
    print(time.time()-time0, "[sec]", float(sum(x-x0)), x.dtype)

として実行すると

21.777454376220703 [sec] 1.5530282978129947e-91 object

mpmath.mp.dps=100
x = np.array(mpmath.arange(0, sample, 1), dtype=object)
papawpowpow(x, imax)

ちなみにarrayを使わずにelementwiseで評価すると

def elewisepapawpowpow(x, imax):
    x0 = x[:]
    time0 = time.time()        
    for n in range(imax):
        xx = [x1**(3) for x1 in x]
        x =  [x1**(mpmath.mpf("1/3")) for x1 in xx]
    diff = [x[i] - x0[i] for i in range(len(x0))]
    print(time.time()-time0, "[sec]", float(sum(diff)), type(x))

となる．

28.9152569770813 [sec] 1.5530282978129947e-91 <class 'list'>

計算速度はnumpy.arrayの場合と比較して30%弱早く，コーディングも使い慣れたnumpy.arrayを使った方がスッキリする．

実行環境

$ lscpu 
アーキテクチャ:                      x86_64
CPU 操作モード:                      32-bit, 64-bit
バイト順序:                          Little Endian
CPU:                                 4
オンラインになっている CPU のリスト: 0-3
コアあたりのスレッド数:              1
ソケットあたりのコア数:              4
ソケット数:                          1
NUMA ノード数:                       1
ベンダー ID:                         AuthenticAMD
CPU ファミリー:                      23
モデル:                              17
モデル名:                            AMD Ryzen 3 2200G with Radeon Vega Graphics
ステッピング:                        0
CPU MHz:                             1479.915
CPU 最大 MHz:                        3500.0000
CPU 最小 MHz:                        1600.0000
BogoMIPS:                            7000.24
仮想化:                              AMD-V
L1d キャッシュ:                      32K
L1i キャッシュ:                      64K
L2 キャッシュ:                       512K
L3 キャッシュ:                       4096K
NUMA ノード 0 CPU:                   0-3
フラグ:                              fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl nonstop_tsc cpuid extd_apicid aperfmperf pni pclmulqdq monitor ssse3 fma cx16 sse4_1 sse4_2 movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb hw_pstate sme ssbd sev ibpb vmmcall fsgsbase bmi1 avx2 smep bmi2 rdseed adx smap clflushopt sha_ni xsaveopt xsavec xgetbv1 xsaves clzero irperf xsaveerptr arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif overflow_recov succor smca

$ cat /proc/cpuinfo | grep MHz
cpu MHz		: 1624.805
cpu MHz		: 1500.378
cpu MHz		: 3700.129
cpu MHz		: 2006.406