3元連立方程式をガウスの消去法で解いてみた

タイトルまんまです（ガウスの消去法についてはwiki参照のこと）。逆行列を求めるのに似たような方法を使うのですが、微妙に忘れかけてたので復習を兼ねて解いてみた。本当は精度の問題とかもあって、もっとちゃんとしたアルゴリズムがあるらしいが、そんなのは知らん。あまりにひどいハードコーディングなので、明日はN元まで拡張するつもり。

view sourceprint?

01	#include <iostream>

02	#include <iomanip>

03

04	void PrintArray(const double *pArray, int col, int row)

05

{

06	for(int i=0; i<col; ++i)

07

{

08	for(int j=0; j<row; ++j)

09

{

10	std::cout << std::setw(10) << pArray[i*row+j] << " ";

11

}

12	std::cout << std::endl;

13

}

14	std::cout << "------------------------------" << std::endl;

15

16

}

17

18

// 入力用データ

19	// x,y,z = 2, 3, 4

20	// 2x + 3y + 2z = 21 (1)

21	// 3x + 4y + 3z = 30 (2)

22	// 5x + 2y + 4z = 32 (3)

23	double mat[4*4] =

24

{

25	2, 3, 2, 21,

26	3, 4, 3, 30,

27	5, 2, 4, 32,

28	1, 1, 1, 1,

29

};

30

31	int main(int argc, char **argv)

32

{

33

//計算用配列

34	double tmp[4*4] = {1};

35

36	//まずは(1)式のxの係数を1にする

37	for(int i=0; i<4; ++i)

38

{

39	tmp[i] = mat[i] / mat[0];

40

}

41	PrintArray(tmp, 4, 4);

42

43	//(2),(3)式のxの係数を１に合わせる

44	for(int i=0; i<4; ++i)

45

{

46	tmp[1*4+i] = mat[1*4+i] / mat[1*4];

47	tmp[2*4+i] = mat[2*4+i] / mat[2*4];

48

}

49	PrintArray(tmp, 4, 4);

50

51	//(1)式から(2),(3)式を引いたものをtmp[2,3][0-3]に格納

52	for(int i=0; i<4; ++i)

53

{

54	tmp[1*4+i] = tmp[0*4+i] - tmp[1*4+i];

55	tmp[2*4+i] = tmp[0*4+i] - tmp[2*4+i];

56

}

57	PrintArray(tmp, 4, 4);

58

59	//次に(2)式のyの係数と(3)式のyの係数を合わせる

60	double yk = tmp[1*4+1] / tmp[2*4+1];

61	for(int i=1; i<4; ++i)

62

{

63	tmp[2*4+i] = tmp[2*4+i] * yk;

64

}

65	PrintArray(tmp, 4, 4);

66

67	//(2)式から(3)式を引く

68	for(int i=1; i<4; ++i)

69

{

70	tmp[2*4+i] = tmp[1*4+i] - tmp[2*4+i];

71

}

72	PrintArray(tmp, 4, 4);

73

74

//z成分を求める

75	double z = tmp[3*4+2] = tmp[2*4+3] / tmp[2*4+2];

76	PrintArray(tmp, 4, 4);

77

78

//y成分を求める

79	double y = tmp[3*4+1] = (tmp[1*4+3] - (z*tmp[1*4+2])) / tmp[1*4+1];

80	PrintArray(tmp, 4, 4);

81

82

//xを求める

83	double x = tmp[3*4+0] = (tmp[0*4+3] - (z*tmp[0*4+2] + y*tmp[0*4+1])) / tmp[0*4+0];

84	PrintArray(tmp, 4, 4);

85

86	std::cout <<"(x,y,z) = (" << x << ", " << y << ", " << z << ")" << std::endl;

87

88

return 0;

89

}

どうでもいいがなんで2次元配列使わなかったんだろう俺…。あと所々でPrintArrayを呼び出してるのはデバッグのためです。あしからず。

linuxからこんにちは

知り合いからいただいたPCを修理することに成功した。電源まわりの調子が悪くて、電源とハードディスクを交換しました。1万弱で新しいPCが手に入ったと思えば安いもんだ。


