おしえて おしえて セジウィック さん
大昔に買ったはいいものの，全然読んでいなかったSegiwickのAlgorithm in Cをぱらぱら見てて思ったんですが，
日本語第1巻で，ユークリッドの互除法に触れてますよね．
ユークリッドの互除法って，gcd(a, b)なら
n_0 = a mod b n_1 = b mod n_0 n_2 = n_0 mod n_1 ...
とやって，n_xが0になるまで続け，n_(x-1)が最大公約数だっていうアルゴリズムだったと思います．とすると，
unsigned long long int gcd(unsigned long long int u, unsigned long long int v) { if(u<v){swap(u, v); } while(v>0) { u%=v; swap(u, v); } return u; }
というコードが自然だと思います．
ところがSegiwick本では
unsigned long long int gcd2(unsigned long long int u, unsigned long long int v) { while(u>0) { if(u<v){swap(u, v); } u-=v; } return v; }
というように，引き算を使っています．
gcd(a, b)において，a = 3 * b + nであるなら，
n = a-b-b-b n = a mod b
どちらも同じことなので，結局は同様の操作をしていることになりますが，
ループ回数のカウント，Cのclock関数で性能を比べてみると，
gcd(123456789012345678, 987654321098765432) Natural (?) way: STEP: 14 Clock: 0 gcd: 2 Segiwick's way: STEP: 13720420 Clock: 125 gcd: 2
というわけで前者の圧勝です．一体なぜSegiwickさんはこのような方式を取ったのか？
いや，全然分からないんですけどね．
昔は除算と剰余が非常に遅く，PentiumIIIで35倍くらいの開きがあったとかなかったとかどこかで目にした覚えがあります．
それが理由でしょうか？それにしたってループのステップ数違いすぎな気もします．
もう少し調べてから加筆するかもしれません．
ソートのアルゴリズムおまけ
まあ要は必要になったから，勉強がてらソートのアルゴリズムを実装しまくっているんですね．
というわけでおまけ．
template<class T> void Sort<T>::insertion_sort(T target[], int length) { T key = 0; for(int i=1; i<length; ++i) { key = target[i]; int j=i-1; while(j>=0 && target[j] > key) { target[j+1] = target[j--]; } target[j+1] = key; } } template<class T> void Sort<T>::shell_sort(T target[], int length, int distance) { T key; while(distance > 0) { for(int i=distance; i<length; ++i) { key = target[i]; int j=i; while(j>=distance && target[j-distance] > target[j]) { swap(target[j], target[j-distance]); j-=distance; } } distance/=2; } } template<class T> void Sort<T>::selection_sort(T target[], int length) { for(int i=0, min; i<length-1; ++i) { min = i; for(int j=i+1; j<length; ++j) { if(target[min] > target[j] ) { min = j; } } swap(target[i], target[min]); } }
何の変哲も無い挿入ソートと選択ソート．
あとシェルソート．シェルソートって挿入ソートの改良版だったんですね．知らなかった．
このアルゴリズムってなんだか直感に反して効率的といわれているので，なんでやねん！とか思ってたのですが，
挿入ソートを効率よくするために…とか言われると分かる気がします．
相変わらず動作保障は無いです．要素20の配列で動いたよ，程度．
ネットを徘徊していると，3-wayのquick sortの実装って余り見当たりませんな．
前にJavaで実装したものがあるので，気が向いたらn-wayにしてここに張ります．
ソートのアルゴリズム
ちょっと分けあってソートについて調べています．
いやはや，奥が深いんですな．
僕は中の下程度の知識しか無いし，今更そんなもん必要じゃないだろと思って大した勉強をしたことが無いんですが．
これまでbubble, insertion, selection, quick, heap, shell, mergeのような，基本的なものしか知りませんでした．
例えばbubble sortひとつとっても，ソート済みのindexを記憶しておくことで探索範囲を徐々に縮めていく方法とか（名前無いの？），
片一方からのみソートしていくんじゃあバランス悪いから交互に攻めようとかいってshaker sortとかいうのがあったり．
しかし僕は頭が豆腐でできているため，本で読んだだけじゃあよく分かりません．
なのでコーディングします．
template<class T> void Sort<T>::bubble_sort(T target[], int length) { for(int i=length-1; i>=1; --i) { for(int j=1; j<=i; ++j) { if(target[j-1] > target[j]) { swap(target[j-1], target[j]); } } } }
これが普通のバブルソート．セジウィック本とほぼ同じようなコードになってしまった.
templateになっているのはご愛嬌．
template<class T> void Sort<T>::bubble_sort2(T target[], int length) { for(int i=length-1; i>=1; --i) { int sorted = 0; for(int j=1; j<=i; ++j) { if(target[j-1] > target[j]) { swap(target[j-1], target[j]); sorted = j; } } if(!sorted){break;} } }
これが改良版バブルソート．おおー．こんないい方法があったのか．考えたことなかった．頭いいなぁ．
template<class T> void Sort<T>::shaker_sort(T target[], int length) { int top = 0, bottom = length-1, sorted = top; while(bottom > top) { for(int j=bottom; j>top; --j) { if(target[j-1] > target[j]) { swap(target[j-1], target[j]); sorted = j-1; } } top = sorted; if(top == bottom){ break;} sorted = top; for(int j=top; j<bottom; ++j) { if(target[j] > target[j+1]) { swap(target[j], target[j+1]); sorted = j+1; } } bottom = sorted; if(top == bottom){ break;} } }
で，これが最終形態のシェーカーソート．前の2つに比べてコードが長いですが，単に左右左右とバブルソート（２つ目のほう）をしているだけです．
いやあ，すごいなぁ．バブルソートを材料に，こんなこと考える人いるのか．感心しますね．世の中にはあたまがいい人であふれてます．
ちなみに僕は情報科学系の学部を出たわけでないんで，アルゴリズムにしてもプログラミングにしても独学なんですね．
だからこういうことに触れる機会が，その道の人に比べて圧倒的に少なかったのです．
いいですよね，情報系の人は．こういうこと勉強すれば単位になるんですから（酷）
自分，何トチ狂ってヘンテコな学科に行っちゃったんだろうとか思います．ちょっと後悔してますが，よかったこともあるのでそれはそれで．
ちなみにコードは写経したわけでないので，うまく動くかは保障外です．あ，swapは単なる交換ですよ．
そんなこんなで貴重な日曜日は過ぎてゆくのだ．
Problem 68, 69
Problem 68は手で解けます．というか手で解きました．
内側の5角形が6〜10の場合と外側の5つが6〜10というケースがありえる
外側が6〜10の方が連結した数字が大きくなるので，後者で考える．
とりあえずこの前提だけでマスは埋められると思います．
そんで，内側の数字の可能性は2通りある（裏表の関係）ということを考慮すれば解けます．
コードを書く方が面倒だったという．
次，記念すべきProblem 69（何が）
問題にある表では1〜10の中で6が最もn/φが大きい．
φがちいさくて，nが大きければ大きくなる．
φはnの倍数が多ければ多いほど小さくなる．
だから・・・ええと，言語不自由で言葉で説明できないのでコードを見てください（汗）
typedef unsigned long long uint64; typedef uint64 natural ; natural phi(natural n) { natural count = 1; for(natural i=2; i<n; ++i) { if(gcd(n,i)==1){ count++;} } return count; } natural f69() { double lgst = 0; natural lgstn = 0; natural diff = 6; for(natural i=6; i<=1000000; i+=diff) { double x = i/(double)phi(i); if(x > lgst) { lgst = x; lgstn = i; diff = lgstn; } } return lgstn; }
一応正しい結果がでたけど，一般性のある方法なのかは謎．そのうち証明を試みてみたいと思います．
くどいけどProject Euler 97
もうちょっとでLV2です．先は長い．
なんだか生活に支障が出てきたので解くスピード下げようかと．．．
とりあえず簡単な問題から解いていく方針に変更．
指標は解いたユーザー数．
97番は90番台なのにやたら解いた人が多かったので見てみたところ，とても簡単．
でっかい素数の下10桁を取得する．
typedef unsigned long long uint64; typedef uint64 natural ; natural f97() { natural n=28433; for(natural i=1; i<=7830457; ++i) { n*=2; n%=10000000000; } ++n; n%=10000000000; return n; }
ちなみに新しい問題が作られて(238番だっけ)，速攻で解かれているみたいですが，
僕はヒープ使い果たしたあたりで一旦やめました．もっとうまい方法があるのかぁ．すごいなぁ．
Project Euler 18と67
三角形の行列を下にたどり，通った数字の最大値を求める．
ご丁寧に，問題18では行列が小さいからbrute forceでもいけるけど，
問題67は同じ問題だけど行列がでかいからキツイぜ，とか書いてある．
これは面白い．当然一般的な解法をここで考えて二つとも解くに限る．
で，僕は動的計画法が得意なのでそれで解きました．
と言っても別に難しいことは無くて，例えば
01
02 03
04 05 06
07 08 09 10
という行列を考える．面倒なので問題のようにピラミッド型じゃなくて三角行列で．
上から初めて，あり得るパスの中で最大のものを保存してゆき，最終的に一番下に出てくるのが最大値候補群となる．
僕は説明が破滅的に苦手なので順を追ってみると
01
03 04
04 05 06
07 08 09 10
青の値を足したものが赤．
これは1番目の行から2番目の行まで行くルート候補の値になる．
次
01
03 04
06 08 09
07 08 09 10
青の値を赤に，緑の値を紫色に足した．
ただし，真ん中の水色の値は両方試して，大きいほうを選択する．
こうすると1番目の行から2番目の行まで行くルートの中で，値が大きいもののみが残る．
最終的に，
01
03 04
06 08 09
15 17 18 19
という行列が出来上がるので，一番下の行から最大のものを選ぶ．
というわけで,同時に2つ解きました．
コーディングは無意味なところにやたら凝ってしまい凄く時間が掛かってしまいました．
typedef unsigned long long uint64;
typedef uint64 natural ;
int index(int x, int y) // index of matrix. we don't use 0 as index
{
if( x <= 0 || y <= 0 || y < x)
{
return -1;
} else
{
return ((y * (y-1) / 2)) + x - 1;
}
}
void createDPMatrix(vector<natural>& s, vector<natural>& d, int x, int y)
{
natural source = d[index(x,y)];
if((index(x,y+1) == -1) || (index(x+1,y+1) == -1)){ return; }
if((index(x+1, y+1) > (int)d.size()-1) || index(x, y+1) > (int)d.size()-1){ return; }
natural left = source + s[index(x, y+1)];
natural right = source + s[index(x+1, y+1)];
if(left > d[index(x,y+1)])
{
d[index(x,y+1)] = left;
createDPMatrix(s, d, x, y+1);
}
if(right > d[index(x+1,y+1)])
{
d[index(x+1,y+1)] = right;
createDPMatrix(s, d, x+1, y+1);
}
}
natural f18()
{
ifstream in("in.txt", ios::binary | ios::in);
if(!in){ cout<<"Can't open file"<<endl; }
vector<natural> inmat;
natural buff;
int x=0, y=1;
natural max=0;
while(in>> buff)
{
inmat.push_back(buff);
if(y == x) { ++y; x=0;}
else { ++x; }
}
vector<natural> dpmat(inmat.size(),0);
dpmat[index(1,1)] = inmat[index(1,1)];
createDPMatrix(inmat, dpmat, 1, 1);
for(int i=1; i<=y; ++i)
{
if(dpmat[index(i, y)] > max)
{
max = dpmat[index(i, y)];
}
}
return max;
}
Project Eulerにて初めて再帰使ったかも．ウソかも．
index()とか，別にいらないんですが，この方が考えやすいと思います．凄く無駄ですが．
三角行列みたいな可変長行列はC++とかだと作りにくくて，
フォーラム見ていると正方行列をそのまま三角行列として扱っている人がちらほらいましたが，これはindexによる擬似三角行列になってます．
偉そうに言いながら，vectorの力借りまくりですが．
あ，フォーラムで下からやっていく方法が．これは気づかなかった．
動的計画法自体の計算量は同じだけど，その方が最後に候補から選ぶ処理が無くなるので賢いですね．精進が足りん．
P_E_17
とか何とか言って（16番の続き），17番．
one, two, three...one thousandまでの文字数をカウント．
うっわ，ひどい．問題も，僕のコードも．
natural f17()
{
natural sum = 0;
map<int, int> lnum;
map<int, int> part;
part.insert(make_pair(1, 3)); //one
part.insert(make_pair(2, 3)); //two
part.insert(make_pair(3, 5)); //three
part.insert(make_pair(4, 4)); //four
part.insert(make_pair(5, 4)); //five
part.insert(make_pair(6, 3)); //six
part.insert(make_pair(7, 5)); //seven
part.insert(make_pair(8, 5)); //eight
part.insert(make_pair(9, 4)); //nine
part.insert(make_pair(10, 3)); //ten
part.insert(make_pair(11, 6)); //eleven
part.insert(make_pair(12, 6)); //twelve
part.insert(make_pair(13, 8)); //thirteen
part.insert(make_pair(14, 8)); //fourteen
part.insert(make_pair(15, 7)); //fifteen
part.insert(make_pair(16, 7)); //sixteen
part.insert(make_pair(17, 9)); //seventeen
part.insert(make_pair(18, 8)); //eighteen
part.insert(make_pair(19, 8)); //nineteen
part.insert(make_pair(20, 6)); //twenty
part.insert(make_pair(30, 6)); //thirty
part.insert(make_pair(40, 5)); //forty
part.insert(make_pair(50, 5)); //fifty
part.insert(make_pair(60, 5)); //sixty
part.insert(make_pair(70, 7)); //seventy
part.insert(make_pair(80, 6)); //eighty
part.insert(make_pair(90, 6)); //ninety
part.insert(make_pair(100, 7)); //hundred
part.insert(make_pair(1000, 8)); //thousand
for(int i=1; i<1000; ++i)
{
int lettercount = 0;
int num = i;
int hundred = num / 100; num%=100;
int ten = num / 10 * 10; num%=10;
int one = num;
if(ten == 10) //10 ~ 19
{
lettercount += part.find(ten + one)->second;
} else if(one)
{
if(ten){ //20 ~ 99
lettercount += part.find(ten)->second;
lettercount += part.find(one)->second;
} else { //01 ~ 09
lettercount += part.find(one)->second;
}
} else if(ten)
{ //10 ~ 90 :step 10
lettercount+= part.find(ten)->second;
}
if(hundred){
if(ten || one)
{ // "3" means the count of "and"
lettercount += part.find(hundred)->second + part.find(100)->second + 3;
} else
{
lettercount += part.find(hundred)->second + part.find(100)->second;
}
}
sum+= lettercount;
}
sum+= part.find(1)->second + part.find(1000)->second;
return sum;
}どうしても，何と言うか一般的なとき方をしにくい．thirteen とか，thirtyをthirとteen / tyに分割するとかはありえるが，
別にそこまでしなくていいか，という感じで適当にやってしまった．もはや指折り数えているのと大差なくて切ない．
と思ったらフォーラムでは大体同じような解き方をしている人が大半みたいなのでまあいいか．
Project Euler 16
Project Euler，一体何のマラソンなんだかよく分からなくなってきました．
16番は2の1000乗の各桁を足すというもの．
Cとかだと組み込み型ではオーバーフローするから工夫が必要だよね，とかそういうことなのかなぁ？
微妙に出題の意図が分からない．
きっともっとHackな方法があるに違いない！とは思うけど，まだマラソンは始まったばかりなのでこの辺で体力を使うのはやめよう．
typedef unsigned long long uint64;
typedef uint64 natural ;
natural f16()
{
const int shoulder = 1000;
const int length = 1000;
int num[length];
int available = 0;
natural sum = 0;
for(int i=1; i<length; ++i)
{
num[i] = 0;
}
num[0] = 1;
for(int i=1; i<=shoulder; ++i)
{
for(int j=available; j>=0; --j)
{
num[j] *=2;
if( num[j] / 10 ){
num[j+1] += num[j] / 10;
num[j] = num[j] % 10;
}
}
if(num[available+1])
{
++available;
}
}
for(int i=0; i<=available; ++i)
{
sum+=num[i];
}
return sum;
}10進数桁ごとの配列にして大きな数を計算する．
前の日記で，typedefの部分をコピペしなかったので今回はコピペ．
この辺，BigDecimalのあるJavaやScript言語で書いたら一瞬なんだろうなぁ．
と思ったのであえてC++です．
Project Eulerとか
Project Eulerというサイトがあるらしい．
要は，数学チックな問題がいっぱいあって，プログラミングして解きましょうという問題集のサイト．
PKU Judge Onlineみたいな感じです．
こういうのハマるので困ります．
自分は低脳なので一つ一つ解くのに時間が掛かるのでさらに困ります．
3桁同士の整数の積となっている，最大の回文数を求める問題について，
フォーラムに投稿しようと思ったのですが，なぜかロックされているのでここで地味に書きます.
んー，こういうのって駄目なのかな？駄目だったら消しますね．
typedef unsigned long long uint64;
typedef uint64 natural ;
natural mirror(natural num)
{
natural mirror = 0;
do
{
mirror *= 10;
mirror += num % 10;
}while(num /= 10);
return mirror;
}
natural f4()
{
natural lterm, rterm;
natural lgst = 0;
natural product = 0;
for(lterm = 100; lterm < 1000; ++lterm)
{
for(rterm=100; rterm <= lterm; ++rterm)
{
product = lterm * rterm;
if(product == mirror(product))
{
if(lgst < product)
{
lgst = product;
}
}
}
}
return lgst;
}フォーラムで見た限りではf4()の中にある2重ループを100〜999と100〜999のレンジで回しているものが多かったけど，
掛け算は交換可能なのだから，こんな風にかけるんじゃないんですか？
って書こうと思ったけどロックされているので…
多分，こういうサイトを小学校くらいでやらせたら受けるんじゃないのかなぁ．
Pythonみたく分かりやすい言語を少し学ばせれば絶対食いつくと思うんですが．
少なくとも僕は食いついて睡眠時間が減っていますよ！
Project EulerのProblem12番
三角数で約数が500より多い（つまり5001以上）ものはいくつ？
まあ，三角数を求めるところは工夫の余地無い気がしますが，
約数を探すところで，ナイーブに1から順番に割れるかどうか？みたくやったらすごく遅くて回答出る前にこっちの我慢の限界．
最適化するな．まだするな．やれ．
というわけで最適化（？）
例えば10を2で割れたなら，その商である5でも割れるんだから，
何も10の約数を探すのに1,2,3,4,5,6,7,8,9,10で割ってみないで，
1で割れる→10でも割れる 10より小さい範囲を探す
2でも割れる→5でも割れる 5より小さい範囲を探す
3,4では割れない．よって4つ．
というコード．5秒．思ったより速くない．
natural f12()
{
natural triangle = 0;
int divnum = 0;
natural upper = 0;
for(natural i=1; ;++i)
{
triangle+=i;
divnum=0;
upper = triangle;
for(natural j=1; j<upper; ++j)
{
if(!(triangle%j))
{
upper = triangle/j;
divnum+=2;
}
}
if(divnum > 500)
{
break;
}
}
return triangle;
}小さい数字でちゃんと動くか知りません．
ちなみに，こういうのをネットで検索してヒントを得るのは主義に反しているので余りやらないんですが，
解いてから色々見ていると，約数じゃなくて素因数分解を前提にしている回答が多い気がしました．
でも，問題には28の場合，1,2,4,7,14,28と明らかに素因数分解ではないですし，どうなんでしょう．
それでもうまくいくものなんですかね.ちゃんと読んで無いですけど．
あと，親切にも日本語に問題を翻訳してくれているページがあるんですね．
僕の場合,英語を読むのに掛かる時間が結構長いので助かります．
いや，これも英語の練習だから！逃げちゃ駄目だ！
うーん，こんなところで詰まるのでは先が思いやられるなぁ．
別に目標とか無いけど．
