Language detail: Mathematica

Coverage: 19.05%
number of '+' ratings
contribution for coverage

Unsolved challenges

codes

Feed

Used modules

next >>

π (Nested Flatten)

せっかくなので数値積分させてみました.

Abs[Pi - npi]を評価すると2x10^(-100)でした.

この値は積分範囲の逆数程度になっている様子.

1
npi = NIntegrate[1/(1 + x^2), {x, -10^100, 10^100}, WorkingPrecision -> 200, MaxRecursion -> 50]

問題文の趣旨をまるで無視。

1
N[Pi,10000]
再帰を用いた迷路探索問題 (Nested Flatten)
再帰のない言語やあっても深さの制限が厳しい言語があるので、
手法は限定しないほうがいいのではないでしょうか
(構造体とかリストとかも)。
要は迷路が解ければいいわけで。

例えばグラフアルゴリズムを標準サポートしている処理系なら、
最短経路も特に難しいということはないかと
(愚直なエージェントを実装したいということがあるかもしれませんが)

num:座標を数字に変換する補助関数
pos:数字を座標に変換する補助関数

実行結果(Mathematicaなのでインデックスがずれます):
{{2,2},{2,3},{2,4},{2,5},{3,5},{4,5},{4,4},{4,3},{5,3},{5,2}}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
<<DiscreteMath`GraphPlot`;
<<DiscreteMath`Combinatorica`;

tmp={
      "******",
      "*8000*",
      "****0*",
      "**000*",
      "*90*0*",
      "******"};
maze=Characters/@tmp;
n=Length@First@maze;
num[i_,j_]:=n (i-1)+j
pos[x_]:={Quotient[x,n,1]+1,Mod[x,n,1]}

start=num@@First@Position[maze,"8"];
goal=num@@First@Position[maze,"9"];
arcs={};
Do[
    If[maze[[i,j]]!="*",
      If[maze[[i,j+1]]!="*",
        AppendTo[arcs,{num[i,j],num[i,j+1]}]];
      If[maze[[i+1,j]]!="*",
        AppendTo[arcs,{num[i,j],num[i+1,j]}]]],
    {i,1,Length@maze-1},{j,1,n-1}];

sol=ShortestPath[FromUnorderedPairs@arcs,start,goal];
If[Head@sol===ShortestPath || Length@sol==1,False,
  pos/@sol]
マルバツゲーム:賢いプレイヤー (Nested Flatten)
#6256(http://ja.doukaku.org/comment/6256/)の続きです

新たに実装したのは23行目以降
nextStates:新しい盤面のリストを作る補助関数
minimaxDecision:ミニマックス戦略
minimaxValue:ミニマックス値を計算する補助関数

マルがミニマックス戦略、バツがランダム・プレーヤーのとき、マルの9954勝0敗46分