OSにはlinux Debianを選択。CPUはAthronX64でグラフィックボードはnvidia Quadro4XLG980と6年ほど前の構成。グラボがDirectX8世代なので、そのうちこいつは交換しようと思いますが、それさえやってしまえばグラフィックプログラミングには困らないスペックになりそうです。

この投稿もlinuxからのものです。ブラウザさえ入っていれば日常の用途にはまったく困らない。いい時代ですね。

これから開発環境を整えようと思います。eclipseも考えたんですが、debianだとちょっと面倒らしいので止めました。なので、ガチのemacsで整えちゃおうかなと。基本的なキーバインドすら覚えていないですが（笑）

OpenGLでテクスチャ表示

楽勝だと思っていたが意外と苦労した。アラインメントや配列について学びなおすいい機会だったと思う。とりあえず、テクスチャを表示するだけコードを載せてみよう。

view sourceprint?

001	#include <gl/glut.h>

002

003

//頂点座標

004	float vertices[][2] =

005

{

006	{-0.5f,  0.5f},

007	{-0.5f, -0.5f},

008	{0.5f, -0.5f},

009	{0.5f,  0.5f},

010

};

011

012	//uv(st)座標

013	float uvs[][2] =

014

{

015	{0.f, 1.f},

016	{0.f, 0.f},

017	{1.f, 0.f},

018	{1.f, 1.f},

019

};

020

021	//テクスチャハンドル

022	GLuint hTextures[1];

023	const unsigned TEXSIZE = 64;

024

025

//画素バッファ

026	GLubyte texImage[TEXSIZE][TEXSIZE][4];

027

028	//アイドル時にも画面は更新する

029	void Idle()

030

{

031	glutPostRedisplay();

032

}

033

034

//描画関数

035	void Display()

036

{

037	glClear(GL_COLOR_BUFFER_BIT);

038	glEnable(GL_TEXTURE_2D);

039	glBindTexture(GL_TEXTURE_2D, hTextures[0]);

040

041	//反時計回りに板ポリを描画

042	glBegin(GL_QUADS);

043

{

044	glVertex2fv(vertices[0]); glTexCoord2fv(uvs[0]);

045	glVertex2fv(vertices[1]); glTexCoord2fv(uvs[1]);

046	glVertex2fv(vertices[2]); glTexCoord2fv(uvs[2]);

047	glVertex2fv(vertices[3]); glTexCoord2fv(uvs[3]);

048

}

049

050

glEnd();

051	glBindTexture(GL_TEXTURE_2D, 0);

052

053	glDisable(GL_TEXTURE_2D);

054	glutSwapBuffers();

055

}

056

057	//テクスチャ初期化

058	void InitTexture()

059

{

060	//テクスチャを１枚生成する

061	glGenTextures(1, hTextures);

062

063	//今回は計算でテクスチャを作る

064	for(unsigned y=0; y<TEXSIZE; ++y)

065

{

066	int yy = y / 8;

067	for(unsigned x=0; x<TEXSIZE; ++x)

068

{

069	int xx = x / 8;

070	int col = ((yy+xx)&1)*255;

071

072	texImage[y][x][0] = col; //r

073	texImage[y][x][1] = 0;   //g

074	texImage[y][x][2] = ~col;//b

075	texImage[y][x][3] = 0xff;//a

076

}

077

}

078

079	glPixelStorei(GL_UNPACK_ALIGNMENT, 4);

080

081	//hTexures[0]の設定をここから行うという宣言

082	glBindTexture(GL_TEXTURE_2D, hTextures[0]);

083

084	//VRAMにイメージを転送

085	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, TEXSIZE, TEXSIZE, 0, GL_RGBA, GL_UNSIGNED_BYTE, texImage);

086

087	//テクスチャマッピング方式。今回は繰り返しなし

088	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);

089	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);

090

091	//補完フィルタ設定。市松模様なのでnearestでおｋ

092	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);

093	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);

094

095	//ポリゴンカラーとの混色設定

096	glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);

097

098	//hTexures[0]の設定終了

099	glBindTexture(GL_TEXTURE_2D, 0);

