Programming in Mathematica

Mathematicaでのプログラミングの例．

FFT

In[130]:=
```FFT[data_, x_, m_] :=
Block[ {n, a, k},
n = Length[data];
a = N[ InverseFourier[data]/Sqrt[n] ];
Re[a[[1]]]+Sum[ 2 Re[a[[k+1]]] Cos[2 Pi k x / n] -
2 Im[a[[k+1]]] Sin[2 Pi k x / n] ,
{k,1,m} ]
];
PlotOver[data_, f_, x_] :=
Plot[ f, {x,0,Length[data]-1},
Epilog -> Map[Point, addX[data]] ];
addX[data_] := Transpose[{Range[0,Length[data]-1],data}];```
In[131]:=
```data = Table[
N[E^(-((2i/128-1)Pi)^2)+0.2Random[]],
{i,0,127}];
ListPlot[addX[data]];```
In[132]:=
`FFT[data, x, 3]`
Out[132]=
```                       Pi x                Pi x
0.37974 - 0.453637 Cos[----] + 0.21434 Cos[----] -
64                  32

3 Pi x                   Pi x
0.0629601 Cos[------] + 0.00224806 Sin[----] -
64                      64

Pi x                   3 Pi x
0.00250614 Sin[----] + 0.00300746 Sin[------]
32                      64```
In[133]:=
`PlotOver[data, %, x]`
Out[133]=
`-Graphics-`

Drawing Binary Trees

In[134]:=
```placeTree[node_] := {node,0} /; ! ListQ[node];
placeTree[{node_,left_,right_}] :=
Block[{lplace,rplace,lpos,rpos,sep},
lplace = placeTree[left];
rplace = placeTree[right];
lpos = Max[rposlist[lplace]];
rpos = Min[lposlist[rplace]];
sep = lpos-rpos+1;
lplace = moveright[lplace,-sep/2.];
rplace = moveright[rplace, sep/2.];
{node,0,lplace,rplace}
]```
In[135]:=
```lposlist[{node_,x_}] := {x};
lposlist[{node_,x_,lplace_,rplace_}] :=
{x,lposlist[lplace]};
rposlist[{node_,x_}] := {x};
rposlist[{node_,x_,lplace_,rplace_}] :=
{x,rposlist[rplace]};```
In[136]:=
```moveright[{node_,x_}, dx_] := {node,x+dx};
moveright[{node_,x_,lplace_,rplace_}, dx_] :=
{node, x+dx, moveright[lplace,dx], moveright[rplace,dx] }```
In[137]:=
```drawPlacedTree[{node_,x_}, y_] := {Point[{x,y}]};
drawPlacedTree[{node_,x_,lplace_,rplace_}, y_] :=
Join[
{ Point[{x,y}],
Line[{{x,y},{lplace[[2]],y-1}}],
Line[{{x,y},{rplace[[2]],y-1}}] },
drawPlacedTree[lplace,y-1],
drawPlacedTree[rplace,y-1]
]```
In[138]:=
```tree = placeTree[
{1,{1,1,
{1,1,1}},
{1,{1,{1,{1,1,1},
1},
1},
{1,1,1}}}]```
Out[138]=
```{1, 0, {1, -2.78125, {1, -3.53125},

{1, -2.03125, {1, -2.53125}, {1, -1.53125}}},

{1, 2.78125, {1, 1.59375,

{1, 0.71875, {1, -0.03125, {1, -0.53125},

{1, 0.46875}}, {1, 1.46875}}, {1, 2.46875}},

{1, 3.96875, {1, 3.46875}, {1, 4.46875}}}}```
In[139]:=
```Show[Graphics[
{PointSize[0.05],drawPlacedTree[tree,0]}
],PlotRange->All]```
Out[139]=
`-Graphics-`

Modifying Built-in Functions

In[140]:=
`E^(a+b Log[x])`
Out[140]=
``` a + b Log[x]
E```
In[141]:=
`FullForm[%]`
Out[141]=
`Power[E, Plus[a, Times[b, Log[x]]]]`

In[142]:=
```Unprotect[Power];
Power[E, a_+b_] := E^a E^b;
Power[E, a_ b_] := (E^b)^a;
Power[E, Log[a_]] := a;
Protect[Power];```
In[143]:=
`??Power`
Out[143]=
```x^y gives x to the power y.

Attributes[Power] = {Listable, OneIdentity, Protected}

0^0 = 1

E^((a_) + (b_)) := E^a*E^b

E^((a_)*(b_)) := (E^b)^a

E^Log[a_] := a

Power/: Default[Power, 2] := 1```
すると簡単化してくれる．
In[144]:=
`E^(a+b Log[x])`
Out[144]=
``` a  b
E  x```
もとに戻しておく．
In[145]:=
`Unprotect[Power]; Clear[Power]; Protect[Power];`

Perfect Hash Function

In[146]:=
`Needs["NumberTheory`NumberTheoryFunctions`"]`
In[147]:=
`pr1[n_] = n^2-n+41;`
In[148]:=
```perfectHash[keys_, pfunc_, x_] :=
Mod[
ChineseRemainderTheorem[
Range[Length[keys]], Map[pfunc, keys] ],
pfunc[x]
]```
2,4,6,8,10,12,19,36を1から8の整数に割り当てる関数を求める．
In[149]:=
```ks = {2,4,6,8,10,12,19,36};
h1[n_] = perfectHash[ks,pr1,n]```
Out[149]=
```                                 2
Mod[16599229625856279, 41 - n + n ]```
In[150]:=
`Map[h1, ks]`
Out[150]=
`{1, 2, 3, 4, 5, 6, 7, 8}`

Dept. CS / Faculty of Eng. / Kobe Univ. / Naoyuki Tamura