回答編集履歴
3
d
test
CHANGED
|
@@ -177,3 +177,59 @@
|
|
|
177
177
|
print(f"1000! = {f}")
|
|
178
178
|
|
|
179
179
|
```
|
|
180
|
+
|
|
181
|
+
|
|
182
|
+
|
|
183
|
+
### 質問の式を mpmath で計算する例
|
|
184
|
+
|
|
185
|
+
|
|
186
|
+
|
|
187
|
+
```python
|
|
188
|
+
|
|
189
|
+
import sys
|
|
190
|
+
|
|
191
|
+
|
|
192
|
+
|
|
193
|
+
import mpmath as mp
|
|
194
|
+
|
|
195
|
+
|
|
196
|
+
|
|
197
|
+
mpmath.mp.dps = 2000
|
|
198
|
+
|
|
199
|
+
|
|
200
|
+
|
|
201
|
+
s = 0
|
|
202
|
+
|
|
203
|
+
l = 0.3
|
|
204
|
+
|
|
205
|
+
n = 1000
|
|
206
|
+
|
|
207
|
+
t = 0.2
|
|
208
|
+
|
|
209
|
+
|
|
210
|
+
|
|
211
|
+
for k in range(n + 1):
|
|
212
|
+
|
|
213
|
+
s += (
|
|
214
|
+
|
|
215
|
+
(-1) ** k
|
|
216
|
+
|
|
217
|
+
/ mp.factorial(1000)
|
|
218
|
+
|
|
219
|
+
* mp.sin(mpmath.pi * l * k)
|
|
220
|
+
|
|
221
|
+
* mp.gamma(l * k + 1)
|
|
222
|
+
|
|
223
|
+
/ (t ** (l * k + 1))
|
|
224
|
+
|
|
225
|
+
)
|
|
226
|
+
|
|
227
|
+
|
|
228
|
+
|
|
229
|
+
s *= -1 / mpmath.pi
|
|
230
|
+
|
|
231
|
+
print(s)
|
|
232
|
+
|
|
233
|
+
# -4.6852969663640794183052611442921398568086793342834616462334696324903867843231395630449988227880810131784478141659219945020203824927023753145340152275411511818373558920285872536714395228721948271826415724657244356679410653902566023614945014042952126178451891146305498687076853150690388109943827589984735451898842139557958621055625632539663091175942064908619685606545219293485810476695033253309084612857919257606784710540821141287519704017988233984166363710615914801845057844847498708240109737744243798479965501611309757531588120324018949625977650842536557835170122470528575372965177074318211049609908479423513835244125701516943697486654794175290551525218964214548834205190130209616768780024865572422562319780254088489151655831942096511712901359285285313156837232190787227348293663807776826088102389004523746146439727884614438552769772052583787278931259025862266267527347055478532365682291188279115875850275919795196649524369131042988665553942619715595168719062131399527261531325133906580769041387916239304444853504881787451700196905154750706687217882182421629624078044956529865543544068656579784986705519747832854686508939351378181527847396122106095160866115322720223585062715943851550247952235652437843682553103198759150039829159286341277372652777332702063251410719991844368057697345823111259489942513149539813031718828345330830157844803109246345868740928996257665817861367467687161910693184414807511579751981740167047404088833058340052830658962523612136648135596475773090234244305096492102145406041366097252467314623275357778005920563778105868762907482707819978490405043196844875157208281722161855859448349229516220667066123504490431317899702446327175215167205030224795654189487061289564293054352001359895654181422873855808898785005875909967936576766340843356222045846500038297018788320027808611619194276127030472392754256051166888548527686436583618833795960420300558879420743599986345799684146483437131704540648998092418292801233242752904307374935321733025575849004051292291211225721865580467562778e-1745
|
|
234
|
+
|
|
235
|
+
```
|
2
d
test
CHANGED
|
@@ -105,3 +105,75 @@
|
|
|
105
105
|
|
|
106
106
|
|
|
107
107
|
|
|
108
|
+
|
|
109
|
+
|
|
110
|
+
|
|
111
|
+
## 追記
|
|
112
|
+
|
|
113
|
+
|
|
114
|
+
|
|
115
|
+
> 1点質問なのですが、どうしても1000項以上計算したい場合は、どのような方法が考えられますでしょうか？アドバイスいただけますと幸いです。
|
|
116
|
+
|
|
117
|
+
|
|
118
|
+
|
|
119
|
+
1000項以上計算したい理由がグラフ描画であるならば、意味がないかもしれません。
|
|
120
|
+
|
|
121
|
+
例えば、頑張って有効数字20桁ぐらい計算したとしても、グラフ上ではその差は表現できないからです。
|
|
122
|
+
|
|
123
|
+
質問の式がなにかの関数のテイラー展開の結果なのだとしたら、最初の10項ぐらいでも有効数字が小数点以下何桁かはあると思いますが、グラフで描画する際の精度としては十分といえます。
|
|
124
|
+
|
|
125
|
+
|
|
126
|
+
|
|
127
|
+
### その上でどうしても1000項ぐらい計算したい場合
|
|
128
|
+
|
|
129
|
+
|
|
130
|
+
|
|
131
|
+
一般的に扱われる倍精度 (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) のライブラリを使う必要があります。
|
|
132
|
+
|
|
133
|
+
任意精度演算ライブラリを使えば、メモリが許す限り、任意の桁数の値を扱うことができます。
|
|
134
|
+
|
|
135
|
+
ただし、欠点として、単精度/倍精度などハードウェアレベルで最適化されている演算ではないので、計算はかなり遅くなるでしょう。
|
|
136
|
+
|
|
137
|
+
|
|
138
|
+
|
|
139
|
+
Python で使える任意精度演算ライブラリとしては sympy についてくる mpmath があります。
|
|
140
|
+
|
|
141
|
+
float では170項ぐらいで階乗はオーバーフローしますが、2000ビットの多倍長演算で1000項まで計算でき、[Wolfram Alpha の結果](https://www.wolframalpha.com/input/?i=1000!) と一致することを確認しました。
|
|
142
|
+
|
|
143
|
+
|
|
144
|
+
|
|
145
|
+
```python
|
|
146
|
+
|
|
147
|
+
import sys
|
|
148
|
+
|
|
149
|
+
|
|
150
|
+
|
|
151
|
+
import mpmath
|
|
152
|
+
|
|
153
|
+
from scipy.special import factorial
|
|
154
|
+
|
|
155
|
+
|
|
156
|
+
|
|
157
|
+
mpmath.mp.dps = 2000 # 1000ビット
|
|
158
|
+
|
|
159
|
+
|
|
160
|
+
|
|
161
|
+
for i in range(1001):
|
|
162
|
+
|
|
163
|
+
f = factorial(i)
|
|
164
|
+
|
|
165
|
+
print(f"{i}! = {f:.2e}")
|
|
166
|
+
|
|
167
|
+
|
|
168
|
+
|
|
169
|
+
if f > sys.float_info.max:
|
|
170
|
+
|
|
171
|
+
break # 171! で倍精度で表される範囲を超える。
|
|
172
|
+
|
|
173
|
+
|
|
174
|
+
|
|
175
|
+
f = mpmath.factorial(mpmath.mpf(1000))
|
|
176
|
+
|
|
177
|
+
print(f"1000! = {f}")
|
|
178
|
+
|
|
179
|
+
```
|
1
d
test
CHANGED
|
@@ -86,7 +86,7 @@
|
|
|
86
86
|
|
|
87
87
|
# 級数を求めるには、行方向に和を取ればよいので、
|
|
88
88
|
|
|
89
|
-
series = 1 / np.pi * seq.sum(axis=1)
|
|
89
|
+
series = -1 / np.pi * seq.sum(axis=1)
|
|
90
90
|
|
|
91
91
|
print(series.shape) # (100,)
|
|
92
92
|
|