100

}

101

102	//エントリポイント

103	int main(int argc, char **argv)

104

{

105

//お決まりの処理

106	glutInit(&argc, argv);

107	glutInitDisplayMode(GLUT_DOUBLE| GLUT_RGBA);

108	glutCreateWindow("DrawTexture");

109	glutDisplayFunc(Display);

110	glutIdleFunc(Idle);

111

112	glClearColor(0,0,0,1);

113	InitTexture();

114

115	//メインループに入る

116	glutMainLoop();

117

118

//終了処理

119	glDeleteTextures(1, hTextures);

120

return 0;

121

}

テクスチャを表示するだけでこの面倒さ。実際にはＢＭＰなどの画像を読み込む処理も要りますから、これからが大変ですな・・・。あと、3次元配列を最初から用意しておくのもどうかと思うので、こいつをクラス化したりもしました。この辺はまたあとで。

VC++でアラインメントを指定する

ビットマップ読み込み関連で失敗したのでメモ。

view sourceprint?

1	struct Hoge

2

{

3

char c;

4

int i;

5

};

上記の構造体サイズは理屈上は5バイトであるが、パディングされてるので実際には8バイトになっている。こいつを確実に5バイトにしたいときは#pragma pack命令を使う。

view sourceprint?

1	#pragma pack(push, 1)

2	struct Hoge

3

{

4

char c;

5

int i;

6

};

7	#pragma pack(pop)

1バイト単位でアラインメントを調整するので、結果としてパディングされないということになる。

サンプルコード

view sourceprint?

01	#include <stdio.h>

02

03

struct A

04

{

05

char c;

06

int i;

07

};

08

09	#pragma pack(push,1)

10

struct B

11

{

12

char c;

13

int i;

14

};

15	#pragma pack(pop)

16

17	int main()

18

{

19	printf("sizeof A: %d\n", sizeof(A));// >> sizeof A: 8

20	printf("sizeof B: %d\n", sizeof(B));// >> sizeof B: 5

21

return 0;

22

}

バイナリでファイルを読み込むときに、使えそうなテクだ。
参考サイト

• Visual C++でアラインメントを扱う3つの機能、#pragma pack, __declspec(align(#)), __alignof演算子についてのメモ

OpenGL+GLSLでポイントライトと距離減衰

距離によって強さが変わるライトの実装に成功した。

ある地点における光の強さは、光源までの距離を d 、一定減衰率を kc 、1時減衰率を k1 、2時減衰率を k2 とすると以下の式で表すことができます。

減衰率kc, k1, k2はそれぞれGL側（ホストプログラム）からシェーダ側に送ることができる。

以下はホスト側のプログラム

view sourceprint?

1

//距離減衰

2	//kc: 一定減衰率 = 0

3	//k1: 1次減衰率 = 0

4	//k2: 2次減衰率 = 1.0 / (lightPos.y * lightPos.y)

5	glLightf(GL_LIGHT0, GL_CONSTANT_ATTENUATION, kc);

6	glLightf(GL_LIGHT0, GL_LINEAR_ATTENUATION, k1);

7	glLightf(GL_LIGHT0, GL_QUADRATIC_ATTENUATION, k2);

んで、フラグメントシェーダ側の記述は以下のとおり

view sourceprint?

01	varying vec3 P;

02	varying vec3 N;

03

04	void main()

05

{

06	//法線やライティング情報を受け取る

07	vec3 L = gl_LightSource[0].position.xyz - P;

08

09	//ライトまでの距離

10	float d = length(L);

11

12	//ライトへの向きベクトルを正規化する

13	L = normalize(L);

14

15

//減衰係数

16	float attenuation = 1.0 / (gl_LightSource[0].constantAttenuation

17	+ gl_LightSource[0].linearAttenuation * d

18	+ gl_LightSource[0].quadraticAttenuation * d * d);

19

20	vec3 N = normalize(N);

21

22	//アンビエント確定

23	vec4 ambient = gl_FrontLightProduct[0].ambient;

24	float dotNL = dot(N, L);

25

26	//ディフューズ確定

27	vec4 diffuse = gl_FrontLightProduct[0].diffuse * max(0.0, dotNL);

28

29	vec3 V = normalize(-P);

30	vec3 H = normalize(L + V);

31	float powNH = pow(max(dot(N, H), 0.0), gl_FrontMaterial.shininess);

32

33	if(dotNL <= 0.0)

34

{

35	powNH = 0.0;

36

}

37

//スペキュラ確定

38	vec4 specular = gl_FrontLightProduct[0].specular * powNH;

39

40

//カラー値統合

41	gl_FragColor = ambient + diffuse + specular;

42

43	//減衰係数を乗じる

44	gl_FragColor *= attenuation;

45

}

