picasaアドレス変更？

いつもbloggerに画像を上げる時にはgoogleのサービスであるpicasaを使っていたんですけど
どうやら現在、一部を除いてブログに画像が表示されていません。

クリックでもリンク先へ飛べないので調べたところ、どうやら画像のアドレスが変わっているようです。
さすがに全部の画像を貼り直すのはとても面倒なので他の方法を考えたんですけど、
残念ながら手作業でリンク先を変更するしかないようです。

調べても引っ掛からなかったから他の人はこんな現象起こってないのかなぁ。
最近Androidケータイに変えてそこから一度更新したので、
それが悪さをしたんじゃないかと思う反面、アドレスの先頭が変わってるので単なる変更の可能性も……。

特にgoogleからのアナウンスはないみたいですね。
と言っても訳後でないと読めませんけど(笑)


すぐに戻すにはちょっと時間が足りないのでぼちぼちって感じですねー。
電験の終わる9/4日以降かな？


夏休みはいろいろする予定だったんだけど遊んで終わっちゃいそうだなぁー……。
そして電験の得点率が6割を越しませんorz
一教科受かればそれでいいんだけどなー。


ここでgalaxy SII tips。
galaxy SIIではホームボタンと電源ボタンでスクリーンショットが撮れます。
これはiPhoneと同じみたいですね。
初期から入っている操作なのでアプリのダウンロードの必要はありません。
操作は、ホームボタンを先に押し、すかさず押したまま電源ボタンを押すことで撮ることができます。

ちょっとしたポップアップが出るので撮れたかどうか確認できます。
画像はギャラリー内に。
でもあんまり使う機会がなさそうってのもあるんだよなぁ。

以上、今回も携帯からの投稿でした。

Android

先日、携帯電話がスマートフォン、サムスン製のgalaxySIIに変わりました。
現在のところ活用できるほど何も入れていませんが、
これからいろいろ入れていく予定です。

iPhoneでなくAndroidを選んだ理由はまあ、
よくある回線の問題や開発の容易さが大きかったかな。

iPhoneのアプリはMacintoshでしか開発できないからなぁ。
その分Androidはプラットホームを選ばないし簡単。
言語もjavaだし。

そんな感じでプログラミングの練習がてら
いつものオセロでも作ってみる感じです。


そのあとはチェスとかかなー？
もちろんそっちの専門ではないので対人打ちだけですけど。
Aiはかなり敷居が……。

他にはARとか触っていけたらいいなーなんて。

今回、ケータイからの投稿なのでこのくらいで。

キーボード掃除

今日から夏休みに入って時間が出来たので前からやろうと思ってたキーボードとの掃除を。

使ってるのがNEC LaVie LLのノートPCなので、 キーを掃除するためには結構面倒な感じで外さなきゃいけないんですよね。

ノートのキーボードはどうやらパンタグラフ（？）みたいな感じになっているようです。
これやる前に必ず写真撮っておきましょう。Insertあたりは結構覚えてないもんです。


壊さないように慎重に取り外し……。





十数分かけて全部取り終えました。
良く見るとスペースキーの場所に何か書いてます。



"H"と"I"かな？もしかしたら組み立て中にイニシャルを書いたとか？笑

ちなみに取り外しはキーの角に指を入れて左右片方ずつ取れば割れずに全部とれました。
シールを引きはがす感じ？

最初の方はストローを加工して水平に上に引っ張って取ってましたけど
時間かかるし手でいったら案外割れずに済んだってのが真実です。

まずはあまり使わないキーで試してみましょう。
（InsertとかPauseとか右Ctrlとか右Altとか……）


実は結構汚いキーボードの内側。



埃やら髪の毛やらがびっしりです。
これを歯ブラシで軽くこすりながら掃除機で吸うこと数分……。

結構きれいになりました。（しかし画像なし）
変な汚れは手元にあったエタノールで処理しました。
さっきのサイン（？）も消えました。

キーの方も同じくエタノール使いながら汚れを綿棒やらで取ります。

そして今までPCを使ってきた勘で組み直し。
ノート特有のONテンキーの数字に助けられつつ……なんとか組み直し完成。




こんな感じになりました。
正解と見比べてみると……。

Enterすぐ左のカッコふたつのキーが逆になってました。笑
はじまりが下にあるって見ただけでなんか違和感あるのに……。

キーを取ってつけ直し。
そんな感じでキーボードの掃除終了でした。

キーの装着の方は、片方のツメに見ながらはめ込むと
反対側は位置が固定されるのでそのまま押すだけで入ります。

あとEnterキーは出来れば外さない方がよかったかな。
補助用（？）にShiftとかEnterなどの大きなキーには金属の棒が入ってるんだけど、
それをつけるのが至難でした。

