ラベル Algorithm の投稿を表示しています。 すべての投稿を表示
ラベル Algorithm の投稿を表示しています。 すべての投稿を表示

2010年6月22日火曜日

線形時間のソート

ソートアルゴリズムにおいて,ソート順が入力要素の比較にのみ基づいて決定されるものを比較ソート(comparison sort)といいます.比較ソートでは,n個の要素をソートするために最悪Ω(nlgn)回の比較を行わなければなりません.ということなので,定数倍を除くとマージソートやヒープソートより早い比較ソートは存在しないのです.

<!-- ここまでコピペ -->

そんな中,O(n)で動作するようなやばいソートもあります.例えば,計数ソート,基数ソート,バケットソートなど….

■計数ソート(counting sort)
このソートでは,配列aがn個の要素a0, a1, …an-1,を持っていた時,以下を仮定します.
0≦ai<k
で,次のように処理を進めます.
1.配列aにおいて,iに等しい要素の個数ciを計算
2.1.の結果を用いて,ciがi以下の要素の個数を含むように計算
3.配列aの各要素aiがどこに位置するかをcaiで調べ,出力配列bに順次代入.

k=O(n)の時,このソートの計算量はO(n)です.

□コード
// 計数ソート
void countingSort(int[] a, int k){
int n=a.length;
int[] b=new int[n];
int[] c=new int[k];
for(int i=0; i<n; i++)
c[a[i]]++;
// c[i]はiに等しい要素の個数を含んでいる
c[0]--;
for(int i=1; i<k; i++)
c[i]+=c[i-1];
// c[i]はi以下の要素の個数-1を含んでいる
for(int i=n-1; i>=0; i--)
b[c[a[i]]--]=a[i];
System.arraycopy(b, 0, a, 0, n);
}


■基数ソート(radix sort)
基数ソートは,i桁目に関するソートを最下位桁から最上位桁まで繰り返します.
ソートは安定ソートなら何でもいいのですが,各桁の値のとる範囲が小さければ計数ソートを使うのが妥当です.
要素数がn,各桁の値の範囲が0からk-1(k種類),桁数がdであるとします.各桁の計数ソートの実行時間は,O(n+k)となります.計数ソートはd回行われるので,基数ソートの実行時間はO(dn+kd)となります.dが定数,k=O(n)ならば,結局O(n)になるので,線形時間で動作するという事になります.

□適用例
ソートする配列
329457657839436720355

1桁目(1の位)に注目してソート
720355436457657329839

2桁目(10の位)に注目してソート.安定ソートを用いているので,下2桁はソートされています.
720329436839355457657

3桁目(100の位)に注目してソート.これで完全にソートされました.
329355436457657720839


□コード
各桁の値は,10進数ではなく2進数で計算しています.2進数では,各桁の値をシフト演算で求める事が出来る,という利点があります.
// 基数ソート
void radixSort(int[] a){
for(int i=0; i<32; i++)
// 安定ソートを用いてi桁目に関する配列をソートする
countingSort(a, i);
}

void countingSort(int[] a, int d){
int n=a.length;
int[] b=new int[n];
int[] c=new int[2];
for(int i=0; i<n; i++)
c[(a[i]>>d)&1]++;
c[0]--;
c[1]+=c[0];
for(int i=n-1; i>=0; i--)
b[c[(a[i]>>d)&1]--]=a[i];
System.arraycopy(b, 0, a, 0, n);
}


■バケットソート(bucket sort)
バケットソートでは,配列の各要素が区間[0, 1)に一様分布する実数だと仮定します.
配列aの要素数をnとします.
このソートでは,区間[0, 1)をn等分し,n個の要素を各バケットに分配します.
例えばn=10なら,
[0,0.1),[0.1,0.2),[0.2,0.3),[0.3,0.4),[0.4,0.5),[0.5,0.6),[0.6,0.7),[0.7,0.8),[0.8,0.9),[0.9,1).
その後,各バケットをソートし,順に繋げていくだけです.
この時のソートアルゴリズムは挿入ソートを用います.
挿入ソートの計算量はO(n2)(n個の要素に対して)ですが,色々あってほぼΘ(1)となります.さらに,バケット数がnですので,結果的にO(n)になるのですよ.