Timing[
  result=Table[game[minimaxDecision,randomDecision],{10000}];
  Count[result,#]&/@{1,-1,0}
  ]

{65.656 Second, {9954, 0, 46}}

マルがランダム・プレーヤー、バツがミニマックス戦略のとき、バツの8024勝0敗1976分

Timing[
  result=Table[game[randomDecision,minimaxDecision],{10000}];
  Count[result,#]&/@{1,-1,0}
  ]

{68.484 Second, {0, 8024, 1976}}

マルバツともにミニマックス戦略のときは引き分け

game[minimaxDecision,minimaxDecision]

0
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
wins={{1,2,3},{4,5,6},{7,8,9},{1,4,7},{2,5,8},{3,6,9},{1,5,9},{3,5,7}};

isWinner[player_]:=MemberQ[Length[Intersection[player,#]]&/@wins,3]

judge[{p1_,p2_}]:=
  If[isWinner@p1,1,
    If[isWinner@p2,-1,
      If[Length@p1==5,0,Null]]]

operators[state_]:=Complement[Range@9,Flatten@state]

game[decision1_,decision2_]:=Module[{state={{},{}},result=Null},
    While[result===Null,
      If[Length@state[[1]]==Length@state[[2]],
        AppendTo[state[[1]],decision1[Sort/@state]],
        AppendTo[state[[2]],decision2[Sort/@state]]];
      result=judge@state];
    result]

randomDecision[state_]:=With[{x=operators@state},
    x[[Random[Integer,{1,Length@x}]]]]

nextStates[{p1_,p2_},ops_]:=
  Map[If[Length@p1==Length@p2,
      {Append[p1,#],p2}&,{p1,Append[p2,#]}&],ops]

Remove[minimaxDecision];
minimaxDecision[s_]:=
  Module[{ops=operators@s,vals,isMax=Length@s[[1]]==Length@s[[2]]},
    vals=minimaxValue[#,Not@isMax]&/@nextStates[s,ops];
    minimaxDecision[s]=ops[[Position[vals,If[isMax,Max,Min]@vals][[1,1]]]]]

minimaxValue[state_,isMax_]:=Module[{result=judge@state},
    If[result=!=Null,result,
      If[isMax,Max,Min][
        minimaxValue[#,Not@isMax]&/@nextStates[state,operators@state]]]]
マルバツゲーム (Nested Flatten)
wins:勝ちのパターン
isWinner:勝者かどうか
judge:ゲームの終了判定
operators:可能な手を返す
game:2つの戦略を対戦させる
randomDecision:ランダム・プレーヤーの戦略

対戦させると、マルの5877勝2835敗1288分

Timing[
  result=Table[game[randomDecision,randomDecision],{10000}];
  Count[result,#]&/@{1,-1,0}
  ]

{7.688 Second, {5877, 2835, 1288}}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
wins={{1,2,3},{4,5,6},{7,8,9},{1,4,7},{2,5,8},{3,6,9},{1,5,9},{3,5,7}};

isWinner[player_]:=MemberQ[Length[Intersection[player,#]]&/@wins,3]

judge[{p1_,p2_}]:=
  If[isWinner@p1,1,
    If[isWinner@p2,-1,
      If[Length@p1==5,0,Null]]]

operators[state_]:=Complement[Range@9,Flatten@state]

game[decision1_,decision2_]:=Module[{state={{},{}},result=Null},
    While[result===Null,
      If[Length@state[[1]]==Length@state[[2]],
        AppendTo[state[[1]],decision1[Sort/@state]],
        AppendTo[state[[2]],decision2[Sort/@state]]];
      result=judge@state];
    result]

randomDecision[state_]:=With[{x=operators@state},
    x[[Random[Integer,{1,Length@x}]]]]
ポリゴンを表示するプログラム (Nested Flatten)
以下の手順で描いています(Mathematica 6)。

1. 頂点を定義
2. ポリゴンを定義
3. 描画

マウスでドラッグすると回転します。

特別なライブラリは使っていません。
(Mathematica 5だと「<< RealTime3D`」が必要です。)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
a = {1, 0, 0};
b = {-1/2, Sqrt@3/2, 0};
c = {-1/2, -Sqrt@3/2, 0};
d = {0, 0, Sqrt@2};

tetra = {
   Polygon[{a, b, c}],
   Polygon[{a, b, d}],
   Polygon[{b, c, d}],
   Polygon[{c, a, d}]
  };

Show[Graphics3D[tetra]]
指定コマンドを別プロセスで起動 (Nested Flatten)
Cygwinシェル上でmath.exeを実行した結果

In[1]:= Run["ps"]
      PID    PPID    PGID     WINPID  TTY  UID    STIME COMMAND
     4876       1    4876       4876  con 1003 15:38:24 /usr/bin/bash
     1348    4876    1348        736  con 1003 17:16:08 /cygdrive/c/Program Files/Wolfram Research/Mathematica/5.2/math
     2368       1    2368       2368  con 1003 17:16:16 /usr/bin/ps
Out[1]= 0

これではリターンコードは戻り値として取得できるが、出力は表示されるだけで再利用できない。

仕方なく、出力を一度ファイルに入れることにする(Windows以外ならもっと簡単なはず)。
結果リストの1番目がリターンコード、2番目が出力

In[3]:= exec@"ps"

Out[3]= {0, {      PID    PPID    PGID     WINPID  TTY  UID    STIME COMMAND,
>          4876       1    4876       4876  con 1003 15:38:24 /usr/bin/bash,
>          1348    4876    1348        736  con 1003 17:16:08\
>      /cygdrive/c/Program Files/Wolfram Research/Mathematica/5.2/math,
>          4768       1    4768       4768  con 1003 17:17:33 /usr/bin/ps}}
1
2
3
4
5
6
7
8
exec[command_] := Module[{
    code,result,
    file = First@ReadList["!mktemp",String]
  },
  code = Run[command <> " > " <> file];
  result = ReadList[file,String];
  Run["rm "<>file];
  {code,result}]
魔方分割数 (Nested Flatten)
グループ未定の数字の最初のものpは、まだ決まってないグループ(X)に。
Xの残りの(n-1)個は、残っている数字から、和が(s/n-p)になるものとする。
最後まで行ったらカウンタに1加算。
初めに戻る。

(ComplementやSubsetsは常にソートされている)

In[3]:= f[3]
Out[3]= 2

In[4]:= f[4]
Out[4]= 392(Core2 6700で0.06秒)

In[5]:= f[5]
Out[5]= 3245664(Core2 6700で690秒)
1
2
3
4
5
6
f[n_] := f[Range[n^2], n, Total@Range[n^2]/n]

f[in_, n_, s_] :=
  If[Length@in == n, 1,
    Total[f[Complement[Rest@in, #], n, s] & /@ (
    Select[Subsets[Rest@in, {n - 1}], First@in + Total@# == s &])]]
最大公約数(除算禁止) (Nested Flatten)

最初の条件は、このコードだけで言えばクリアしています。

最後の条件への答えは、このコードだけで言えば0です。

多言語クックブックにこういうお題を載せる意図がよくわかりません。

1
GCD[Fibonacci@1999, Fibonacci@2000]
正整数のゲーデル数化? (Nested Flatten)

この手の問題はMathematicaだと簡単に書けますね。

1
goedel[n_] := Times @@ MapIndexed[Prime[#2[[1]]]^#1 &, IntegerDigits[n]]
PageRankの計算 (Nested Flatten)
グラフが疎な場合はSparseArrayを使うと劇的に速くなります。

変更点は2点です。
・行列をSparseArrayにする
・固有ベクトルを一つだけ求める
(絶対値最大固有値の固有ベクトルが返る)

使い方は変わりません。

ノード数1000、エッジ数10000のグラフで試したところ、
2桁ほど速くなりました。
(私のMathematica 5.2で30秒から0.38秒に)
1
2
3
4
5
6
7
8
pageRank2[graph_Graph] := Module[{
  aMatrix = SparseArray@ToAdjacencyMatrix@graph,
  eigenVector
},
  aMatrix = SparseArray@Transpose[(#/Total@#) & /@ aMatrix];
  eigenVector = First@Eigenvectors[N@aMatrix, 1];
  eigenVector/Total@eigenVector
]
行列やグラフの基本的な操作が用意されているので、
それらを使って実装しました。

隣接リストからそのまま計算してもいいのですが、
グラフを引数にする関数にしておいたほうが、
汎用性は高いでしょう。

aList = {{2, 3, 4, 5, 7}, {1}, {1, 2}, {2, 3, 5}, {1, 3, 4, 6}, {1, 5}, {5}};
graph = FromAdjacencyLists[aList, Type -> Directed];

このようにグラフを作ってから、

pageRank@graph

として計算します。結果は

{0.303514, 0.166134, 0.140575, 0.105431, 0.178914, 0.0447284, 0.0607029}

(主記憶で処理できないようなものには対応していません。)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
Needs@"DiscreteMath`Combinatorica`";

pageRank[graph_Graph] := Module[{
  aMatrix = ToAdjacencyMatrix@graph,
  eigenVector
  },
  aMatrix = Transpose[(#/Total@#)&/@aMatrix];
  eigenVector = First@Eigenvectors@N@aMatrix;
  eigenVector/Total@eigenVector
]
n日後を返す関数を返す関数 (Nested Flatten)
日本限定、nは0以上、西暦限定、暦はユリウス暦とグレゴリオ暦のみ、妥当な日付が入力されることする

・入力が1873年1月1日より前ならユリウス暦で計算
・・結果が1892年12月2日より後ならグレゴリオ暦に変換
・・そうでなければそのまま
・そうでなければグレゴリオ暦で計算

Mathematicaでは、1752年9月14日にユリウス暦からグレゴリオ暦に切り替わったことになっていることに注意

> nDaysLater[1]@{1752, 9, 2}

{1752, 9, 3} (異常なし)

> nDaysLater[1]@{1872, 12, 2}

{1873, 1, 1} (暦が切り替わった)

> nDaysLater[5]@{2007, 7, 20}

{2007, 7, 25} (異常なし)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
<<Miscellaneous`Calendar` (Version 6では「<<Calendar`」)

nDaysLater[n_] := Function[{d},
    Module[{tmp},
      If[{10000, 100, 1}. d < 18730101,
        tmp = DaysPlus[d, n, Calendar -> Julian];
        If[18721202 < {10000, 100, 1}. tmp,
          DaysPlus[{1873, 1, 1}, DaysBetween[{1872, 12, 3}, tmp, Calendar -> Julian]],
          tmp],
        DaysPlus[d, n]]]]
exp(pi * sqrt(n))が整数に近くなるnを探す (Nested Flatten)
1から199まで(Range)の整数から、選ぶ(Select)、
基準は、
「x=Exp[Pi Sqrt[n]]として(With)、Abs[x-Round@x]が10^-4以下のもの」

「10^-4を0.0001と書いてはいけない」という問題かな
1
2
3
4
Select[Range@199,
  Function[{n}, 
    With[{x = Exp[Pi Sqrt@n]},
      Abs[x - Round@x] <= 10^-4]]]
アレイのuniq (Nested Flatten)
マニュアルの方法をちょっとだけ改良しました
1
2
3
uniq[x_] := Reap[Sow[True, x], _, #1 &][[2]]

uniq@{3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9}
重複無し乱数 (Nested Flatten)
1~nまでのリストをRangeで生成して,
RandomSampleでシャッフルしました.

実行例:
In[2]:= bingo[10]
Out[2]= {6, 3, 10, 8, 1, 2, 9, 5, 4, 7}

RandomSampleで再利用無しの要素取得,
RandomChoiceで再利用OKの要素取得が
できるようになってます.
1
bingo[n_] := RandomSample[Range[n]];
全ての組み合わせ (Nested Flatten)
Tuplesというそのままの関数があります.
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
In[1]:= Tuples[{{1, 2, 3, 4}, {"a", "b", "c"}}]

Out[1]= {{1, "a"}, {1, "b"}, {1, "c"}, {2, "a"}, {2, "b"}, {2, 
  "c"}, {3, "a"}, {3, "b"}, {3, "c"}, {4, "a"}, {4, "b"}, {4, "c"}}

In[2]:= Tuples[{{0, 1}, {"a", "b"}, {"foo", "bar"}}]

Out[2]= {{0, "a", "foo"}, {0, "a", "bar"}, {0, "b", "foo"}, {0, "b", 
  "bar"}, {1, "a", "foo"}, {1, "a", "bar"}, {1, "b", "foo"}, {1, "b", 
  "bar"}}
/*コメント*/を取り除く (Nested Flatten)
普通に正規表現で...
1
2
removecomment[s_] := Module[{},
   StringReplace[s, RegularExpression["/\\*(.+?)(\\*/|$)"] -> ""]];
議席数をドント方式で (Nested Flatten)
ytakenakaさんのコメントの
http://homepage2.nifty.com/PAF00305/math/apportionment/node4.html
を参考に書いてみました.

Floor[定数k×得票数] の合計 == 議席数

を満たす最大のkを求めます.

実行例:
In[2]:=
  seki = 100;
  hyou = {123, 4, 56, 78};
  don[seki, hyou]
Out[2]:= {48, 1, 21, 30}
1
2
3
4
5
6
don[seki_, hyou_] := Module[{k},
   k = NMaximize[{
       k, Total[Floor[k hyou]] == seki
       }, k][[1]];
   Floor[k hyou]
   ];
ピラミッドを作る (Nested Flatten)
実行例:

In[2]:= pyramid[4]
   *
  ***
 *****
*******
1
2
3
4
5
pyramid[n_] := 
  Module[{}, 
   Print[StringJoin[Table[" ", {n - #}], Table["*", {2 # - 1}]]] & /@ 
     Range[n];
   ];
next >>

Index

Feed

Other

Link

Pathtraq

loading...