まわりからも掃除機で吸えるのでEnterだけ付けたままでもよかったなー。

それと掃除機を近づけ過ぎると固定用のプラスティックも外れたのでそこも要注意かな。


夏休み初日に何やってんだろ（笑）

Windowsプログラミングによるオセロ

学校の課題の休憩中にWinAPIの練習がてらいつものオセロを作ってみました。

Drawing関数はInvalidateRect()とかでWM_PAINT内で処理させた方がよかったかな。

最初に何よりも迷ったのがVisualC++でのWinプロジェクトの作成。
コンソールだとウィンドウとか出せないので、Winプロジェクトを作る必要があるんだけれど、
わざわざそんな初歩的なこと書いてないし、
なにより参考HPが基本的に以前のVerのVisualStudioだったりとか。

新規プロジェクト→Win32→Win32プロジェクトから作れます。
コンソールアプリケーションじゃだめ。

中身自体は空のプロジェクトでOK。そのあといつものように.cppを追加しとけば大丈夫。
とりあえずAPIの基本的な動き方はわかった（つもり）ので、
これからチェスとか作れればいいなぁなんて。

ちょっと前にもあったけどネタ。
http://www.gizmodo.jp/2011/06/667.html

これに対応する予定は今のところありません（笑）


以下リスト。

#include<windows.h>
#include <tchar.h>
#include<stdlib.h>
#include<stdio.h>


//const
#define WINDOW_WIDTH (430)
#define WINDOW_HEIGHT (340)


//grobal
RECT g_windowPos;
int board[10][10] = {0};
int buf_board[10][10] = {0};
int rev_board[10][10] = {0};
static int act = 1,wait = 2;//1 = player1,2 = player2
int rev_act = 2,rev_wait = 1;


//Prototype
HWND Create(HINSTANCE hInst);
LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wp,LPARAM lp);
void OthelloInit(HWND hWnd);
void Drawing(HWND hWnd);
int Algo(int x,int y,HWND hWnd);
void Check(int x,int y,int i,int j,int *flag);
void REV(int flag,HWND hWnd);
int NUMSTONE(int color);


//Start
int WINAPI WinMain(HINSTANCE hInst,
HINSTANCE hPrevInst,
LPSTR lpCmdLine,
int showCmd)
{
HWND hWnd;
MSG msg;


//Create Window
hWnd = Create(hInst);


if ( hWnd == NULL){
MessageBox(NULL,_T("Creating Window is failed"),_T("ERROR"),MB_OK);
return 1;
}


//Show Window
ShowWindow( hWnd, SW_SHOW);
UpdateWindow(hWnd);


OthelloInit(hWnd);


//Message Loop
while( 1 ){
BOOL ret = GetMessage( &msg,NULL,0,0);//Get Message
if( ret == 0 || ret == -1){
//if get Message of Shutdown,
//or failed GetMessage() (return -1),EndLoop
break;
} else {
//Message
TranslateMessage( &msg );
DispatchMessage( &msg );
}
}
return 0;
}


HWND Create(HINSTANCE hInst){
WNDCLASSEX wc;


// Window Class
wc.cbSize = sizeof(wc); //const size
wc.style = CS_HREDRAW | CS_VREDRAW; //style
wc.lpfnWndProc = WndProc; //window procedure
wc.cbClsExtra = 0; //Extra Info1
wc.cbWndExtra = 0; //Extra Info2
wc.hInstance = hInst; //Instance handle
wc.hIcon = (HICON)LoadImage( //Icon
NULL,MAKEINTRESOURCE(IDI_APPLICATION),IMAGE_ICON,
0,0,LR_DEFAULTSIZE | LR_SHARED
);
wc.hIconSm = wc.hIcon; //child Icon
wc.hCursor = (HCURSOR)LoadImage( //Mouse Cursor
NULL,MAKEINTRESOURCE(IDC_ARROW),IMAGE_CURSOR,
0,0,LR_DEFAULTSIZE | LR_SHARED
);
wc.hbrBackground = (HBRUSH)GetStockObject(WHITE_BRUSH);//Window Background
wc.lpszMenuName = NULL; //Menu Name
wc.lpszClassName = _T("Default Class Name");//Window Class Name


//Submit Window Class
if ( RegisterClassEx( &wc ) == 0){
return NULL;
}


//Window's location
g_windowPos.left = (GetSystemMetrics( SM_CXSCREEN ) - WINDOW_WIDTH )/2;
g_windowPos.top = (GetSystemMetrics( SM_CYSCREEN ) - WINDOW_HEIGHT )/2;
g_windowPos.right = g_windowPos.left + WINDOW_WIDTH;
g_windowPos.bottom = g_windowPos.top + WINDOW_HEIGHT;


//Create Window
return CreateWindow(
wc.lpszClassName, //Window Class Name
_T("Sample Program"), //Title bar
WS_OVERLAPPEDWINDOW & ~WS_THICKFRAME, //Window type
g_windowPos.left, //X
g_windowPos.top, //Y
WINDOW_WIDTH, //Window Width
WINDOW_HEIGHT, //Window Height
NULL, //Parent Window handle
NULL, //Menu Handle
hInst, //Instance Handle
NULL //Extra Data
);
}


