C言語関係掲示板

過去ログ

No.260.3乗根を解くプログラム


No.1564

3乗根を解くプログラムに関して
投稿者---NEPIA(2002/05/21 22:21:01)


C言語を勉強している学生です。
 3乗根をニュートンラプソン法を使って解くプログラムを作成する問題で、いろいろ調べてみたのですがどうやっても3乗根を解いているプログラムの意味が分かりません。どなたかよろしければ教えてください

No.1571

Re:3乗根を解くプログラムに関して
投稿者---かずま(2002/05/22 03:00:34)


> C言語を勉強している学生です。

数学は勉強していますか。

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

double cubic_root(double a)
{
    double x1, x2 = a;

    do {
        x1 = x2;
        x2 = (a - x1*x1*x1) / (3*x1*x1) + x1;
    } while (fabs((x2 - x1) / x1) > 1e-15);
    return x2;
}
		
int main()
{
    double a;

    while (scanf("%lf", &a) == 1)
        printf("%.16g\n", cubic_root(a));
    return 0;
}
^ をべき乗の演算子とします。
a の 3乗根を求めたい場合、y = x^3 - a のグラフを描いて
y = 0 すなわち x軸との交点を求めればいいことは分かりますね。
このグラフ上の任意の点 (x1, x1^3 - a) における接線はどうなるでしょう。
微分して、y' = 3 * x^2 だから、この点における接線の傾きは 3 * x1^2
接線の方程式は、y - (x1^3 - a) = (3 * x1^2) (x - x1)
これが x軸と交わる点 は、y = 0 を代入して、
x = -(x1^3 - a) / (3 * x1^2) + x1
この値を x2 とすると、x2 は x1 より、求める解にずっと近づいていますね。
ちゃんと図を描かないと理解できませんよ。
さて、これを繰り返すと、x2 は、どんどん真の値に近づいて、一つ前の値 x1
との差が 小さくなります。このとき、繰り返しを終了し、x2 を解とします。

No.1572

Re:3乗根を解くプログラムに関して
投稿者---ジャスミン茶(2002/05/22 03:16:25)


はじめまして、ジャスミン茶です。

せっかくなのでニュートン・ラプソン法というものを調べてみました。
私も完全には理解できていないのですが、一応簡単に説明してみます。
ただし、微分の知識が必要になります。

nの3乗根をaとします。
するとaを3乗したものがxになるわけです。
式であらわすと「 a*a*a == n 」です。
さて、これを数学的に表現すると「 a~3 = n 」になります。(以下、べき乗の表現は「~」です。)
xy座標(空間)を使って関数として表現すると、右肩上がりの三次方程式が得られます。

「 y = x~3 - N 」(あえて、nをNとしています。)

この関数においてx座標をaと置くとy座標はa~3 - Nとなります。
つまりこの三次方程式の解として (a, a~3 - N)***(1)*** が得られます。
この時「 y = a~3 - N = 0 」と置く事でNの3乗根であるaが求められるわけです。

これが素直に求められれば何の問題も無いのですが、そう簡単には行きません。
そこで、aを超えない最大の整数を考えます。(便宜上aのままです。理解し辛い場合bにするなどして下さい。)

for (a=1; a*a*a < N; a++) {
} // aの3乗がNを超えるまでa = a + 1を繰り返す。
a--; // 超えてしまっているので1を引いておく。

(あ、言い忘れてしまいましたが、今回の説明ではN >= 1の場合しか適用できません。
Nが負、1以下の少数、0の場合は0で除算が入ってエラーが起きますので場合分けが必要です。多分。)

さて、これでaは求める値に大分近づきました。
問題はここからです。
aは確かに求める真の値に近づきましたが、求める真の値をAとするならば、依然 a<A です。
というわけで A<a となるように範囲を絞っていきます。

ある関数(2次以上)において微分を行うと、その関数と接する直線の傾きが得られます。
今回は
y = x~3 - N より、
y' = 3x~2
となります。

先ほどの解(a, a~3 - N)より「 x = a 」における「 y = x~3 - N 」と接する接線は
y' = 3x~2にx=aを代入する事で
「 y = (3a~2)x - b 」***(2)***
となります。(bは任意です。後で代入で消してしまいます。)
***(2)***が x=a の時、y座標は「 y = 3a~3 - b 」***(3)***です。

