回答編集履歴
3
d
answer
CHANGED
|
@@ -87,4 +87,32 @@
|
|
|
87
87
|
|
|
88
88
|
f = mpmath.factorial(mpmath.mpf(1000))
|
|
89
89
|
print(f"1000! = {f}")
|
|
90
|
+
```
|
|
91
|
+
|
|
92
|
+
### 質問の式を mpmath で計算する例
|
|
93
|
+
|
|
94
|
+
```python
|
|
95
|
+
import sys
|
|
96
|
+
|
|
97
|
+
import mpmath as mp
|
|
98
|
+
|
|
99
|
+
mpmath.mp.dps = 2000
|
|
100
|
+
|
|
101
|
+
s = 0
|
|
102
|
+
l = 0.3
|
|
103
|
+
n = 1000
|
|
104
|
+
t = 0.2
|
|
105
|
+
|
|
106
|
+
for k in range(n + 1):
|
|
107
|
+
s += (
|
|
108
|
+
(-1) ** k
|
|
109
|
+
/ mp.factorial(1000)
|
|
110
|
+
* mp.sin(mpmath.pi * l * k)
|
|
111
|
+
* mp.gamma(l * k + 1)
|
|
112
|
+
/ (t ** (l * k + 1))
|
|
113
|
+
)
|
|
114
|
+
|
|
115
|
+
s *= -1 / mpmath.pi
|
|
116
|
+
print(s)
|
|
117
|
+
# -4.6852969663640794183052611442921398568086793342834616462334696324903867843231395630449988227880810131784478141659219945020203824927023753145340152275411511818373558920285872536714395228721948271826415724657244356679410653902566023614945014042952126178451891146305498687076853150690388109943827589984735451898842139557958621055625632539663091175942064908619685606545219293485810476695033253309084612857919257606784710540821141287519704017988233984166363710615914801845057844847498708240109737744243798479965501611309757531588120324018949625977650842536557835170122470528575372965177074318211049609908479423513835244125701516943697486654794175290551525218964214548834205190130209616768780024865572422562319780254088489151655831942096511712901359285285313156837232190787227348293663807776826088102389004523746146439727884614438552769772052583787278931259025862266267527347055478532365682291188279115875850275919795196649524369131042988665553942619715595168719062131399527261531325133906580769041387916239304444853504881787451700196905154750706687217882182421629624078044956529865543544068656579784986705519747832854686508939351378181527847396122106095160866115322720223585062715943851550247952235652437843682553103198759150039829159286341277372652777332702063251410719991844368057697345823111259489942513149539813031718828345330830157844803109246345868740928996257665817861367467687161910693184414807511579751981740167047404088833058340052830658962523612136648135596475773090234244305096492102145406041366097252467314623275357778005920563778105868762907482707819978490405043196844875157208281722161855859448349229516220667066123504490431317899702446327175215167205030224795654189487061289564293054352001359895654181422873855808898785005875909967936576766340843356222045846500038297018788320027808611619194276127030472392754256051166888548527686436583618833795960420300558879420743599986345799684146483437131704540648998092418292801233242752904307374935321733025575849004051292291211225721865580467562778e-1745
|
|
90
118
|
```
|
2
d
answer
CHANGED
|
@@ -51,4 +51,40 @@
|
|
|
51
51
|
plt.show()
|
|
52
52
|
```
|
|
53
53
|
|
|
54
|
-
|
|
54
|
+
|
|
55
|
+
|
|
56
|
+
## 追記
|
|
57
|
+
|
|
58
|
+
> 1点質問なのですが、どうしても1000項以上計算したい場合は、どのような方法が考えられますでしょうか？アドバイスいただけますと幸いです。
|
|
59
|
+
|
|
60
|
+
1000項以上計算したい理由がグラフ描画であるならば、意味がないかもしれません。
|
|
61
|
+
例えば、頑張って有効数字20桁ぐらい計算したとしても、グラフ上ではその差は表現できないからです。
|
|
62
|
+
質問の式がなにかの関数のテイラー展開の結果なのだとしたら、最初の10項ぐらいでも有効数字が小数点以下何桁かはあると思いますが、グラフで描画する際の精度としては十分といえます。
|
|
63
|
+
|
|
64
|
+
### その上でどうしても1000項ぐらい計算したい場合
|
|
65
|
+
|
|
66
|
+
一般的に扱われる倍精度 (Python の float 型) の範囲では、1000!は表現できる範囲を超えてしまうので、[任意精度演算](https://ja.wikipedia.org/wiki/%E4%BB%BB%E6%84%8F%E7%B2%BE%E5%BA%A6%E6%BC%94%E7%AE%97) のライブラリを使う必要があります。
|
|
67
|
+
任意精度演算ライブラリを使えば、メモリが許す限り、任意の桁数の値を扱うことができます。
|
|
68
|
+
ただし、欠点として、単精度/倍精度などハードウェアレベルで最適化されている演算ではないので、計算はかなり遅くなるでしょう。
|
|
69
|
+
|
|
70
|
+
Python で使える任意精度演算ライブラリとしては sympy についてくる mpmath があります。
|
|
71
|
+
float では170項ぐらいで階乗はオーバーフローしますが、2000ビットの多倍長演算で1000項まで計算でき、[Wolfram Alpha の結果](https://www.wolframalpha.com/input/?i=1000!) と一致することを確認しました。
|
|
72
|
+
|
|
73
|
+
```python
|
|
74
|
+
import sys
|
|
75
|
+
|
|
76
|
+
import mpmath
|
|
77
|
+
from scipy.special import factorial
|
|
78
|
+
|
|
79
|
+
mpmath.mp.dps = 2000 # 1000ビット
|
|
80
|
+
|
|
81
|
+
for i in range(1001):
|
|
82
|
+
f = factorial(i)
|
|
83
|
+
print(f"{i}! = {f:.2e}")
|
|
84
|
+
|
|
85
|
+
if f > sys.float_info.max:
|
|
86
|
+
break # 171! で倍精度で表される範囲を超える。
|
|
87
|
+
|
|
88
|
+
f = mpmath.factorial(mpmath.mpf(1000))
|
|
89
|
+
print(f"1000! = {f}")
|
|
90
|
+
```
|
1
d
answer
CHANGED
|
@@ -42,7 +42,7 @@
|
|
|
42
42
|
print(seq.shape) # (t, n) の配列
|
|
43
43
|
|
|
44
44
|
# 級数を求めるには、行方向に和を取ればよいので、
|
|
45
|
-
series = 1 / np.pi * seq.sum(axis=1)
|
|
45
|
+
series = -1 / np.pi * seq.sum(axis=1)
|
|
46
46
|
print(series.shape) # (100,)
|
|
47
47
|
|
|
48
48
|
fig, ax = plt.subplots()
|