//Window Procedure
LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wp,LPARAM lp)
{
static int x,y;
static BOOL draw = FALSE;
HWND hButtonWnd;
HINSTANCE hInst;
hInst = (HINSTANCE)GetWindowLong(hWnd,GWL_HINSTANCE);


switch( msg ){
case WM_LBUTTONUP:
x = LOWORD(lp);
y = HIWORD(lp);
if(x >= 20 && x <= 55+34*7 && y >= 20 && y <= 55+34*7){
Algo(x,y,hWnd);
Drawing(hWnd);
}
return 0;
case WM_RBUTTONUP:
if(act == 1){
act = 2;
wait = 1;
}else if(act == 2){
act = 1;
wait = 2;
}
Drawing(hWnd);
return 0;
case WM_CREATE:
hButtonWnd = CreateWindow(
_T("BUTTON"),_T("UNDO"),
WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
310,20,100,30,hWnd,(HMENU)1000,hInst,NULL);
hButtonWnd = CreateWindow(
_T("BUTTON"),_T("RESTART"),
WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
310,260,100,30,hWnd,(HMENU)2000,hInst,NULL);
break;
case WM_COMMAND:
switch(LOWORD(wp)){
case 1000:
REV(2,hWnd);
Drawing(hWnd);
break;
case 2000:
OthelloInit(hWnd);
return 0;
}
break;
case WM_CLOSE:
if( MessageBox(hWnd,_T("shutdown?"),_T("checking"),MB_YESNO ) == IDNO){
return 0;
}
break;
case WM_DESTROY://when Window Destroy
PostQuitMessage(0);
return 0;
}


//other Messages is DEFAULT
return DefWindowProc(hWnd,msg,wp,lp);


}


void OthelloInit(HWND hWnd){
int i,j;


act = 1;
wait = 2;
for (i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
board[i][j] = 0;
}
}


for (i = 0;i <; 10;i++){
board[i][0] = 8;
board[i][9] = 8;
board[0][i] = 8;
board[9][i] = 8;
rev_board[i][0] = 8;
rev_board[i][9] = 8;
rev_board[0][i] = 8;
rev_board[9][i] = 8;


}
board[4][5] = 1;
board[5][4] = 1;
board[4][4] = 2;
board[5][5] = 2;
rev_board[4][5] = 1;
rev_board[5][4] = 1;
rev_board[4][4] = 2;
rev_board[5][5] = 2;


EnableWindow(GetDlgItem(hWnd,1000),FALSE);


Drawing(hWnd);
}


void Drawing(HWND hWnd){
int i,j;
char word[50];
RECT rc;
HDC hDC;
HBRUSH hBrush;
HPEN hPen;
WCHAR prev_word[50];


hDC = GetDC(hWnd);
/*
GetClientRect(hWnd,&rc);
InvalidateRect(hWnd,&rc,FALSE);
*/
hPen = CreatePen(PS_NULL,0,0);
hBrush = CreateSolidBrush(RGB(255,255,255));
SelectObject(hDC,hPen);

SelectObject(hDC,hBrush);
Rectangle(hDC,310,160,400,200);
DeleteObject(hPen);
DeleteObject(hBrush);

//GetStockObject(BLACK_PEN);
SelectObject(hDC,(HPEN)GetStockObject(BLACK_PEN));
for (i = 0;i < 8;i++){
for (j = 0;j < 8;j++){
hBrush = CreateSolidBrush(RGB(0,255,0));
SelectObject(hDC,hBrush);
Rectangle(hDC,20+34*i,20+34*j,55+34*i,55+34*j);
DeleteObject(hBrush);
switch(board[i+1][j+1]){
case 1:
hBrush = CreateSolidBrush(RGB(0,0,0));
SelectObject(hDC,hBrush);
Ellipse(hDC,22+34*i,22+34*j,53+34*i,53+34*j);
DeleteObject(hBrush);
break;
case 2:
hBrush = CreateSolidBrush(RGB(255,255,255));
SelectObject(hDC,hBrush);
Ellipse(hDC,22+34*i,22+34*j,53+34*i,53+34*j);
DeleteObject(hBrush);
break;
}
}
}
SetRect(&rc,310,100,400,150);
switch(act){
case 1:
DrawText(hDC,_T("Black turn"),-1,&rc,0);
break;
case 2:
DrawText(hDC,_T("White turn"),-1,&rc,0);
break;
}

sprintf_s(word,30,"black = %d\nwhite = %d",NUMSTONE(1),NUMSTONE(2));
mbstowcs_s(0,prev_word,30,word,_TRUNCATE);
SetRect(&rc,310,160,400,200);
DrawText(hDC,prev_word,-1,&rc,DT_WORDBREAK);
ReleaseDC(hWnd,hDC);
}


