テーラー展開

[beginner]

テーラー展開を使ってsin(x)の値を求めるプログラムを作ったのですが値が正しくでてきません。どこが間違っているか教えてもらえると有り難いです。
pow関数を使わない方法でお願いします。
#include<stdio.h>
#include<math.h>
double mysin(double,double);
int main(void) {
double x,y,PI=atan(1.0)*4.0,k=1.0;
printf("角度を度数法で入力してください\n");
scanf("%lf",&y);
x=(2.0*PI/360.0)*y;
printf("sin(%lf)=%lf",x,1+mysin(x,k));
return 0;
}
double mysin(double x,double k) {
double my_power(double);
double kaijou(double);
return (my_power(x)/kaijou(k));
}
double my_power(double x) {
double a=1.0,b=a,c=0.0;
int i,j;
for(j=1;j<=8;++j) {
for(i=20-j;1<=i;--i)
a*=x;
if(j%2==0)
b=-b;
c+=b;
}
return c;
}
double kaijou(double k) {
int i,j;
double a=1.0,b=0.0;
for(j=1;j<=8;++j) {
for(i=20-j;1<=i;--i)
a*=i;
b+=a;
}
return b;
}

[bitter_fox]

beginnerさん、フォーラムルールは読まれましたでしょうか？？

フォーラムルールに従って、どのような値を入力したらどのような値が出力されたのか、またどの値が出力されたら成功なのかを明確にしてください。
また、コードを載せる際はcodeタグを使っていただきますようにお願いします。

コード:

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

double mysin(double,double);

int main(void) {
	double x,y,PI=atan(1.0)*4.0,k=1.0;

	printf("角度を度数法で入力してください\n");
	scanf("%lf",&y);
	x=(2.0*PI/360.0)*y;

	printf("sin(%lf)=%lf",x,1+mysin(x,k));

	return 0;
}

double mysin(double x,double k) {
	double my_power(double);
	double kaijou(double);

	return (my_power(x)/kaijou(k));
}

double my_power(double x) {
	double a=1.0,b=a,c=0.0;
	int i,j;

	for(j=1;j<=8;++j) {
		for(i=20-j;1<=i;--i)
			a*=x;
		if(j%2==0)
			b=-b;
		c+=b;
	}

	return c;
}

double kaijou(double k) {
	int i,j;
	double a=1.0,b=0.0;

	for(j=1;j<=8;++j) {
		for(i=20-j;1<=i;--i)
			a*=i;
		b+=a;
	}

	return b;
}

う～ん、my_power関数内で、a *= xとしてるが、それ以外でaを使って無い上にcを返してるので、無意味では無いでしょうか？そもそも、cって何の値なのでしょう？
あと、kaijouの引数kがひとつも関数内で使われて無いのですが、大丈夫でしょうか？

[hr][追記]これって、どういうアルゴリズムなのでしょうか？
割り算がひとつしかないと言うことは、通分をされていると言うことでしょうか？

コメントをつけてもう一度、コードをあげてください。
特に、my_powerとkaijouについては細かくお願いします。(ここには、いくつかバグがあります・・・)

ちなみに、sinのテイラー展開で求めるのであれば、以下の要領の方が簡単なような気がします。。。
あくまでヒントです。(sinのテイラー展開をそのまま書いたものです)

コード:

int main()
{
	弧度を取得、

以下を複数回行う(nは試行毎に2ずつ増えていく)
	バッファに、(x^n/n!)を代入する(nは一以上の奇数)

	もし、試行回数が2で割り切れるときは
	答えに、バッファを加算、

	割り切れないときは
	答えに、バッファを減算。
}

[hr][修正]
階乗の!が抜けていたので修正しました。

[たいちう]

人によっては数式の方が理解しやすいのではないかと、Wikipediaへのリンクを。
http://ja.wikipedia.org/wiki/%E3%83%86% ... 5%E9%96%8B

ここにsin xの式があるので、そのままプログラムを作ればよいです。

[みけCAT]

自分で
http://assam.cims.hokudai.ac.jp/~josch/ ... aurin2.htm
を参考に作ったコードもあげておきます。

► スポイラーを表示

コード:

double sin(double angle) {
	double result;
	double temp;
	int i,j;
	double fugou=-1;
	result=angle;
	for(i=3;i<=31;i+=2) {
		temp=fugou;
		for(j=0;j<i;j++)temp*=angle;
		for(j=1;j<=i;j++)temp/=j;
		result+=temp;
		fugou=-fugou;
	}
	return result;
}

[beginner]

bitter_foxさん、たいちうさん、みけCATさんコメントありがとうございます。
＞＞bitter_foxさん
　　　フォーラムルールを読んでいませんでした。すみません。分子のｘと分母の階乗を別々に足すということをしていました。
＞＞みけCATさん
　　　できれば作っていただいたコードの説明をお願いします。

[bitter_fox]

フォーラムルールを読んでいませんでした。すみません。

大丈夫ですよ、次回からお願いしますね。