□コード
// バケットソート
void bucketSort(double[] a){
int n=a.length;
LinkedList<Double>[] list=new LinkedList[n];
for(int i=0; i<n; i++)
list[i]=new LinkedList<Double>();
for(int i=0; i<n; i++)
list[(int)(n*a[i])].add(a[i]);
for(int i=0, k=0; i<n; i++){
Double[] b=list[i].toArray(new Double[0]);
// リストlist[i]を挿入ソートでソートする
insertionSort(b);
for(double d : b)
a[k++]=d;
}
}

// 挿入ソート
void insertionSort(Double[] a){
int n=a.length;
for(int j=1; j<n; j++){
double key=a[j];
int i=j-1;
while(i>=0&&a[i]>key){
a[i+1]=a[i];
i--;
}
a[i+1]=key;
}
}


元ネタは,アルゴリズムイントロダクション第1巻より第9章 線形時間のソーティングでした.

2010年6月20日日曜日

順列生成

順列を生成.

再帰による辞書式順でない順列生成
public void generate(int[] a){
recursive(a,0);
}

private void recursive(int[] a,int j){
if(j==a.length){
System.out.println(Arrays.toString(a));
return;
}
for(int i=j;i<a.length;i++){
// a[i]とa[j]を入れ替え
int t=a[j]; a[j]=a[i]; a[i]=t;
recursive(a,j+1);
// a[i]とa[j]を入れ替え
t=a[j]; a[j]=a[i]; a[i]=t;
}
}

再帰による辞書式順順列生成

public void generate(int[] a){
recursive(a,0);
}

private void recursive(int[] a,int j){
if(j==a.length){
System.out.println(Arrays.toString(a));
return;
}
for(int i=j;i<a.length;i++){
// j~iを循環右ローテート
int t=a[i];
for(int k=i;k>j;k--)
a[k]=a[k-1];
a[j]=t;
recursive(a,j+1);
// j~iを循環左ローテート
for(int k=j;k<i;k++)
a[k]=a[k+1];
a[i]=t;
}
}

辞書式順でない順列生成では,単純に2つの要素を交換していますが,辞書式順の順列生成では,その操作を右ローテートに置き換えています.
右ローテートだと辞書式順になる理由は以下の通り.

a1a2…an
と要素が並んでいて,
ai>ai+1
とします.(色々説明が面倒なので要素aiを数とします)
aiを先頭に置いた時,辞書式順で初めに来る列(=最小値)は何かというと,
a1…ai-1aiai+1…an
のaiを先頭に持ってきた
ai[a1…ai-1ai+1…an]
になります.ちょうどこれは,元の列のa1…aiを右にローテートしたものです.
というわけで,i=1~nについてこの操作を実行すると,
a1[a2a3…an-1an]
a2[a1a3…an-1an]

an-1[a1a2…an-2an]
an[a1a2…an-2an-1]
とn個の列が出来ます.
それぞれの列について再帰させ,[]内でも同様の操作を行うと,1つの列あたりn-1個の列が出来ます.

a1[a2[a3a4…an-1an]]
a1[a3[a2a4…an-1an]]

a1[an-1[a2a3…an-2an]]
a1[an[a2a3…an-2an-1]]

[[]]の一番内側の列は,2項目までを決定した時の最小値です.
というわけで,n回再帰することで,辞書式順の順列が作れます.

2010年6月13日日曜日

ヒープソート

ヒープソート

■コード
int heapSize;
int parent(int i){ return (i-1)/2; }
int left(int i){ return i*2+1; }
int right(int i){ return i*2+2; }
void swap(int[] a, int i, int j){
int t=a[i]; a[i]=a[j]; a[j]=t;
}