int Algo(int x,int y,HWND hWnd){
int locx,locy,i,j,flag;
for(i = 0;i < 8;i++){
if (x >= 20+34*i && x <= 55+34*i){
locx = i+1;
}
if (y >= 20+34*i && y <= 55+34*i){
locy = i+1;
}
}


if(board[locx][locy] != 0){
return -1;
}


if(locx >= 1 && locx <= 8 && locy >= 1 && locy <= 8){


/*置き石判定*/
flag = 0;
REV(flag,hWnd);
for (i = -1;i <= 1;i++){
for(j = -1;j <= 1;j++){
Check(locx,locy,i,j,&flag);
}
}
if(flag != 0){
REV(flag,hWnd);
board[locx][locy] = act;
if(act == 2){
rev_act = act;
rev_wait = wait;
wait = act;
act = 1;
} else {
rev_act = act;
rev_wait = wait;
wait = act;
act = 2;
}
}
}
return 0;
}


void Check(int x,int y,int i,int j,int *flag){
int k,num;
if(board[x+i][y+j] == wait){
num = 0;
while(1){
if(board[x+i][y+j] == wait){
num++;
} else if (board[x+i][y+j] == act){
for(k = 0;k < num;k++){
board[x-i*k][y-j*k] = act;
}
*flag = 1;
break;
} else if(board[x+i][y+j] == 8 || board[x+i][y+j] == 0){
break;
}
x += i;
y += j;
}
}
}


void REV(int flag,HWND hWnd){
int i,j;
if(flag == 0){
for(i = 0;i < 10;i++){
for(j = 0;j &t; 10;j++){
buf_board[i][j] = board[i][j];
}
}
} else if(flag == 1){
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
rev_board[i][j] = buf_board[i][j];
}
}
EnableWindow(GetDlgItem(hWnd,1000),TRUE);
} else if(flag == 2){
act = rev_act;
wait = rev_wait;
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
board[i][j] = rev_board[i][j];
}
}
EnableWindow(GetDlgItem(hWnd,1000),FALSE);
}
}

int NUMSTONE(int color){
int i,j,num;
num = 0;
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
if(board[i][j] == color){
num++;
}
}
}
return num;
}

課題・ガウスの消去法3（ピボット）

ガウスの消去法1
ガウスの消去法2

引き続いてガウスの消去法。

前回はいったん答えが出たように思いましたが、このような

場合には正常な値が出ませんでした。

PRINTを前進消去の直後に入れると解るんですけど、
2行目以下の値がERRORになっています。

これは手計算してみると解るんですけど、
2行2列が0になってしまっています。

これによって2行目を2行2列で割って、1にするという計算が誤作動を起こします。
（定数割る0なので）

これを回避するために、ピボット（1にしたい対角要素）に合わせて
行を入れ替えるという方法があります。


↓



まず、このピボットの選択のためには今見ている行以下の行で、
見たい列の中の最大値がある行を現在の行と入れ替える方法を使います。

今の場合、見ている行は2行2列で、2行目以下（つまり2行目と3行目）で
2列目が最大になる行を選べばいいことになります。

ややこしいですけど2列目が一番大きい3行目を2行目と入れ替えるということです。


前進消去のループのしょっぱなには今見ている行列をbufに入れて、
bufで割ってしまっているのでピボットの選択はこれより前、
かつfor(i=0~)よりあとに入れます。


この式を例にとった場合、列の最大値をa_maxとでも置くと、

if(a_max > a_matrix[0][0]){
a_max = a_matrix[0][0];
}

if(a_max > a_matrix[1][0]){
a_max = a_matrix[1][0];
}
if(a_max > a_matrix[2][0]){
a_max = a_matrix[2][0];
}
//ここまで1行目(i=0)の場合

if(a_max > a_matrix[1][1]){
a_max = a_matrix[1][1];
}
if(a_max > a_matrix[2][1]){
a_max = a_matrix[2][1];
}
//ここまで2行目(i=1)の場合