＞＞みけCATさん
　　　できれば作っていただいたコードの説明をお願いします。

みけCATさんではないですが、軽く説明します。。

コード:

double sin(double angle) {           // angle 角度
    double result;                   // result 結果が代入される変数
    double temp;                     // temp 一時的な値を保存しておく変数
    int i,j;                         // i, j ループ要のカウンタ
    double fugou=-1;                 // fugou 符号の情報が入っています。この変数のおかげて非常にスマートなプログラムになってます。
 
    result=angle;                    // sinのテイラー展開の第1項目(x)のところです。
 
    for(i=3;i<=31;i+=2) {            // 第2項目(x^3/3!)以降です。
        temp=fugou;                  // tempに符号情報を与えてあげます。(1)
 
        for(j=0;j<i;j++)temp*=angle; // べき乗をします。分子(x^i)の部分です。(2)
 
        for(j=1;j<=i;j++)temp/=j;    // 階乗で先の値を割ります。分母(i!)の部分です(3)
 
        result+=temp;                // 結果に、tempの値を加算します。(4)
 
        fugou=-fugou;                // fugouの正負を逆転させます。(5)
    }

    return result;
}

みけCATさんのプログラムは、非常にスマートになってます。その鍵となってるのは、fugouです。
fugouには、符号情報が入っており、tempに対して代入してあげることで、互い違いに加算したり減算したりするsinのテイラー展開をうまく処理してます。(tempに対しては乗算と除算のみを行っているのでうまくいっています)

[hr]
まずforの試行回数が奇数回目のときを見ていきましょう。
(1)で、tempに対して、-1を代入してあげます。
temp = -1

(2)では、angleをi乗してtempに代入してあげます。(pow関数にあたります)(iは何回目かによって変動します(3, 7, 11, 15, 19, ...))
temp = -(angle^i)

(3)では、iの階乗をtempに代入してあげます。
temp = -((angle^i) / i!)

(4)で、結果にtempを加算してあげます。
result = result - ((angle^i) / i!)

(5)で、fugouの正負を逆転します。つぎのループに対しての準備になります。
fugou = -(-1) = 1

[hr]
次にforの試行回数が偶数回目のときを見ていきます。
(1)で、tempに対して、1を代入してあげます。
temp = 1

(2)では、angleをi乗してtempに代入してあげます。(pow関数にあたります)(iは何回目かによって変動します(5, 9, 13, 17, 21, ...))
temp = (angle^i)

(3)では、iの階乗をtempに代入してあげます。
temp = ((angle^i) / i!)

(4)で、結果にtempを加算してあげます。
result = result + ((angle^i) / i!)

(5)で、fugouの正負を逆転します。つぎのループに対しての準備になります。
fugou = -1

このように、奇数回目と偶数回目を行ったりきたりしてsin(angle)に収束します。

下の、青で囲まれているのが、奇数回のとき、赤で囲まれているが、偶数回のときです。
[hr][修正]
値の動きの説明が、見にくかったので修正しました。
［修正］
fugouの所を修正しました。
別の変数に代入する直前にまで、fugouを掛けるのを遅らせると、いかなるアルゴリズムに対しても同様のことが出来ます。

[みけCAT]

みけCATさんではないですが、軽く説明します。。

ありがとうございます。
その通りです。

もう少し効率のいいコードも書いてみました。
前のコードはangle^iとi!をいちいち計算しなおしているので、効率が悪いです。

コード:

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

double mysin(double);

int main(void) {
	double result,resultreal;
	double angle;
	/*角度の入力*/
	printf("角度を度数法で入力してください。\n>");
	scanf("%lf",&angle);
	/*弧度法に変換*/
	angle=angle*3.14159265358979323846/180.0;
	/*検証用に本物のsinで計算*/
	resultreal=sin(angle);
	/*自作関数で計算*/
	result=mysin(angle);
	/*結果の表示*/
	printf("本物:%1.15f\n",resultreal);
	printf("自作:%1.15f\n",result);
	return 0;
}

double mysin(double angle) {		/*angleは弧度法*/
	double result;					/*結果を代入する変数*/
	double temp;					/*現在の項の値の変数*/
	double i;						/*ループ変数*/
	double maeresult;				/*前回の結果（終了判定用）*/
	result=temp=angle;				/*一番目の項*/
	for(i=3;;i+=2) {				/*ループ*/
		/*この時点でtempはangle^(i-2)/(i-2)!のはず（符号は除く）*/
		temp=-temp;					/*tempの符号を反転する*/
		temp=temp*angle*angle;		/*tempにangle^2を掛ける*/
		temp=temp/(i-1)/i;			/*tempの分母をi!にする*/
		maeresult=result;			/*前回の結果を保存する*/
		result+=temp;				/*現在の項を加算する*/
		if(result==maeresult)break;	/*結果が変わらなくなったら抜ける*/
	}
	return result;					/*結果を返す*/
}

[beginner]

bitter_foxさん、みけCATさん大変わかりやすい説明です。
ありがとうございます。
また質問すると思うのでその時はよろしくお願いします。