この直線と元の三次関数は「 x=a 」の時接するのですから、
お互いのyの値も等しいという事になります。
この事から(1)と(3)より(yを比較して)、
「 a~3 - N = 3a~3 - b 」が成り立ちます。
これをbで整理すると、「 b = 2a~3 + N 」***(4)***
と、なります。

これによって直線(2)に(4)を代入する事でbを「2a~3 + N」で表現できます。
(つまりNの値が変わると(2)のy切片の値が変わる…! ここが重要です。)

さて代入後の式は、
「y = (3a~2)x - 2a~3 - N」
です。
(これは当然直線です。ここも注意してください。すなわち、(3a~2)が傾きで、残りの「-2a~3 - N」がy切片です。)

この式は直線なので「 y=0 」とする事で 「 x=? 」(解の座標)を求める事が出来ます。
省略して解を書くと「 x = 1/3(n/a~2 + 2a)」***(5)***となります。

さて、この時、aが求める値よりも少し小さいとしてグラフを書いてみてください。
(想像するだけでも良いんですが、書くとよく分かります。)
真の値をAとするならば「a < A < (5)」になります。(確認しておきますが、aはx座標の値です。)

ここでNを10としてaを2と置くならば、(5)は2.16666…になります。
ぉお!かなり近いですね。!
そう、接線のy軸(y=0)との交点のx座標は元の関数のそれと非常に近いんです。
という事は得られた「x=a=2.16666…」をもう一度(5)に代入すると……?!

はい、もうお分かりですね。^^;
そう、a に 1/3(n/a~2 + 2a)を代入すればいいのです。
これを何回か繰り返せばかなり近似値に近づきます。

というわけで、
実際のプログラムはこのようになります。

#include <stdio.h>
#define KAZU 10

double sanjokon(double x)
{
double a;
int i,n;

for (a=1; a*a*a<x; a++){
}
a--;

for (i=0; i<10; i++) {
a = (x/(a*a)+2*a)/3;
}

return a;
}

main()
{
double x=KAZU;

printf("%f", sanjokon(x));
}

ふぅ…

No.1574

Re:3乗根を解くプログラムに関して
投稿者---ジャスミン茶(2002/05/22 03:58:58)


0、負、小数点以下の数字にも対応しました。

#include <stdio.h>
#define KAZU -0.1

double sanjokon(double x)
{
	double a;
	int i,n;

	if (x==0)
		return 0;
	if (x>0) {
		for (a=1; a*a*a<x; a=a+0.1) {
		}
		a=a-0.5;
	} else {
		for (a=-1; a*a*a>x; a=a-0.1) {
		}
		a=a+0.5;
	}

	for (i=0; i<10; i++) {
		a = (x/(a*a)+2*a)/3;
	}

	return a;
}

main()
{
	double x=KAZU;

	printf("%f", sanjokon(x));
}


やっと<pre>の使い方が分かった…

No.1575

Re:3乗根を解くプログラムに関して
投稿者---ジャスミン茶(2002/05/22 04:22:20)


すいません、再び修正します。
なぜか使わない int n がありました。
さらに最初のforの後のa=a +- xxxも不要です。
(微調整しようと思ったのですが、forループの収束の方が圧倒的に効果がありました。)

#include <stdio.h>
#define KAZU -0.1

double sanjokon(double x)
{
	double a;
	int i;

	if (x==0)
		return 0;
	if (x>0) {
		for (a=0.5; a*a*a<x; a+=1) {}
	} else {
		for (a=-0.5; a*a*a>x; a-=1) {}
	}

	for (i=0; i<10; i++) {
		a = (x/(a*a)+2*a)/3;
	}

	return a;
}

main()
{
	double x=KAZU;

	printf("%f", sanjokon(x));
}


スレ汚してすいません。

No.1576

Re:3乗根を解くプログラムに関して
投稿者---ジャスミン茶(2002/05/22 06:32:36)


度々本当にすいません。
研究の結果、aの初期値が0でなければループ回数を大きく取ることで、
かなり真の値に近づく事が分かりました。
前のaの調整は引数の絶対値が大きい時に数千、数万回のループが発生してしまい効率が良くないです。
というわけでaを1に決め打ちしてループ回数を増やしました。
これで実行速度が上がると思います。^^;
(引数があまりにも極端に小さい場合aが0になってやばいかな?、
と思ったのですがとりあえず大丈夫でした。
あまりに小さいとx==0としてはじかれますし、そうでないときはaがおかしくなるだけでした。
逆に大きすぎる場合ですと、環境によってはエラーになります。
共に9~1000とか1~-1000ぐらいで試しました。
これはscanfが悪いんだと思います。)