/**
* iを根とする2分木がヒープとなるように,a[i]をヒープの中に“滑り落とす”
* 但し,left(i),right(i)を根とする2分木がそれぞれヒープであるとする
*/
void heapify(int[] a, int i){
int left=left(i);
int right=right(i);
int largest=i;
if(left<heapSize && a[left]>a[largest])
largest=left;
if(right<heapSize && a[right]>a[largest])
largest=right;
if(largest!=i){
swap(a, i, largest);
heapify(a, largest);
}
}

/** ヒープを構築 */
void buildHeap(int[] a){
for(int i=(a.length-1)/2; i>=0; i--)
heapify(a, i);
}

/** Heap sort(ヒープソート) */
public void heapSort(int[] a){
heapSize=a.length;
buildHeap(a);
for(int i=a.length-1;i>=1;i--){
swap(a, 0, i);
heapSize--;
heapify(a, 0);
}
}

■特徴
・ヒープ構造を用いる
・最悪計算量O(nlogn)
・データによる計算量の変化が小さい

2010年5月15日土曜日

マージソート

何となく作りたくなったので,マージソート.

■コード
/** Marge sort(マージソート) */
public void margeSort(int[] a){
int n=a.length;
if(n==1) return;
// b.length+c.length=a.length
int[] b=new int[n/2];
int[] c=new int[n-n/2];
System.arraycopy(a, 0, b, 0, b.length);
System.arraycopy(a, n/2, c, 0, c.length);
margeSort(b);
margeSort(c);
// ソートされたbとcをaにマージ
for(int i=0, j=0, k=0;i<n;i++)
if(k==c.length)
a[i]=b[j++];
else if(j==b.length)
a[i]=c[k++];
else if(b[j]<c[k])
a[i]=b[j++];
else
a[i]=c[k++];
}


■特徴
・分割統治法を用いたソートアルゴリズム
・最悪計算量O(nlogn)
・安定ソートである

■原理
ある配列(数列)aとbが以下であるとします.
a1, …, am
b1, …, bn
この時,
ai<=ai+1
bj<=bj+1
とすれば,つまり,aもbもソートされているとすれば,それらをマージした以下の配列(数列)cが作れる.
c1, …, cm+n
これを再帰的に適用します.

■適用例
8,3,4,6,2,5,7,1

まずは,"分割"します.

83462571

83462571

83462571

83462571

配列の要素が一つになるまで"分割"しました.要素数が1の配列が"ソートされている事"は,明らかですので,ここからはソートされている配列を"統治"していきます.

38462517

ここで,要素数が2の配列が4つ出来ましたが,それぞれソートされています.

34681257

12345678

というわけでゴール.正しくソートされました.

2010年4月17日土曜日

プログラム・プロムナード 最大長方形の面積

下記の連載「プログラム・プロムナード」には,有益なアルゴリズムが色々と載っているので実装してみました.

プログラム・プロムナード
http://www.ipsj.or.jp/07editj/promenade/index.html

2002年4月号 最大長方形の面積(和田英一)

import java.util.*;
import java.awt.*;

