Benchmarkingパッケージを使ったMathematicaのベンチマーク
私物laptop
Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz
# Core: 2
# Threading: 4
{"MachineName" -> "--", "System" -> "Linux x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.3.1",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 1.392,
"TotalTime" -> 9.945, "Results" -> {{"Data Fitting", 0.537},
{"Digits of Pi", 0.341}, {"Discrete Fourier Transform", 0.68},
{"Eigenvalues of a Matrix", 0.608}, {"Elementary Functions", 0.567},
{"Gamma Function", 0.467}, {"Large Integer Multiplication", 0.467},
{"Matrix Arithmetic", 0.334}, {"Matrix Multiplication", 0.754},
{"Matrix Transpose", 1.61}, {"Numerical Integration", 0.676},
{"Polynomial Expansion", 0.092}, {"Random Number Sort", 1.16},
{"Singular Value Decomposition", 0.765}, {"Solving a Linear System", 0.887}}}
------------------------------------------------------------------------------------------------
とあるサーバー１:
Intel(R) Xeon(R) CPU X5550 @ 2.67GHz
# Core: 4
# Threading: 8
{"MachineName" -> "--", "System" -> "Linux x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 1.332,
"TotalTime" -> 10.388, "Results" -> {{"Data Fitting", 1.059},
{"Digits of Pi", 0.507}, {"Discrete Fourier Transform", 0.545},
{"Eigenvalues of a Matrix", 1.001}, {"Elementary Functions", 0.389},
{"Gamma Function", 0.66}, {"Large Integer Multiplication", 0.668},
{"Matrix Arithmetic", 0.536}, {"Matrix Multiplication", 0.408},
{"Matrix Transpose", 1.022}, {"Numerical Integration", 1.002},
{"Polynomial Expansion", 0.136}, {"Random Number Sort", 1.159},
{"Singular Value Decomposition", 0.651}, {"Solving a Linear System", 0.645}}}
------------------------------------------------------------------------------------------------
とあるサーバー２:
Intel(R) Xeon(R) CPU E5607 @ 2.27GHz
# Core: 4?8?
# Threading: 8
{"MachineName" -> "--", "System" -> "Linux x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 1.056,
"TotalTime" -> 13.103, "Results" -> {{"Data Fitting", 0.962},
{"Digits of Pi", 0.648}, {"Discrete Fourier Transform", 0.6},
{"Eigenvalues of a Matrix", 1.076}, {"Elementary Functions", 0.481},
{"Gamma Function", 0.874}, {"Large Integer Multiplication", 0.916},
{"Matrix Arithmetic", 0.706}, {"Matrix Multiplication", 0.541},
{"Matrix Transpose", 1.485}, {"Numerical Integration", 1.342},
{"Polynomial Expansion", 0.14}, {"Random Number Sort", 1.542},
{"Singular Value Decomposition", 0.878}, {"Solving a Linear System", 0.912}}}
------------------------------------------------------------------------------------------------
Model Name: iMac
Processor Name: Intel Core i7
Processor Speed: 3.1 GHz
# Core: 4
# Threading: 8
{"MachineName" -> "--", "System" -> "Mac OS X x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.3.1",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 1.985,
"TotalTime" -> 6.975, "Results" -> {{"Data Fitting", 0.689},
{"Digits of Pi", 0.295}, {"Discrete Fourier Transform", 0.369},
{"Eigenvalues of a Matrix", 0.474}, {"Elementary Functions", 0.506},
{"Gamma Function", 0.353}, {"Large Integer Multiplication", 0.33},
{"Matrix Arithmetic", 0.632}, {"Matrix Multiplication", 0.259},
{"Matrix Transpose", 0.726}, {"Numerical Integration", 0.62},
{"Polynomial Expansion", 0.086}, {"Random Number Sort", 0.819},
{"Singular Value Decomposition", 0.435}, {"Solving a Linear System", 0.382}}}
------------------------------------------------------------------------------------------------
Model Name: iMac
Processor Name: Intel Core 2 Duo
Processor Speed: 3.06 GHz
# Core: 2
# Threading: 2(?)
{"MachineName" -> "--", "System" -> "Mac OS X x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 0.613,
"TotalTime" -> 22.586, "Results" -> {{"Data Fitting", 1.713},
{"Digits of Pi", 0.67}, {"Discrete Fourier Transform", 1.548},
{"Eigenvalues of a Matrix", 1.356}, {"Elementary Functions", 2.113},
{"Gamma Function", 0.768}, {"Large Integer Multiplication", 0.769},
{"Matrix Arithmetic", 2.632}, {"Matrix Multiplication", 1.646},
{"Matrix Transpose", 1.804}, {"Numerical Integration", 1.949},
{"Polynomial Expansion", 0.267}, {"Random Number Sort", 1.669},
{"Singular Value Decomposition", 1.992}, {"Solving a Linear System", 1.69}}}
-----------------------------------------------------------------------------------------------
raspberry pi 2
Coretex-A7
# Core: 4
# threading: 4(?)
{"MachineName" -> "--", "System" -> "Linux ARM (32-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 24, 2016", "BenchmarkResult" -> 0.029,
"TotalTime" -> 475.891,
"Results" -> {{"Data Fitting", 10.904}, {"Digits of Pi", 4.805}, {
"Discrete Fourier Transform", 30.977}, {
"Eigenvalues of a Matrix", 35.914}, {
"Elementary Functions", 25.479}, {"Gamma Function", 6.437}, {
"Large Integer Multiplication", 6.735}, {
"Matrix Arithmetic", 6.07}, {"Matrix Multiplication", 151.389}, {
"Matrix Transpose", 7.202}, {"Numerical Integration", 11.684}, {
"Polynomial Expansion", 1.086}, {"Random Number Sort", 8.376}, {
"Singular Value Decomposition", 75.031}, {
"Solving a Linear System", 93.802}}}
----------------------------------------------------------------------------------------------------
とあるサーバー
Intel(R) Xeon(R) CPU E5-1680 v3 @ 3.20GHz
{"MachineName" -> "--", "System" -> "Linux x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.4.1",
"Date" -> "June 1, 2016", "BenchmarkResult" -> 2.06, "TotalTime" -> 6.721,
"Results" -> {{"Data Fitting", 0.568}, {"Digits of Pi", 0.293},
{"Discrete Fourier Transform", 0.274}, {"Eigenvalues of a Matrix",
0.78}, {"Elementary Functions", 0.209}, {"Gamma Function", 0.395},
{"Large Integer Multiplication", 0.389}, {"Matrix Arithmetic", 0.214},
{"Matrix Multiplication", 0.282}, {"Matrix Transpose", 0.611},
{"Numerical Integration", 0.705}, {"Polynomial Expansion", 0.098},
{"Random Number Sort", 0.984}, {"Singular Value Decomposition", 0.589},
{"Solving a Linear System", 0.33}}}
2016年5月7日土曜日
2016年1月24日日曜日
Mathematica ベンチマーク @ Raspberry pi 2
Raspberry pi 2 を買ったのでMathematica(無償)のベンチマークを測ってみる。
mathematicaはnoobs から rasbian (jessei) installでデフォルトで導入されていた。
mathematicaはnoobs から rasbian (jessei) installでデフォルトで導入されていた。
--------------- 参考 (手元のノート core-i7)-----------------------
In[1]:= Needs["Benchmarking`"]
In[2]:= Timing@Benchmark[]
Out[2]= {14.596, InputForm[{
"MachineName" -> "----", "System" -> "Linux x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.0",
"Date" -> "January 24, 2016", "BenchmarkResult" -> 1.197,
"TotalTime" -> 11.567,
"Results" -> {{"Data Fitting", 0.629}, {"Digits of Pi", 0.536}, {
"Discrete Fourier Transform", 0.728}, {
"Eigenvalues of a Matrix", 0.935}, {
"Elementary Functions", 0.774}, {"Gamma Function", 0.456}, {
"Large Integer Multiplication", 0.46}, {
"Matrix Arithmetic", 0.416}, {"Matrix Multiplication", 0.983}, {
"Matrix Transpose", 1.468}, {"Numerical Integration", 0.62}, {
"Polynomial Expansion", 0.092}, {"Random Number Sort", 1.104}, {
"Singular Value Decomposition", 1.268}, {
"Solving a Linear System", 1.098}}}]}
-------------- Raspberry pi 2 (オーバクロック無し) ----------------------
In[1]:= Needs["Benchmarking`"]
In[2]:= Timing@Benchmark[]
Out[2]= {551.620000, InputForm[{
"MachineName" -> "----", "System" -> "Linux ARM (32-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 24, 2016", "BenchmarkResult" -> 0.029,
"TotalTime" -> 475.891,
"Results" -> {{"Data Fitting", 10.904}, {"Digits of Pi", 4.805}, {
"Discrete Fourier Transform", 30.977}, {
"Eigenvalues of a Matrix", 35.914}, {
"Elementary Functions", 25.479}, {"Gamma Function", 6.437}, {
"Large Integer Multiplication", 6.735}, {
"Matrix Arithmetic", 6.07}, {"Matrix Multiplication", 151.389}, {
"Matrix Transpose", 7.202}, {"Numerical Integration", 11.684}, {
"Polynomial Expansion", 1.086}, {"Random Number Sort", 8.376}, {
"Singular Value Decomposition", 75.031}, {
"Solving a Linear System", 93.802}}}]}
線形代数の差がすごいな。topで観察する限りElementary Functions の計算でcpu使用率が400%になる。
Raspberry pi 2 は４コアなので、フル性能使えているんだと思うけど、並列計算用karnelの起動できない。(LaunchKarnels[]の実行結果が空)
(ParallelDoかなんか実行た際に並列計算用のパッケージ?みたいなのをダウンロード後使えるようになった)
一方面白いことに、有料版Mathematicaはライセンスの制約上フロントエンドカーネルは２つまで(mathematicaを独立に２つまで)しか起動できないが、
Raspberry pi の場合はその制限が無いようだ。
他に気づいたことは環境設定が無いといこと、MathematicaScript、MatheKarnelが無いこと・・・（残念）
追記:代わりにwolframというコマンドがあるっぽい
Raspberry 版のMathematicaは何かしら機能制限がありそうだと思ってはいたけどちゃんと確かめてなかったので他にも制約がありそうです。
本当はラズパイ並列して、mathematicaタスク分散とか・・・、もう少し遊べそうかと思ったけど無理そうだな。
まー遊び道具としてもうちょっといじってみるか。
速度の問題はあるけど5000円程度でmathematicaが使えると思うと非常に安いのでは？
Mathematicaは日本人特別価格の設定でHome edition で7万くらい(日本以外は$300，日本から海外版購入不可)だからね．
-----------
raspberry piを安価pcとして使用した際に気づいたこと。
microsdの相性の問題で起動しなかったので、調べて買うこと。
wikiにまとまっているのでそれを買えば問題ないはず。
os がarmなのでubuntuのようにパッケージが豊富でないこと
例えばfirefox、やchromeがaptになく代わりにiceweaseelやmidoriがある。
chromium-browserは一つ前のversionのosではあったようだけど、最新版ではaptに無いため、debを拾ってきて インストールする必要がある。
Adobe flashが無いので困る人には困るかもしれない。
操作性は思ったより遅くはないけど機敏ではないかな
(El captanよりは応答が良いと思うのは私だけだろうか)
動画を見るのは結構大変かも？
オーバクロックのturboモードはめちゃくちゃ不安定だったけどpi2は割と安定してそう。
-------------- Raspberry pi 2 (オーバクロック無し) ----------------------
In[1]:= Needs["Benchmarking`"]
In[2]:= Timing@Benchmark[]
Out[2]= {551.620000, InputForm[{
"MachineName" -> "----", "System" -> "Linux ARM (32-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 24, 2016", "BenchmarkResult" -> 0.029,
"TotalTime" -> 475.891,
"Results" -> {{"Data Fitting", 10.904}, {"Digits of Pi", 4.805}, {
"Discrete Fourier Transform", 30.977}, {
"Eigenvalues of a Matrix", 35.914}, {
"Elementary Functions", 25.479}, {"Gamma Function", 6.437}, {
"Large Integer Multiplication", 6.735}, {
"Matrix Arithmetic", 6.07}, {"Matrix Multiplication", 151.389}, {
"Matrix Transpose", 7.202}, {"Numerical Integration", 11.684}, {
"Polynomial Expansion", 1.086}, {"Random Number Sort", 8.376}, {
"Singular Value Decomposition", 75.031}, {
"Solving a Linear System", 93.802}}}]}
線形代数の差がすごいな。topで観察する限りElementary Functions の計算でcpu使用率が400%になる。
(ParallelDoかなんか実行た際に並列計算用のパッケージ?みたいなのをダウンロード後使えるようになった)
一方面白いことに、有料版Mathematicaはライセンスの制約上フロントエンドカーネルは２つまで(mathematicaを独立に２つまで)しか起動できないが、
Raspberry pi の場合はその制限が無いようだ。
速度の問題はあるけど5000円程度でmathematicaが使えると思うと非常に安いのでは？
Mathematicaは日本人特別価格の設定でHome edition で7万くらい(日本以外は$300，日本から海外版購入不可)だからね．
-----------
raspberry piを安価pcとして使用した際に気づいたこと。
microsdの相性の問題で起動しなかったので、調べて買うこと。
wikiにまとまっているのでそれを買えば問題ないはず。
os がarmなのでubuntuのようにパッケージが豊富でないこと
例えばfirefox、やchromeがaptになく代わりにiceweaseelやmidoriがある。
chromium-browserは一つ前のversionのosではあったようだけど、最新版ではaptに無いため、debを拾ってきて インストールする必要がある。
Adobe flashが無いので困る人には困るかもしれない。
操作性は思ったより遅くはないけど機敏ではないかな
動画を見るのは結構大変かも？
オーバクロックのturboモードはめちゃくちゃ不安定だったけどpi2は割と安定してそう。
2014年5月27日火曜日
MathemathicaScript の引数
MathematicaScriptの引数の引き渡しをいつも忘れてしまうので覚え書き．
MathematicaScriptの引数は
とでもすればよい．
MathematicaScriptの引数は
$ScriptCommandLine渡され，一番目はScriptのファイル名，2番目から引数となる．引数の制御は
If [Length[$ScriptCommandLine] !=2, Print["usage: MathematicaScript -script" <> $ScriptCommandLine[[1]] <> " value1"]; Exit[1] ]とでもしておけば良い．引数はstringで受け取るので，実数とかにしたければ
ToExpression@$ScriptCommandLine[[2]]
とでもすればよい．
2013年12月1日日曜日
Mathmaticaを用いた数値対角化の計算時間
Mathematicaを用いた対角化の計算速度をplotしてみました．
対角化した行列のclassは
ランダムではなくスパーズではない一般のエルミート行列(赤ドット)と
ランダムではなくスパースではない一般のユニタリー行列(青ドット)．
また，行列は200桁精度で近似評価されています．
厳密評価のまま計算するとキットモット遅い．
行列の構成方法は気にしないと言う事でご勘弁．
実線は計算時間を３次の多項式でfittingしたもの．
最後に描画とfittingのプログラムを置いときます．
(もちろん計算時間のデータファイルがないと動きませんが・・・
欲しいなら上げますが，欲しいか？)
対角化した行列のclassは
ランダムではなくスパーズではない一般のエルミート行列(赤ドット)と
ランダムではなくスパースではない一般のユニタリー行列(青ドット)．
また，行列は200桁精度で近似評価されています．
厳密評価のまま計算するとキットモット遅い．
行列の構成方法は気にしないと言う事でご勘弁．
実線は計算時間を３次の多項式でfittingしたもの．
最後に描画とfittingのプログラムを置いときます．
(もちろん計算時間のデータファイルがないと動きませんが・・・
欲しいなら上げますが，欲しいか？)
import numpy as np
import matplotlib.pyplot as plt
def f(x,para):
return para[0]*x**3 + para[1]*x**2 + para[2]*x + para[3]
data = np.loadtxt("time_prec200.dat").transpose()
fig = plt.figure(figsize=(8,4))
ax0 = fig.add_subplot(1,2,1)
ax1 = fig.add_subplot(1,2,2)
colors=["red", "blue"]
for i, ax in enumerate([ax0,ax1]):
ax.plot(data[0],data[1]/3600,'o', label='Hermitian matrix',markersize=10, color=colors[0])
ax.plot(data[0],data[2]/3600,'o', label='Unitary matrix',markersize=10,color=colors[1])
# ax.plot(data[0],data[3]/3600,'og')
for j in range(1,3):
x = np.arange(100, 10000,1)
para = np.polyfit(data[0],data[j],3)
y = f(x,para)
ax.plot(x,y/3600,'-',linewidth=3,color=colors[j-1])
if i==1:
ax.loglog()
ax.set_xlim(50,10000)
else:
ax.set_ylim(0,5)
ax.set_xlim(0,1000)
ax.set_ylabel("calculation time (hours)",fontsize=16)
le = ax.legend(loc=2,numpoints=1)
le.draw_frame(False)
ax.grid()
ax.set_xlabel("matrix dimension",fontsize=16)
sup = plt.suptitle("diagonalization time with Mathematica (precision=200)",fontsize=15)
sup.set_position((0.5,1.0))
plt.tight_layout()
plt.savefig("calc_time.png")
plt.show()
2013年11月16日土曜日
Mathematicaの糖衣構文
Mathematicaの関数へのアクセスは糖衣構文によって色々な人に受け入れやすい事になっています．
多分，同じ数学概念でも棲んでいる世界によって解釈が異なったりするからだと思います．ただ，私には甘すぎては苦手です．
可読性って概念はMathematicaにはないのかね…
他人のプログラム読む時，絶望を覚えるんですが・・・
例えば以下の例は全てSqrt[3]を返します．
最近Mathematicaに関する事ばかりだな・・・
多分，同じ数学概念でも棲んでいる世界によって解釈が異なったりするからだと思います．ただ，私には甘すぎては苦手です．
可読性って概念はMathematicaにはないのかね…
他人のプログラム読む時，絶望を覚えるんですが・・・
例えば以下の例は全てSqrt[3]を返します．
In[1]:= f[x_, y_: 2] := x^(1/y)
In[2]:= f[3]
f[3, 2]
3~f~2
3 // f
3 // f[#, 2] &
f @ 3
f @@ {3}
f @@ {3, 2}
f[#, 2] & @ 3
Apply[f, {3}]
Apply[f, {3, 2}]
Apply[f, {##}] & @@ {3, 2}
Out[2]= Sqrt[3]
Out[3]= Sqrt[3]
Out[4]= Sqrt[3]
Out[5]= Sqrt[3]
Out[6]= Sqrt[3]
Out[7]= Sqrt[3]
Out[8]= Sqrt[3]
Out[9]= Sqrt[3]
Out[10]= Sqrt[3]
Out[11]= Sqrt[3]
Out[12]= Sqrt[3]
Out[13]= Sqrt[3]
最近Mathematicaに関する事ばかりだな・・・
Mathematica Fourier変換について
バグを疑い始めたので，プログラムを純化させてみる．
In[1]:= IdentityTr[x_] := InverseFourier[Fourier[x]]
In[2]:= a = {0, 0, 0.33333333333333333, 0}
b = {0, 0, 0.333333333333333333, 0}
Precision[a[[3]]]
Precision[b[[3]]]
Out[2]= {0, 0, 0.333333, 0}
Out[3]= {0, 0, 0.33333333333333333, 0}
Out[4]= MachinePrecision
Out[5]= 17.5229
In[8]:= NestList[IdentityTr, a, 20]
NestList[IdentityTr, b, 20]
Out[8]= {{0, 0, 0.333333, 0}, {0., 0., 0.333333, 0.}, {0., 0.,
0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0.,
0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333,
0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0.,
0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0.,
0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333,
0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0.,
0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}}
Out[9]= {{0, 0, 0.33333333333333333, 0}, {0.*10^-18,
0.*10^-18 + 0.*10^-18 I, 0.33333333333333333,
0.*10^-18 + 0.*10^-18 I}, {0.*10^-17, 0.*10^-18 + 0.*10^-18 I,
0.3333333333333333, 0.*10^-18 + 0.*10^-18 I}, {0.*10^-17,
0.*10^-17 + 0.*10^-17 I, 0.333333333333333,
0.*10^-17 + 0.*10^-17 I}, {0.*10^-16, 0.*10^-16 + 0.*10^-16 I,
0.333333333333333, 0.*10^-16 + 0.*10^-16 I}, {0.*10^-15,
0.*10^-15 + 0.*10^-15 I, 0.33333333333333,
0.*10^-15 + 0.*10^-15 I}, {0.*10^-14, 0.*10^-15 + 0.*10^-15 I,
0.3333333333333, 0.*10^-15 + 0.*10^-15 I}, {0.*10^-14,
0.*10^-14 + 0.*10^-14 I, 0.333333333333,
0.*10^-14 + 0.*10^-14 I}, {0.*10^-13, 0.*10^-13 + 0.*10^-13 I,
0.33333333333, 0.*10^-13 + 0.*10^-13 I}, {0.*10^-12,
0.*10^-12 + 0.*10^-12 I, 0.33333333333,
0.*10^-12 + 0.*10^-12 I}, {0.*10^-11, 0.*10^-11 + 0.*10^-11 I,
0.3333333333, 0.*10^-11 + 0.*10^-11 I}, {0.*10^-10,
0.*10^-11 + 0.*10^-11 I, 0.333333333,
0.*10^-11 + 0.*10^-11 I}, {0.*10^-10, 0.*10^-10 + 0.*10^-10 I,
0.33333333, 0.*10^-10 + 0.*10^-10 I}, {0.*10^-9,
0.*10^-9 + 0.*10^-9 I, 0.33333333,
0.*10^-9 + 0.*10^-9 I}, {0.*10^-8, 0.*10^-8 + 0.*10^-8 I, 0.3333333,
0.*10^-8 + 0.*10^-8 I}, {0.*10^-7, 0.*10^-8 + 0.*10^-8 I, 0.333333,
0.*10^-8 + 0.*10^-8 I}, {0.*10^-7, 0.*10^-7 + 0.*10^-7 I, 0.33333,
0.*10^-7 + 0.*10^-7 I}, {0.*10^-6, 0.*10^-6 + 0.*10^-6 I, 0.3333,
0.*10^-6 + 0.*10^-6 I}, {0.*10^-5, 0.*10^-5 + 0.*10^-5 I, 0.3333,
0.*10^-5 + 0.*10^-5 I}, {0.*10^-4, 0.*10^-4 + 0.*10^-4 I, 0.333,
0.*10^-4 + 0.*10^-4 I}, {0.*10^-3, 0.*10^-4 + 0.*10^-4 I, 0.33,
0.*10^-4 + 0.*10^-4 I}}
2013年11月9日土曜日
Mathematicaでの任意制度計算
Mathematicaが手に入ったので，いろいろテストしてみました．
その結果， 任意精度計算を行うにあたって精度の保持はユーザーが監視しないと行けない事が分りました．
また，監視を怠ると正しい結果が得られない事も分りました．
簡単な例として離散Fourier変換(FFTか？)と逆変換を繰り返す演算をさせます．
Wolfram Mathematica Document Foureirの説明にある通り，
Inputのlistが厳密な数値である場合，
N関数を適用し精度はMachinePrecisionに落とされます．
そのためFoureir関数は厳密精度計算を行ってくれません．
しかしながらユーザーがInputのlistにN関数を適用し精度をユーザー側でコントロールしたい状況が発生しますが，安易な方法では正しい結果が得られません．
その結果が In[3]です．
一方listにN関数を作用させて作った初期条件の場合，全ての結果のPrecisionが0になっています．
これでは使い物になりません．１回のFourier関数の演算でどうやら10^2-10^10程度の桁落ちが発生し最終的にゼロが帰ってくる様です．
上記問題を回避する為にはSetPrecision関数を用いて，強制的に精度を維持する事です
(func1関数を参照)．
ただしこの方法が正当な回避方法かどうか私には分りませんが，少なくとも正しい答えが望んだ精度で帰ってくる様です(In[10]以降)．
その結果， 任意精度計算を行うにあたって精度の保持はユーザーが監視しないと行けない事が分りました．
また，監視を怠ると正しい結果が得られない事も分りました．
簡単な例として離散Fourier変換(FFTか？)と逆変換を繰り返す演算をさせます．
Wolfram Mathematica Document Foureirの説明にある通り，
Inputのlistが厳密な数値である場合，
N関数を適用し精度はMachinePrecisionに落とされます．
そのためFoureir関数は厳密精度計算を行ってくれません．
しかしながらユーザーがInputのlistにN関数を適用し精度をユーザー側でコントロールしたい状況が発生しますが，安易な方法では正しい結果が得られません．
その結果が In[3]です．
In[1]:= func1[y_, iter_, dps_: 10] := Module[{i, x},
x = y;
For[i = 0, i < iter, i++;
If[ Mod[i, 2] != 0, x = SetPrecision[Fourier[x], dps],
x = SetPrecision[InverseFourier[x], dps]];
];
Return[x]]
func0[y_, iter_] := Module[{i, x},
x = y;
For[i = 0, i < iter, i++;
If[ Mod[i, 2] != 0, x = Fourier[x], x = InverseFourier[x]];
];
Return[x]]
In[3]:= a = {0, 0, 1/3, 0}
res = func0[a, 100]
Precision[res[[1]]]
$MachinePrecision
b = N[{0, 0, 1/3, 0}, 10]
res = func0[b, 100]
Precision[res[[1]]]
Out[3]= {0, 0, 1/3, 0}
Out[4]= {0., 0., 0.333333, 0.}
Out[5]= MachinePrecision
Out[6]= 15.9546
Out[7]= {0, 0, 0.3333333333, 0}
Out[8]= {0.*10^5, 0.*10^4 + 0.*10^4 I, 0.*10^5, 0.*10^4 + 0.*10^4 I}
Out[9]= 0.
listに厳密な数値を代入した場合，正しい結果がMachinPrecisionで帰ってきます．一方listにN関数を作用させて作った初期条件の場合，全ての結果のPrecisionが0になっています．
これでは使い物になりません．１回のFourier関数の演算でどうやら10^2-10^10程度の桁落ちが発生し最終的にゼロが帰ってくる様です．
上記問題を回避する為にはSetPrecision関数を用いて，強制的に精度を維持する事です
(func1関数を参照)．
ただしこの方法が正当な回避方法かどうか私には分りませんが，少なくとも正しい答えが望んだ精度で帰ってくる様です(In[10]以降)．
In[10]:= a = {0, 0, 1/3, 0}
res = func1[a, 100, 10]
Precision[res[[3]]]
b = N[{0, 0, 1/3, 0}, 10]
res = func1[b, 100, 10]
Precision[res[[3]]]
Out[10]= {0, 0, 1/3, 0}
Out[11]= {0, 0.*10^-21 + 0.*10^-21 I, 0.3333333333,
0.*10^-21 + 0.*10^-21 I}
Out[12]= 10.
Out[13]= {0, 0, 0.3333333333, 0}
Out[14]= {0, 0.*10^-21 + 0.*10^-21 I, 0.3333333333,
0.*10^-21 + 0.*10^-21 I}
Out[15]= 10.
In[34]:= Exit[]
登録:
投稿 (Atom)