ふーん

メモメモ.

Sympyは純粋なPythonパッケージなので,高度に最適化された(CやFORTRANで書かれた)Numpy, SciPyモジュールの様に高速であるって事はまず無い.

sympy.prodよりsympy.Mulの方がかなり速い – Qiita

Pythonオブジェクト(リスト型とか)が欲しい場合は,randomモジュール.
numpy.ndarrayが欲しくて,数値計算をしたい場合はそのままnumpyの世界へ.

import random
import numpy as np
from functools import reduce
from operator import mul
import sympy


#nums_ = [int(i) for i in np.random.randint(1, 100, 100)]
#np.random.randint(1, 100, 100).tolist()
nums = [random.randrange(1, 100) for _ in range(100)]

res1 = nums[0]
for x in nums[1:]:
    res1 *= x

res2 = reduce(lambda x, y: x*y, nums)
res3 = reduce(mul, nums)
res4 = np.prod(nums, dtype=np.float64)
res5 = np.prod(np.array(nums), dtype=np.float64)
res6 = sympy.prod(nums)
res7 = sympy.Mul(*nums)

FireShot Capture 257 - JupyterLab Alpha Preview - http___localhost_8888_lab

Pythonは任意精度なので,メモリが許す限り計算ができる.
numpyはCとかに近くて,与えられた型の精度に依存する.
intを与えれば,プラットフォームに応じてint32 or int64,
dtypeを指定すれば,それに準じた精度になる.

今の場合,int64ではサチって0.floatでもサチっている分が誤差っている.

この手の場合,素直にピュアに書くのがリーズナブルチョイスじゃないだろうか.

%%timeit
res1 = nums[0]
for x in nums[1:]:
    res1 *= x

13.4 µs ± 298 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit reduce(lambda x, y: x*y, nums)
%timeit reduce(mul, nums)
%timeit np.prod(nums, dtype=np.float64)
%timeit np.prod(np.array(nums), dtype=np.float64)
%timeit sympy.prod(nums)
%timeit sympy.Mul(*nums)

24 µs ± 1.92 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
13 µs ± 156 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
22.6 µs ± 471 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
20.5 µs ± 1.79 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
13.4 µs ± 566 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
6.66 µs ± 44.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%%timeit
nums = [random.randrange(1, 100) for _ in range(100)]
res = nums[0]
for x in nums[1:]:
    res *= x

289 µs ± 18.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit reduce(lambda x, y: x*y, [random.randrange(1, 100) for _ in range(100)])
%timeit reduce(mul, [random.randrange(1, 100) for _ in range(100)])
%timeit np.prod([random.randrange(1, 100) for _ in range(100)], dtype=np.float64)
%timeit np.prod(np.array([random.randrange(1, 100) for _ in range(100)]), dtype=np.float64)
%timeit sympy.prod([random.randrange(1, 100) for _ in range(100)])
%timeit sympy.Mul(*[random.randrange(1, 100) for _ in range(100)])

305 µs ± 20.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
279 µs ± 2.79 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
303 µs ± 15.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
298 µs ± 12.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
297 µs ± 32.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
1.25 ms ± 55.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
 
 
 

「sympy.Mul」は知らなかった.
これって,「sympy.core.mul.Mul」を呼んでいるのか.

ちなみに,「sympy.prod」は中身をみると,

Source code for sympy.core.mul – SymPy 1.1.1 documentation


def prod(a, start=1):
    """Return product of elements of a. Start with int 1 so if only
       ints are included then an int result is returned.

    Examples
    ========

    >>> from sympy import prod, S
    >>> prod(range(3))
    0
    >>> type(_) is int
    True
    >>> prod([S(2), 3])
    6
    >>> _.is_Integer
    True

    You can start the product at something other than 1:

    >>> prod([1, 2], 3)
    6

    """
    return reduce(operator.mul, a, start)

「reduce(mul, nums)」と全く同じ.

近似値を求めるなら,mathモジュールを使うという手もある.

import math
math.exp(sum(map(math.log, nums)))

9.50370395101778e+166

%timeit math.exp(sum(map(math.log, nums)))

11.6 µs ± 75.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit math.exp(sum(map(math.log, [random.randrange(1, 100) for _ in range(100)])))

289 µs ± 24.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

思ったよりは遅いけど.

現実的な問題(有効数字を気にする様な)ではかなり有効な方法.
 
 
 

%timeit reduce(mul, range(1, 10))
%timeit math.exp(sum(map(math.log, range(1, 10))))
%timeit sympy.Mul(*[i for i in range(1, 10)])

1.33 µs ± 77.5 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
2.47 µs ± 204 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
2.56 µs ± 51.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit reduce(mul, range(1, 100))
%timeit math.exp(sum(map(math.log, range(1, 100))))
%timeit sympy.Mul(*[i for i in range(1, 100)])

12.8 µs ± 341 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
12 µs ± 125 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
12.8 µs ± 407 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
 
 
 

print(reduce(mul, range(1, 4)))
print(math.exp(sum(map(math.log, range(1, 4)))))
print(10**(sum(map(math.log10, range(1, 4)))))
%timeit reduce(mul, range(1, 4))
%timeit math.exp(sum(map(math.log, range(1, 4))))
%timeit 10**(sum(map(math.log10, range(1, 4))))