public class Main{

// 単純なループによる実装 O(n^6)
private int count1(int[][] a){
int res=0;
int m=a[0].length;
int n=a.length;
for(int y1=0;y1<n;y1++){
for(int x1=0;x1<m;x1++){
// 始点-(x1, y1)
for(int y2=y1;y2<n;y2++){
for(int x2=x1;x2<m;x2++){
// 終点-(x2, y2)
int c=0;
// (x1, y1)-(x2, y2)から構成される長方形内の1の個数を数える
for(int j=y1;j<=y2;j++)
for(int i=x1;i<=x2;i++)
if(a[j][i]==1)
c++;
// 長方形内のマスが全て1であり
// 長方形の大きさ(=c)が今までの最大値より大きい場合は更新
if(c==(x2-x1+1)*(y2-y1+1))
if(c>res)
res=c;
}
}
}
}
return res;
}

// 動的計画法による実装
private int count2(int[][] a){
int m=a[0].length;
int n=a.length;
int max=0;
LinkedList<Point>[][] lists=(LinkedList<Point>[][])new LinkedList[n][m];;
LinkedList<Point> temp=new LinkedList<Point>();

for(int j=0;j<n;j++)
for(int i=0;i<m;i++)
lists[j][i]=new LinkedList<Point>();

for(int j=0;j<n;j++){
for(int i=0;i<m;i++){
if(a[j][i]==0)
continue;
int up, left;
for(up=0;j-up>=0&&a[j-up][i]==1;up++)
;
for(left=0;i-left>=0&&a[j][i-left]==1;left++)
;
// 上隣=0,左隣=0
if((i==0||a[j][i-1]==0)&&(j==0||a[j-1][i]==0)){
// (1, 1)
lists[j][i].add(new Point(1, 1));
}
// 上隣=0,左隣=1
else if(j==0||a[j-1][i]==0){
// (左方に続く1の個数, 1)
lists[j][i].add(new Point(left, 1));
}
// 上隣=0,左隣=1
else if(i==0||a[j][i-1]==0){
// (1, 上方に続く1の個数)
lists[j][i].add(new Point(1, up));
}
// 上隣=1,左隣=1
else{
for(Point p:lists[j][i-1])
lists[j][i].add(new Point(p.x+1, Math.min(p.y, up)));
for(Point p:lists[j-1][i])
lists[j][i].add(new Point(Math.min(p.x, left), p.y+1));
}

// 重複要素の削除
for(Point p:lists[j][i])
if(!temp.contains(p))
temp.add(p);
lists[j][i].clear();

// qより大きい長方形が存在したらqをリストに追加しない
for(Point q:temp){
for(Point p:temp){
if(q!=p&&q.x<=p.x&&q.y<=p.y){
q=null;
break;
}
}
if(q!=null)
lists[j][i].add(q);
}
temp.clear();

for(Point p:lists[j][i])
max=Math.max(max, p.x*p.y);
}
}
return max;
}

public static void main(String[] args){
int[][] a;
a=new int[50][50];
for(int j=0;j<a.length;j++)
for(int i=0;i<a[0].length;i++)
a[j][i]=(int)(Math.random()*2);
System.out.println("1:"+new Main().count1(a));
System.out.println("2:"+new Main().count2(a));
}
}


マージしないでも動作するかな?と思って,マージの実装を外して検証してみましたが,リストに保存する長方形が爆発的に増えた為,メモリが足りませんでした.

2010年3月11日木曜日

深さ優先探索,動的計画法,メモ化再帰

ITmediaにて,chokudai(高橋直大)さん連載の「最強最速アルゴリズマー養成講座」の最新記事が公開されました.

最強最速アルゴリズマー養成講座:アルゴリズマーの登竜門、「動的計画法・メモ化再帰」はこんなに簡単だった (1/5) - ITmedia エンタープライズ
http://www.itmedia.co.jp/enterprise/articles/1003/06/news002.html

というわけで,格子の問題を,解説を元にサクッと組んでみました.


