【minjijo.cのソースコード】
/*
* 最小二乗法の計算を行うプログラム
* minjijo.c
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h> /* strcspn, strcmp, strcat */
#include <ctype.h> /* isdigit */
#include <math.h> /* exp */
#define N 22 /* 入力できるデータ点の数 + 2 */
#define BUF_SIZE 256
/*
* 自作のbool型の定義
*/
typedef enum {
false,
true
} bool;
/*
* ユーザー定義関数の関数原型宣言
*/
void initialize(void);
bool menu(void);
void select_func(void);
bool is_valid_integer(char buf[]);
bool input_data(void);
bool is_valid_real_number(char buf[]);
void ffv(int p, double a, double *b);
void calc(void);
void output(void);
void output_result(void);
void foutput_result(void);
void gen_filename(char filename[]);
bool retry(void);
/*
* グローバル変数
*/
int f, g, n;
double xx, yy, p, q, h, s, fx, gx, c[4][4], d[4],
x[N], y[N], x2[N], x3[N], a[N][4], b[4][N], xi;
/*
* メイン関数
*/
int main(void)
{
/* 初期化 */
initialize();
while (1) {
/* メニュー */
if (!menu())
break;
/* 関数の種類の選択 */
select_func();
/* データの入力 */
if (!input_data()) {
printf("\n入力が中断されました。最初からやり直します。\n\n");
continue;
}
/* 計算 */
calc();
/* 出力 */
output();
/* 操作を繰り返すかどうかの選択 */
if (!retry())
break;
}
printf("プログラムを終了します。\n");
return EXIT_SUCCESS;
}
/*
* 初期化をする関数
*/
void initialize(void)
{
/* 乱数系列の設定 */
srand((unsigned int)time(NULL));
}
/*
* メニューの関数
* true: 計算を行う、false: 終了する。
*/
bool menu(void)
{
bool flag = true;
char buf[BUF_SIZE], ch, zz;
/* 標題の表示 */
printf("【最小二乗法の計算を行います】\n\n"
"次のいずれかを選択してください。\n\n");
/* 処理の選択 */
while (1) {
printf("1. 計算を行う。\n2. 終了する。\n\n1、2のどちらですか? ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%c%c", &ch, &zz);
if ((ch != '1' && ch != '2') || zz != '\n') {
printf("無効な入力です。\n1または2を入力してください。\n\n");
continue;
} else if (ch == '2') {
flag = false;
break;
} else
break;
}
return flag;
}
/*
* 関数の種類を選択する関数
* is_valid_integer()関数を使用する。
*/
void select_func(void)
{
char buf[BUF_SIZE], zz;
printf("\n最小二乗法によって\n y = a*f(x) + b*g(x)\n"
"の形の曲線[または直線]を当てはめます。\n\n"
"基本関数f(x)、g(x)を1～4の番号で選択してください。\n");
while(1) {
printf("f(x) = [1:(x)、2:(1/x)、3:(e^x)]--> ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%d%c", &f, &zz);
buf[strcspn(buf, "\n")] = '\0';
if (!is_valid_integer(buf)) {
printf("無効な入力です。\n半角整数値を入力してください。\n");
continue;
} else if ((f < 1 || 3 < f) || zz != '\n') {
printf("無効な入力です。\n1から3の整数を入力してください。\n");
continue;
} else
break;
}
while(1) {
printf("g(x) = [1:(x)、2:(1/x)、3:(e^x)、4:(定数)]--> ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%d%c", &g, &zz);
buf[strcspn(buf, "\n")] = '\0';
if (!is_valid_integer(buf)) {
printf("無効な入力です。\n半角整数値を入力してください。\n");
continue;
} else if ((g < 1 || 4 < g) || zz != '\n') {
printf("無効な入力です。\n1から4の整数を入力してください。\n");
continue;
} else
break;
}
}
/*
* 文字列が正しい形式の整数値かどうかを判定する関数
* true: 正しい形式、false: 無効な形式。
*/
bool is_valid_integer(char buf[])
{
int i = 0;
/* 空文字列ならば偽を返す */
if (buf[0] == '\0')
return false;
/* 最初は符号(省略可) */
if (buf[i] == '+' || buf[i] == '-')
i++;
/* 次に数字(省略不可) */
if (!isdigit((unsigned char)buf[i]))
return false;
i++;
/* 次に数字列(省略可) */
while (isdigit((unsigned char)buf[i]))
i++;
/* 他に余計な文字がなければ真、あれば偽を返す */
return buf[i] == '\0';
}
/*
* データを入力する関数
* true: 次に進む、false: 中断してメニューに戻る。
* ffv()関数、is_valid_integer()関数およびis_valid_real_number関数を使用する。
*/
bool input_data(void)
{
char buf[BUF_SIZE], ch, zz;
int i;
double prev_x = 0.0;
while(1) {
printf("(x, y)のデータの数nは何点ですか？"
"(ただし、2 ≦ n ≦ %dとする)\nn = ? ", N - 2);
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%d%c", &n, &zz);
buf[strcspn(buf, "\n")] = '\0';
if (!is_valid_integer(buf)) {
printf("無効な入力です。\n半角整数値を入力してください。\n");
continue;
} else if ((n < 2 || N-2 < n) || zz != '\n') {
printf("無効な入力です。\n2から%dの整数を入力してください。\n",
N - 2);
continue;
}
printf("データ x の値は「小から大」の順に入力する。\n"
"途中で中断したい場合は、qを入力してください。\n");
for (i = 1; i <= n; i++) {
while (1) {
printf("x = ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%lf%c", &x[i], &zz);
buf[strcspn(buf, "\n")] = '\0';
/* 中断判定 */
if (strcmp(buf, "q") == 0 || strcmp(buf, "Q") ==0)
return false;
if (!is_valid_real_number(buf)) {
printf("無効な入力です。\n半角実数値を入力してください\n");
continue;
} else if (i > 1 && x[i] < prev_x) {
printf("データ x の値が「小から大」の順ではありません。\n"
"%d個目のデータから入力をやり直してください。\n"
"入力を中断したい場合は、qを入力してください。\n", i);
prev_x = x[i];
continue;
} else {
prev_x = x[i];
break;
}
}
while (1) {
printf("y = ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%lf%c", &y[i], &zz);
buf[strcspn(buf, "\n")] = '\0';
/* 中断判定 */
if (strcmp(buf, "q") == 0 || strcmp(buf, "Q") ==0)
return false;
if (!is_valid_real_number(buf)) {
printf("無効な入力です。\n半角実数値を入力してください\n");
continue;
} else
break;
}
/* 関数を呼び出す */
xi = x[i];
ffv(f, xi, &fx); ffv(g, xi, &gx);
x2[i] = fx; x3[i] = gx;
a[i][1] = x2[i]; a[i][2] = x3[i]; a[i][3] = y[i];
b[1][i] = a[i][1]; b[2][i] = a[i][2];
}
while (1) {
printf("\n正しく入力しましたか？(y/n) ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%c%c", &ch, &zz);
if ((ch == 'y' || ch == 'Y') && zz == '\n')
goto do_exit;
else if ((ch == 'n' || ch == 'N') && zz == '\n') {
printf("入力をやり直してください。\n\n");
break;
} else {
printf("無効な入力です。\nyまたはnを入力してください。\n");
continue;
}
}
}
do_exit:
return true;
}
/*
* 文字列が正しい形式の実数値かどうかを判定する関数
* true: 正しい形式、false: 無効な形式。
*/
bool is_valid_real_number(char buf[])
{
int i = 0;
/* 空文字列ならば偽を返す */
if (buf[0] == '\0')
return false;
/* 最初は符号(省略可) */
if (buf[i] == '+' || buf[i] == '-')
i++;
/* 次に数字(省略不可) */
if (!isdigit((unsigned char)buf[i]))
return false;
i++;
/* 次に数字列(省略可) */
while (isdigit((unsigned char)buf[i]))
i++;
/* 次に小数点(省略可) */
if (buf[i] == '.') {
i++;
/* 次に数字(小数点がある場合は省略不可、なければ省略可) */
if (!isdigit((unsigned char)buf[i]))
return false;
i++;
/* 次に数字列(省略可) */
while (isdigit((unsigned char)buf[i]))
i++;
}
/* 次にeまたはE(省略可) */
if (buf[i] == 'e' || buf[i] == 'E') {
i++;
/* 次に符号(省略可) */
if (buf[i] == '+' || buf[i] == '-')
i++;
/* 次に数字(eまたはEがある場合は省略不可、なければ省略可) */
if (!isdigit((unsigned char)buf[i]))
return false;
i++;
/* 次に数字列(省略可) */
while (isdigit((unsigned char)buf[i]))
i++;
}
/* 他に余計な文字がなければ真、あれば偽を返す */
return buf[i] == '\0';
}
/*
* 基本関数の関数値を求める関数
*/
void ffv(int p, double a, double *b)
{
switch (p) {
case 1:
*b = a;
break;
case 2:
*b = 1.0 / a;
break;
case 3:
*b = exp(a);
break;
case 4:
*b = 1.0;
break;
default:
*b = a;
break;
}
}
/*
* 計算をする関数
*/
void calc(void)
{
int i, j, l;
/* tA・A を計算して配列 c[2][3] に入れる*/
for (i = 1; i <= 2; i++)
for (j = 1; j <= 3; j++) {
c[i][j] = 0.0;
for (l = 1; l <= n; l++)
c[i][j] += b[i][l] * a[l][j];
}
/* 正規方程式をガウス・ジョルダン法で解く */
for (i = 1; i <= 2; i++) {
p = c[i][i];
for (j = i; j <= 3; j++)
c[i][j] /= p;
for (l = 1; l <= 2; l++) {
if (l != i) {
q = c[l][i];
for (j = i; j <= 3; j++)
c[l][j] -= q * c[i][j];
}
}
}
/* 答を配列d[1]、d[2]に入れる */
for (i = 1; i <= 2; i++)
d[i] = c[i][3];
printf("\n求めた基本関数の係数の出力\n"
" a = d[1] = %.15e\n b = d[2] = %.15e\n", d[1], d[2]);
}
/*
* 出力をする関数
* output_result関数およびfoutput_result関数を使用する。
*/
void output(void)
{
char buf[BUF_SIZE], ch, zz;
while (1) {
while (1) {
printf("数表を出力しますか？(y/n) ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%c%c", &ch, &zz);
if ((ch != 'y' && ch != 'Y' && ch != 'n' && ch != 'N') ||
zz != '\n') {
printf("無効な入力です。\nyまたはnを入力してください。\n\n");
continue;
} else
break;
}
/* 数表を出力する場合 */
if ((ch == 'y' || ch == 'Y') && zz == '\n') {
while(1) {
printf("\nエンターキーを押せば数表を出力します。\n");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%c", &zz);
if (zz != '\n') {
printf("無効な入力です。\n");
printf("エンターキーを押してください。\n");
continue;
} else
break;
}
/* グラフを描くための準備（数表を出力） */
h = (x[n] - x[1]) / 50.0; xx = x[1];
while (1) {
printf("数表を画面に出力しますか、ファイルに出力しますか？\n"
"画面なら1を、ファイルなら2を入力してください: ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%c%c", &ch, &zz);
/* 画面出力の場合 */
if (ch == '1' && zz == '\n') {
output_result();
goto do_exit;
}
/* ファイル出力の場合 */
else if (ch == '2' && zz == '\n') {
foutput_result();
goto do_exit;
}
/* 無効な入力の場合 */
else {
printf("無効な入力です。\n1または2を入力してください: ");
continue;
}
}
}
/* 数表を出力しない場合 */
else if ((ch == 'n' || ch == 'N') && zz == '\n')
break;
/* 無効な入力の場合 */
else
printf("無効な入力です。\nyまたはnを入力してください: ");
}
do_exit:
;
}
/*
* 数表を画面出力する関数
* ffv関数を使用する。
*/
void output_result(void)
{
int i;
for (i = 0; i <= 50; i++) {
ffv(f, xx, &fx);
ffv(g, xx, &gx);
yy = d[1] * fx + d[2] * gx;
printf("%.15e,%.15e\n", xx, yy);
xx += h;
}
}
/*
* 数表をファイル出力する関数
* gen_filename関数およびffv関数を使用する。
*/
void foutput_result(void)
{
char filename[BUF_SIZE];
int i;
FILE *fout;
gen_filename(filename);
fout = fopen(filename, "w");
if(fout == NULL) {
printf("ファイルを開けません。ファイル出力を中断します。\n");
return;
}
for (i = 0; i <= 50; i++) {
ffv(f, xx, &fx);
ffv(g, xx, &gx);
yy = d[1] * fx + d[2] * gx;
fprintf(fout, "%.15e,%.15e\n", xx, yy);
xx += h;
}
fclose(fout);
}
/*
* 現在日付時刻からファイル名を生成する関数
*/
void gen_filename(char filename[])
{
char strnumrand[10];
time_t tt = time(NULL);
struct tm *t = localtime(&tt);
sprintf(filename, "minjijo%04d%02d%02d_%02d%02d%02d_",
1900+(t->tm_year), 1+(t->tm_mon), t->tm_mday,
t->tm_hour, t->tm_min, t->tm_sec);
sprintf(strnumrand, "%04d", rand() % 10000);
strcat(filename, strnumrand);
strcat(filename, ".csv");
}
/*
* 操作を繰り返すかどうかを選択する関数
* true: 繰り返す、false: 終了する。
*/
bool retry(void)
{
bool flag;
char buf[BUF_SIZE], ch, zz;
while (1) {
printf("\n繰り返しますか？(y/n) ");
fgets(buf, sizeof buf, stdin);
sscanf(buf, "%c%c", &ch, &zz);
if ((ch == 'y' || ch == 'Y') && zz == '\n') {
printf("\n+++++\n\n");
flag = true;
break;
} else if ((ch == 'n' || ch == 'N') && zz == '\n') {
flag = false;
break;
}
printf("無効な入力です。\nyまたはnを入力してください。\n");
}
return flag;
}
C:\Users\skonishi\Documents>minjijo
【最小二乗法の計算を行います】
次のいずれかを選択してください。
1. 計算を行う。
2. 終了する。
1、2のどちらですか? 1
最小二乗法によって
y = a*f(x) + b*g(x)
の形の曲線[または直線]を当てはめます。
基本関数f(x)、g(x)を1～4の番号で選択してください。
f(x) = [1:(x)、2:(1/x)、3:(e^x)]--> 1
g(x) = [1:(x)、2:(1/x)、3:(e^x)、4:(定数)]--> 2
(x, y)のデータの数nは何点ですか？(ただし、2 ≦ n ≦ 20とする)
n = ? 7
データ x の値は「小から大」の順に入力する。
途中で中断したい場合は、qを入力してください。
x = 0.2
y = 12.1
x = 0.5
y = 4.9
x = 1.0
y = 2.9
x = 0.5
データ x の値が「小から大」の順ではありません。 ←xの小さい値から順番に入力されるように促される。
4個目のデータから入力をやり直してください。
入力を中断したい場合は、qを入力してください。
x = q ←qを入力すると中断される。
入力が中断されました。最初からやり直します。
【最小二乗法の計算を行います】
次のいずれかを選択してください。
1. 計算を行う。
2. 終了する。
1、2のどちらですか? 1
最小二乗法によって
y = a*f(x) + b*g(x)
の形の曲線[または直線]を当てはめます。
基本関数f(x)、g(x)を1～4の番号で選択してください。
f(x) = [1:(x)、2:(1/x)、3:(e^x)]--> 1
g(x) = [1:(x)、2:(1/x)、3:(e^x)、4:(定数)]--> 2
(x, y)のデータの数nは何点ですか？(ただし、2 ≦ n ≦ 20とする)
n = ? 7
データ x の値は「小から大」の順に入力する。
途中で中断したい場合は、qを入力してください。
x = 0.2
y = 12.1
x = 0.5
y = 4.9
x = 1.0
y = 2.9
x = 2.0
y = 2.1
x = 4.0
y = 2.1
x = 8.0
y = 3.4
x = 10.0
y = 4.3
正しく入力しましたか？(y/n) y
求めた基本関数の係数の出力
a = d[1] = 3.980929968509960e-01
b = d[2] = 2.401049801925565e+00
数表を出力しますか？(y/n) y
エンターキーを押せば数表を出力します。
[エンター]キー
数表を画面に出力しますか、ファイルに出力しますか？
画面なら1を、ファイルなら2を入力してください: 2
繰り返しますか？(y/n) y
+++++
【最小二乗法の計算を行います】
次のいずれかを選択してください。
1. 計算を行う。
2. 終了する。
1、2のどちらですか? 1
最小二乗法によって
y = a*f(x) + b*g(x)
の形の曲線[または直線]を当てはめます。
基本関数f(x)、g(x)を1～4の番号で選択してください。
f(x) = [1:(x)、2:(1/x)、3:(e^x)]--> 1
g(x) = [1:(x)、2:(1/x)、3:(e^x)、4:(定数)]--> 4
(x, y)のデータの数nは何点ですか？(ただし、2 ≦ n ≦ 20とする)
n = ? 13
データ x の値は「小から大」の順に入力する。
途中で中断したい場合は、qを入力してください。
x = 0.898
y = -291.1
x = 0.937
y = -288.4
x = 0.987
y = -289.1
x = 1.046
y = -285.7
x = 1.111
y = -283.5
x = 1.188
y = -280.0
x = 1.297
y = -277.2
x = 1.414
y = -273.9
x = 1.548
y = -268.6
x = 1.692
y = -261.5
x = 1.976
y = -254.0
x = 2.257
y = -243.1
x = 2.591
y = -231.0
正しく入力しましたか？(y/n) y
求めた基本関数の係数の出力
a = d[1] = 3.524999140730006e+01
b = d[2] = -3.226773336336213e+02
数表を出力しますか？(y/n) y
エンターキーを押せば数表を出力します。
[エンター]キー
数表を画面に出力しますか、ファイルに出力しますか？
画面なら1を、ファイルなら2を入力してください: 2
繰り返しますか？(y/n) y
+++++
【最小二乗法の計算を行います】
次のいずれかを選択してください。
1. 計算を行う。
2. 終了する。
1、2のどちらですか? 2
プログラムを終了します。
2.000000000000000e-01,1.208486760899803e+01
3.960000000000000e-01,6.220901902322604e+00
5.920000000000001e-01,4.291498422253299e+00
7.880000000000000e-01,3.360714796652551e+00
9.840000000000000e-01,2.831814771020857e+00
1.180000000000000e+00,2.504537704017705e+00
1.376000000000000e+00,2.292724947624504e+00
1.572000000000000e+00,2.153187561231423e+00
1.768000000000000e+00,2.061888261150641e+00
1.964000000000000e+00,2.004385094860960e+00
2.160000000000000e+00,1.971478003719247e+00
2.356000000000000e+00,1.957028408698759e+00
2.552000000000000e+00,1.956783563827992e+00
2.748000000000000e+00,1.967704024751765e+00
2.944000000000000e+00,1.987559764361657e+00
3.140000000000001e+00,2.014677552126638e+00
3.336000000000001e+00,2.047777386753186e+00
3.532000000000001e+00,2.085863389545205e+00
3.728000000000001e+00,2.128149115523809e+00
3.924000000000001e+00,2.174005248370517e+00
4.120000000000001e+00,2.222922225163377e+00
4.316000000000001e+00,2.274483044919919e+00
4.512000000000000e+00,2.328343164275197e+00
4.708000000000000e+00,2.384215414970064e+00
4.904000000000000e+00,2.441858538189741e+00
5.100000000000000e+00,2.501068362749014e+00
5.295999999999999e+00,2.561670942200058e+00
5.491999999999999e+00,2.623517161488912e+00
5.687999999999999e+00,2.686478458691412e+00
5.883999999999999e+00,2.750443401820268e+00
6.079999999999998e+00,2.815314927749707e+00
6.275999999999998e+00,2.881008098511797e+00
6.471999999999998e+00,2.947448265286760e+00
6.667999999999997e+00,3.014569556200635e+00
6.863999999999997e+00,3.082313623206559e+00
7.059999999999997e+00,3.150628597700830e+00
7.255999999999997e+00,3.219468215405176e+00
7.451999999999996e+00,3.288791079351330e+00
7.647999999999996e+00,3.358560036193294e+00
7.843999999999996e+00,3.428741646025061e+00
8.039999999999996e+00,3.499305729747375e+00
8.235999999999995e+00,3.570224981064992e+00
8.431999999999995e+00,3.641474632598162e+00
8.627999999999995e+00,3.713032167503266e+00
8.823999999999995e+00,3.784877069526600e+00
9.019999999999994e+00,3.856990605645379e+00
9.215999999999994e+00,3.929355636444658e+00
9.411999999999994e+00,4.001956450187068e+00
9.607999999999993e+00,4.074778617192074e+00
9.803999999999993e+00,4.147808861682605e+00
9.999999999999993e+00,4.221034948702513e+00
8.980000000000000e-01,-2.910228413498659e+02
9.318600000000000e-01,-2.898292766408147e+02
9.657200000000000e-01,-2.886357119317635e+02
9.995800000000000e-01,-2.874421472227124e+02
1.033440000000000e+00,-2.862485825136612e+02
1.067300000000000e+00,-2.850550178046100e+02
1.101160000000000e+00,-2.838614530955588e+02
1.135020000000000e+00,-2.826678883865076e+02
1.168880000000000e+00,-2.814743236774564e+02
1.202740000000000e+00,-2.802807589684053e+02
1.236600000000000e+00,-2.790871942593541e+02
1.270460000000000e+00,-2.778936295503029e+02
1.304320000000000e+00,-2.767000648412517e+02
1.338180000000000e+00,-2.755065001322005e+02
1.372040000000000e+00,-2.743129354231494e+02
1.405900000000000e+00,-2.731193707140982e+02
1.439760000000000e+00,-2.719258060050470e+02
1.473620000000000e+00,-2.707322412959958e+02
1.507480000000000e+00,-2.695386765869446e+02
1.541340000000000e+00,-2.683451118778935e+02
1.575200000000000e+00,-2.671515471688423e+02
1.609060000000000e+00,-2.659579824597911e+02
1.642920000000000e+00,-2.647644177507399e+02
1.676780000000000e+00,-2.635708530416887e+02
1.710640000000000e+00,-2.623772883326376e+02
1.744500000000000e+00,-2.611837236235864e+02
1.778360000000000e+00,-2.599901589145352e+02
1.812220000000000e+00,-2.587965942054840e+02
1.846080000000000e+00,-2.576030294964328e+02
1.879940000000000e+00,-2.564094647873817e+02
1.913800000000000e+00,-2.552159000783305e+02
1.947660000000000e+00,-2.540223353692793e+02
1.981520000000000e+00,-2.528287706602281e+02
2.015380000000000e+00,-2.516352059511769e+02
2.049240000000000e+00,-2.504416412421258e+02
2.083100000000000e+00,-2.492480765330746e+02
2.116960000000000e+00,-2.480545118240234e+02
2.150820000000000e+00,-2.468609471149722e+02
2.184679999999999e+00,-2.456673824059211e+02
2.218539999999999e+00,-2.444738176968699e+02
2.252399999999999e+00,-2.432802529878187e+02
2.286259999999999e+00,-2.420866882787676e+02
2.320119999999998e+00,-2.408931235697164e+02
2.353979999999998e+00,-2.396995588606652e+02
2.387839999999998e+00,-2.385059941516141e+02
2.421699999999998e+00,-2.373124294425629e+02
2.455559999999998e+00,-2.361188647335117e+02
2.489419999999997e+00,-2.349253000244605e+02
2.523279999999997e+00,-2.337317353154094e+02
2.557139999999997e+00,-2.325381706063582e+02
2.590999999999997e+00,-2.313446058973070e+02
©2017 KONISHI, Shoichiro.