まあほとんど教本のまんまなんですが・・・。

gl_LightSource[0].constantAttenuation
gl_LightSource[0].linearAttenuation
gl_LightSource[0].quadraticAttenuation

がそれぞれホスト側で設定した減衰係数です。

画像の式に当たるのがコメントで減衰係数と書いてあるところ。んで、最後にカラー値を統合したgl_FragColorに減衰係数をかけて明るさを確定してるわけですね。このあたりはレイトレースでやったことがかなり活きてるかんじ。

OpenGLで簡易シャドウ

OpenGLでのお話。今まで敬遠してきたシャドウイングにチャレンジ。簡易シャドウというらしい。原理はまだよくわからないが、使用している行列や平面のパラメトリック表現から察するに、ライトの位置から平面にびたーっと張り付くように頂点変形をさせるのが影の変換行列ってやつらしい。

影の変換行列
平面の方程式を ax + by + cz + d = 0 、光源ベクトルをe(e_x, e_y, e_z) とすると、影の変換行列(Ms)は

となります。行列の意味についていろいろ考えたいが、とりあえず後回し。

しかし、よくよく見るとポリゴンが重なってる部分は影が濃く出てますね。自分の影を見ればわかりますが、透明なガラスならともかく、この影は不自然ですね。おそらくこれが「簡易」シャドウといわれる所以でしょう。

上記のマトリクスを最初はHTMLに書き込もうと思ったが、画像化してくれる超便利なツールがオンラインで使用できるのでそちらを使った。LATEX形式だが、知らなくても適当にボタンを押してるだけでそれっぽく仕上がるナイスな出来！

VC++のキーワードハイライト機能を拡張する

VC++2008のお話。

"C:\Program Files\Microsoft Visual Studio 9.0\Common7\IDE"フォルダ下にある"usertype.dat"を開く。なければ作成。

登録したいキーワードを追加する。記述の仕方はいたって簡単で、ハイライトしたいキーワードを１行ずつ追加するだけ。以下の例はGLSLの組み込み型。

次にVCのメニューバーから"ツール"->"オプション"->"テキストエディタ"->"ファイル拡張子"を選択。追加したい拡張子と、編集用のエディタを選択する。

私の場合はバーテクスシェーダファイルは"vert"、ピクセルシェーダファイルは"frag"で登録しました。

最後にVCを再起動して終了。

これだけだとほかのブログにもあるので、さらに小ネタを追加。先ほど追加した"usertype.dat"ですが、新しい型が出てくるたびにフォルダから探して編集するのが面倒なので、ショートカットを作成しちゃいましょう。

最後にショートカットの"プロパティ"->"リンク先"をテキストエディタに変更して、もともとのリンク先はそのまま引数として使う。

画像の1番が追加する使用したいエディタへのパス。2番がもともと記述されていた"usertype.dat"へのパスです。ショートカットのアイコンがエディタのものに変わっていれば成功。このショートカットをデスクトップや、プロジェクトのディレクトリに放り込んでおけばいつでも1クリックで追加したい型を編集可能になります。

OpenSSH と Squidでプロキシサーバー

学校のwebフィルターがあまりにも鬱陶しかったので家に鯖を立てて串を通すことにした。
Linuxはdebian, opensshとsquidを入れて、学校のwinからputtyというソフトを使ってキャッシュをコピーしてくる。
ほとんど担任の入れ知恵だが、これで調べ物が楽にできる。


