今回作成するツールでは、制御文を組み込んでみたいので、始めからウェブページで紹介されていた 前橋氏の電卓をベースに手を加えていくことにする。
まずは、逆ポーランド記法の例を実行してみました。定義も短く簡単に理解できました。 次に、中間記法のの電卓の例を見てみました。中間記法だと、演算子の優先順序や括弧の ことを考慮しなくてはいけません。例をみてみると、どうやら演算子の優先順序は、 後から定義された演算子が優先度が高くなるように、 bisonが自動的に行ってくれるらしいです。 また、構文規則も逆ポーランド記法のときとほとんど変わりません。 全てのプログラムが50行ほどですんでしまいました。自分で使ってみて、 bisonはとても強力なツールだとやっと理解できました。
//complex.h
/* 複素数 */
struct complex
{
double re;
double im;
};
typedef struct complex complex;
void equals_complex(complex*,complex);
void re_equals_complex(complex*,double);
void im_equals_complex(complex*,double);
void complex_plus_complex(complex*,complex,complex);
void complex_minas_complex(complex*,complex,complex);
void complex_multi_complex(complex*,complex,complex);
void complex_divide_complex(complex*,complex,complex);
//complex.c
#include "complex.h"
void equals_complex(complex* c,complex a)
{
c->re = a.re;
c->im = a.im;
}
void re_equals_complex(complex* c,double a)
{
c->re = a;
c->im = 0;
}
void im_equals_complex(complex* c,double a)
{
c->re = 0;
c->im = a;
}
void complex_plus_complex(complex* c,complex a,complex b)
{
c->re = a.re + b.re;
c->im = a.im + b.im;
}
void complex_minas_complex(complex* c,complex a,complex b)
{
c->re = a.re - b.re;
c->im = a.im - b.im;
}
void complex_multi_complex(complex* c,complex a,complex b)
{
c->re = a.re * b.re - a.im * b.im;
c->im = a.re * b.im + a.im * b.re;
}
void complex_divide_complex(complex* c,complex a,complex b)
{
double alpha;
alpha = b.re * b.re + b.im * b.im;
c->re = (a.re * b.re + a.im * b.im) / alpha;
c->im = (a.im * b.re - a.re * b.im) / alpha;
}
次に、bisonの規則を変更します。 出てくる値としては、doubleのほかに、complexがありますから、YYSTYPE宣言から、union宣言に変更 します。また、構文規則も複素数用に変更します。
%union {
double val; /* 数値を返すため */
complex comp; /* 複素数を扱う */
}
%token <val> NUM /* 単純な倍精度数値 */
%token <comp> COMP /* 倍制度複素数 */
%type <comp> cexp
/* BISON宣言 */
%token NUM
%left '-' '+'
%left '*' '/'
%left NEG /* negation--単項マイナス */
/* 文法規則が続く */
%%
input: /* 空文字列 */
| input line
;
line: '\n'
| cexp '\n' { printf ("re:%.5lf im:%.5lf\n", $1.re, $1.im); }
| error '\n' { yyerrok; }
;
cexp: NUM { re_equals_complex(&$$,$1); }
| 'i' NUM { im_equals_complex(&$$,$2); }
| NUM 'i' { im_equals_complex(&$$,$1); }
| 'i' { im_equals_complex(&$$,1); }
| cexp '+' cexp { complex_plus_complex(&$$,$1,$3); }
| cexp '-' cexp { complex_minas_complex(&$$,$1,$3); }
| cexp '*' cexp { complex_multi_complex(&$$,$1,$3); }
| cexp '/' cexp { complex_divide_complex(&$$,$1,$3); }
| '-' cexp %prec NEG { $$.re = -$2.re; $$.im = -$2.im; }
| '(' cexp ')' { equals_complex(&$$,$2); }
;
実行結果は次のようになります。
(1 + i)/(2 + 2i) re:0.50000 im:0.00000 1 + i re:1.00000 im:1.00000 1 + 2i re:1.00000 im:2.00000 1 + i2 re:1.00000 im:2.00000 1 + parse error 1 + 2i + 3 + 4i re:4.00000 im:6.00000 (1 + i)*(1+ i) (1 + i)*(1 + i) re:0.00000 im:2.00000
//matrix.h
/* Use BLAS */
#include "f2c.h"
#include "blaswrap.h"
#define true 1
#define false 0
typedef int bool;
/* double行列 */
struct matrix
{
double** data;
int row;
int column;
};
typedef struct matrix matrix;
bool alloc_matrix(matrix*,int,int);
void free_matrix(matrix*);
void show_matrix(matrix);
bool zero_matrix(matrix);
bool is_matrix(matrix);
bool copy_matrix(matrix*,matrix);
bool matrix_plus_matrix(matrix*,matrix,matrix);
bool matrix_minus_matrix(matrix*,matrix,matrix);
bool matrix_multi_matrix(matrix*,matrix,matrix);
bool double_multi_matrix(matrix*,double,matrix);
bool minus_matrix(matrix*,matrix);
//matrix.c
#include "matrix.h"
#include "math.h"
#ifndef NULL
#define NULL 0
#endif /* !define NULL */
/* Matrix Construct */
bool alloc_matrix(matrix* m,int row,int col)
{
int i;
m->row = row;
m->column = col;
m->data = NULL;
if ((row <= 0) || (col <= 0)) return;
m->data = (double**)malloc(sizeof(double*)*row);
m->data[0] = (double*)malloc(sizeof(double)*row*col);
for (i = 1;i < row;i++) m->data[i] = m->data[i-1] + col;
// printf("Alloc Matrix (%d,%d)\n",row,col);
return true;
}
/* Matrix Destruct */
void free_matrix(matrix* m)
{
if (m->data != NULL) {
free(m->data[0]);
free(m->data);
// printf("Free Matrix (%d,%d)\n",m->row,m->column);
m->row = -1; m->column = -1;
}
}
/* Matrix Show */
void show_matrix(matrix m)
{
int i,j;
if (!is_matrix(m)) return;
for (i=0;i<m.row;i++) {
for (j=0;j<m.column;j++) {
printf("%.5lf ",m.data[i][j]);
}
printf("\n");
}
printf("\n");
}
/* Matrix initilize zero */
bool zero_matrix(matrix m)
{
int i,j;
if (!is_matrix(m)) return false;
for (i=0;i<m.row;i++) {
for (j=0;j<m.column;j++) {
m.data[i][j] = 0;
}
}
return true;
}
/* Matrix initilize random */
bool random_matrix(matrix m)
{
int i,j;
int max = 10;
if (!is_matrix(m)) return false;
for (i=0;i<m.row;i++) {
for (j=0;j<m.column;j++) {
m.data[i][j] = rand()%(max) - max/2;
}
}
return true;
}
/* initilize rand() */
void init_random()
{
int x;
x = (srand((unsigned)time(NULL)), rand() % 100);
}
/* m is Matrix? */
bool is_matrix(matrix m)
{
if (m.row <= 0 || m.column <= 0 || m.data == NULL) return false;
else return true;
}
/* Copy Matrix x -> y */
bool copy_matrix(matrix *y,matrix x)
{
int inc,N;
int i,j;
if (!is_matrix(x)) return false;
inc = 1;
N = (x.row)*(x.column);
// free_matrix(y);
alloc_matrix(y,x.row,x.column);
dcopy_(&N,x.data[0],&inc,y->data[0],&inc);
return true;
}
/* Matrix Plus C = A + B */
bool matrix_plus_matrix(matrix* c,matrix a,matrix b)
{
int inc,N;
double alpha;
if (!is_matrix(a) | !is_matrix(b) | (a.row != b.row) | (a.column != b.column)) return false;
copy_matrix(c,a);
inc = 1;
alpha = 1;
N = (c->row)*(c->column);
daxpy_(&N,&alpha,b.data[0],&inc,c->data[0],&inc);
return true;
}
/* Matrix Minas C = A - B */
bool matrix_minus_matrix(matrix* c,matrix a,matrix b)
{
int inc,N;
double alpha;
if (!is_matrix(a) | !is_matrix(b) | (a.row != b.row) | (a.column != b.column)) return false;
copy_matrix(c,a);
inc = 1;
alpha = -1;
N = (c->row)*(c->column);
daxpy_(&N,&alpha,b.data[0],&inc,c->data[0],&inc);
return true;
}
/* Matrix Multi C = A * B */
bool matrix_multi_matrix(matrix* c,matrix a,matrix b)
{
char transa[1];
char transb[1];
int M,N,K;
int lda,ldb,ldc;
double alpha,beta;
double* ta;
double* tb;
double* tc;
int i,j;
if (!is_matrix(a) | !is_matrix(b) | (a.column != b.row)) return false;
free_matrix(c);
alloc_matrix(c,a.row,b.column);
transa[0] = 'N'; transb[0] = 'N';
M = a.row; N = b.column; K = a.column;
lda = M; ldb = K; ldc = M;
alpha = 1.0; beta = 0.0;
ta = (double*)malloc(sizeof(double)*M*K);
tb = (double*)malloc(sizeof(double)*K*N);
tc = (double*)malloc(sizeof(double)*M*N);
for (i=0;i<M;i++) {
for (j=0;j<K;j++) {
ta[i+j*M] = a.data[i][j];
}
}
for (i=0;i<K;i++) {
for (j=0;j<N;j++) {
tb[i+j*K] = b.data[i][j];
}
}
dgemm_(transa,transb,&M,&N,&K,&alpha,ta,&lda,tb,&ldb,&beta,tc,&ldc);
for (i=0;i<M;i++) {
for (j=0;j<N;j++) {
c->data[i][j] = tc[i+j*M];
}
}
free(ta);
free(tb);
free(tc);
return true;
}
BLASを呼ぶことにより、高速な計算ができているか確かめてみます。 1000 ×1000行列の積の速度を測ってみました。(Athlon900 + Windows2000 + Cygwin)
| BLAS | 速度[ms] |
| 標準BLAS | 46247 |
| 最適化BLAS | 2523 |
リンクするものを、blas.a(標準BLAS)からlibatlas.a(最適化BLAS)に変えただけで高速化 されていますので、BLASの機能がうまく使えていることが分かります。
これを、先ほどの電卓に組み込もうかと思いましたが、変数が扱えたほうが便利なので、 後から前橋氏の電卓に組み込んでみます。
最初に定義されているものは定義、初期Cコードで、%{から}%までの 間に書きます。ここでは、YY_INPUTが再定義されています。これはlexの 入力ルーチンを置き換えるためにします。このプログラムでは、コマンド行編集機能を 実装するためにここを置き換えているようです。もし、GUIツールから呼び出したいときも、 ここの処理を変更すれば、テキストフィールドから受け取った文字列を渡すことも 可能です。
次に、スタート状態を追加しています。これは読み込みモードのようなもので、 例えば、COMMENTというスタート状態になると、このCOMMENTがついた規則としか マッチしません。また、ここでINITIALというものがありますが、これはいつでも 利用可能でありスタート状態が1つも存在しないときの初期状態です。
次に、切り出してくる字句を設定します。 < INITIAL > "define" return DEFINE から < INITIAL > "%" return MOD までは、字句を見てその種類を返しているだけです。 もし切り出してくる文字の種類を増やしたかったらここに、 < INITIAL > ßin" return SIN の ように定義を加えてあげます。lexの字句解析は自動的にできるだけ長い文字列を字句として 切り出してくれるため便利です。
次は、名前の規則です。このプログラムでは、Cなどと同様に、英字とハイフンから始まって 英数字かハイフンが続くものを名前として受け取ります。読み込まれた字句は、グローバル 変数のyytextに保持されています。これをbisonに渡すために、yylval共用体にコピーします。
次の規則は、整数と小数点の規則です。整数の場合、0以外では、先頭に0が来ることはないため 定義が2つに分かれています。
最後にコメントの規則があり、これでlexの定義部分は終了です。 このコメントをC言語のような、//や/* 〜 */にも対応させてみたいです。
calc.lの最後は、YYINPUTの置き換えです。
この次に、%tokenで終端記号を定義しています。 定義されているものは、表2「終端記号の種類」にあるものです。
ここで定義されているものが，lexの定義の中で使われます。 字句を追加する場合には、実際にはここに追加します。
また、%token < expression > INT_LITERAL というのは、 この字句のデータ型を指定しています。指定できるものは、 unionで定義したものだけです。
次に、%typeで非終端記号を定義します。 上から規則を確認していきます。 bisonでは、繰り返しを再帰で表します。
<解釈単位> := <定義または式> | <解釈単位> <定義または式>
<定義または式> := <関数定義> | <式> ";" | <エラー>
<関数定義> := "define" <名前> "(" <パラメータ列> ")" "{" <式列> "}"
| "define" <名前> "(" ")" "{" <式列> "}"
<パラメータ列> := <名前> | <パラメータ列> "," <名前>
<式列> := <式> | <式列> "," <式>
<式> := <同値式> | <名前> "=" <同値式>
<同値式> := <関係式> | <同値式> "==" <関係式> | <同値式> "!=" <関係式>
<関係式> := <加算式> | <関係式> ">" <加算式> | <関係式> ">=" <加算式>
| <関係式> "<" <加算式> | <関係式> "<=" <加算式>
<加算式> := <掛算式> | <加算式> "+" <掛算式> | <加算式> "-" <掛算式>
<掛算式> := <関連性のない式> | <掛算式> "*" <関連性のない式>
| <掛算式> "/" <関連性のない式> | <掛算式> "%" <関連性のない式>
<関連性のない式> := <後置式> | "-" <関連性のない式>
<後置式> := <一次式>
<一次式> := <名前> "(" 式列 ")" | <名前> "(" ")"
| <if文> | <while文> | "(" <式> ")" | <名前>
| <整数リテラル> | <小数リテラル>
<if文> := "if" <式> "{" <式列> "}"
| "if" <式> "{" <式列> "}" "else" "{" <式列> "}"
<while文> := "while" <式> "{" <式列> "}"
<名前> := [a-zA-Z_][a-zA-Z0-9_]*
<整数リテラル> := "0" | [1-9][0-9]*
<小数リテラル> := [0-9]*"."[0-9]*
このプログラムでは、構文規則により演算子の優先順位と、結合規則が定義されています。 例えば、加算式では、加算式に掛算式を足したものが加算式になると定義されていることから、 掛算をし終わったものを足すので、必然的に掛算のほうが優先順位は高くなります。 この例では、優先度が高いものが下にくるようになっています。また演算子の定義は 全て左再帰となっていることから、左から結合していくことが分かります。 これは、%leftで字句を定義したのと同じ意味となります。 もし、べき乗を定義するとすれば、これは右から結合していくので、右再帰となります。
この構文規則にかかれているCのコードは、解析木を作るものです。 中間記法の電卓では、構文解析をしながら計算をしてましたが、 このプログラムでは、構文解析木を作ってから最後に実行します。 clc_eval_expressionが生成された構文木を実行させる関数です。
CLC.h
電卓のインタプリタ関連。
create.c
解析木の作成。計算できるものはあらかじめ計算。
DBG.h
デバッグ用関数。
error.c
エラー関数。yyerrorも含む。
eval.c
解析木の評価。
interface.c
電卓のインタプリタ。
MEM.h
メモリ確保関数定義。
util.c メモリ確保関数、変数や関数の追加や検索。
まずは、flexの定義にこれらの字句を切り取ってくるように登録します。
<INITIAL>"sin" return SIN;
<INITIAL>"cos" return COS;
<INITIAL>"tan" return TAN;
<INITIAL>"^" return POW;
次に、bisonの構文規則にこれらの規則を追加します。 これらの優先度は、掛算よりも高くなります。 まずは、トークンとして登録しておきます。
%token SIN COS TAN POW
次に、構文規則を追加します。 今回追加する、sin,cos,tanは、関数の名前が前にあるので、前置式という字句ということに しておきます。また、sin sin 1のようにつながらず、sin(sin 1)のようにするために、 これらの関数は後置式をとります。また、べき乗は右から結合するため、右再帰となっています。
<関連性のない式> := <後置式> | <前置式> <前置式> := "sin" <後置式> | "cos" <後置式> | "tan" <後置式> | "-" <後置式> <後置式> := <一次式> | <一次式> "^" <後置式>
これをbisonの定義に直すと次のようになります。
%type <expression> prefix_expression
unary_expression
: postfix_expression
| prefix_expression
;
prefix_expression
: SIN postfix_expression
{
$$ = clc_create_sin_expression($2);
}
| COS postfix_expression
{
$$ = clc_create_cos_expression($2);
}
| TAN postfix_expression
{
$$ = clc_create_tan_expression($2);
}
| SUB unary_expression
{
$$ = clc_create_minus_expression($2);
}
;
postfix_expression
: primary_expression
| primary_expression POW postfix_expression
{
$$ = clc_create_binary_expression(POW_EXPRESSION, $1, $3);
}
;
次に、ここで使われている関数をCのプログラムに追加します。 最初に、これらの関数を新しい式のタイプとして登録します。 また、sinなどの構文木や、評価をするための関数の定義を加えます。 これらは、マイナスの式の関数を参考にして作ります。
//calc.hへの変更
typedef enum {
...
POW_EXPRESSION, /* '^' 追加 */
SIN_EXPRESSION, /* 'sin' 追加 */
COS_EXPRESSION, /* 'cos' 追加 */
TAN_EXPRESSION, /* 'tan' 追加 */
...
} ExpressionType;
Expression *clc_create_sin_expression(Expression *operand); /* sin の構文木追加 */
Expression *clc_create_cos_expression(Expression *operand); /* cos の構文木追加 */
Expression *clc_create_tan_expression(Expression *operand); /* tan の構文木追加 */
Value clc_eval_sin_expression(LocalEnvironment *env, Expression *operand); /* sin の評価追加 */
Value clc_eval_cos_expression(LocalEnvironment *env, Expression *operand); /* cos の評価追加 */
Value clc_eval_tan_expression(LocalEnvironment *env, Expression *operand); /* tan の評価追加 */
次に、解析木を作る関数を作ります。 ここで加える関数は、sin,cos,tanだけです。 べき乗は加える必要はありません。実際に計算するほうにだけ加えればいいだけです。
//create.cへの変更
/* sin の構文木追加 */
Expression *
clc_create_sin_expression(Expression *operand)
{
if (operand->type == INT_EXPRESSION
|| operand->type == DOUBLE_EXPRESSION) {
Value v;
v = clc_eval_sin_expression(NULL, operand);
/* Notice! Overwriting operand expression. */
*operand = convert_value_to_expression(&v);
return operand;
} else {
Expression *exp;
exp = clc_alloc_expression(SIN_EXPRESSION);
exp->u.minus_expression = operand;
return exp;
}
}
/* cos の構文木追加 */
Expression *
clc_create_cos_expression(Expression *operand)
{
if (operand->type == INT_EXPRESSION
|| operand->type == DOUBLE_EXPRESSION) {
Value v;
v = clc_eval_cos_expression(NULL, operand);
/* Notice! Overwriting operand expression. */
*operand = convert_value_to_expression(&v);
return operand;
} else {
Expression *exp;
exp = clc_alloc_expression(COS_EXPRESSION);
exp->u.minus_expression = operand;
return exp;
}
}
/* tan の構文木追加 */
Expression *
clc_create_tan_expression(Expression *operand)
{
if (operand->type == INT_EXPRESSION
|| operand->type == DOUBLE_EXPRESSION) {
Value v;
v = clc_eval_tan_expression(NULL, operand);
/* Notice! Overwriting operand expression. */
*operand = convert_value_to_expression(&v);
return operand;
} else {
Expression *exp;
exp = clc_alloc_expression(TAN_EXPRESSION);
exp->u.minus_expression = operand;
return exp;
}
}
次に、式の評価をする関数を加えます。 intとdoubleの評価に、今回追加した式の計算を追加します。 ここに追加するものは、べき乗だけです。
//eval.cへの変更
static int
eval_binary_int(ExpressionType operator, int left, int right)
{
...
case POW_EXPRESSION: /* '^'追加 */
result = (int)pow(left,right);
break;
...
case MINUS_EXPRESSION: /* FALLTHRU */
case SIN_EXPRESSION: /* FALLTHRU */ /* sin追加 */
case COS_EXPRESSION: /* FALLTHRU */ /* cos追加 */
case TAN_EXPRESSION: /* FALLTHRU */ /* tan追加 */
...
}
static void
eval_binary_double(ExpressionType operator, double left, double right,
Value *result)
{
if (operator == ADD_EXPRESSION || operator == SUB_EXPRESSION
|| operator == MUL_EXPRESSION || operator == DIV_EXPRESSION
|| operator == MOD_EXPRESSION || operator == POW_EXPRESSION) { /* '~'追加 */
result->type = DOUBLE_VALUE;
}
...
case POW_EXPRESSION: /* '~'追加 */
result->u.double_value = pow(left, right);
break;
...
case MINUS_EXPRESSION: /* FALLTHRU */
case SIN_EXPRESSION: /* FALLTHRU */ /* sin追加 */
case COS_EXPRESSION: /* FALLTHRU */ /* cos追加 */
case TAN_EXPRESSION: /* FALLTHRU */ /* tan追加 */
...
}
また、マイナスの評価関数を参考にして、sin,cos,tanの 評価関数を作ります。
/* sinの評価追加 */
Value
clc_eval_sin_expression(LocalEnvironment *env, Expression *operand)
{
Value operand_val;
Value result;
operand_val = eval_expression(env, operand);
if (operand_val.type == INT_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = sin(operand_val.u.int_value);
} else {
DBG_assert(operand_val.type == DOUBLE_VALUE,
("operand_val.type..%d", operand_val.type));
result.type = DOUBLE_VALUE;
result.u.double_value = sin(operand_val.u.double_value);
}
return result;
}
/* cosの評価追加 */
Value
clc_eval_cos_expression(LocalEnvironment *env, Expression *operand)
{
Value operand_val;
Value result;
operand_val = eval_expression(env, operand);
if (operand_val.type == INT_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = cos(operand_val.u.int_value);
} else {
DBG_assert(operand_val.type == DOUBLE_VALUE,
("operand_val.type..%d", operand_val.type));
result.type = DOUBLE_VALUE;
result.u.double_value = cos(operand_val.u.double_value);
}
return result;
}
/* tanの評価追加 */
Value
clc_eval_tan_expression(LocalEnvironment *env, Expression *operand)
{
Value operand_val;
Value result;
operand_val = eval_expression(env, operand);
if (operand_val.type == INT_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = tan(operand_val.u.int_value);
} else {
DBG_assert(operand_val.type == DOUBLE_VALUE,
("operand_val.type..%d", operand_val.type));
result.type = DOUBLE_VALUE;
result.u.double_value = tan(operand_val.u.double_value);
}
return result;
}
static Value
eval_expression(LocalEnvironment *env, Expression *expr)
{
...
case POW_EXPRESSION: /* FALLTHRU */ /* '~'追加 */
...
case MINUS_EXPRESSION:
v = clc_eval_minus_expression(env, expr->u.minus_expression);
break;
case SIN_EXPRESSION: /* sin追加 */
v = clc_eval_sin_expression(env, expr->u.minus_expression);
break;
case COS_EXPRESSION: /* cos追加 */
v = clc_eval_cos_expression(env, expr->u.minus_expression);
break;
case TAN_EXPRESSION: /* tan追加 */
v = clc_eval_tan_expression(env, expr->u.minus_expression);
break;
...
}
これを実行した結果は以下のようになります。
>4^0.5; >>2.000000 >2^3; >>8 >2.0^3.0; >>8.000000 >2^-3; parse error. near token - >>3 >2^(-3); >>0 >2^(-3.0); >>0.125000 >pi = 3.14159; >>3.141590 >sin pi; >>0.000003 >sin(pi/2); >>1.000000 >cos pi; >>-1.000000 >cos(pi/2); >>0.000001 >tan pi; >>-0.000003 >tan(pi/4); >>0.999999
どうやら思い通りに動作しているようです。
<一次式> := <名前> "(" 式列 ")" | <名前> "(" ")"
| <if文> | <while文> | <return文> | "(" <式> ")" | <名前>
| <整数リテラル> | <小数リテラル>
<return文> := "return" <式>
まずは、flexとbisonの設定を変更します。
字句として新たに加えるものは、returnだけです。
終端記号としては、return文を加えます。
//flexの設定
<INITIAL>"return" return RETURN;
//bisonの設定
%token RETURN
%type return_expression
primary_expression
...
| return_expression
...
;
return_expression
: RETURN expression
{
$$ = clc_create_return_expression($2);
}
;
次に、ここで使っている関数をCのプログラムに加えます。 式の種類として、return式を追加します。
//calc.hの変更点
typedef enum {
...
RETURN_EXPRESSION, /* return 追加 */
...
} ExpressionType;
typedef struct {
Expression *return_expression;
} ReturnExpression;
struct Expression_tag {
...
ReturnExpression return_expression;
...
}
Expression *clc_create_return_expression(Expression *expression);
//create.cの変更点
Expression *
clc_create_return_expression(Expression *expression)
{
Expression *exp;
exp = clc_alloc_expression(RETURN_EXPRESSION);
exp->u.return_expression.return_expression = expression;
return exp;
}
//eval.cの変更点
static int
eval_binary_int(ExpressionType operator, int left, int right)
{
...
case RETURN_EXPRESSION: /* FALLTHRU */ /* return追加 */
...
}
static void
eval_binary_double(ExpressionType operator, double left, double right,
Value *result)
{
...
case RETURN_EXPRESSION: /* FALLTHRU */ /* return追加 */
...
}
static Value
eval_return_expression(LocalEnvironment *env,Expression *expression)
{
Value result;
result = eval_expression(env, expression);
return result;
}
static Value
eval_expression(LocalEnvironment *env, Expression *expr)
{
...
case RETURN_EXPRESSION:
v = eval_return_expression(env,expr->u.return_expression.return_expression);
break;
...
}
次に、式列中でreturn式がでたら、それより後ろの式の評価をしないで、値を 返すようにします。返す値に、評価した式にreturnが含まれていたかを 設定する変数を加えます。それを見て、残りの式を実行するのかを決定します。
//calc.hの変更点
/* return値の設定 */
typedef enum {
RETURN_VALUE = 1;
NON_RETURN_VALUE
} ReturnType;
typedef struct {
ValueType type;
ReturnType return_type; /* return値の追加 */
union {
int int_value;
double double_value;
} u;
} Value;
//eval.cの変更点
static Value
eval_return_expression(LocalEnvironment *env,Expression *expression)
{
Value result;
result = eval_expression(env, expression);
result.return_type = RETURN_VALUE;
return result;
}
static Value
eval_expression_list_expression(LocalEnvironment *env,
Expression *expression, Expression *next)
{
Value v;
v = eval_expression(env, expression);
if (v.return_type != RETURN_VALUE) {
if (next) {
v = eval_expression(env, next);
}
}
return v;
}
実行結果は次のようになります。
%>define add(a,b) {
%>a + b,return a,0
%>}
%>add(1,2);
1
%>add(2,3);
2
%>define add2(a,b) {
%>a + b,a,0
%>}
%>add2(1,2);
0
%>add2(2,3);
0
思い通りに動作しています。 return式を組み込んだほうは、途中で式の評価が終了していることが 確認できます。このことより、ユーザー定義関数の返り値を自由に設定できます。
//flexの変更点
<INITIAL>"i" {
Expression *expression = clc_alloc_expression(COMPLEX_EXPRESSION);
expression->u.complex_value.r = 0; expression->u.complex_value.i = 1;
yylval.expression = expression;
return COMPLEX_LITERAL;
}
primary_expression
...
| COMPLEX_LITERAL
;
//bisonの変更点
%token <expression> COMPLEX_LITERAL
次に、複素数用のCのプログラムを追加します。 複素数式を追加して、式の共用体の値にも複素数を追加します。
//calc.hの変更点
typedef enum {
...
COMPLEX_EXPRESSION, /* complex追加 */
...
} ExpressionType;
typedef enum {
...
COMPLEX_VALUE, /* complex追加 */
...
} ValueType;
struct Expression_tag {
...
Complex complex_value; /* complex追加 */
...
};
typedef struct {
...
Complex complex_value;
...
};
次に、複素数の解析木を作るための関数を作ります。 intやdoubleの解析木を作っている関数を参考にして作ります。
//create.cの変更点
static Expression
convert_value_to_expression(Value *v)
{
Expression expr;
if (v->type == INT_VALUE) {
expr.type = INT_EXPRESSION;
expr.u.int_value = v->u.int_value;
} else if (v->type == DOUBLE_VALUE) {
expr.type = DOUBLE_EXPRESSION;
expr.u.double_value = v->u.double_value;
} else { /* complex追加 */
DBG_assert(v->type == COMPLEX_VALUE, ("v->type..%d\n", v->type));
expr.type = COMPLEX_EXPRESSION;
expr.u.complex_value = v->u.complex_value;
}
return expr;
}
Expression *
clc_create_binary_expression(ExpressionType operator,
Expression *left,
Expression *right)
{
if ((left->type == INT_EXPRESSION
|| left->type == DOUBLE_EXPRESSION
|| left->type == COMPLEX_EXPRESSION) /* complex追加 */
&& (right->type == INT_EXPRESSION
|| right->type == DOUBLE_EXPRESSION
|| right->type == COMPLEX_EXPRESSION)) { /* complex追加 */
...
}
Expression *
clc_create_sin_expression(Expression *operand)
{
if (operand->type == INT_EXPRESSION
|| operand->type == DOUBLE_EXPRESSION
|| operand->type == COMPLEX_EXPRESSION) { /* complex追加 */
...
}
Expression *
clc_create_cos_expression(Expression *operand)
{
if (operand->type == INT_EXPRESSION
|| operand->type == DOUBLE_EXPRESSION
|| operand->type == COMPLEX_EXPRESSION) { /* complex追加 */
...
}
Expression *
clc_create_tan_expression(Expression *operand)
{
if (operand->type == INT_EXPRESSION
|| operand->type == DOUBLE_EXPRESSION
|| operand->type == COMPLEX_EXPRESSION) { /* complex追加 */
...
}
次に、複素数式がきたときの処理を追加します。
//eval.cの変更点
/* complex追加 */
static Value
eval_complex_expression(Complex complex_value)
{
Value v;
v.type = COMPLEX_VALUE;
v.u.complex_value = complex_value;
return v;
}
static int
eval_binary_int(ExpressionType operator, int left, int right)
{
...
case COMPLEX_EXPRESSION: /* FALLTHRU *//* complex追加 */
...
}
static void
eval_binary_double(ExpressionType operator, double left, double right,
Value *result)
{
...
case COMPLEX_EXPRESSION: /* FALLTHRU *//* complex追加 */
...
}
次に、複素数の解析木を計算する関数を作ります。
ここで使われている複素数の関数は、あらかじめ作っておいた関数で、
掛け算や割り算、sinやcosの複素数用の関数です。
関数の中は小数点の計算をする関数を参考にしました。
/* complex追加 */
static void
eval_binary_complex(ExpressionType operator, Complex left, Complex right,
Value *result)
{
if (operator == ADD_EXPRESSION || operator == SUB_EXPRESSION
|| operator == MUL_EXPRESSION || operator == DIV_EXPRESSION
|| operator == MOD_EXPRESSION || operator == POW_EXPRESSION) { /* '^'追加 */
result->type = COMPLEX_VALUE;
} else {
DBG_assert(operator == EQ_EXPRESSION || operator == NE_EXPRESSION
|| operator == GT_EXPRESSION || operator == GE_EXPRESSION
|| operator == LT_EXPRESSION || operator == LE_EXPRESSION,
("operator..%d\n", operator));
result->type = INT_VALUE;
}
switch (operator) {
case INT_EXPRESSION: /* FALLTHRU */
case DOUBLE_EXPRESSION: /* FALLTHRU */
case COMPLEX_EXPRESSION: /* FALLTHRU *//* complex追加 */
case IDENTIFIER_EXPRESSION: /* FALLTHRU */
case EXPRESSION_LIST_EXPRESSION: /* FALLTHRU */
case ASSIGN_EXPRESSION:
DBG_assert(0, ("bad case...%d", operator));
break;
case ADD_EXPRESSION:
result->u.complex_value.r = left.r + right.r;
result->u.complex_value.i = left.i + right.i;
break;
case SUB_EXPRESSION:
result->u.complex_value.r = left.r - right.r;
result->u.complex_value.i = left.i - right.i;
break;
case MUL_EXPRESSION:
result->u.complex_value = complex_Mul(left, right);
break;
case DIV_EXPRESSION:
result->u.complex_value = complex_Div(left, right);
break;
case MOD_EXPRESSION:
result->u.complex_value = complex_Mod(left, right);
break;
case POW_EXPRESSION: /* '~'追加 */
result->u.complex_value = complex_Pow(left, right);
break;
case EQ_EXPRESSION:
result->u.int_value = (left.r == right.r) && (left.i == right.i);
break;
case NE_EXPRESSION:
result->u.int_value = (left.r != right.r) || (left.i != right.i);
break;
case GT_EXPRESSION:
result->u.int_value = (int)complex_gt(left.r, left.i, right.r, right.i);
break;
case GE_EXPRESSION:
result->u.int_value = (int)complex_ge(left.r, left.i, right.r, right.i);
break;
case LT_EXPRESSION:
result->u.int_value = (int)complex_lt(left.r, left.i, right.r, right.i);
break;
case LE_EXPRESSION:
result->u.int_value = (int)complex_le(left.r, left.i, right.r, right.i);
break;
case MINUS_EXPRESSION: /* FALLTHRU */
case SIN_EXPRESSION: /* FALLTHRU */ /* sin追加 */
case COS_EXPRESSION: /* FALLTHRU */ /* cos追加 */
case TAN_EXPRESSION: /* FALLTHRU */ /* tan追加 */
case IF_EXPRESSION: /* FALLTHRU */
case WHILE_EXPRESSION: /* FALLTHRU */
case RETURN_EXPRESSION: /* FALLTHRU */ /* return追加 */
case FUNCTION_CALL_EXPRESSION: /* FALLTHRU */
case EXPRESSION_TYPE_NUM: /* FALLTHRU */
default:
DBG_assert(0, ("bad default...%d", operator));
}
}
次に、演算子の左右をみて、計算する関数を呼び出す関数に、 複素数を計算する部分を加えます。左右とも複素数の場合は 複素数の計算をして、片方が複素数なら複素数でないほうを複素数に変換してから 計算します。
Value
clc_eval_binary_expression(LocalEnvironment *env,
ExpressionType operator,
Expression *left, Expression *right)
{
Value left_val;
Value right_val;
Value result;
left_val = eval_expression(env, left);
right_val = eval_expression(env, right);
if (left_val.type == INT_VALUE
&& right_val.type == INT_VALUE) {
result.type = INT_VALUE;
result.u.int_value
= eval_binary_int(operator,
left_val.u.int_value, right_val.u.int_value);
} else if (left_val.type == DOUBLE_VALUE
&& right_val.type == DOUBLE_VALUE) {
eval_binary_double(operator,
left_val.u.double_value, right_val.u.double_value,
&result);
} else if (left_val.type == COMPLEX_VALUE
&& right_val.type == COMPLEX_VALUE) {
eval_binary_complex(operator,
left_val.u.complex_value, right_val.u.complex_value,
&result);
} else {
/* cast to complex */
if (left_val.type == COMPLEX_VALUE) {
if (right_val.type == DOUBLE_VALUE) {
right_val.u.complex_value.r = right_val.u.double_value;
right_val.u.complex_value.i = 0;
} else {
right_val.u.complex_value.r = right_val.u.int_value;
right_val.u.complex_value.i = 0;
}
right_val.type = COMPLEX_VALUE;
eval_binary_complex(operator,
left_val.u.complex_value, right_val.u.complex_value,
&result);
} else if (right_val.type == COMPLEX_VALUE) {
if (left_val.type == DOUBLE_VALUE) {
left_val.u.complex_value.r = left_val.u.double_value;
left_val.u.complex_value.i = 0;
} else {
left_val.u.complex_value.r = left_val.u.int_value;
left_val.u.complex_value.i = 0;
}
left_val.type = COMPLEX_VALUE;
eval_binary_complex(operator,
left_val.u.complex_value, right_val.u.complex_value,
&result);
/* cast to double */
} else if (left_val.type == DOUBLE_VALUE) {
right_val.u.double_value = right_val.u.int_value;
right_val.type = DOUBLE_VALUE;
eval_binary_double(operator,
left_val.u.double_value, right_val.u.double_value,
&result);
} else {
left_val.u.double_value = left_val.u.int_value;
left_val.type = DOUBLE_VALUE;
eval_binary_double(operator,
left_val.u.double_value, right_val.u.double_value,
&result);
}
}
return result;
}
最後に、マイナスやsinなどの式の評価や、評価が終わった複素数の式を表示するような プログラムを追加します。
Value
clc_eval_minus_expression(LocalEnvironment *env, Expression *operand)
{
Value operand_val;
Value result;
operand_val = eval_expression(env, operand);
if (operand_val.type == INT_VALUE) {
result.type = INT_VALUE;
result.u.int_value = -operand_val.u.int_value;
} else if (operand_val.type == DOUBLE_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = -operand_val.u.double_value;
} else if (operand_val.type == COMPLEX_VALUE) { /* complex追加 */
result.type = COMPLEX_VALUE;
result.u.complex_value.r = -operand_val.u.complex_value.r;
result.u.complex_value.i = -operand_val.u.complex_value.i;
} else {
DBG_assert(0,
("operand_val.type..%d", operand_val.type));
}
return result;
}
/* sinの評価追加 */
Value
clc_eval_sin_expression(LocalEnvironment *env, Expression *operand)
{
Value operand_val;
Value result;
operand_val = eval_expression(env, operand);
if (operand_val.type == INT_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = sin(operand_val.u.int_value);
} else if (operand_val.type == DOUBLE_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = sin(operand_val.u.double_value);
} else if (operand_val.type == COMPLEX_VALUE) { /* complex追加 */
result.type = COMPLEX_VALUE;
result.u.complex_value = complex_Sin(operand_val.u.complex_value);
} else {
DBG_assert(0,
("operand_val.type..%d", operand_val.type));
}
return result;
}
/* cosの評価追加 */
Value
clc_eval_cos_expression(LocalEnvironment *env, Expression *operand)
{
Value operand_val;
Value result;
operand_val = eval_expression(env, operand);
if (operand_val.type == INT_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = cos(operand_val.u.int_value);
} else if (operand_val.type == DOUBLE_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = cos(operand_val.u.double_value);
} else if (operand_val.type == COMPLEX_VALUE) { /* complex追加 */
result.type = COMPLEX_VALUE;
result.u.complex_value = complex_Cos(operand_val.u.complex_value);
} else {
DBG_assert(0,
("operand_val.type..%d", operand_val.type));
}
return result;
}
/* tanの評価追加 */
Value
clc_eval_tan_expression(LocalEnvironment *env, Expression *operand)
{
Value operand_val;
Value result;
operand_val = eval_expression(env, operand);
if (operand_val.type == INT_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = tan(operand_val.u.int_value);
} else if (operand_val.type == DOUBLE_VALUE) {
result.type = DOUBLE_VALUE;
result.u.double_value = tan(operand_val.u.double_value);
} else if (operand_val.type == COMPLEX_VALUE) { /* complex追加 */
result.type = COMPLEX_VALUE;
result.u.complex_value = complex_Tan(operand_val.u.complex_value);
} else {
DBG_assert(0,
("operand_val.type..%d", operand_val.type));
}
return result;
}
static Value
eval_expression(LocalEnvironment *env, Expression *expr)
{
...
case COMPLEX_EXPRESSION:
v = eval_complex_expression(expr->u.complex_value);
break;
...
}
void
clc_eval_expression(Expression *expression)
{
...
} else if (v.type == COMPLEX_VALUE) {
printf("\t"); complex_print(v.u.complex_value); printf("\n");
} else {
...
}
これを実行した結果は次のようになります。
%>1 + 2;
3
%>1 + i;
1.000000+1.000000 i
%>(1 + i) + (2 + 2*i);
3.000000+3.000000 i
%>(1 + i) - (2 + 2*i);
-1.000000-1.000000 i
%>(1 + i)*(1 + i);
0.000000+2.000000 i
%>(1 + i)^2;
0.000000+2.000000 i
%>(1 + i)^(1 + i);
0.273957+0.583701 i
%>sin(1 + i);
1.298458+0.634964 i
この計算結果は、ちょっと見ただけではあっているか分かりませんが、
前期で履修した複素解析の知識を利用して計算してみたところ、
計算結果は正しかったので、プログラムは正しく動いていると
思われます。
そこで、これらの関数は"{"から"}"までの間に限り複数行に渡って書くことができるとします。 そして、"{"を読み込んだあとは、改行を空白として処理して、"}"が出たら改行として 処理するようにします。
また、最後にセミコロンがきたときだけ実行結果を表示するようにします。
そのため、次のような変更を加えました。
//flexの変更点
static int read_type = 0;
<INITIAL>"{" {
read_type++; return LC;
}
<INITIAL>"}" {
read_type--; if (read_type < 0) read_type = 0; return RC;
}
<INITIAL>[ \t] ;
<INITIAL>"\n" {
if (read_type == 0) {
return EOL;
}
}
static int
tty_input(char *buf, int max_size)
{
int len;
if (st_readline_buffer == NULL) {
st_readline_used_len = 0;
st_readline_buffer = readline("%>");
if (st_readline_buffer == NULL) {
return 0;
}
len = strlen(st_readline_buffer);
st_readline_buffer = realloc(st_readline_buffer,sizeof(char)*(len+2));
st_readline_buffer[len] = '\n';
st_readline_buffer[len+1] = '\0';
}
len = smaller(strlen(st_readline_buffer) - st_readline_used_len,
max_size);
strncpy(buf, &st_readline_buffer[st_readline_used_len], len);
st_readline_used_len += len;
if (st_readline_buffer[st_readline_used_len] == '\0') {
free(st_readline_buffer);
st_readline_buffer = NULL;
}
return len;
}
//bisonの変更点
%token EOL
translation_unit
: definition_or_expression EOL
| translation_unit definition_or_expression EOL
;
definition_or_expression
: function_definition
| expression SEMICOLON
{
clc_eval_expression($1,NOT_SHOW_VALUE);
}
| expression
{
clc_eval_expression($1,SHOW_VALUE);
}
| error
{
yyclearin;
clc_reopen_current_storage();
}
;
"{"がきたら、関数の中ということで、read_typeの値を1増やします。 これが整数なのは、括弧が入れ子になっている場合に対応するためです。 例えば、関数定義の中にif式があるときなどです。 もし、read_typeが0ならば、それは関数の中ではないので、改行を 字句として切り出し、それ以外なら空白として処理します。
また、readlineは改行を読み取ってきませんので、改行を最後に加えるように しています。
bisonの定義では、受理する規則の最後に改行を追加しました。 また、式のあとにセミコロンがあるときは、評価結果を表示せず、 ないときにだけ評価結果を表示するように変更しています。
これにあわせて、Cのプログラムにも変更を加えます。
//calc.hの変更点
typedef enum {
SHOW_VALUE = 1,
NOT_SHOW_VALUE
} ShowType;
void clc_eval_expression(Expression *expr,ShowType show_type); /* 結果を表示するオプション追加 */
//eval.cの変更点
void
clc_eval_expression(Expression *expression,ShowType show_type) /* 結果を表示するオプション追加 */
{
...
if (show_type == SHOW_VALUE && clc_current_interpreter->input_mode == CLC_TTY_INPUT_MODE) {
if (v.type == INT_VALUE) {
...
}
これを実行した結果は次のようになります。
%>a = 1 + i
1.000000+1.000000 i
%>define add(a,b) {
%>a + b
%>}
%>add(a,1)
2.000000+1.000000 i
%>add(2,1)
3
%>add(2,1);
<行列式> := <ランダム行列宣言> | <行列宣言>
<ランダム行列宣言> := "rand" "(" <式> "," <式> ")"
<行列宣言> := "[" <行列要素> "]"
<行列要素> := <式列> | <行列要素> ";" <式列>
制約：式列の数は同じでなくてはならない。
行列の値には変数は使えない。
この定義によると、例えば行列としては、rand(3,4)や、[1,2,3;4,5,6;7,8,9]のようなものがあることになります。 また、行列を宣言するときに式列を受け取りますので、1 + 2やsin(pi)といったものも 使えることになります。ただし、この式は行列の値として代入するときに、値として 計算できる必要があります。また、それぞれの行の要素の数は同じでなくてはいけません。 ここらの構文規則では表しにくいところは、制約として日本語で書いてあります。 実際には、式列の数が同じではない場合には、最も数が多いものに列数をあわせ、 足りない部分には0を入れておくことにします。
flexとbisonの定義をまずは変更します。
//flexの変更点
<INITIAL>"rand" return RAND;
<INITIAL>"[" return LB;
<INITIAL>"]" return RB;
//bisonの変更点
%union {
char *identifier;
ParameterList *parameter_list;
MatrixUnitList *matrix_unit_list;
Expression *expression;
}
%token RAND LB RB
%type <expression> matrix_expression random_matrix matrix_define
%type <matrix_unit_list> matrix_unit
primary_expression
...
| matrix_expression
...
/* <行列式> */
matrix_expression
: random_matrix
| matrix_define
;
/* <ランダム行列宣言> */
random_matrix
: RAND LP expression COMMA expression RP
{
$$ = clc_create_random_matrix($3, $5);
}
;
/* <行列宣言> */
matrix_define
: LB matrix_unit RB
{
$$ = clc_create_matrix($2);
}
;
/* <行列要素> */
matrix_unit
: expression_list
{
$$ = clc_create_matrix_unit($1);
}
| matrix_unit SEMICOLON expresion_list
{
$$ = clc_chain_matrix_unit($1,$3);
}
;
次に、これにあわせてCのプログラムも変更します。 まずは行列式を宣言しておきます。 また、行列の要素となる構造体を宣言します。 この作り方は、式列の場合と似ています。 最後に、bisonで使用している関数を定義します。
//calc.hの変更点
typedef enum {
...
MATRIX_EXPRESSION, /* matrix追加 */
...
} ExpressionType;
typedef struct {
ExpressionList *expression_list;
ExpressionList *next_list;
} MatrixUnitList;
struct Expression_tag {
...
Matrix matrix_value; /* matrix追加 */
...
};
typedef enum {
...
MATRIX_VALUE, /* matrix追加 */
...
} ValueType;
typedef struct {
...
Matrix matrix_value; /* matrix追加 */
...
} Value;
Expression *clc_create_random_matrix(Expression* row,Expression* col); /* ランダムな要素が格納された行列 */
Expression *clc_create_matrix(MatrixUnitList *list); /* 指定された行列作成 */
MatrixUnitList *clc_create_matrix_unit(Expression* unit); /* 行列の要素 */
MatrixUnitList *clc_chain_matrix_unit(MatrixUnitList* list,Expression* list); /* 行列の要素 */
次に、解析木を作る関数を追加します。
//create.cの変更点
static Expression
convert_value_to_expression(Value *v)
{
...
} else if (v->type == MATRIX_VALUE) { /* matrix追加 */
expr.type = MATRIX_EXPRESSION;
expr.u.matrix_value = v->u.matrix_value;
}
...
}
Expression *
clc_create_binary_expression(ExpressionType operator,
Expression *left,
Expression *right)
{
if ((left->type == INT_EXPRESSION
|| left->type == DOUBLE_EXPRESSION
|| left->type == COMPLEX_EXPRESSION /* complex追加 */
|| left->type == MATRIX_EXPRESSION) /* matrix追加 */
&& (right->type == INT_EXPRESSION
|| right->type == DOUBLE_EXPRESSION
|| right->type == COMPLEX_EXPRESSION /* complex追加 */
|| right->type == MATRIX_EXPRESSION)) { /* matrix追加 */
...
}
Expression *
clc_create_minus_expression(Expression *operand)
{
if (operand->type == INT_EXPRESSION
|| operand->type == DOUBLE_EXPRESSION
|| operand->type == COMPLEX_EXPRESSION /* complex追加 */
|| operand->type == MATRIX_EXPRESSION) { /* matrix追加 */
...
}
ここまでは、複素数を追加したときと同じようなことです。 次に行列を宣言する関数を作ります。 最初に、縦横の大きさを決めて、値はランダムな行列を作る関数です。 これは、与えられた式が整数か小数のときにそれらを、配列の大きさとして 利用します。それ以外の場合には作成が失敗します。
Expression *
clc_create_random_matrix(Expression* row,Expression* col) /* ランダムな要素が格納された行列 */
{
Expression *expression;
int i_col = -1,i_row = -1;
if (row->type == INT_EXPRESSION
|| row->type == DOUBLE_EXPRESSION) {
if (row->type == INT_EXPRESSION) {
i_row = row->u.int_value;
} else if (row->type == DOUBLE_EXPRESSION) {
i_row = (int)(row->u.double_value);
} else {
clc_runtime_error(BAD_EXPRESSION_ERR, "can not use this expression to matrix size. ");
i_row = -1;
}
}
if (col->type == INT_EXPRESSION
|| col->type == DOUBLE_EXPRESSION) {
if (col->type == INT_EXPRESSION) {
i_col = col->u.int_value;
} else if (col->type == DOUBLE_EXPRESSION) {
i_col = (int)(col->u.double_value);
} else {
clc_runtime_error(BAD_EXPRESSION_ERR, "can not use this expression to matrix size. ");
i_col = -1;
}
}
if (i_row <= 0 || i_row <= 0)
clc_runtime_error(BAD_EXPRESSION_ERR, "can not use this expression to matrix size. ");
expression = clc_alloc_expression(MATRIX_EXPRESSION);
alloc_matrix(&(expression->u.matrix_value),i_row,i_col);
init_random();
random_matrix((expression->u.matrix_value));
return expression;
}
次に、行列の要素が指定された行列を作成します。 これは、最初に渡された値を一度全部確認し、 縦横の大きさを決定してから、配列に値を代入します。
また、行列に値を代入する場合には、与えられた式を評価して得られた
値を代入します。
Expression
*clc_create_matrix(MatrixUnitList *list) /* 指定された行列作成 */
{
Expression *expression;
ExpressionList unit;
MatrixUnitList* next_list;
Value v;
int i_row = 0,i_col = 0,t = 0;
next_list = list;
//サイズの確認
while (next_list) {
i_row++;
unit = *(next_list->expression_list);
t = 0;
while(true) {
t++;
if (unit.next) {
if (unit.next->type == EXPRESSION_LIST_EXPRESSION) {
unit = unit.next->u.expression_list;
} else {
break;
}
} else {
break;
}
}
i_col = larger(i_col,t);
next_list = next_list->next_list;
}
//値を代入
expression = clc_alloc_expression(MATRIX_EXPRESSION);
alloc_matrix(&(expression->u.matrix_value),i_row,i_col);
zero_matrix((expression->u.matrix_value));
next_list = list;
i_row = 0;
while (next_list) {
unit = *(next_list->expression_list);
i_col = 0;
while(true) {
v = clc_eval_expression_value(NULL, unit.expression);
if (v.type == INT_VALUE) {
(expression->u.matrix_value).data[i_row][i_col] = v.u.int_value;
} else if (v.type == DOUBLE_VALUE) {
(expression->u.matrix_value).data[i_row][i_col] = v.u.double_value;
} else {
clc_runtime_error(BAD_EXPRESSION_ERR,
"can not use this expression to matrix value. ");
}
if (unit.next) {
if (unit.next->type == EXPRESSION_LIST_EXPRESSION) {
unit = unit.next->u.expression_list;
} else {
break;
}
} else {
break;
}
i_col++;
}
next_list = next_list->next_list;
i_row++;
}
return expression;
}
最後に、列式を連結して行列の要素を作成する関数を作ります。
MatrixUnitList *clc_create_matrix_unit(Expression* unit) /* 行列の要素 */
{
MatrixUnitList* list = clc_malloc(sizeof(MatrixUnitList));
if (unit->type == EXPRESSION_LIST_EXPRESSION) {
list->expression_list = &unit->u.expression_list;
list->next_list = NULL;
} else {
DBG_assert(0, ("unit->type..%d\n", unit->type));
}
return list;
}
MatrixUnitList *clc_chain_matrix_unit(MatrixUnitList* list,Expression* unit) /* 行列の要素 */
{
MatrixUnitList* next_list = list;
MatrixUnitList* parent_list = list;
if (unit->type == EXPRESSION_LIST_EXPRESSION) {
while (next_list) {
parent_list = next_list;
next_list = next_list->next_list;
}
parent_list->next_list = clc_create_matrix_unit(unit);
} else {
DBG_assert(0, ("unit->type..%d\n", unit->type));
}
return list;
}
次に、行列を評価する関数を作ります。 ここで問題になるのは、評価し終わって用済みになった 行列をいつ開放するかということです。 なぜなら、計算をする関数では、与えられた行列が新しく作成されたもの なのか、変数からのポインタなのか分からないため、むやみに開放できないのです。 メモリ管理はCプログラムの最大の難点だったりもします。
そこで、行列のデータは自前で開放しないようにします。
今回はガーベッジコレクション(GC)を利用します。
この説明によると、一般的には、mallocをGC_MALLOCに変換するだけでよいよう
なことが書かれていました。そこで行列の確保の関数を次のように変更しました。
bool alloc_matrix(Matrix* m,int row,int col)
{
int i;
m->row = row;
m->column = col;
m->data = NULL;
if ((row <= 0) || (col <= 0)) return false;
m->data = (double**)GC_MALLOC(sizeof(double*)*row);
if (m->data == 0) {
printf("Out of Memory Error\n");
return false;
}
m->data[0] = (double*)GC_MALLOC(sizeof(double)*row*col);
if (m->data[0] == 0) {
printf("Out of Memory Error\n");
return false;
}
for (i = 1;i < row;i++) m->data[i] = m->data[i-1] + col;
return true;
}
また、ちゃんとメモリを開放しているのか確かめてみました。 その結果、ちゃんとメモリを開放しているのかは確認できませんでした。 なぜならば、僕のパソコンの環境では、mallocしてfreeしなくても、 確保したメモリ領域が自動的に開放されているらしいからです。 そのため、GCの効果は確認できませんでした。
しかし、開放されるだろうということを期待して、GCを使いメモリの開放は しないようにします。
では、行列の評価をする関数を作ります。
//eval.cの変更点
/* matrix追加 */
static Value
eval_matrix_expression(Matrix matrix_value)
{
Value v;
v.type = MATRIX_VALUE;
copy_matrix(&(v.u.matrix_value),matrix_value);
return v;
}
static int
eval_binary_int(ExpressionType operator, int left, int right)
{
...
case MATRIX_EXPRESSION: /* FALLTHRU *//* matrix追加 */
}
static void
eval_binary_double(ExpressionType operator, double left, double right,
Value *result)
{
...
case MATRIX_EXPRESSION: /* FALLTHRU *//* matrix追加 */
...
}
static void
eval_binary_complex(ExpressionType operator, Complex left, Complex right,
Value *result)
{
...
case MATRIX_EXPRESSION: /* FALLTHRU *//* matrix追加 */
}
ここまでは、複素数のときと同じです。 次に、行列の計算をする関数を作ります。
/* matrix追加 */
static void
eval_binary_matrix(ExpressionType operator, Matrix left, Matrix right,
Value *result)
{
if (operator == ADD_EXPRESSION || operator == SUB_EXPRESSION
|| operator == MUL_EXPRESSION) {
result->type = MATRIX_VALUE;
} else {
DBG_assert(operator == EQ_EXPRESSION || operator == NE_EXPRESSION,
("operator..%d\n", operator));
result->type = INT_VALUE;
}
switch (operator) {
case INT_EXPRESSION: /* FALLTHRU */
case DOUBLE_EXPRESSION: /* FALLTHRU */
case COMPLEX_EXPRESSION: /* FALLTHRU *//* complex追加 */
case MATRIX_EXPRESSION: /* FALLTHRU *//* matrix追加 */
case IDENTIFIER_EXPRESSION: /* FALLTHRU */
case EXPRESSION_LIST_EXPRESSION: /* FALLTHRU */
case ASSIGN_EXPRESSION:
DBG_assert(0, ("bad case...%d", operator));
break;
case ADD_EXPRESSION:
matrix_plus_matrix(&(result->u.matrix_value),left,right);
break;
case SUB_EXPRESSION:
matrix_minus_matrix(&(result->u.matrix_value),left,right);
break;
case MUL_EXPRESSION:
matrix_multi_matrix(&(result->u.matrix_value),left,right);
break;
case DIV_EXPRESSION:
case MOD_EXPRESSION:
case POW_EXPRESSION: /* '~'追加 */
DBG_assert(0, ("bad case...%d", operator));
break;
case EQ_EXPRESSION:
result->u.int_value = (matrix_equal_matrix(left,right) == true);
break;
case NE_EXPRESSION:
result->u.int_value = (matrix_equal_matrix(left,right) != true);
break;
case GT_EXPRESSION:
case GE_EXPRESSION:
case LT_EXPRESSION:
case LE_EXPRESSION:
case MINUS_EXPRESSION: /* FALLTHRU */
case SIN_EXPRESSION: /* FALLTHRU */ /* sin追加 */
case COS_EXPRESSION: /* FALLTHRU */ /* cos追加 */
case TAN_EXPRESSION: /* FALLTHRU */ /* tan追加 */
case IF_EXPRESSION: /* FALLTHRU */
case WHILE_EXPRESSION: /* FALLTHRU */
case RETURN_EXPRESSION: /* FALLTHRU */ /* return追加 */
case FUNCTION_CALL_EXPRESSION: /* FALLTHRU */
case EXPRESSION_TYPE_NUM: /* FALLTHRU */
default:
DBG_assert(0, ("bad default...%d", operator));
}
}
行列が計算できる演算子の場合にのみ計算をしています。 行列にある演算子は、+,-,*,==,!=です。
次に、左右のオペランドをみて、使う計算を選択する関数に、 行列への分岐を追加します。行列は、行列同士でしか 計算できないこととします。
Value
clc_eval_binary_expression(LocalEnvironment *env,
ExpressionType operator,
Expression *left, Expression *right)
{
...
} else if (left_val.type == MATRIX_VALUE
&& right_val.type == MATRIX_VALUE) {
eval_binary_matrix(operator,
left_val.u.matrix_value, right_val.u.matrix_value,
&result);
}
...
}
Value
clc_eval_minus_expression(LocalEnvironment *env, Expression *operand)
{
...
} else if (operand_val.type == MATRIX_VALUE) { /* matrix追加 */
result.type = MATRIX_VALUE;
minus_matrix(result.u.matrix_value,operand_val.u.matrix_value);
}
...
}
static Value
eval_expression(LocalEnvironment *env, Expression *expr)
{
...
case MATRIX_EXPRESSION: /* matrix追加 */
v = eval_matrix_expression(expr->u.matrix_value);
break;
...
}
void
clc_eval_expression(Expression *expression,ShowType show_type)
{
...
} else if (v.type == MATRIX_VALUE) {
show_matrix(v.u.matrix_value); printf("\n");
}
...
}
実行結果は次のようになりました。
%>a = rand(3,3)
198.00000 -1.00000 413.00000
-180.00000 174.00000 411.00000
-433.00000 473.00000 50.00000
%>x = 1
1
%>pi = 3.14159
3.141590
%>b = [1,2,3;x+1,x*2,x+3;sin(pi/2),cos(pi/4),tan(pi/4)]
1.00000 2.00000 3.00000
2.00000 2.00000 4.00000
1.00000 0.70711 1.00000
%>c = [1,0,0;0,1,0;0,0,1]
1.00000 0.00000 0.00000
0.00000 1.00000 0.00000
0.00000 0.00000 1.00000
%>b + c
2.00000 2.00000 3.00000
2.00000 3.00000 4.00000
1.00000 0.70711 2.00000
%>c = [1,1,1;1,1,1;1,1,1]
1.00000 1.00000 1.00000
1.00000 1.00000 1.00000
1.00000 1.00000 1.00000
%>c = [1,1,1;1,1,1;1,1,1]
1.00000 1.00000 1.00000
1.00000 1.00000 1.00000
1.00000 1.00000 1.00000
%>b * c
6.00000 6.00000 6.00000
8.00000 8.00000 8.00000
2.70711 2.70711 2.70711
その前に、この行列計算をする関数を作っておきます。 作り方としては、lapackのdgesv_を呼ぶようにするだけです。
bool matrix_dgesv(Matrix* x,Matrix a,Matrix b)
{
int N;
int NHRS;
double* A;
int LDA;
int* IPIV;
double* B;
int LDB;
int INFO;
int i,j;
if (!is_matrix(a) || !is_matrix(b)) {
printf("計算できる行列ではありません。");
return false;
}
if ((a.row != a.column) || (a.row != b.row)) {
printf("計算できる行列ではありません。");
return false;
}
N = a.row;
NHRS = b.column;
A = (double*)GC_MALLOC(sizeof(double)*N*N);
for (i = 0; i < N;i++) {
for (j = 0; j < N;j++) {
A[i+j*N] = a.data[i][j];
}
}
LDA = N;
IPIV = (int*)GC_MALLOC(sizeof(int)*N);
B = (double*)GC_MALLOC(sizeof(double)*N);
for (i = 0; i < N;i++) {
for (j = 0; j < NHRS;j++) {
B[i+j*N] = b.data[i][j];
}
}
LDB = N;
dgesv_(&N,&NHRS,A,&LDA,IPIV,B,&LDB,&INFO);
alloc_matrix(x,N,NHRS);
for (i = 0; i < N; i++) {
for (j = 0; j < NHRS;j++) {
x->data[i][j] = B[i+j*N];
}
}
GC_FREE(A); GC_FREE(B); GC_FREE(IPIV);
if (INFO != 0) {
printf("計算失敗... %d\n",INFO);
return false;
}
return true;
}
これをテストデータで実行させてみます。
Matrix m1,m2,m3;
alloc_matrix(&m1, 3, 3);
alloc_matrix(&m2, 3, 3);
alloc_matrix(&m3, 3, 3);
init_random();
zero_matrix(m1);
zero_matrix(m2);
zero_matrix(m3);
m1.data[0][0] = 1; m1.data[0][1] = 2; m1.data[0][2] = 3;
m1.data[1][0] = 2; m1.data[1][1] = 8; m1.data[1][2] = 11;
m1.data[2][0] = 3; m1.data[2][1] = 22; m1.data[2][2] = 35;
m2.data[0][0] = 1; m2.data[0][1] = 1; m2.data[0][2] = -1;
m2.data[1][0] = 0; m2.data[1][1] = 2; m2.data[1][2] = 0;
m2.data[2][0] = 0; m2.data[2][1] = 3; m2.data[2][2] = 1;
show_matrix(m1);
show_matrix(m2);
matrix_dgesv(&m3,m1,m2);
show_matrix(m3);
free_matrix(&m1);
free_matrix(&m2);
free_matrix(&m3);
1.00000 2.00000 3.00000 2.00000 8.00000 11.00000 3.00000 22.00000 35.00000 1.00000 1.00000 -1.00000 0.00000 2.00000 0.00000 0.00000 3.00000 1.00000 1.58333 1.00000 -1.66667 -1.54167 0.00000 1.33333 0.83333 -0.00000 -0.66667
思い通りに計算できています。
では、まず新しい演算子をflexに追加します。
//flexの変更
<INITIAL>"\" return BS;
//bisonの変更
%token BS
postfix_expression
: primary_expression
| primary_expression POW postfix_expression
{
$$ = clc_create_binary_expression(POW_EXPRESSION, $1, $3);
}
| postfix_expression BS primary_expression
{
$$ = clc_create_binary_expression(BS_EXPRESSION, $1, $3);
}
;
また、これに対するCのプログラムを追加します。
まずは定義を追加します。
//calc.cの変更点
typedef enum {
...
BS_EXPRESSION, /* '\' 追加 */
...
} ExpressionType;
次に、評価関数に変更を加えます。この演算子を使えるのは、 行列計算のみです。
static int
eval_binary_int(ExpressionType operator, int left, int right)
{
...
case BS_EXPRESSION: /* '\'追加 */
DBG_assert(0, ("bad case...%d", operator));
break;
...
}
static void
eval_binary_double(ExpressionType operator, double left, double right,
Value *result)
{
...
case BS_EXPRESSION: /* '\'追加 */
DBG_assert(0, ("bad case...%d", operator));
break;
...
}
static void
eval_binary_complex(ExpressionType operator, Complex left, Complex right,
Value *result)
{
...
case BS_EXPRESSION: /* '\'追加 */
DBG_assert(0, ("bad case...%d", operator));
break;
...
}
static void
eval_binary_matrix(ExpressionType operator, Matrix left, Matrix right,
Value *result)
{
static void
eval_binary_matrix(ExpressionType operator, Matrix left, Matrix right,
Value *result)
{
if (operator == ADD_EXPRESSION || operator == SUB_EXPRESSION
|| operator == MUL_EXPRESSION || operator == BS_EXPRESSION) { /* '\'追加 */
result->type = MATRIX_VALUE;
...
case BS_EXPRESSION: /* '\'追加 */
matrix_dgesv(&(result->u.matrix_value),left,right);
break;
...
}
static Value
eval_expression(LocalEnvironment *env, Expression *expr)
{
...
case BS_EXPRESSION: /* FALLTHRU */ /* '\'追加 */
...
}
この実行結果は次のようになりました。
%>A = [1,2,3;2,8,11;3,22,35] 1.00000 2.00000 3.00000 2.00000 8.00000 11.00000 3.00000 22.00000 35.00000 %>b = [1,1,-1;0,2,0;0,3,1] 1.00000 1.00000 -1.00000 0.00000 2.00000 0.00000 0.00000 3.00000 1.00000 %>x = A\b 1.58333 1.00000 -1.66667 -1.54167 0.00000 1.33333 0.83333 -0.00000 -0.66667
struct _calc_complex_matrix
{
double** data_r;
double** data_i;
int row;
int column;
};
typedef struct _calc_complex_matrix ComplexMatrix;
bool alloc_complex_matrix(ComplexMatrix*,int,int);
void free_complex_matrix(ComplexMatrix*);
void show_complex_matrix(ComplexMatrix);
bool zero_complex_matrix(ComplexMatrix);
bool random_complex_matrix(ComplexMatrix m);
bool is_complex_matrix(ComplexMatrix);
bool copy_complex_matrix(ComplexMatrix*,ComplexMatrix);
bool complex_matrix_equal_complex_matrix(ComplexMatrix,ComplexMatrix);
bool complex_matrix_plus_complex_matrix(ComplexMatrix*,ComplexMatrix,ComplexMatrix);
bool complex_matrix_minus_complex_matrix(ComplexMatrix*,ComplexMatrix,ComplexMatrix);
bool complex_matrix_multi_complex_matrix(ComplexMatrix*,ComplexMatrix,ComplexMatrix);
bool double_multi_complex_matrix(ComplexMatrix*,double,ComplexMatrix);
bool minus_complex_matrix(ComplexMatrix*,ComplexMatrix);
bool complex_matrix_dgesv(ComplexMatrix*,ComplexMatrix,ComplexMatrix);
/* ComplexMatrix Construct */
bool alloc_complex_matrix(ComplexMatrix* m,int row,int col)
{
int i;
m->row = row;
m->column = col;
m->data_r = NULL;
m->data_i = NULL;
if ((row <= 0) || (col <= 0)) return false;
m->data_r = (double**)GC_MALLOC(sizeof(double*)*row);
m->data_i = (double**)GC_MALLOC(sizeof(double*)*row);
if ((m->data_r == 0) || (m->data_i == 0)) {
printf("Out of Memory Error\n");
return false;
}
m->data_r[0] = (double*)GC_MALLOC(sizeof(double)*row*col);
m->data_i[0] = (double*)GC_MALLOC(sizeof(double)*row*col);
if ((m->data_r[0] == 0) || (m->data_i[0] == 0)) {
printf("Out of Memory Error\n");
return false;
}
for (i = 1;i < row;i++) m->data_r[i] = m->data_r[i-1] + col;
for (i = 1;i < row;i++) m->data_i[i] = m->data_i[i-1] + col;
return true;
}
/* ComplexMatrix Destruct */
void free_complex_matrix(ComplexMatrix* m)
{
if (m->data_r != NULL) {
GC_FREE(m->data_r[0]);
GC_FREE(m->data_r);
m->row = -1; m->column = -1;
}
if (m->data_i != NULL) {
GC_FREE(m->data_i[0]);
GC_FREE(m->data_i);
m->row = -1; m->column = -1;
}
}
/* ComplexMatrix Show */
void show_complex_matrix(ComplexMatrix m)
{
int i,j;
if (!is_complex_matrix(m)) return;
for (i=0;i<m.row;i++) {
for (j=0;j<m.column;j++) {
printf("%.5f",m.data_r[i][j]);
if (m.data_i[i][j] >= 0) printf("+");
printf("%.5f",m.data_i[i][j]);
printf("i ");
}
printf("\n");
}
printf("\n");
}
/* ComplexMatrix initilize zero */
bool zero_complex_matrix(ComplexMatrix m)
{
int i,j;
if (!is_complex_matrix(m)) return false;
for (i=0;i<m.row;i++) {
for (j=0;j<m.column;j++) {
m.data_r[i][j] = 0;
m.data_i[i][j] = 0;
}
}
return true;
}
/* ComplexMatrix initilize random */
bool random_complex_matrix(ComplexMatrix m)
{
int i,j;
int max = 10;
if (!is_complex_matrix(m)) return false;
for (i=0;i<m.row;i++) {
for (j=0;j<m.column;j++) {
m.data_r[i][j] = rand()%(max) - max/2;
m.data_i[i][j] = rand()%(max) - max/2;
}
}
return true;
}
/* m is ComplexMatrix? */
bool is_complex_matrix(ComplexMatrix m)
{
if (m.row <= 0 || m.column <= 0 || m.data_r == NULL || m.data_i == NULL) return false;
else return true;
}
/* a == b */
bool complex_matrix_equal_complex_matrix(ComplexMatrix a,ComplexMatrix b)
{
int i,j;
if (!is_complex_matrix(a) || !is_complex_matrix(b)) return false;
if ((a.row != b.row) || (a.column != b.column)) return false;
for (i = 0; i < a.row; i++) {
for (j = 0; j < a.column; j++) {
if ((d_abs(a.data_r[i][j] - b.data_r[i][j]) > DBL_MIN)
|| (d_abs(a.data_i[i][j] - b.data_i[i][j]) > DBL_MIN)) return false;
}
}
return true;
}
/* Copy ComplexMatrix x -> y */
bool copy_complex_matrix(ComplexMatrix *y,ComplexMatrix x)
{
int inc,N;
if (!is_complex_matrix(x)) return false;
inc = 1;
N = (x.row)*(x.column);
alloc_complex_matrix(y,x.row,x.column);
dcopy_(&N,x.data_r[0],&inc,y->data_r[0],&inc);
dcopy_(&N,x.data_i[0],&inc,y->data_i[0],&inc);
return true;
}
/* ComplexMatrix Plus C = A + B */
bool complex_matrix_plus_complex_matrix(ComplexMatrix* c,ComplexMatrix a,ComplexMatrix b)
{
int inc,N;
double alpha;
if (!is_complex_matrix(a) | !is_complex_matrix(b) | (a.row != b.row) | (a.column != b.column))
return false;
copy_complex_matrix(c,a);
inc = 1;
alpha = 1;
N = (c->row)*(c->column);
daxpy_(&N,&alpha,b.data_r[0],&inc,c->data_r[0],&inc);
daxpy_(&N,&alpha,b.data_i[0],&inc,c->data_i[0],&inc);
return true;
}
/* ComplexMatrix Minas C = A - B */
bool complex_matrix_minus_complex_matrix(ComplexMatrix* c,ComplexMatrix a,ComplexMatrix b)
{
int inc,N;
double alpha;
if (!is_complex_matrix(a) | !is_complex_matrix(b) | (a.row != b.row) | (a.column != b.column))
return false;
copy_complex_matrix(c,a);
inc = 1;
alpha = -1;
N = (c->row)*(c->column);
daxpy_(&N,&alpha,b.data_r[0],&inc,c->data_r[0],&inc);
daxpy_(&N,&alpha,b.data_i[0],&inc,c->data_i[0],&inc);
return true;
}
/* ComplexMatrix Multi C = A * B */
bool complex_matrix_multi_complex_matrix(ComplexMatrix* c,ComplexMatrix a,ComplexMatrix b)
{
char transa[1];
char transb[1];
int M,N,K;
int lda,ldb,ldc;
double alpha,beta;
double* ta_r,*ta_i;
double* tb_r,*tb_i;
double* tc_r,*tc_i;
int i,j;
if (!is_complex_matrix(a) | !is_complex_matrix(b) | (a.column != b.row))
return false;
alloc_complex_matrix(c,a.row,b.column);
transa[0] = 'N'; transb[0] = 'N';
M = a.row; N = b.column; K = a.column;
lda = M; ldb = K; ldc = M;
alpha = 1.0; beta = 0.0;
ta_r = (double*)malloc(sizeof(double)*M*K);
tb_r = (double*)malloc(sizeof(double)*K*N);
tc_r = (double*)malloc(sizeof(double)*M*N);
ta_i = (double*)malloc(sizeof(double)*M*K);
tb_i = (double*)malloc(sizeof(double)*K*N);
tc_i = (double*)malloc(sizeof(double)*M*N);
for (i=0;i<M;i++) {
for (j=0;j<K;j++) {
ta_r[i+j*M] = a.data_r[i][j];
ta_i[i+j*M] = a.data_i[i][j];
}
}
for (i=0;i<K;i++) {
for (j=0;j<N;j++) {
tb_r[i+j*K] = b.data_r[i][j];
tb_i[i+j*K] = b.data_i[i][j];
}
}
alpha = 1.0; beta = 0.0;
dgemm_(transa,transb,&M,&N,&K,&alpha,ta_r,&lda,tb_r,&ldb,&beta,tc_r,&ldc);
alpha = -1.0; beta = 1.0;
dgemm_(transa,transb,&M,&N,&K,&alpha,ta_i,&lda,tb_i,&ldb,&beta,tc_r,&ldc);
alpha = 1.0; beta = 0.0;
dgemm_(transa,transb,&M,&N,&K,&alpha,ta_r,&lda,tb_i,&ldb,&beta,tc_i,&ldc);
alpha = 1.0; beta = 1.0;
dgemm_(transa,transb,&M,&N,&K,&alpha,ta_i,&lda,tb_r,&ldb,&beta,tc_i,&ldc);
for (i=0;i<M;i++) {
for (j=0;j<N;j++) {
c->data_r[i][j] = tc_r[i+j*M];
c->data_i[i][j] = tc_i[i+j*M];
}
}
free(ta_r);
free(tb_r);
free(tc_r);
free(ta_i);
free(tb_i);
free(tc_i);
return true;
}
bool minus_complex_matrix(ComplexMatrix* a,ComplexMatrix b)
{
int i,j;
if (!copy_complex_matrix(a,b)) return false;
for (i = 0; i < a->row;i++) {
for (j = 0; j < a->column;j++) {
a->data_r[i][j] = -a->data_r[i][j];
a->data_i[i][j] = -a->data_i[i][j];
}
}
return true;
}
では、これを今回のプログラムに組み込みます。
まずは、flexとbisonの定義を追加します。
//flexの変更点
<INITIAL>"crand" return CRAND;
//bisonの変更点
%token CRAND
random_matrix
: RAND LP expression COMMA expression RP
{
$$ = clc_create_random_matrix($3, $5);
}
| CRAND LP expression COMMA expression RP
{
$$ = clc_create_random_complex_matrix($3, $5);
}
;
ランダムな値の複素行列を作る場合の新しい命令を追加しました。 値を指定する場合には、行列のものと共用します。与えられた数に 複素数があったら複素行列となり、整数や小数だけなら行列となるようにします。
//calc.hの変更点
typedef enum {
...
COMPLEX_MATRIX_EXPRESSION, /* complex_matrix追加 */
...
} ExpressionType;
struct Expression_tag {
...
ComplexMatrix complex_matrix_value; /* complex_matrix追加 */
...
};
typedef struct {
...
ComplexMatrix complex_matrix_value; /* complex_matrix追加 */
...
} Value;
typedef enum {
...
COMPLEX_MATRIX_VALUE, /* complex matrix追加 */
...
} ValueType;
Expression *clc_create_random_complex_matrix(Expression* row,Expression* col); /* ランダムな要素が格納された複素行列 */
//create.cの変更点
Expression *
clc_create_random_complex_matrix(Expression* row,Expression* col) /* ランダムな要素が格納された行列 */
{
Expression *expression;
int i_col = -1,i_row = -1;
if (row->type == INT_EXPRESSION
|| row->type == DOUBLE_EXPRESSION) {
if (row->type == INT_EXPRESSION) {
i_row = row->u.int_value;
} else if (row->type == DOUBLE_EXPRESSION) {
i_row = (int)(row->u.double_value);
} else {
clc_runtime_error(BAD_EXPRESSION_ERR, "can not use this expression to matrix size. ");
i_row = -1;
}
}
if (col->type == INT_EXPRESSION
|| col->type == DOUBLE_EXPRESSION) {
if (col->type == INT_EXPRESSION) {
i_col = col->u.int_value;
} else if (col->type == DOUBLE_EXPRESSION) {
i_col = (int)(col->u.double_value);
} else {
clc_runtime_error(BAD_EXPRESSION_ERR, "can not use this expression to matrix size. ");
i_col = -1;
}
}
if (i_row <= 0 || i_row <= 0)
clc_runtime_error(BAD_EXPRESSION_ERR, "can not use this expression to matrix size. ");
expression = clc_alloc_expression(COMPLEX_MATRIX_EXPRESSION);
alloc_complex_matrix(&(expression->u.complex_matrix_value),i_row,i_col);
init_random();
random_complex_matrix((expression->u.complex_matrix_value));
return expression;
}
Expression
*clc_create_matrix(MatrixUnitList *list) /* 指定された行列作成 */
{
Expression *expression;
ExpressionList unit;
MatrixUnitList* next_list;
Expression exp;
Value v;
int i_row = 0,i_col = 0,t = 0;
int create_complex = 0;
next_list = list;
//サイズの確認
while (next_list) {
i_row++;
unit = *(next_list->expression_list);
t = 0;
while(true) {
t++;
exp = *unit.expression;
v = clc_eval_expression_value(NULL, unit.expression);
if (v.type == COMPLEX_VALUE) create_complex = 1;
if (unit.next) {
if (unit.next->type == EXPRESSION_LIST_EXPRESSION) {
unit = unit.next->u.expression_list;
} else {
break;
}
} else {
break;
}
}
i_col = larger(i_col,t);
next_list = next_list->next_list;
}
if (create_complex == 0) { //matrix
//値を代入
expression = clc_alloc_expression(MATRIX_EXPRESSION);
alloc_matrix(&(expression->u.matrix_value),i_row,i_col);
zero_matrix((expression->u.matrix_value));
next_list = list;
i_row = 0;
while (next_list) {
unit = *(next_list->expression_list);
i_col = 0;
while(true) {
v = clc_eval_expression_value(NULL, unit.expression);
if (v.type == INT_VALUE) {
(expression->u.matrix_value).data[i_row][i_col] = v.u.int_value;
} else if (v.type == DOUBLE_VALUE) {
(expression->u.matrix_value).data[i_row][i_col] = v.u.double_value;
} else {
clc_runtime_error(BAD_EXPRESSION_ERR,
"can not use this expression to matrix value. ");
}
if (unit.next) {
if (unit.next->type == EXPRESSION_LIST_EXPRESSION) {
unit = unit.next->u.expression_list;
} else {
break;
}
} else {
break;
}
i_col++;
}
next_list = next_list->next_list;
i_row++;
}
} else { //complex_matrix
//値を代入
expression = clc_alloc_expression(COMPLEX_MATRIX_EXPRESSION);
alloc_complex_matrix(&(expression->u.complex_matrix_value),i_row,i_col);
zero_complex_matrix((expression->u.complex_matrix_value));
next_list = list;
i_row = 0;
while (next_list) {
unit = *(next_list->expression_list);
i_col = 0;
while(true) {
v = clc_eval_expression_value(NULL, unit.expression);
if (v.type == INT_VALUE) {
(expression->u.complex_matrix_value).data_r[i_row][i_col]
= v.u.int_value;
(expression->u.complex_matrix_value).data_i[i_row][i_col] =0;
} else if (v.type == DOUBLE_VALUE) {
(expression->u.complex_matrix_value).data_r[i_row][i_col]
= v.u.double_value;
(expression->u.complex_matrix_value).data_i[i_row][i_col] = 0;
} else if (v.type == COMPLEX_VALUE) {
(expression->u.complex_matrix_value).data_r[i_row][i_col]
= v.u.complex_value.r;
(expression->u.complex_matrix_value).data_i[i_row][i_col]
= v.u.complex_value.i;
} else {
clc_runtime_error(BAD_EXPRESSION_ERR,
"can not use this expression to matrix value. ");
}
if (unit.next) {
if (unit.next->type == EXPRESSION_LIST_EXPRESSION) {
unit = unit.next->u.expression_list;
} else {
break;
}
} else {
break;
}
i_col++;
}
next_list = next_list->next_list;
i_row++;
}
}
return expression;
}
//eval.cの変更点
static Value
eval_complex_matrix_expression(ComplexMatrix complex_matrix_value)
{
Value v;
v.type = COMPLEX_MATRIX_VALUE;
copy_complex_matrix(&(v.u.complex_matrix_value),complex_matrix_value);
return v;
}
static int
eval_binary_int(ExpressionType operator, int left, int right)
{
...
case COMPLEX_MATRIX_EXPRESSION: /* FALLTHRU *//* complex_matrix追加 */
...
}
static void
eval_binary_double(ExpressionType operator, double left, double right,
Value *result)
{
...
case COMPLEX_MATRIX_EXPRESSION: /* FALLTHRU *//* complex_matrix追加 */
...
}
static void
eval_binary_complex(ExpressionType operator, Complex left, Complex right,
Value *result)
{
...
case COMPLEX_MATRIX_EXPRESSION: /* FALLTHRU *//* matrix追加 */
...
}
static void
eval_binary_matrix(ExpressionType operator, Matrix left, Matrix right,
Value *result)
{
...
case COMPLEX_MATRIX_EXPRESSION: /* FALLTHRU *//* complex_matrix追加 */
...
}
static void
eval_binary_complex_matrix(ExpressionType operator, ComplexMatrix left, ComplexMatrix right,
Value *result)
{
if (operator == ADD_EXPRESSION || operator == SUB_EXPRESSION
|| operator == MUL_EXPRESSION) {
result->type = COMPLEX_MATRIX_VALUE;
} else {
DBG_assert(operator == EQ_EXPRESSION || operator == NE_EXPRESSION,
("operator..%d\n", operator));
result->type = INT_VALUE;
}
switch (operator) {
case INT_EXPRESSION: /* FALLTHRU */
case DOUBLE_EXPRESSION: /* FALLTHRU */
case COMPLEX_EXPRESSION: /* FALLTHRU *//* complex追加 */
case MATRIX_EXPRESSION: /* FALLTHRU *//* matrix追加 */
case COMPLEX_MATRIX_EXPRESSION: /* FALLTHRU *//* complex_matrix追加 */
case IDENTIFIER_EXPRESSION: /* FALLTHRU */
case EXPRESSION_LIST_EXPRESSION: /* FALLTHRU */
case ASSIGN_EXPRESSION:
DBG_assert(0, ("bad case...%d", operator));
break;
case ADD_EXPRESSION:
complex_matrix_plus_complex_matrix(&(result->u.complex_matrix_value),left,right);
break;
case SUB_EXPRESSION:
complex_matrix_minus_complex_matrix(&(result->u.complex_matrix_value),left,right);
break;
case MUL_EXPRESSION:
complex_matrix_multi_complex_matrix(&(result->u.complex_matrix_value),left,right);
break;
case DIV_EXPRESSION:
case MOD_EXPRESSION:
case POW_EXPRESSION: /* '~'追加 */
case BS_EXPRESSION: /* '\'追加 */
DBG_assert(0, ("bad case...%d", operator));
break;
case EQ_EXPRESSION:
result->u.int_value = (complex_matrix_equal_complex_matrix(left,right) == true);
break;
case NE_EXPRESSION:
result->u.int_value = (complex_matrix_equal_complex_matrix(left,right) != true);
break;
case GT_EXPRESSION:
case GE_EXPRESSION:
case LT_EXPRESSION:
case LE_EXPRESSION:
case MINUS_EXPRESSION: /* FALLTHRU */
case SIN_EXPRESSION: /* FALLTHRU */ /* sin追加 */
case COS_EXPRESSION: /* FALLTHRU */ /* cos追加 */
case TAN_EXPRESSION: /* FALLTHRU */ /* tan追加 */
case IF_EXPRESSION: /* FALLTHRU */
case WHILE_EXPRESSION: /* FALLTHRU */
case RETURN_EXPRESSION: /* FALLTHRU */ /* return追加 */
case FUNCTION_CALL_EXPRESSION: /* FALLTHRU */
case EXPRESSION_TYPE_NUM: /* FALLTHRU */
default:
DBG_assert(0, ("bad default...%d", operator));
}
}
Value
clc_eval_binary_expression(LocalEnvironment *env,
ExpressionType operator,
Expression *left, Expression *right)
{
...
} else if (left_val.type == COMPLEX_MATRIX_VALUE
&& right_val.type == COMPLEX_MATRIX_VALUE) {
eval_binary_complex_matrix(operator,
left_val.u.complex_matrix_value, right_val.u.complex_matrix_value,
&result);
} else {
...
if (left_val.type == COMPLEX_MATRIX_VALUE || right_val.type == COMPLEX_MATRIX_VALUE) {
if (left_val.type == MATRIX_VALUE) {
cast_to_complex_matrix(&(left_val.u.complex_matrix_value),
left_val.u.matrix_value);
left_val.type = COMPLEX_MATRIX_VALUE;
eval_binary_complex_matrix(operator,
left_val.u.complex_matrix_value, right_val.u.complex_matrix_value,
&result);
} else if (right_val.type == MATRIX_VALUE) {
cast_to_complex_matrix(&(right_val.u.complex_matrix_value),
right_val.u.matrix_value);
right_val.type = COMPLEX_MATRIX_VALUE;
eval_binary_complex_matrix(operator,
left_val.u.complex_matrix_value, right_val.u.complex_matrix_value,
&result);
} else {
DBG_assert(0, ("bad calc left...%d right...%d", left_val.type,right_val.type));
}
}
...
}
Value
clc_eval_minus_expression(LocalEnvironment *env, Expression *operand)
{
...
} else if (operand_val.type == COMPLEX_MATRIX_VALUE) { /* complex matrix追加 */
result.type = COMPLEX_MATRIX_VALUE;
minus_complex_matrix(&(result.u.complex_matrix_value),operand_val.u.complex_matrix_value);
}
...
}
static Value
eval_expression(LocalEnvironment *env, Expression *expr)
{
...
case COMPLEX_MATRIX_EXPRESSION: /* complex matrix追加 */
v = eval_complex_matrix_expression(expr->u.complex_matrix_value);
break;
...
}
void
clc_eval_expression(Expression *expression,ShowType show_type)
{
...
} else if (v.type == COMPLEX_MATRIX_VALUE) {
show_complex_matrix(v.u.complex_matrix_value); printf("\n");
}
...
}
実行結果は次のようになりました。
%>a = [1,2;3,4] 1.00000 2.00000 3.00000 4.00000 %>a = [1,2;3,4+2*i] 1.00000+0.00000i 2.00000+0.00000i 3.00000+0.00000i 4.00000+2.00000i %>b = [i,sin(1+i);cos(i),(i+i)^(1+i)] 0.00000+1.00000i 1.29846+0.63496i 1.54308+-0.00000i -0.26565+0.31982i %>a + b 1.00000+1.00000i 3.29846+0.63496i 4.54308+0.00000i 3.73435+2.31982i %>c = crand(3,3) 4.00000+3.00000i 4.00000-1.00000i 3.00000-3.00000i 1.00000-4.00000i 4.00000-1.00000i 0.00000+2.00000i -3.00000+2.00000i 3.00000+0.00000i -5.00000-4.00000i
calc.h いろいろな関数定義 calc.l 字句解析(flex定義) calc.y 構文解析(bison定義) CLC.h インターフェース定義ファイル complex.c 複素数計算ルーチン complex.h 複素数計算関数定義 create.c 解析木作成 DBG.h デバッグ処理 error.c エラー処理 eval.c 解析木の評価 interface.c インターフェース main.c Makefile matrix.c 行列計算ルーチン matrix.h 行列計算関数定義 MEM.h util.c
プログラムをコンパイルするには、まずdebugフォルダとmemoryフォルダでmakeをします。 次に、トップフォルダでmakeをします。そのときに、BLASとLAPACKとGCのライブラリを リンクしますので、それらを別に用意する必要があります。
必要に応じてmakeファイルの以下の部分を 変更して下さい。
BLAS = \ ./blas/cblaswr.o\ ./blas/libatlas.a LAPACK = \ ./blas/libF77.a\ ./blas/libI77.a\ ./blas/lapack.a\ ./blas/libg2c.a GC = \ ./gc/gc.a
今回のプログラムはまだ個人的には満足はしていません。
なぜなら、 bisonの定義でまだコンフリクトがたくさん発生しているからです。
時間があったらここらへんの修正をしたいと思っています。
また、当初の計画では、GUIベースのアプリケーションにしたかったのですが、
大変そうだったので今回はCUIで我慢しておくことにしました。
また心配なのは、メモリ管理のあたりです。ガーベッジコレクションが
思い通りに動いているか謎です。