if(a_max > a_matrix[2][2]){
a_max = a_matrix[2][2];
}
//ここまで3行目(i=2)の場合

次にこれを簡略化すると、jはこの時使用されていないのでjを使うと、

for(j = i;j < 3;j++){
if(a_max > a_matrix[j][i]){
a_max = a_matrix[j][i];
}
}

となります。これでその列の最大値が求まりましたが、
このままでは何行目にその最大値があったのかわかりません。
（最大値があった行と、もとの行を入れ替えるために必要）

なので適当な変数にその「最大値があった行」をいれておきます。
ここでは使われない(予定）のkをつかいます。
（できるだけ宣言する変数は少ない方がいいので、いままで宣言した変数の中で
その時値を変えても大丈夫な変数があればそれを使います。）

for(j = i;j < 3;j++){
if(a_max > a_matrix[j][i]){
a_max = a_matrix[j][i];
k = j;
}
}

これで最大値のある行が把握できました。

これからこのk行とi行とを入れ替えます。
入れ替えには1つ目を適当な変数に一度格納しておき、2つ目を1つ目に入れた後に
適当な変数から2つ目にいれるという方法を使います。

もとの行はi行、入れ替える行はk行なので

for(j = 0;j < 3;j++){
buf = a_matrix[i][j];
a_matrix[i][j] = a_matrix[k][j];
a_matrix[k][j] = buf;
}
buf = b_matrix[i];
b_matrix[i] = b_matrix[k];
b_matrix[k] = buf;

となります。

これらをすべてリストに直したものが以下にあります。
a_maxはbufで代用できるのでbufに変更、
for(i=0~)の最初に最大値の初期化としてbuf = 0;が入っています。

そして最大値buf（元a_max）は、マイナスであろうと0でなければいいので
絶対値をfabs()でとっています。fabs()はmath.hを#includeすれば使えます。

他にbuf(a_max)が0の時にERRORとして終了させると、
連立方程式が不成立の時に強制的に終了させることが出来ます。

以上、おつさまですた。

#include<stdio.h>
#include<math.h>

/*prototype*/
int PRINT(double a_matrix[3][3],double b_matrix[3]);

int main(void){
int i,j,k;
double a_matrix[3][3] = {{1,1,1},{1,1,3},{2,4,5}};
double b_matrix[3] = {6,8,19};
double buf;

PRINT(a_matrix,b_matrix);

//前進消去
for(i = 0;i < 3;i++){
buf = 0;
for(j = i;j < 3;j++){
if(buf < fabs(a_matrix[j][i])){
buf = a_matrix[j][i];
k = j;
}
}

for(j = 0;j < 3;j++){
buf = a_matrix[i][j];
a_matrix[i][j] = a_matrix[k][j];
a_matrix[k][j] = buf;
}
buf = b_matrix[i];
b_matrix[i] = b_matrix[k];
b_matrix[k] = buf;

buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] / buf;

for(j = i + 1;j < 3;j++){
buf = a_matrix[j][i];
for(k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

//後退代入
for (i = 2;i >= 1;i--){
for (j = i-1;j >= 0;j--){
buf = a_matrix[j][i];
for (k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

PRINT(a_matrix,b_matrix);
return 0;
}

int PRINT(double a_matrix[3][3],double b_matrix[3]){
int i,j;

for(i = 0;i < 3;i++){
for(j = 0;j < 3;j++){
printf("%6.4f ",a_matrix[i][j]);
}
printf("%6.4f\n",b_matrix[i]);
}
printf("\n");
return 0;
}

課題・ガウスの消去法2（後退代入）

ガウスの消去法1
ガウスの消去法3

前回に引き続き、ガウスの消去法。
前進消去が終わって上三角行列になったので、ここからは答えを求める後退代入。
用語はあっているか知りませんけど。

現時点でここまで求めることが出来ています。


ここから3行目を利用して2行目と1行目を消していきます。

1. 2行目3列が0になるように3行目を定数倍して2行目から引く
   （2行目マイナス、3行目3列かける2行目）
2. 1行目3列が0になるように3行目を定数倍して2行目から引く
   （1行目マイナス、3行目3列かける1行目）
3. 1行目2列が0になるように2行目を定数倍して1行目から引く
   （1行目マイナス、2行目2列かける1行目）

上の式が、
1.

2行目3列が0になって、

2.

1行目3列が0になって、

3.

1行目2列が0になり、答えが求まります。


前進消去と似たような感じなので手早くプログラムを羅列していきます。

/*1*/
buf = a_matrix[1][2];
a_matrix[1][0] = a_matrix[1][0] - a_matrix[2][0] * buf;
a_matrix[1][1] = a_matrix[1][1] - a_matrix[2][1] * buf;
a_matrix[1][2] = a_matrix[1][2] - a_matrix[2][2] * buf;
b_matrix[1] = b_matrix[1] - b_matrix[2]* buf;

/*2*/
buf = a_matrix[0][2];
a_matrix[0][0] = a_matrix[0][0] - a_matrix[2][0] * buf;
a_matrix[0][1] = a_matrix[0][1] - a_matrix[2][1] * buf;
a_matrix[0][2] = a_matrix[0][2] - a_matrix[2][2] * buf;
b_matrix[0] = b_matrix[0] - b_matrix[2] * buf;

/*3*/
buf = a_matrix[0][1];
a_matrix[0][0] = a_matrix[0][0] - a_matrix[1][0] * buf;
a_matrix[0][1] = a_matrix[0][1] - a_matrix[1][1] * buf;
a_matrix[0][2] = a_matrix[0][2] - a_matrix[1][2] * buf;
b_matrix[0] = b_matrix[0] - b_matrix[1] * buf;

これを前までのプログラムの下に入れてあげると（すべてのforループの外）、

x1 = 3,
x2 = 2,
x3 = 1

という答えになります。

ここから、今まで書いたプログラムの中に上手く入るようにしていきます。
このくらい短いプログラムの省略だとまず恒常的に、一番回るのが遅いループを探します。

この場合、buf = a_matrix[1][2]の右側の配列番号が、手順3まで変化しません。
反対に左側の部分は、手順2で0に変化しています。
この一番外のループをiにします。

次に、2番目に遅く回るループを探します。
さっきの左側の部分は手順が次に行くまで変化しません。これをjとします。

残りは手順内で0~2と変化するのでこれをkとします。

するとこのプログラムは、

for (i = 2;i >= 1;i--){
for (j = i-1;j >= 0;j--){
buf = a_matrix[j][i];
for (k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}

となります。

前回を含めたプログラムリストに直すと、

#include

/*prototype*/
int PRINT(double a_matrix[3][3],double b_matrix[3]);

int main(void){
int i,j,k;
double a_matrix[3][3] = {{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3] = {9,17,8};
double buf;

PRINT(a_matrix,b_matrix);

/*前進消去*/
for(i = 0;i < 3;i++){
buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] / buf;

for(j = i + 1;j < 3;j++){
buf = a_matrix[j][i];
for(k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

PRINT(a_matrix,b_matrix);

/*後退代入*/
for (i = 2;i >= 1;i--){
for (j = i-1;j >= 0;j--){
buf = a_matrix[j][i];
for (k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

PRINT(a_matrix,b_matrix);
return 0;
}

int PRINT(double a_matrix[3][3],double b_matrix[3]){
int i,j;

for(i = 0;i < 3;i++){
for(j = 0;j < 3;j++){
printf("%6.4f ",a_matrix[i][j]);
}
printf("%6.4f\n",b_matrix[i]);
}
printf("\n");
return 0;
}

これで連立一次方程式のガウス消去法による解が求まりました。

しかしこれだけでは

のような連立一次方程式の解は得られません。
次はピボットについて書く予定です。

おつさまですた。

課題・ガウスの消去法1（前進消去法）

ガウスの消去法2
ガウスの消去法3

前回から授業が1回空いて、今回はガウスの消去法です。
2つ前はシンプソン公式による近似で、1つ前は改良オイラー法による常微分方程式でした。
ガウスの消去法は前進消去とか後退代入とかのやつです。
今回は前進消去だけを扱って、後退代入は次の記事になるとおもいます。

今回はこの式を題材に使います。

まあ簡単な連立一次方程式なんですけど、項の数が多くなってくると手計算は大変ですよね、
ということでプログラム化です。

行列についてはあまり言及しません。
学校で行くと去年の前期のことなんで忘れてるってのもありますけど＾＾；
計算自体は線形代数です。

行列式の計算は

となることを利用しています。

まずは1行目を左端の2で割って、

とします。（左端を1にしたい）

次にその1行目に、「1行目の左の数が2行目の左の数になるように」定数nをかけます。
そして2行目から、n倍した1行目を引きます。するとこのように、

2行目の左端が0になりました。
（同じように3行目も1行目をn倍して引いてます。）


1行目の左端が1に、それ以外の左端が0になったので、
次は2行目2列が1になるように同じように定数（この場合では2）で割ります。



そしてその2行目を使って3行目の真ん中を0にします。（1/2をかけている）


再び、3行目の右が1になるように定数で割ります。



ここまでが前進消去です。左下に0が三角形に並んでいます。

今回はこの「前進消去」をプログラム化していきます。



これまでの流れとしては、行列をi行j列とすると、

1. 1行目の1列が1になるように1行目を割る。
   （1行目1列で1行目を割る）
2. 2行目の1列が0になるように1行目を定数倍して2行目から引く。
   （2行目マイナス、2行目1列かける1行目）
3. 3行目1列が0になるように、1行目を定数倍して3行目から引く。
   （3行目マイナス、3行目1列かける1行目）
4. 2行目2列が1になるように2行目を割る。
   （2行目2列で2行目を割る）
5. 3行目2列が0になるように、2行目を定数倍して3行目から引く。
   （3行目マイナス、3行目2列かける2行目）
6. 3行目3列が1になるように3行目を割る。

という感じです。


1~3をプログラム化すると、

int i,j,k;
double a_matrix[3][3]={{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3]={9,17,8};

/*1*/
a_matrix[0][0] = a_matrix[0][0] / a_matrix[0][0];//ここでa_matrix[0][0]が変わる
a_matrix[0][1] = a_matrix[0][1] / a_matrix[0][0];
a_matrix[0][2] = a_matrix[0][2] / a_matrix[0][0];
b_matrix[0] = b_matrix[0] /a_matrix[0][0];

/*2*/
a_matrix[1][0] = a_matrix[1][0] - a_matrix[0][0] * a_matrix[1][0];//ここでa_matrix[1][0]が変わる
a_matrix[1][1] = a_matrix[1][1] - a_matrix[0][1] * a_matrix[1][0];
a_matrix[1][2] = a_matrix[1][2] - a_matrix[0][2] * a_matrix[1][0];
b_matrix[1] = b_matrix[1] - b_matrix[0] * a_matrix[1][0];

/*3*/
a_matrix[2][0] = a_matrix[2][0] - a_matrix[0][0] * a_matrix[2][0];//ここでa_matrix[1][0]が変わる
a_matrix[2][1] = a_matrix[2][1] - a_matrix[0][1] * a_matrix[2][0];
a_matrix[2][2] = a_matrix[2][2] - a_matrix[0][2] * a_matrix[2][0];
b_matrix[2] = b_matrix[2] - b_matrix[0] * a_matrix[2][0];

と書きたいところですけど、コメントのようにその時点で使いたい数字が変更されてしまうので、
違う数値に格納しておきます。これをbufとでも付けておきます。
修正して4~6も入れたのが以下です。

修正後

int i,j,k;
double a_matrix[3][3]={{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3]={9,17,8};
double buf;

/*1*/
buf = a_matrix[0][0];
a_matrix[0][0] = a_matrix[0][0] / buf;
a_matrix[0][1] = a_matrix[0][1] / buf;
a_matrix[0][2] = a_matrix[0][2] / buf;
b_matrix[0] = b_matrix[0] /buf;

/*2*/
buf = a_matrix[1][0];
a_matrix[1][0] = a_matrix[1][0] - a_matrix[0][0] * buf;
a_matrix[1][1] = a_matrix[1][1] - a_matrix[0][1] * buf;
a_matrix[1][2] = a_matrix[1][2] - a_matrix[0][2] * buf;
b_matrix[1] = b_matrix[1] - b_matrix[0] * buf;

/*3*/
buf = a_matrix[2][0];
a_matrix[2][0] = a_matrix[2][0] - a_matrix[0][0] * buf;
a_matrix[2][1] = a_matrix[2][1] - a_matrix[0][1] * buf;
a_matrix[2][2] = a_matrix[2][2] - a_matrix[0][2] * buf;
b_matrix[2] = b_matrix[2] - b_matrix[0] * buf;

/*4*/
buf = a_matrix[1][1];
a_matrix[1][0] = a_matrix[1][0] / buf;
a_matrix[1][1] = a_matrix[1][1] / buf;
a_matrix[1][2] = a_matrix[1][2] / buf;
b_matrix[1] = b_matrix[1] /buf;

/*5*/
buf = a_matrix[2][1];
a_matrix[2][0] = a_matrix[2][0] - a_matrix[1][0] * buf;
a_matrix[2][1] = a_matrix[1][1] - a_matrix[1][1] * buf;
a_matrix[2][2] = a_matrix[1][2] - a_matrix[1][2] * buf;
b_matrix[2] = b_matrix[2] - b_matrix[1] * buf;

/*6*/
buf = a_matrix[2][2];
a_matrix[2][0] = a_matrix[2][0] / buf;
a_matrix[2][1] = a_matrix[2][1] / buf;
a_matrix[2][2] = a_matrix[2][2] / buf;
b_matrix[2] = b_matrix[2] /buf;

これで一応の答えは出たかなー？
まあご周知の通りこんなプログラム汎用性がないわけです。
行列増やしたらまた打ち直しとか。

そんなわけでここからループ化していきます。
実際、自分がプログラムに書くのってここからなんですけど＾＾；

まず、感じとしては、同じ一連の数字のところに同じ変数を入れます。
その前に、プログラムとして同じ処理を見極めます。

手順1と4,6はbufで割る計算、2,3,5はaにbufをかけてaから引く計算です。
なので順番は1→2,3→4→5→6とします。

手順1のbuf = a[0][0];を手順4,6と比べてみると、[0][0]と[1][1],[2][2]であることがわかります。
なので両方同じように増えるのでa[i][i]と置いてみます。
ループとして書くと、

for(i = 0;i < 3;i++){
/*処理*/
}

そして、aをbufで割る計算のところは、aの左側がiと同じ、右側は0~2と増えていきます。
右側をjと置きます。右側をjと置くと、jは0,1,2なので

for(j = 0;j < 3;j++){
/*処理*/
}

とすればいいように思います。
以下は実行しないでください。

int i,j,k;
double a_matrix[3][3]={{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3]={9,17,8};
double buf;

/*1*/
for(i = 0;i < 3;i++){
buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] /buf;/*何度も割ってはいけないのでjの外に出す*/
}

/*2*/
buf = a_matrix[1][0];
a_matrix[1][0] = a_matrix[1][0] - a_matrix[0][0] * buf;
a_matrix[1][1] = a_matrix[1][1] - a_matrix[0][1] * buf;
a_matrix[1][2] = a_matrix[1][2] - a_matrix[0][2] * buf;
b_matrix[1] = b_matrix[1] - b_matrix[0] * buf;

/*3*/
buf = a_matrix[2][0];
a_matrix[2][0] = a_matrix[2][0] - a_matrix[0][0] * buf;
a_matrix[2][1] = a_matrix[2][1] - a_matrix[0][1] * buf;
a_matrix[2][2] = a_matrix[2][2] - a_matrix[0][2] * buf;
b_matrix[2] = b_matrix[2] - b_matrix[0] * buf;

/*5*/
buf = a_matrix[2][1];
a_matrix[2][0] = a_matrix[2][0] - a_matrix[1][0] * buf;
a_matrix[2][1] = a_matrix[1][1] - a_matrix[1][1] * buf;
a_matrix[2][2] = a_matrix[1][2] - a_matrix[1][2] * buf;
b_matrix[2] = b_matrix[2] - b_matrix[1] * buf;

次に、手順2,3,5を考えます。
2,3,5ともに1と同じくaの右は0~2です。
2,3ともに左側は1~2です。
しかし5では左側が2のみです。

この場合では手順2,3をペアに、5をひと固まりとして考えます。

手順2,3ではiの値は0で、aの左は1~2です。
手順5ではiの値は1で、aの左は2です。

こうすると、aの左はi+1から始まり、2で終われば良いように思います。
これはb_matrix[i] = b_matrix[i] /buf;の下に書かなければいけないため、
さっきのjループの外側です。なのでもう一度、変数jを利用します。

その中にまたループを作ってaの右を0~2しなくてはいけないので、
同じようにkでループを作ります。

for(j = i+1;j < 3;j++){
for(k = 0;k < 3;k++){
/*処理*/
}
}

これらすべてをプログラムとして完成させると、

#include<stdio.h>

/*prototype*/
int PRINT(double a_matrix[3][3],double b_matrix[3]);

int main(void){
int i,j,k;
double a_matrix[3][3] = {{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3] = {9,17,8};
double buf;

PRINT(a_matrix,b_matrix);

for(i = 0;i < 3;i++){
buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] / buf;

for(j = i + 1;j < 3;j++){
buf = a_matrix[j][i];
for(k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}
PRINT(a_matrix,b_matrix);
return 0;
}

int PRINT(double a_matrix[3][3],double b_matrix[3]){
int i,j;
for(i = 0;i < 3;i++){
for(j = 0;j < 3;j++){
printf("%6.4f ",a_matrix[i][j]);
}
printf("%6.4f\n",b_matrix[i]);
}
printf("\n");
return 0;
}

となります。
これで前進消去は完成です。
PRINTは行列を表示するためのサブルーチンです。

後退代入も同じような手段なのでここまで書ければ簡単だと思います。
それとループの打ち切り数3は変更がきくように変数に書き換えておくと便利かと思います。
同じく行列も50x50くらいで作っておくと、項の数が多くなっても融通がききます。

次回は後退代入とピボット（の予定）です。

おつさまですた。