　　　　　　　 ／￣￣￣＼
　 　 　 　 ／　　　　　　　 ＼
　　　　　/　　 ─　　　─　　ヽ
　　　 　 |　　 （●） 　（●）　　|
　　　　　＼　　 （__人__）　__,／
　　　　　／　　　｀ ⌒´　　　＼
　　　＿/((┃))＿＿＿＿_＿i　|　ｷｭｯｷｭｯ
..　／ /ヽ,,⌒) ￣￣￣￣￣（,,ノ ＼
／　 /＿＿＿＿＿＿＿＿＿ヽ.. 　＼
.￣￣￣￣￣￣￣￣￣￣￣￣￣￣
　　　　 　　　　＿＿＿_
　　　　 　 　／⌒ 　 ー ＼
　 　 　 　／ （●） 　（●） ＼　　＋
　　　 　/ :::::⌒（__人__）⌒:::::ヽ
　 　 　 |　　　　 |r┬-|　　　　|　　+
.　 　 　 ＼_　　　`ー'´ 　 _,／
　 　 　　／　 　 　 　 　 　 ＼　　　　　+
　 　 　 |￣￣￣￣￣￣￣￣ |　　ﾄﾝ
　　　_（,,) 　　 自　由　 　　 (,,）_
　 ／ 　|　　 　　　　　　　 　 　|　 ＼
／　 　 |＿＿＿＿＿＿＿＿_|　　　＼

それにしてもクリエイティブ系の学科なのにゲーム会社のサイトすら見れないってどういうことよ・・・。

PythonでFizzBuzz

4月→新入社員→研修→FizzBuzz!!

ということでやってみた。社会人じゃないけど。
一番苦労したのがprintだと改行が上手くいかないので最後にカンマをつけるってところ。

※下記のプログラムには誤りがあります！

view sourceprint?

1	for i in range(1, 100):

2	if(i % 3 == 0) or (i % 5 == 0):

3	print 'fizz',

4	if(i % 5 == 0):

5	print 'buzz',

6

else:

7

print i,

8	print '\n',

この手のアルゴリズムネタはスクリプト言語でやるのが楽しいね。

（追記）
とか格好つけて嘘っぱちなコードをネット上に晒すとかアホもいいとこだ。改めて実行すると上記のアルゴリズムじゃダメダメですね。即効で直したものを上げておこう。あー恥ずかしい！

view sourceprint?

1	for i in range(1, 100):

2	if(i % 3 == 0) or (i % 5 == 0):

3	if(i % 3 == 0):

4	print 'fizz',

5	if(i % 5 == 0):

6	print 'buzz',

7

else:

8

print i,

9	print '\n',

参考サイト

• とりあ えず: [Python][お勉強] Python入門(16) - printステートメント

VC++のビルド時に作られた中間ファイルを削除する

久々にプログラマらしいお仕事をしたので書き残そう。
詳細はあとで編集。とりあえず体裁整える前のコードだけ上げておく。

view sourceprint?

01	# -*- coding: mbcs -*-

02	from fnmatch import fnmatch

03

import os

04

05	root = os.getcwd()

06	print root

07

08	for dirs in os.walk(root):

09	path = dirs[0] + '\\'

10	for f in dirs[2]:

11	if fnmatch(f, '*.obj') or fnmatch(f, '*.pch'):

12

print f

13	os.remove(path+f)

FBX SDKのメモリリークに対応する

Autodesk FBX SDK 2010.2を使っていると単純なプログラムでもメモリリークを起こすことがわかった。
メモリリークを起こす例

view sourceprint?

1	//マネージャの生成

2	mySdkManager = KFbxSdkManager::Create();

3

4	//マネージャの破棄（多分ここでリークしている）

5	mySdkManager->Destroy();

こいつに対応するにはIOSREF.FreeIOSettings()という命令をマネージャを破棄する直前に追加すればよい。

view sourceprint?

1	//マネージャの生成

2	mySdkManager = KFbxSdkManager::Create();

3

4	//マネージャの破棄

5	IOSREF.FreeIOSettings();//<-追加

6	mySdkManager->Destroy();

参考リンク

• Memory Leaks in 2010.0
• FBX SDKでのメモリリーク