print(reduce(mul, range(1, 5)))
print(math.exp(sum(map(math.log, range(1, 5)))))
print(10**(sum(map(math.log10, range(1, 5)))))
%timeit reduce(mul, range(1, 5))
%timeit math.exp(sum(map(math.log, range(1, 5))))
%timeit 10**(sum(map(math.log10, range(1, 5))))

print(reduce(mul, range(1, 6)))
print(math.exp(sum(map(math.log, range(1, 6)))))
print(10**(sum(map(math.log10, range(1, 6)))))
%timeit reduce(mul, range(1, 6))
%timeit math.exp(sum(map(math.log, range(1, 6))))
%timeit 10**(sum(map(math.log10, range(1, 6))))

print(reduce(mul, range(1, 10)))
print(math.exp(sum(map(math.log, range(1, 10)))))
print(10**(sum(map(math.log10, range(1, 10)))))
%timeit reduce(mul, range(1, 10))
%timeit math.exp(sum(map(math.log, range(1, 10))))
%timeit 10**(sum(map(math.log10, range(1, 10))))

print(reduce(mul, range(1, 100)))
print(math.exp(sum(map(math.log, range(1, 100)))))
print(10**(sum(map(math.log10, range(1, 100)))))
%timeit reduce(mul, range(1, 100))
%timeit math.exp(sum(map(math.log, range(1, 100))))
%timeit 10**(sum(map(math.log10, range(1, 100))))

6
6.0
6.0
864 ns ± 78.9 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
2.13 µs ± 342 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
1.8 µs ± 134 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
24
23.999999999999993
23.999999999999993
816 ns ± 54.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
1.99 µs ± 219 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
1.99 µs ± 254 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
120
119.99999999999997
119.99999999999996
989 ns ± 108 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
2.01 µs ± 111 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
1.78 µs ± 44.9 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
362880
362880.00000000047
362880.0000000002
1.27 µs ± 46.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
2.5 µs ± 267 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
2.12 µs ± 130 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000
9.332621544395465e+155
9.332621544393672e+155
12.7 µs ± 484 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
12.1 µs ± 282 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
8.16 µs ± 123 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

for i in range(4, 31):
    print('-'*30)
    print(f'prod(range(1, {i})):')
    print(reduce(mul, range(1, i)))
    print(math.exp(sum(map(math.log, range(1, i)))))
    print(10**(sum(map(math.log10, range(1, i)))))
    print('-'*30)

prod(range(1, 4)):
6
6.0

6.0


prod(range(1, 5)):
24
23.999999999999993

23.999999999999993


prod(range(1, 6)):
120
119.99999999999997

119.99999999999996


prod(range(1, 7)):
720
720.0000000000001

719.9999999999993


prod(range(1, 8)):
5040
5040.000000000002

5039.999999999996


prod(range(1, 9)):
40320
40320.00000000002

40319.99999999999


prod(range(1, 10)):
362880
362880.00000000047

362880.0000000002


prod(range(1, 11)):
3628800
3628800.0000000084

3628800.000000002


prod(range(1, 12)):
39916800
39916800.00000003

39916800.00000005


prod(range(1, 13)):
479001600
479001600.000001

479001599.99999994


prod(range(1, 14)):
6227020800
6227020800.0000105

6227020800.000002


prod(range(1, 15)):
87178291200
87178291200.00017

87178291200.00005


prod(range(1, 16)):
1307674368000
1307674368000.003

1307674368000.001


prod(range(1, 17)):
20922789888000
20922789888000.055

20922789888000.027


prod(range(1, 18)):
355687428096000
355687428096000.7

355687428096000.6


prod(range(1, 19)):
6402373705728000
6402373705727994.0

6402373705728017.0


prod(range(1, 20)):
121645100408832000
1.2164510040883208e+17

1.216451004088323e+17


prod(range(1, 21)):
2432902008176640000
2.43290200817664e+18

2.432902008176646e+18


prod(range(1, 22)):
51090942171709440000
5.109094217170942e+19

5.10909421717096e+19


prod(range(1, 23)):
1124000727777607680000
1.1240007277776115e+21

1.1240007277776123e+21


prod(range(1, 24)):
25852016738884976640000
2.5852016738885062e+22

2.585201673888509e+22


prod(range(1, 25)):
620448401733239439360000
6.204484017332436e+23

6.204484017332426e+23


prod(range(1, 26)):
15511210043330985984000000
1.5511210043331066e+25

1.5511210043331063e+25


prod(range(1, 27)):
403291461126605635584000000
4.0329146112660826e+26

4.0329146112660785e+26


prod(range(1, 28)):
10888869450418352160768000000
1.0888869450418425e+28

1.0888869450418425e+28


prod(range(1, 29)):
304888344611713860501504000000
3.0488834461171542e+29

3.0488834461171602e+29


prod(range(1, 30)):
8841761993739701954543616000000
8.841761993739751e+30

8.841761993739792e+30

広告
カテゴリー: 未分類 パーマリンク

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中