public class Main{
private int w, h;
private int[][] memo;

private int recursive1(int[][] grid, int x, int y){
if(y==0&&x==w-1){
return 1;
}
else if(y<0||x>=w||grid[y][x]==1){
return 0;
}
else {
return recursive1(grid, x+1, y)+recursive1(grid, x, y-1);
}
}

// 単純な深さ優先探索(DFS) O(2^(h+m))
private int count1(int[][] grid){
w=grid[0].length;
h=grid.length;
return recursive1(grid,0,h-1);
}

// 動的計画法(DP) O(h*w)
private int count2(int[][] grid){
w=grid[0].length;
h=grid.length;
int[] a=new int[w];

for(int i=0;i<w&&grid[h-1][i]==0;i++){
a[i]=1;
}

for(int j=h-2;j>=0;j--){
for(int i=1;i<w;i++){
if(grid[j][i]==1)
a[i]=0;
else
a[i]+=a[i-1];
}
}
return a[w-1];
}

private int recursive3(int[][] grid, int x, int y){
if(y==0&&x==w-1){
return 1;
}
else if(y<0||x>=w||grid[y][x]==1){
return 0;
}
else {
if(memo[y][x]!=0)
return memo[y][x];
else
return memo[y][x]=recursive3(grid, x+1, y)+recursive3(grid, x, y-1);
}
}

// メモ化再帰
private int count3(int[][] grid){
w=grid[0].length;
h=grid.length;
memo=new int[h][w];
return recursive3(grid,0,h-1);
}

public static void main(String[] args){
// 0:通過可能
// 1:通過不可能
int[][] grid=new int[][]
{
{0,0,0,0,0,0},
{0,1,0,0,0,0},
{0,0,0,0,1,0},
{0,0,1,0,0,0},
{0,0,0,0,0,0},
};
System.out.println("1:"+new Main().count1(grid));
System.out.println("2:"+new Main().count2(grid));
System.out.println("3:"+new Main().count3(grid));
}
}


深さ優先探索,動的計画法,メモ化再帰の3つです.意外なほど簡単に組めたので驚き.

「最強最速アルゴリズマー養成講座」の記事一覧は下記から.

「最強最速アルゴリズマー養成講座」最新記事一覧 - ITmedia Keywords
http://www.itmedia.co.jp/keywords/algorithmer.html

2010年3月6日土曜日

三分探索

三分探索.二分探索じゃなくて三分探索.

二分探索は,単調関数の零点を求めるのに使われますが,
三分探索は,凸関数の極値を求めるのに使われます.

参考

実数探索三種類解説 - nodchipの日記
http://d.hatena.ne.jp/nodchip/20090303/1236058357

3分探索 - ICPC突破専用ザク
http://d.hatena.ne.jp/ir5/20090630/1246349028

fを上に凸な関数とします.
区間[a, b]を3等分して[a, c1], [c1, c2], [c2, b]の3つの区間をつくります.
この時,極大値は,上の3区間の内のどれかにあり,以下が成立します.

・極大値が[a, c1]にある(Fig.1)
⇒f(c1)>f(c2)が成立.b=c2とする.
・極大値が[c1, c2]にある(Fig.2)
⇒a=c1となってもb=c2となってもよい .
・極大値が[c2, b]にある(Fig.3)
⇒f(c1)<f(c2)が成立.a=c1とする.



Fig.1


Fig.2


Fig.3

言い換えれば,
f(c1)>f(c2)⇒極大値は[a, c2]に存在
f(c1)<f(c2)⇒極大値は[c1, b]に存在

よって,各ステップ毎に探索空間を2/3にすることが出来ます.

これを実装すると下のようになります.


double search(double left, double right, int iteration){
for(int i=0;i if(f((left*2+right)/3)>f((left+right*2)/3))
right=(left+right*2)/3;
else
left=(left*2+right)/3;
}
return (left+right)/2;
}


・例1
f(x)=x-x
df/dx=(-ln(x)-1)x-x
ですので,
df/dx=0
とおけば,
x=1/e≒0.367879441
となります.

結果:
初期値
[0, 4]
100回反復計算させた時の解
x=0.3678794519063162

・例2
鋸波
x-floor(x)
連続関数ですが,x∈Zの時は微分不可能であり極大値をとります.

結果:
初期値
[0, 4]
100回反復計算させた時の解
x=1.9999999999999996
上手く計算できたようです.

・例3
δ関数っぽいもの
δ関数は,
δ(x)=0 (x≠0)
δ(0)=∞
ですが,それっぽいものとして,

f(x)=0 (x≠0)
f(0)=1
に三分法を適用.

結果:
初期値
[-2, 2]
100回反復計算させた時の解
x=1.9999999999999996

駄目でした.連続で無いとやはり無理ですね.