>>かずまさんへ

はじめまして。
失礼ながらいくつかかずまさんのコードの表現を使わせてもらいました。
あと、知らない関数や表現がありましたので、参考にさせてもらいます。
(それと、かずまさんのですと0の時にエラーが出ますのでreturn 0を返した方がよろしいかと思います。生意気言ってスミマセン。)

#include <stdio.h>

double sanjokon(double x)
{
	double a=1;
	int i;

	if (x==0)
		return 0;
	for (i=0; i<200; i++) {
		a = (x/(a*a)+2*a)/3;
	}
	return a;
}

main()
{
	double x;

	puts("数字を入力してください。\n数字以外を入力した時は終了します。\n");
	while (scanf("%lf", &x) == 1) {
		printf("%.30g の三乗根は %.30g です。\n\n", x, sanjokon(x));
	}
	puts("*** End ***");
	return 0;
}


これで追加修正は最後にします。(多分)
…初心者の成長の様子が手にとるように分かったんじゃないでしょうか。(汗
みんなこうやってうまくなるんでしょう。きっと。(恥
ぁぅー…

No.1611

Re:3乗根を解くプログラムに関して
投稿者---かずま(2002/05/28 20:31:56)


> 研究の結果、aの初期値が0でなければループ回数を大きく取ることで、
> かなり真の値に近づく事が分かりました。

> for (i=0; i<200; i++) {

ループ回数は、なぜ 200 なのでしょうか。

cubic_root(1000) は 15回しかループを回りません。
sanjokon(1000) は 200回ループを回ります。

cubic_root(1e54) は 209回ループを回って、1e+18 という答えを出します。
sanjokon(1e54) は 200回ループを回って、3.04787781156801e+18 という答えを出します。
double cubic_root(double a)
{
    double x, x2 = a;

    if (a == 0)
        return 0;
    do {
        x = x2;
        x2 = (a/(x*x) + x + x) * (1.0/3.0);
    } while (fabs((x2 - x) / x) > 1e-15);
    return x2;
}
fabs((x2 - x) / x) > 1e-15 で、10進15桁の精度を保証しています。
初期値を a でなく、1 にすると、1、2回余計にループを回ります。
ループ内の式は、簡略化して、乗除算の回数を減らしました。

No.1615

Re:3乗根を解くプログラムに関して
投稿者---ジャスミン茶(2002/05/29 09:28:50)


>> for (i=0; i<200; i++) {
>
>ループ回数は、なぜ 200 なのでしょうか。

えと、このスレッドの目的は三乗根を求めるプログラムを解説することなので、
できる限り理解しやすい簡潔なコードを目指したからです。
そのために変数も1つしか使っていませんし、その初期値も1です。
(実際はカウンタ用にint型をもう一つ使っていますが…)
最初はaにxを代入しようと思ったのですが、分かり難いと思って止めました。
ループ回数が大きくなってしまいましたが、かなり理解しやすいコードになったのではないでしょうか?
あれ以上に短くするのは私にはもう無理です。
また、標準ライブラリの関数も使用していません。
(その方が理解が早いと思います。分からない事を調べていたら分からない事が出て来て、またそれを調べるような事態を避けたかったのです。)


>cubic_root(1000) は 15回しかループを回りません。
>sanjokon(1000) は 200回ループを回ります。

うーん、比較されても困るのですが、これは私のコードが良くないという事でしょうか?
あくまでも解説用ですので…(まだ初心者ですし^^;)
ただ、かずまさんの場合はループの回数は引数次第ですよね。


>cubic_root(1e54) は 209回ループを回って、1e+18 という答えを出します。
>sanjokon(1e54) は 200回ループを回って、3.04787781156801e+18 という答えを出します。
>
double cubic_root(double a)
{
    double x, x2 = a;

    if (a == 0)
        return 0;
    do {
        x = x2;
        x2 = (a/(x*x) + x + x) * (1.0/3.0);
    } while (fabs((x2 - x) / x) > 1e-15);
    return x2;
}
fabs((x2 - x) / x) > 1e-15 で、10進15桁の精度を保証しています。
>初期値を a でなく、1 にすると、1、2回余計にループを回ります。
>ループ内の式は、簡略化して、乗除算の回数を減らしました。

えーと、本気で実行速度を求めるのであれば、テーブル方式にして引数の値で分岐して、
初期値を近似値にセットしてしまえば、わずか数回のループですんでしまうと思います。
(ただ、この辺は実行ファイルの大きさ、メモリ消費量とのトレードオフになります。
さらに言えば戻り値はdoubleなので、引数の大きさによっては桁の精度を緩くしてもいいと思います。
場合分けを行い、回数を決め打ちすれば更に速くなるでしょう。)
えと、それからwhileの条件式もループに含まれるので、ループ回数が少なくても実際の演算量は増えますから、ループが少なければ無条件に速いということはないと思います。
また、外部関数を呼び出していますのでメモリも多少食うかと。


>(fabs((x2 - x) / x) > 1e-15)
この書き方初めて知りました。
参考になります。




No.1617

Re:3乗根を解くプログラムに関して
投稿者---かずま(2002/05/29 23:01:13)


>>ループ回数は、なぜ 200 なのでしょうか。
>
>えと、このスレッドの目的は三乗根を求めるプログラムを解説することなので、
>できる限り理解しやすい簡潔なコードを目指したからです。

ループ回数が 200だと、理解しやすい簡潔なコードになるのでしょうか。


>ループ回数が大きくなってしまいましたが、かなり理解しやすいコードになったのではないでしょうか?

200 は大きいのでしょうか。これでは、求められない値もありますよね。
なぜ 200 という値をループ回数に選んだのかが知りたかっただけなんですが、
まあ、いいでしょう。とりあえず、今は、整数部が 50桁程度の値に対しては、
正しい解を求めることができるということで納得しておきます。失礼しました。



No.1622

Re:3乗根を解くプログラムに関して
投稿者---ジャスミン茶(2002/05/30 04:54:41)


>ループ回数が 200だと、理解しやすい簡潔なコードになるのでしょうか。

いえ、そういうわけではありません。
a = (x/(a*a)+2*a)/3;
という三乗根の近似値を求める単純な計算式があり、それを用いているという事を示しただけです。
ループ回数には特に意味はありません。
私のコードですと、引数とaの初期値、そしてループの回数が結果の精度に影響します。
私が伝えたかったのは、なぜ三乗根の計算が上記の簡単な計算式で求められるのかということなんです。
ループが200回なのは、それだけあればある程度の引数までは精度が期待出来ると思ったからです。
この回数は基本的に任意であって、いくつでも良いんです。
引数がaの初期値に近ければ10回位のループでいいと思います。


>200 は大きいのでしょうか。これでは、求められない値もありますよね。
>なぜ 200 という値をループ回数に選んだのかが知りたかっただけなんですが、
>まあ、いいでしょう。とりあえず、今は、整数部が 50桁程度の値に対しては、
>正しい解を求めることができるということで納得しておきます。失礼しました。

うーん、なぜそんなに突っかかるのでしょう?^^;
確かに実用性のないコードだと思いますが、実際にこれを使ってくれと言っている訳ではないですし、
あくまでも原理を理解してもらうためのコードなんですが…
200という数字は適当です。100でも良いです。
ループの回数と引数と初期値が全てですから。
実際のところ、関数側は引数は分からないですから初期値も決めようがありません。
引数にあわせて初期値を決定するのは三乗根の説明には不要です。
そのためにループを大きく取って精度を高くしたのです。
(実を言えば、200回位で50桁程の精度が期待できるため、実際はそのように調整しました。
しかし、三乗根を求める計算の説明とは本質的に関係は無いと思います。)

かずまさんのコードは誤差が一定以下になるまでループを繰り返すので、
ループの回数は不定ですよね。
それに対して私のは誤差に関係なく一定の回数だけループを繰り返します。
要するに、私のは誤差が出てしまうから良くない、という事なのでしょうか?
それならば、かずまさんのようにすれば問題はなくなりますが…
いずれにせよ、「三乗根を求める計算を解説する」という事とは関係がないように思われます。

No.1661

Re:すいませんありがとうございました
投稿者---NEPIA(2002/06/04 20:52:14)


教えてもらっておきながら、返事遅れてすいません。パソコンが壊れてしまって。この質問に対していろいろ議論していただいて大変勉強になりました。

戻る


「初心者のためのポイント学習C言語」 Last modified:2002.06.28
Copyright(c) 2000-2002 TOMOJI All Rights Reserved